Буду очень благодарен тому, кто поможет исправить ошибки в моей попытке организовать умножение длинных чисел с Быстрым Преобразованием Фурье. Само преобразование понимаю лишь отчасти, хотя перечитал литературы... От Википедии до Кнута) В итоге пробовал перевести этот алгоритм на С++ на Делфи. Подумал, что среди людей знающих плюсы больше тех кто знает паскаль, чем наоборот) Поэтому тема здесь.
Заранее благодарен!!
Код:
int rev[MAXN];
base wlen_pw[MAXN];
void fft (base a[], int n, bool invert) {
for (int i=0; i<n; ++i)
if (i < rev[i])
swap (a[i], a[rev[i]]);
for (int len=2; len<=n; len<<=1) {
double ang = 2*PI/len * (invert?-1:+1);
int len2 = len>>1;
base wlen (cos(ang), sin(ang));
wlen_pw[0] = base (1, 0);
for (int i=1; i<len2; ++i)
wlen_pw[i] = wlen_pw[i-1] * wlen;
for (int i=0; i<n; i+=len) {
base t,
*pu = a+i,
*pv = a+i+len2,
*pu_end = a+i+len2,
*pw = wlen_pw;
for (; pu!=pu_end; ++pu, ++pv, ++pw) {
t = *pv * *pw;
*pv = *pu - t;
*pu += t;
}
}
}
if (invert)
for (int i=0; i<n; ++i)
a[i] /= n;
}
void calc_rev (int n, int log_n) {
for (int i=0; i<n; ++i) {
rev[i] = 0;
for (int j=0; j<log_n; ++j)
if (i & (1<<j))
rev[i] |= 1<<(log_n-1-j);
}
}
void multiply (const vector<int> & a, const vector<int> & b, vector<int> & res) {
vector<base> fa (a.begin(), a.end()), fb (b.begin(), b.end());
size_t n = 1;
while (n < max (a.size(), b.size())) n <<= 1;
n <<= 1;
fa.resize (n), fb.resize (n);
fft (fa, false), fft (fb, false);
for (size_t i=0; i<n; ++i)
fa[i] *= fb[i];
fft (fa, true);
res.resize (n);
for (size_t i=0; i<n; ++i)
res[i] = int (fa[i].real() + 0.5);int carry = 0;
for (size_t i=0; i<n; ++i) {
res[i] += carry;
carry = res[i] / 10;
res[i] %= 10;
}
}
(брал с сайта:
http://e-maxx.ru/algo/fft_multiply)
Мой выводящий нули вариант на паскале... а иногда вообще вылетающий в рантайм(
Код:
type complex = record
re,im:double;
end;
type arcomp = array of complex;
type arint = array of integer;
var a,b:arcomp;
r:arint;
s:string;
i,n1,n2,t:integer;
rev: arint;
wlen_pw:arcomp;
function umnoj(a,b:complex):complex;
begin
umnoj.re:=a.re*b.re-a.im*b.im;
umnoj.im:=a.re*b.im+a.im*b.re;
end;
function summ(a,b:complex):complex;
begin
summ.re:=a.re+b.re;
summ.im:=a.im+b.im;
end;
function diff(a,b:complex):complex;
begin
diff.re:=a.re-b.re;
diff.im:=a.im-b.im;
end;
function base(a,b:double):complex;
begin
base.re:=a;
base.im:=b;
end;
procedure fft (a:arcomp;invert:boolean; n:integer );
var len,len2,k,j,d:integer;
ang:double;
wlen,t:complex;
wlen_pw:arcomp;
procedure swap(var a,b:complex);
var c:complex;
begin
c:=a;
a:=b;
b:=c;
end;
begin
for i:=0 to n do
if i<rev[i] then
swap (a[i], a[rev[i]]);
len:=2;
setlength(wlen_pw,n);
while (len<=n)do begin
if invert then
ang := (-2)*PI/len
else
ang := 2*PI/len;
len2 := len shr 1;
wlen.re:=cos(ang);
wlen.im:=sin(ang);
wlen_pw[0]:= base(1, 0);
for i:=1 to len2-1 do
wlen_pw[i] := umnoj(wlen_pw[i-1],wlen);
i:=0;
while i<n do begin
k:=0; j:=0; d:=0;
while k<>len2 do begin
t:=umnoj(a[i+len2+j],wlen_pw[d]);
a[i+len2+j]:=diff(a[i+k],t);
inc(k); inc(j); inc(d);
end;
i:=i+len;
end;
len:=len shl 1;
end;
if invert then
for i:=0 to n do begin
a[i].re := a[i].re/n;
a[i].im := a[i].im/n;
end;
end;
procedure calc_rev (n,log_n:integer);
var i,j:integer;
begin
for i:=0 to n-1 do begin
rev[i]:=0;
for j:=0 to log_n-1 do
if (i and (1 shl j))<>0 then
rev[i]:=rev[i]or (1 shl (log_n-1-j));
end;
end;
procedure multiply (a,b:arcomp; var res: arint;len1,len2:integer);
var n,_as,bs:integer;
i,carry,k:integer;
fa,fb:arcomp;
function max(a,b:integer):integer;
begin
if a>b then max:=a else max:=b;
end;
begin
n:=1;
setlength(fa,len1);
setlength(fb,len2);
for i:=0 to len1-1 do
fa[i]:=a[i];
for i:=0 to len2-1 do
fb[i]:=b[i];
_as:=len1;
bs:=len2;
while (n<max(_as,bs)) do
n:=n*2;
n:=n*2;
setlength(fa,n);
setlength(fb,n);
setlength(rev,n);
calc_rev(n,round(ln(n)));
fft(fb,false,n);
fft(fa,false,n);
for i:= 0 to n-1 do
fa[i]:= umnoj(fa[i],fb[i]);
fft(fa,true,n);
setlength(res,n);
for i:= 0 to n-1 do
res[i] := round(fa[i].re + 0.5);
carry:=0;
for i:= 0 to n-1 do begin
res[i]:=res[i]+carry;
carry := round(res[i] / 10);
res[i] :=res[i] div 10;
end;
k:=n-1;
while (res[k] = 0) do
dec(k);
setlength(res,k+1);
end;