Sz\303\241m\303\255t\303\263g\303\251pes sz\303\241melm\303\251let J\303\241rai Antal Ezek a programok csak szeml\303\251ltet\303\251sre szolg\303\241lnak

LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEhRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0Yn

restart;

# # This procedure do the bit reversation of the # number x which is k bit long. # reverse:=proc(x,k) local xx,i,y; xx:=x; y:=0; for i to k do if type(xx,odd) then y:=2*y+1; xx:=(xx-1)/2; else y:=2*y; xx:=xx/2; fi; od; y; end;

Zio2JEkieEc2Ikkia0dGJTYlSSN4eEdGJUkiaUdGJUkieUdGJUYlRiVDJj5GKEYkPkYqIiIhPyhGKSIiIkYwRiZJJXRydWVHJSpwcm90ZWN0ZWRHQCUtSSV0eXBlR0YyNiRGKEkkb2RkR0YyQyQ+RiosJiomIiIjRjBGKkYwRjBGMEYwPkYoLCYqJiNGMEY8RjBGKEYwRjBGQCEiIkMkPkYqLCRGO0YwPkYoLCRGP0YwRipGJUYlRiU=

# # This procedure do the complex butterfly operation # between the two terms A[i] and A[j] of the # table A. The multiplier is w. # cbutterfly:=proc(A,i,j,w) local X,Y; X:=A[i]; Y:=A[j]*w; A[i]:=X+Y; A[j]:=X-Y; end;

Zio2JkkiQUc2IkkiaUdGJUkiakdGJUkid0dGJTYkSSJYR0YlSSJZR0YlRiVGJUMmPkYqJkYkNiNGJj5GKyomJkYkNiNGJyIiIkYoRjQ+Ri4sJkYqRjRGK0Y0PkYyLCZGKkY0RishIiJGJUYlRiU=

# # This is a complex FFT procedure. # It use the butterfly procedure to operate on the vector A. # The number of rounds is k in the FFT. # T is a table of the powers of primitive root of unity. # cfft:=proc(A,T,k) local l,s,w,t; for l from 0 to k-1 do for s from 0 to 2^l-1 do w:=T[s]; for t from 0 to 2^(k-l-1)-1 do cbutterfly(A,2^(k-l)*s+t,2^(k-l)*s+t+2^(k-l-1),w); od; od; od; end;

Zio2JUkiQUc2IkkiVEdGJUkia0dGJTYmSSJsR0YlSSJzR0YlSSJ3R0YlSSJ0R0YlRiVGJT8oRikiIiEiIiIsJkYnRi9GLyEiIkkldHJ1ZUclKnByb3RlY3RlZEc/KEYqRi5GLywmKSIiI0YpRi9GL0YxRjJDJD5GKyZGJjYjRio/KEYsRi5GLywmKUY3LChGJ0YvRilGMUYvRjFGL0YvRjFGMi1JK2NidXR0ZXJmbHlHRiU2JkYkLCYqJilGNywmRidGL0YpRjFGL0YqRi9GL0YsRi8sKEZERi9GLEYvRj5GL0YrRiVGJUYl

# # The pre procedure do the preparation for the # table T[0..2^k-1)] of powers of the primitive # root of unity. # pre:=proc(T,k) local i,xx; Digits:=2*Digits; for i from 0 to 2^k-1 do T[i]:=evalf(exp(-2*Pi*I*reverse(i,k)/2^(k+1))); od; Digits:=Digits/2; end;

Zio2JEkiVEc2Ikkia0dGJTYkSSJpR0YlSSN4eEdGJUYlRiVDJT5JJ0RpZ2l0c0dGJSwkKiYiIiMiIiJGLEYwRjA/KEYoIiIhRjAsJilGL0YmRjBGMCEiIkkldHJ1ZUclKnByb3RlY3RlZEc+JkYkNiNGKC1JJmV2YWxmR0Y3NiMtSSRleHBHRiU2IywkKiosJComRi9GMF4jRjBGMEYwRjBJI1BpR0Y3RjAtSShyZXZlcnNlR0YlNiRGKEYmRjApRi8sJkYmRjBGMEYwRjVGNT5GLCwkKiYjRjBGL0YwRixGMEYwRiVGJUYl

k:=5; A:=array(0..2^k-1); for i from 0 to 2^k-1 do A[i]:=0 od: T:=array(0..2^(k-1)-1); pre(T,k-1):

IiIm

PTYiNiM7IiIhIiNKRVxbbCE=

PTYiNiM7IiIhIiM6RVxbbCE=

A[1]:=1+I; cfft(A,T,5): print(A);

LCYiIiJGI14jRiNGIw==

SSJBRzYi

for i from 0 to 2^k-1 do j:=reverse(i,k); if i<j then x:=A[i]; A[i]:=A[j]; A[j]:=x; fi; od: print(A);

SSJBRzYi

cfft(A,T,5): print(A);

SSJBRzYi

for i from 0 to 2^k-1 do j:=reverse(i,k); if i<j then x:=A[i]; A[i]:=A[j]; A[j]:=x; fi; od: print(A);

SSJBRzYi

# # This procedure do the inverse complex butterfly operation # between the two terms A[i] and A[j] of the # table A. The power of primitive unity is w. # icbutterfly:=proc(A,i,j,w) local X,Y,e; X:=A[i]+A[j]; Y:=A[i]-A[j]; A[i]:=X; A[j]:=Y/w; end;

Zio2JkkiQUc2IkkiaUdGJUkiakdGJUkid0dGJTYlSSJYR0YlSSJZR0YlSSJlR0YlRiVGJUMmPkYqLCYmRiQ2I0YmIiIiJkYkNiNGJ0YyPkYrLCZGMEYyRjMhIiI+RjBGKj5GMyomRitGMkYoRjdGJUYlRiU=

# # This procedure is a complex IFFT procedure. # It use the ibutterfly procedure to operate on the vector A. # The number of round in the FFT is k and T is a table of the # powers of primitive root of unity. # icfft:=proc(A,T,k) local l,s,w,t; for l from k-1 to 0 by -1 do for s from 0 to 2^l-1 do w:=T[s]; for t from 0 to 2^(k-l-1)-1 do icbutterfly(A,2^(k-l)*s+t,2^(k-l)*s+t+2^(k-l-1),w); od; od; od; end;

Zio2JUkiQUc2IkkiVEdGJUkia0dGJTYmSSJsR0YlSSJzR0YlSSJ3R0YlSSJ0R0YlRiVGJT8oRiksJkYnIiIiRi8hIiJGMCIiIUkldHJ1ZUclKnByb3RlY3RlZEc/KEYqRjFGLywmKSIiI0YpRi9GL0YwRjJDJD5GKyZGJjYjRio/KEYsRjFGLywmKUY3LChGJ0YvRilGMEYvRjBGL0YvRjBGMi1JLGljYnV0dGVyZmx5R0YlNiZGJCwmKiYpRjcsJkYnRi9GKUYwRi9GKkYvRi9GLEYvLChGREYvRixGL0Y+Ri9GK0YlRiVGJQ==

for i from 0 to 2^k-1 do A[i]:=0 od: A[1]:=1+I; cfft(A,T,5): print(A);

LCYiIiJGI14jRiNGIw==

SSJBRzYi

icfft(A,T,k): print(A);

SSJBRzYi

# # The cdigmul procedure do the digit-by-digit # multiplication of the two numbers after the # cfft's. The result will be in the first table. # cdigmul:=proc(T,S,k) local i; for i from 0 to 2^k-1 do T[i]:=T[i]*S[i]; od; end;

Zio2JUkiVEc2IkkiU0dGJUkia0dGJTYjSSJpR0YlRiVGJT8oRikiIiEiIiIsJikiIiNGJ0YsRiwhIiJJJXRydWVHJSpwcm90ZWN0ZWRHPiZGJEYoKiZGNEYsJkYmRihGLEYlRiVGJQ==

A:=array(0..2^k-1); for i from 0 to 2^k -1 do A[i]:=0 od: B:=array(0..2^k-1); for i from 0 to 2^k -1 do B[i]:=0 od:

PTYiNiM7IiIhIiNKRVxbbCE=

PTYiNiM7IiIhIiNKRVxbbCE=

A[0]:=1: A[1]:=1: A[2]:=1: print(A);

SSJBRzYi

# # This is the polynom multiplication, do multiplication or squaring. # A and B are the two polynomials, the FFT and IFFT use k rounds. # polmulcfft:=proc(A,B,k) global T; if A=B then cfft(A,T,k); cdigmul(A,A,k); icfft(A,T,k); else cfft(A,T,k); cfft(B,T,k); cdigmul(A,B,k); icfft(A,T,k); fi; end;

Zio2JUkiQUc2IkkiQkdGJUkia0dGJUYlRiVGJUAlL0YkRiZDJS1JJWNmZnRHRiU2JUYkSSJUR0YlRictSShjZGlnbXVsR0YlNiVGJEYkRictSSZpY2ZmdEdGJUYtQyZGKy1GLDYlRiZGLkYnLUYwRiNGMkYlNiNGLkYl

polmulcfft(A,A,k): print(map(x->x/2^k,A));

PTYiNiM7IiIhIiNKRVxbbEFGJiwmJCIrJSoqKioqKioqKiEjNSIiIiomJEYmRiZGLV4jRi1GLUYtRi0sJiQiKyoqKioqKioqPiEiKkYtRi5GLSIiIywmJCIrKioqKioqKipIRjRGLUYuRi0iIiQsJiQiKysrKys/RjRGLUYuRi0iIiYsJiQiKz5FZSEzKSEjP0YtRi5GLSIiJSwmJCIrKysrKzVGNEYtRi5GLSIiKCwmJCIrSXpzXkEhIz5GLUYuRi0iIicsJiQiK19dTGg1RkohIiJGLkYtIiM1LCYkIis+O0I+OUZKRk9GLkYtIiM2LCYkIitBTTRWaUZKRk9GLkYtIiIpLCZGL0YtRi5GLSIiKiwmJCIrKysrdj1GSkZPRi5GLSIjOiwmJCIrUT1FKHkmRkpGLUYuRi0iIzksJiQiK1FpaCQqb0ZKRi1GLkYtIiM4LCYkIis6XyNlSyJGSkYtRi5GLSIjNywmJCIrd00pKWVGRkpGT0YuRi0iI0AsJiQiK1E8JT5wIkZKRi1GLkYtIiM/LCYkIisrKytdaUZKRi1GLkYtIiNCLCYkIitELVJnd0ZBRi1GLkYtIiNBLCYkIitHJj1fJT5GSkYtRi5GLSIjO0ZdcCIjPEZdcCIjPUZdcCIjPiwmJCIrKysrREpGSkYtRi5GLUYnLCYkIispeUddSSdGSkZPRi5GTyIjSSwmJCIrNygqXHh4RkpGT0YuRk8iI0gsJkZib0ZPRi5GLSIjRywmJCIrUV87NioqRkFGT0YuRi0iI0UsJkZSRi1GLkYtIiNGLCYkIisrdmQxcCEjQUZPRi5GLSIjQ0ZZIiNELCZGZm5GLUYuRi0=

# # The cfftpre procedure do the preparation for the # table for cfft, where x is the number, A is the # table, m the modulus, and k is the number of # rounds for the cfft, hence the table is 2^k long. # cfftpre:=proc(x,A,m,k) local i,xx; xx:=x; for i from 0 to 2^k-1 do A[i]:=irem(xx,m,'xx'); od; end;

Zio2JkkieEc2IkkiQUdGJUkibUdGJUkia0dGJTYkSSJpR0YlSSN4eEdGJUYlRiVDJD5GK0YkPyhGKiIiISIiIiwmKSIiI0YoRjBGMCEiIkkldHJ1ZUclKnByb3RlY3RlZEc+JkYmNiNGKi1JJWlyZW1HRjY2JUYrRicuRitGJUYlRiU=

# # The cnorm procedure do the normalization # after the icfft; A is the table, m the modulus, # and k is the number of rounds for the cfft, # hence the table is 2^k long. The fraction parts are # left in the table A. # cnorm:=proc(A,m,k) local i,x; x:=0; for i from 2^k-1 to 0 by -1 do A[i]:=A[i]/2^k; x:=m*x+round(A[i]); A[i]:=A[i]-round(A[i]); od; x; end;

Zio2JUkiQUc2IkkibUdGJUkia0dGJTYkSSJpR0YlSSJ4R0YlRiVGJUMlPkYqIiIhPyhGKSwmKSIiI0YnIiIiRjIhIiJGM0YtSSV0cnVlRyUqcHJvdGVjdGVkR0MlPiZGJDYjRikqJkY4RjJGMEYzPkYqLCYqJkYmRjJGKkYyRjItSSZyb3VuZEc2JEY1SShfc3lzbGliR0YlNiNGOEYyPkY4LCZGOEYyRj5GM0YqRiVGJUYl

# # This is the main procedure, do the multiplication or the squaring. # a and b are the two numbers, m is the modulus for the preparation. # The complex FFT and IFFT use k rounds. # mulcfft:=proc(a,b,m,k) global A,B,T; if a=b then cfftpre(a,A,m,k); cfft(A,T,k); cdigmul(A,A,k); icfft(A,T,k); cnorm(A,m,k); else cfftpre(a,A,m,k); cfft(A,T,k); cfftpre(b,B,m,k); cfft(B,T,k); cdigmul(A,B,k); icfft(A,T,k); cnorm(A,m,k); fi; end;

Zio2JkkiYUc2IkkiYkdGJUkibUdGJUkia0dGJUYlRiVGJUAlL0YkRiZDJy1JKGNmZnRwcmVHRiU2JkYkSSJBR0YlRidGKC1JJWNmZnRHRiU2JUYvSSJUR0YlRigtSShjZGlnbXVsR0YlNiVGL0YvRigtSSZpY2ZmdEdGJUYyLUkmY25vcm1HRiU2JUYvRidGKEMpRixGMC1GLTYmRiZJIkJHRiVGJ0YoLUYxNiVGP0YzRigtRjU2JUYvRj9GKEY3RjlGJTYlRi9GP0YzRiU=

mulcfft(123456789,987654321,20,5);

IjNwX2o3NmpLPjc=

123456789*987654321;

IjNwX2o3NmpLPjc=

print(A);

SSJBRzYi

# # The CtoR procedure do the conversion from complex # representation to real representation. # The result will be in the same table. # CtoR:=proc(A,T,k) local i,x,y; x:=Re(A[0]); y:=Im(A[0]); A[0]:=2*(x+y)+I*2*(x-y); x:=Re(A[1]); y:=Im(A[1]); A[1]:=2*x-I*2*y; for i to k-1 do CtoRsteps(A,T,2^i); od; end;

Zio2JUkiQUc2IkkiVEdGJUkia0dGJTYlSSJpR0YlSSJ4R0YlSSJ5R0YlRiVGJUMpPkYqLUkjUmVHJSpwcm90ZWN0ZWRHNiMmRiQ2IyIiIT5GKy1JI0ltR0YwRjE+RjIsKComIiIjIiIiRipGPEY8KiZGO0Y8RitGPEY8KiYqJkY7RjxeI0Y8RjxGPCwmRipGPEYrISIiRjxGPD5GKi1GLzYjJkYkNiNGPD5GKy1GN0ZFPkZGLCZGOkY8KiYsJEY/RjxGPEYrRjxGQj8oRilGPEY8LCZGJ0Y8RjxGQkkldHJ1ZUdGMC1JKkN0b1JzdGVwc0dGJTYlRiRGJilGO0YpRiVGJUYl

# # The CtoRstep procedure do the conversion from complex # representation to real representation for one pair # ll, uu with weight w. # CtoRstep:=proc(ll,uu,w) local al,be,ga,de,xi,et,a,b,x,y,u,v; xi:=Re(w); et:=Im(w); al:=Re(ll); be:=Im(ll); ga:=Re(uu); de:=Im(uu); x:=al+ga; y:=be-de; a:=al-ga; b:=be+de; u:=et*b-xi*a; v:=xi*b+et*a; [x+v+I*(y+u),x-v+I*(u-y)] end;

Zio2JUkjbGxHNiJJI3V1R0YlSSJ3R0YlNi5JI2FsR0YlSSNiZUdGJUkjZ2FHRiVJI2RlR0YlSSN4aUdGJUkjZXRHRiVJImFHRiVJImJHRiVJInhHRiVJInlHRiVJInVHRiVJInZHRiVGJUYlQy8+Ri0tSSNSZUclKnByb3RlY3RlZEc2I0YnPkYuLUkjSW1HRjlGOj5GKS1GODYjRiQ+RiotRj1GQD5GKy1GODYjRiY+RiwtRj1GRT5GMSwmRikiIiJGK0ZKPkYyLCZGKkZKRiwhIiI+Ri8sJkYpRkpGK0ZNPkYwLCZGKkZKRixGSj5GMywmKiZGLkZKRjBGSkZKKiZGLUZKRi9GSkZNPkY0LCYqJkYtRkpGMEZKRkoqJkYuRkpGL0ZKRko3JCwoRjFGSkY0RkoqJl4jRkpGSiwmRjJGSkYzRkpGSkZKLChGMUZKRjRGTSomRmduRkosJkYzRkpGMkZNRkpGSkYlRiVGJQ==

# # The CtoRsteps procedure do the conversion from complex # representation to real representation for one series. # The index i is the lower index for the first pair. # The result will be in the same table. # CtoRsteps:=proc(A,T,i) local k,j,z; k:=i; j:=2*i-1; while j>k do z:=CtoRstep(A[k],A[j],T[k]); A[k]:=z[1]; A[j]:=z[2]; k:=k+1; j:=j-1; od; end;

Zio2JUkiQUc2IkkiVEdGJUkiaUdGJTYlSSJrR0YlSSJqR0YlSSJ6R0YlRiVGJUMlPkYpRic+RiosJiomIiIjIiIiRidGMkYyRjIhIiI/KEYlRjJGMkYlMkYpRipDJz5GKy1JKUN0b1JzdGVwR0YlNiUmRiQ2I0YpJkYkNiNGKiZGJkY8PkY7JkYrNiNGMj5GPSZGKzYjRjE+RiksJkYpRjJGMkYyPkYqLCZGKkYyRjJGM0YlRiVGJQ==

# # The RtoC procedure do the conversion from real # representation to complex representation. # The result will be in the same table. # RtoC:=proc(A,T,k) local i,x,y; x:=Re(A[0]); y:=Im(A[0]); A[0]:=x+y+I*(x-y); x:=Re(A[1]); y:=Im(A[1]); A[1]:=2*x-I*2*y; for i to k-1 do RtoCsteps(A,T,2^i); od; end;

Zio2JUkiQUc2IkkiVEdGJUkia0dGJTYlSSJpR0YlSSJ4R0YlSSJ5R0YlRiVGJUMpPkYqLUkjUmVHJSpwcm90ZWN0ZWRHNiMmRiQ2IyIiIT5GKy1JI0ltR0YwRjE+RjIsKEYqIiIiRitGOiomXiNGOkY6LCZGKkY6RishIiJGOkY6PkYqLUYvNiMmRiQ2I0Y6PkYrLUY3RkE+RkIsJiomIiIjRjpGKkY6RjoqJiwkKiZGSUY6RjxGOkY6RjpGK0Y6Rj4/KEYpRjpGOiwmRidGOkY6Rj5JJXRydWVHRjAtSSpSdG9Dc3RlcHNHRiU2JUYkRiYpRklGKUYlRiVGJQ==

# # The RtoCstep procedure do the conversion from real representation # to complex representation for one pair ll, uu using weight w. # RtoCstep:=proc(ll,uu,w) local al,be,ga,de,xi,et,a,b,x,y,u,v; xi:=Re(w); et:=Im(w); al:=Re(ll); be:=Im(ll); ga:=Re(uu); de:=Im(uu); x:=al+ga; y:=be-de; a:=al-ga; b:=be+de; u:=xi*a+et*b; v:=xi*b-et*a; [x-v+I*(u+y),x+v+I*(u-y)] end;

Zio2JUkjbGxHNiJJI3V1R0YlSSJ3R0YlNi5JI2FsR0YlSSNiZUdGJUkjZ2FHRiVJI2RlR0YlSSN4aUdGJUkjZXRHRiVJImFHRiVJImJHRiVJInhHRiVJInlHRiVJInVHRiVJInZHRiVGJUYlQy8+Ri0tSSNSZUclKnByb3RlY3RlZEc2I0YnPkYuLUkjSW1HRjlGOj5GKS1GODYjRiQ+RiotRj1GQD5GKy1GODYjRiY+RiwtRj1GRT5GMSwmRikiIiJGK0ZKPkYyLCZGKkZKRiwhIiI+Ri8sJkYpRkpGK0ZNPkYwLCZGKkZKRixGSj5GMywmKiZGLUZKRi9GSkZKKiZGLkZKRjBGSkZKPkY0LCYqJkYtRkpGMEZKRkoqJkYuRkpGL0ZKRk03JCwoRjFGSkY0Rk0qJl4jRkpGSiwmRjJGSkYzRkpGSkZKLChGMUZKRjRGSiomRmduRkosJkYzRkpGMkZNRkpGSkYlRiVGJQ==

# # The RtoCsteps procedure do the conversion from real # representation to complex representation for one series. # The index i is the lower index for the first pair. # The result will be in the same table. # RtoCsteps:=proc(A,T,i) local k,j,z; k:=i; j:=2*i-1; while j>k do z:=RtoCstep(A[k],A[j],T[k]); A[k]:=z[1]; A[j]:=z[2]; k:=k+1; j:=j-1; od; end;

Zio2JUkiQUc2IkkiVEdGJUkiaUdGJTYlSSJrR0YlSSJqR0YlSSJ6R0YlRiVGJUMlPkYpRic+RiosJiomIiIjIiIiRidGMkYyRjIhIiI/KEYlRjJGMkYlMkYpRipDJz5GKy1JKVJ0b0NzdGVwR0YlNiUmRiQ2I0YpJkYkNiNGKiZGJkY8PkY7JkYrNiNGMj5GPSZGKzYjRjE+RiksJkYpRjJGMkYyPkYqLCZGKkYyRjJGM0YlRiVGJQ==

# # The srfftpre procedure do the preparation for the # signed table for rfft; x is the number, A is the # table, m the even positive modulus, and k is the # number of rounds for the rfft, hence the table is 2^k long. # All number is multipled by the factor f. # srfftpre:=proc(x,A,m,k,f) local i,c,d,xx; xx:=x; c:=0; for i from 0 to 2^k-1 do d:=irem(xx,m,'xx'); if d>=m/2 then d:=d-m+c; c:=1; else d:=d+c; c:=0; fi; A[i]:=evalf(d*f); d:=irem(xx,m,'xx'); if d>=m/2 then d:=d-m+c; c:=1; else d:=d+c; c:=0; fi; A[i]:=A[i]+I*evalf(d*f); od; end;

Zio2J0kieEc2IkkiQUdGJUkibUdGJUkia0dGJUkiZkdGJTYmSSJpR0YlSSJjR0YlSSJkR0YlSSN4eEdGJUYlRiVDJT5GLkYkPkYsIiIhPyhGK0YyIiIiLCYpIiIjRihGNEY0ISIiSSV0cnVlRyUqcHJvdGVjdGVkR0MoPkYtLUklaXJlbUdGOjYlRi5GJy5GLkAlMSwkKiYjRjRGN0Y0RidGNEY0Ri1DJD5GLSwoRi1GNEYnRjhGLEY0PkYsRjRDJD5GLSwmRixGNEYtRjRGMT4mRiY2I0YrLUkmZXZhbGZHRjo2IyomRi1GNEYpRjRGPEZBPkZOLCZGTkY0KiZeI0Y0RjRGUEY0RjRGJUYlRiU=

# # The rnorm procedure do the normalization after the irfft. # A is the table, m the modulus, and k is the number of # rounds for the rfft, hence the table is 2^k long. # Before conversion, all entry in A multiplied by the factor f. # The fraction parts are left in the table A. # rnorm:=proc(A,m,k,f) local i,x; x:=0; for i from 2^k-1 to 0 by -1 do A[i]:=evalf(A[i]*f); x:=m*x+round(Im(A[i])); A[i]:=A[i]-I*round(Im(A[i])); x:=m*x+round(Re(A[i])); A[i]:=A[i]-round(Re(A[i])); od; x; end;

Zio2JkkiQUc2IkkibUdGJUkia0dGJUkiZkdGJTYkSSJpR0YlSSJ4R0YlRiVGJUMlPkYrIiIhPyhGKiwmKSIiI0YnIiIiRjMhIiJGNEYuSSV0cnVlRyUqcHJvdGVjdGVkR0MnPiZGJDYjRiotSSZldmFsZkdGNjYjKiZGOUYzRihGMz5GKywmKiZGJkYzRitGM0YzLUkmcm91bmRHNiRGNkkoX3N5c2xpYkdGJTYjLUkjSW1HRjY2I0Y5RjM+RjksJkY5RjMqJiwkXiNGM0YzRjNGQkYzRjQ+RissJkZBRjMtRkM2Iy1JI1JlR0Y2RklGMz5GOSwmRjlGM0ZRRjRGK0YlRiVGJQ==

# # The rdigmul procedure do the digit-by-digit # multiplication of the two numbers after the # rfft's. The result will be in the first table. # rdigmul:=proc(A,B,k) local i; A[0]:=Re(A[0])*Re(B[0])+I*Im(A[0])*Im(B[0]); for i from 1 to 2^k-1 do A[i]:=A[i]*B[i]; od; end;

Zio2JUkiQUc2IkkiQkdGJUkia0dGJTYjSSJpR0YlRiVGJUMkPiZGJDYjIiIhLCYqJi1JI1JlRyUqcHJvdGVjdGVkRzYjRiwiIiItRjI2IyZGJkYtRjVGNSooXiNGNUY1LUkjSW1HRjNGNEY1LUY8RjdGNUY1PyhGKUY1RjUsJikiIiNGJ0Y1RjUhIiJJJXRydWVHRjM+JkYkRigqJkZFRjUmRiZGKEY1RiVGJUYl

# # This is the main procedure, do the multiplication # or the squaring using real FFT and IFFF. # a and b are the two numbers, m is the modulus for # the preparation. The FFT and IFFT use k rounds. # mulrfft:=proc(a,b,m,k) local f1,f2,r; global A,B,T; if type(k,odd) then f1:=2^(-(k+1)/2-HWS); f2:=2^(2*HWS-2); else f1:=2^(-k/2-HWS); f2:=2^(2*HWS-3); fi; if a=b then srfftpre(a,A,m,k,f1); print(A); cfft(A,T,k); print(A); CtoR(A,T,k); print(A); rdigmul(A,A,k); print(A); RtoC(A,T,k); print(A); icfft(A,T,k); print(A); rnorm(A,m,k,f2); else srfftpre(a,A,m,k,f1); cfft(A,T,k); CtoR(A,T,k); srfftpre(b,B,m,k,f1); cfft(B,T,k); CtoR(B,T,k); rdigmul(A,B,k); RtoC(A,T,k); icfft(A,T,k); rnorm(A,m,k,f2); fi; end;

Zio2JkkiYUc2IkkiYkdGJUkibUdGJUkia0dGJTYlSSNmMUdGJUkjZjJHRiVJInJHRiVGJUYlQyRAJS1JJXR5cGVHJSpwcm90ZWN0ZWRHNiRGKEkkb2RkR0YxQyQ+RiopIiIjLCgqJiMiIiJGN0Y7RihGOyEiIkY6RjxJJEhXU0dGJUY8PkYrKUY3LCYqJkY3RjtGPUY7RjtGN0Y8QyQ+RiopRjcsJkY5RjxGPUY8PkYrKUY3LCZGQUY7IiIkRjxAJS9GJEYmQy8tSSlzcmZmdHByZUdGJTYnRiRJIkFHRiVGJ0YoRiotSSZwcmludEdGMTYjRlAtSSVjZmZ0R0YlNiVGUEkiVEdGJUYoRlEtSSVDdG9SR0YlRlZGUS1JKHJkaWdtdWxHRiU2JUZQRlBGKEZRLUklUnRvQ0dGJUZWRlEtSSZpY2ZmdEdGJUZWRlEtSSZybm9ybUdGJTYmRlBGJ0YoRitDLEZNRlRGWC1GTjYnRiZJIkJHRiVGJ0YoRiotRlU2JUZhb0ZXRigtRllGY28tRmVuNiVGUEZhb0YoRmduRmluRltvRiU2JUZQRmFvRldGJQ==

HWS:=16; k:=5; A:=array(0..2^k-1); B:=array(0..2^k-1); T:=array(0..2^k-1); pre(T,k):

IiM7

IiIm

PTYiNiM7IiIhIiNKRVxbbCE=

PTYiNiM7IiIhIiNKRVxbbCE=

PTYiNiM7IiIhIiNKRVxbbCE=

mulrfft(123456789,987654321,18,5);

IjNwX2o3NmpLPjc=

123456789*987654321;

IjNwX2o3NmpLPjc=

print(A);

SSJBRzYi

HWS:=16; k:=15; A:=array(0..2^k-1); B:=array(0..2^k-1); T:=array(0..2^k-1); pre(T,k):

IiM7

IiM6

PTYiNiM7IiIhIiZuRiRFXFtsIQ==

PTYiNiM7IiIhIiZuRiRFXFtsIQ==

PTYiNiM7IiIhIiZuRiRFXFtsIQ==

mulrfft(123456789,987654321,16,5);

IjNwX2o3NmpLPjc=

123456789*987654321;

IjNwX2o3NmpLPjc=

# # This procedure adds 1 to the number represented by the # bit vector b (used in reverse bit order) on length k. # incrementreverse:=proc(b,k) local i,c; c:=b; for i to k do if c[i]=0 then c[i]:=1; return c; else c[i]:=0 fi; od; c; end;

Zio2JEkiYkc2Ikkia0dGJTYkSSJpR0YlSSJjR0YlRiVGJUMlPkYpRiQ/KEYoIiIiRi1GJkkldHJ1ZUclKnByb3RlY3RlZEdAJS8mRik2I0YoIiIhQyQ+RjJGLU9GKT5GMkY0RilGJUYlRiU=

# # This procedure convert the bits in interval II of the # bit vector b. The bit b[1] represents the highest bit, # 2^(m-2+op(1,II)), usually 2^(m-1). # convertreverse:=proc(b,m,II) local i,lw; lw:=0; for i from op(1,II) to op(2,II) do lw:=lw+b[i]*2^(m-1-i+op(1,II)) od; lw; end;

Zio2JUkiYkc2IkkibUdGJUkjSUlHRiU2JEkiaUdGJUkjbHdHRiVGJUYlQyU+RioiIiE/KEYpLUkjb3BHJSpwcm90ZWN0ZWRHNiQiIiJGJ0YzLUYwNiQiIiNGJ0kldHJ1ZUdGMT5GKiwmRipGMyomJkYkNiNGKUYzKUY2LCpGJkYzRjMhIiJGKUY/Ri9GM0YzRjNGKkYlRiVGJQ==

# # This procedure do a Fermat FFT round modulo 2^(2^m)+1 on # array A having 2^n elements. The siccor size is 2^k and # the weights are get from the conversion of bits in interval II # from a counter starting with zero and incremented after # each butterfly sequence. # ffftround:=proc(A,n,m,k,II) local i,j,x,y,w,b; j:=2^k; b:=[0$i=1..m+1]; while j<2^n do w:=2^convertreverse(b,m,II); b:=incrementreverse(b,m+1); for i from j-2^k while i<j do x:=A[i]; y:=A[i+2^k]; A[i]:=x+w*y mod (2^(2^m)+1); A[i+2^k]:=x-w*y mod (2^(2^m)+1); od; j:=j+2^(k+1); od; end;

Zio2J0kiQUc2IkkibkdGJUkibUdGJUkia0dGJUkjSUlHRiU2KEkiaUdGJUkiakdGJUkieEdGJUkieUdGJUkid0dGJUkiYkdGJUYlRiVDJT5GLCkiIiNGKD5GMDcjLUkiJEclKnByb3RlY3RlZEc2JCIiIS9GKzsiIiIsJkYnRj5GPkY+PyhGJUY+Rj5GJTJGLClGNEYmQyY+Ri8pRjQtSS9jb252ZXJ0cmV2ZXJzZUdGJTYlRjBGJ0YpPkYwLUkxaW5jcmVtZW50cmV2ZXJzZUdGJTYkRjBGPz8oRissJkYsRj5GMyEiIkY+RiUyRitGLEMmPkYtJkYkNiNGKz5GLiZGJDYjLCZGK0Y+RjNGPj5GUy1JJG1vZEdGJTYkLCZGLUY+KiZGL0Y+Ri5GPkY+LCYpRjQpRjRGJ0Y+Rj5GPj5GVi1GZW42JCwmRi1GPkZobkZPRmluPkYsLCZGLEY+KUY0LCZGKEY+Rj5GPkY+RiVGJUYl

# # This procedure do a Fermat IFFT round modulo 2^(2^m)+1 on # array A having 2^n elements. The siccor size is 2^k and # the weights are get from the conversion of bits in interval II # from a counter starting with zero and incremented after # each butterfly sequence. # iffftround:=proc(A,n,m,k,II) local i,j,x,y,w,b; j:=2^k; b:=[0$i=1..m+1]; while j<2^n do w:=2^convertreverse(b,m,II); b:=incrementreverse(b,m+1); for i from j-2^k while i<j do x:=A[i]; y:=A[i+2^k]; A[i]:=x+y mod (2^(2^m)+1); A[i+2^k]:=(x-y)/w mod (2^(2^m)+1); od; j:=j+2^(k+1); od; end;

Zio2J0kiQUc2IkkibkdGJUkibUdGJUkia0dGJUkjSUlHRiU2KEkiaUdGJUkiakdGJUkieEdGJUkieUdGJUkid0dGJUkiYkdGJUYlRiVDJT5GLCkiIiNGKD5GMDcjLUkiJEclKnByb3RlY3RlZEc2JCIiIS9GKzsiIiIsJkYnRj5GPkY+PyhGJUY+Rj5GJTJGLClGNEYmQyY+Ri8pRjQtSS9jb252ZXJ0cmV2ZXJzZUdGJTYlRjBGJ0YpPkYwLUkxaW5jcmVtZW50cmV2ZXJzZUdGJTYkRjBGPz8oRissJkYsRj5GMyEiIkY+RiUyRitGLEMmPkYtJkYkNiNGKz5GLiZGJDYjLCZGK0Y+RjNGPj5GUy1JJG1vZEdGJTYkLCZGLUY+Ri5GPiwmKUY0KUY0RidGPkY+Rj4+RlYtRmVuNiQqJiwmRi1GPkYuRk9GPkYvRk9GaG4+RiwsJkYsRj4pRjQsJkYoRj5GPkY+Rj5GJUYlRiU=

# # This is the digit-by-digit multiplication modulo 2^(2^m)+1 # of arrays A and B having 2^n elements. # fdigmul:=proc(A,B,n,m) local i; for i from 0 to 2^n-1 do A[i]:=A[i]*B[i] mod (2^(2^m)+1) od; end;

Zio2JkkiQUc2IkkiQkdGJUkibkdGJUkibUdGJTYjSSJpR0YlRiVGJT8oRioiIiEiIiIsJikiIiNGJ0YtRi0hIiJJJXRydWVHJSpwcm90ZWN0ZWRHPiZGJEYpLUkkbW9kR0YlNiQqJkY1Ri0mRiZGKUYtLCYpRjApRjBGKEYtRi1GLUYlRiVGJQ==

# # This procedure divides by 2^k modulo 2^(2^m)+1 # all elements of array A having 2^n elements. # fshift:=proc(A,n,m,k) local i,r; r:=1/2^k mod (2^(2^m)+1); for i from 0 to 2^n-1 do A[i]:=A[i]*r mod (2^(2^m)+1) od; end;

Zio2JkkiQUc2IkkibkdGJUkibUdGJUkia0dGJTYkSSJpR0YlSSJyR0YlRiVGJUMkPkYrLUkkbW9kR0YlNiQqJCkiIiNGKCEiIiwmKUYzKUYzRiciIiJGOEY4PyhGKiIiIUY4LCYpRjNGJkY4RjhGNEkldHJ1ZUclKnByb3RlY3RlZEc+JkYkNiNGKi1GLzYkKiZGQEY4RitGOEY1RiVGJUYl

# # This procedure multiplies by sqrt(2) modulo 2^(2^m)+1 # all odd-indexed elements in the upper half of # array A having 2^n elements. # mulsqrt:=proc(A,n,m) local i,sqrttwo; sqrttwo:=2^(3*2^(m-2))-2^(2^(m-2)); for i from 2^(n-1)+1 to 2^n-1 by 2 do A[i]:=A[i]*sqrttwo mod (2^(2^m)+1) od; end;

Zio2JUkiQUc2IkkibkdGJUkibUdGJTYkSSJpR0YlSShzcXJ0dHdvR0YlRiVGJUMkPkYqLCYpIiIjLCQqJiIiJCIiIilGLywmRidGM0YvISIiRjNGM0YzKUYvRjRGNj8oRiksJilGLywmRiZGM0YzRjZGM0YzRjNGLywmKUYvRiZGM0YzRjZJJXRydWVHJSpwcm90ZWN0ZWRHPiZGJDYjRiktSSRtb2RHRiU2JComRkFGM0YqRjMsJilGLylGL0YnRjNGM0YzRiVGJUYl

# # This procedure divides by sqrt(2) modulo 2^(2^m)+1 # all odd-indexed elements in the upper half of # array A having 2^n elements. # divsqrt:=proc(A,n,m) local i,sqrthalf; sqrthalf:=1/(2^(3*2^(m-2))-2^(2^(m-2))) mod (2^(2^m)+1); for i from 2^(n-1)+1 to 2^n-1 by 2 do A[i]:=A[i]*sqrthalf mod (2^(2^m)+1) od; end;

Zio2JUkiQUc2IkkibkdGJUkibUdGJTYkSSJpR0YlSSlzcXJ0aGFsZkdGJUYlRiVDJD5GKi1JJG1vZEdGJTYkKiQsJikiIiMsJComIiIkIiIiKUYzLCZGJ0Y3RjMhIiJGN0Y3RjcpRjNGOEY6RjosJilGMylGM0YnRjdGN0Y3PyhGKSwmKUYzLCZGJkY3RjdGOkY3RjdGN0YzLCYpRjNGJkY3RjdGOkkldHJ1ZUclKnByb3RlY3RlZEc+JkYkNiNGKS1GLjYkKiZGSEY3RipGN0Y8RiVGJUYl

# # This procedure do Fermat FFT modulo 2^(2^m)+1 # for array A having 2^n elements. # ffft:=proc(A,n,m) local i; for i from 0 to n-2 do ffftround(A,n,m,n-1-i,1..i) od; if m>=n-1 then ffftround(A,n,m,0,1..n-1) else mulsqrt(A,n,m); ffftround(A,n,m,0,1..n-2) fi; end;

Zio2JUkiQUc2IkkibkdGJUkibUdGJTYjSSJpR0YlRiVGJUMkPyhGKSIiISIiIiwmRiZGLSIiIyEiIkkldHJ1ZUclKnByb3RlY3RlZEctSSpmZmZ0cm91bmRHRiU2J0YkRiZGJywoRiZGLUYtRjBGKUYwO0YtRilAJTEsJkYmRi1GLUYwRictRjQ2J0YkRiZGJ0YsO0YtRjpDJC1JKG11bHNxcnRHRiVGIy1GNDYnRiRGJkYnRiw7Ri1GLkYlRiVGJQ==

# # This procedure do Fermat IFFT modulo 2^(2^m)+1 # for array A having 2^n elements. # iffft:=proc(A,n,m) local i; if m>=n-1 then iffftround(A,n,m,0,1..n-1) else iffftround(A,n,m,0,1..n-2); divsqrt(A,n,m) fi; for i from n-2 to 0 by -1 do iffftround(A,n,m,n-1-i,1..i) od; end;

Zio2JUkiQUc2IkkibkdGJUkibUdGJTYjSSJpR0YlRiVGJUMkQCUxLCZGJiIiIkYuISIiRictSStpZmZmdHJvdW5kR0YlNidGJEYmRiciIiE7Ri5GLUMkLUYxNidGJEYmRidGMztGLiwmRiZGLiIiI0YvLUkoZGl2c3FydEdGJUYjPyhGKUY5Ri9GM0kldHJ1ZUclKnByb3RlY3RlZEctRjE2J0YkRiZGJywoRiZGLkYuRi9GKUYvO0YuRilGJUYlRiU=

# # This procedure do multiplication or squaring using Fermat FFT # and IFFT modulo 2^(2^m)+1 for array having 2^n elements. # 2^(m-1)-k bits is in one array element. # ffftmul:=proc(a,b,n,m) local i,c,A,B; if a=b then A:=Array(0..2^n-1,convert(a,base,2^(2^(m-1)-k))); ffft(A,n,m); fdigmul(A,A,n,m); iffft(A,n,m); fshift(A,n,m,n); c:=0; for i from 0 to 2^n-1 do c:=c+(2^(2^(m-1)-k))^i*A[i] od; else A:=Array(0..2^n-1,convert(a,base,2^(2^(m-1)-k))); ffft(A,n,m); B:=Array(0..2^n-1,convert(b,base,2^(2^(m-1)-k))); ffft(B,n,m); fdigmul(A,B,n,m); iffft(A,n,m); fshift(A,n,m,n); c:=0; for i from 0 to 2^n-1 do c:=c+(2^(2^(m-1)-k))^i*A[i] od; fi; c; end;

Zio2JkkiYUc2IkkiYkdGJUkibkdGJUkibUdGJTYmSSJpR0YlSSJjR0YlSSJBR0YlSSJCR0YlRiVGJUMkQCUvRiRGJkMpPkYsLUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiQ7IiIhLCYpIiIjRiciIiJGPCEiIi1JKGNvbnZlcnRHRjU2JUYkSSViYXNlR0YlKUY7LCYpRjssJkYoRjxGPEY9RjxJImtHRiVGPS1JJWZmZnRHRiU2JUYsRidGKC1JKGZkaWdtdWxHRiU2JkYsRixGJ0YoLUkmaWZmZnRHRiVGSS1JJ2ZzaGlmdEdGJTYmRixGJ0YoRic+RitGOD8oRipGOEY8RjlJJXRydWVHRjU+RissJkYrRjwqJilGQkYqRjwmRiw2I0YqRjxGPEMrRjJGRz5GLS1GNDYkRjctRj82JUYmRkFGQi1GSDYlRi1GJ0YoLUZLNiZGLEYtRidGKEZNRk9GUkZTRitGJUYlRiU=

n:=6; m:=4; k:=4; ffftmul(123456789,987654321,n,m,k);

IiIn

IiIl

IiIl

IjNwX2o3NmpLPjc=

123456789*987654321;

IjNwX2o3NmpLPjc=

LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEhRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0Yn