Arndt - Algorithms for Programmers (523138), страница 14
Текст из файла (страница 14)
, ϕ(pa ), ϕ(pb ), . . . ) is at least by a factor of two smaller than ϕ(m) = . . . · ϕ(pa ) · ϕ(pb ) · . . .because both ϕ(pa ) and ϕ(pb ) are even, so Z/mZ∗ cannot be cyclic in that case. The same argumentholds for m = 2k0 · pk if k0 > 1. For m = 2k Z/mZ∗ is cyclic only for k = 1 and k = 2 because of theabove mentioned irregularity of R(2k ).Pseudo code (following [19]) for a function that returns the order of some element x in Z/mZ:Code 4.6 (Order of an element in Z/mZ) Return the order of an element x in Z/mZfunction order(x,m){if gcd(x,m)!=1 then return 0 // x not a unith := phi(m) // number of elements of ring of unitse := h{n, p[], k[]} := factorization(h) // h==product(i=0..n-1,p[i]**k[i])for i:=0 to n-1{f := p[i]**k[i]e := e / fg1 := x**e mod mwhile g1!=1{g1 := g1**p[i] mod me := e * p[i]p[i] := p[i] - 1}}return e}Pseudo code for a function that returns some element x in Z/mZ of maximal order:Code 4.7 (Element of maximal order in Z/mZ) Return an element that has maximal order inZ/mZfunction maxorder_element(m){R := maxorder(m)for x:=1 to m-1{if order(x,m)==R then}// never reached}return xFor prime m the function returns a primitive root.
It is a good idea to have a table of small primes stored(which will also be useful in the factorization routine) and restrict the search to small primes and only ifthe modulus is greater than the largest prime of the table proceed with a loop as above:Code 4.8 (Element of maximal order in Z/mZ) Return an element that has maximal order inZ/mZ, use a precomputed table of primesfunction maxorder_element(m, pt[], np)// pt[0..np-1] = 2,3,5,7,11,13,17,...{if m==2 then return 1R := maxorder(m)for i:=0 to np-1{if order(pt[i],m)==R then return xCHAPTER 4.
NUMBER THEORETIC TRANSFORMS (NTTS)}73}// hardly ever reachedfor x:=pt[np-1] to m-1 step 2{if order(x,m)==R then return x}// never reached[FXT: maxorder element mod in mod/maxorder.cc]There is no problem if the prime table contains primes ≥ m: The first loop will finish before order() iscalled with an element ≥ m, because before that can happen, the element of maximal order is found.4.3Pseudocode for NTTsTo implement NTTs (‘mod-m FFTs’) one basically must implement modulo-arithmetics and replacee± 2 π i/n by an n-th root of unity in Z/mZ in the code.
[FXT: class mod in mod/mod.h]For the back-transform one uses the (mod m) inverse r̄ of r (an element of order n) that was used forthe forward transform. To check whether r̄ exists one tests whether gcd(r, m) = 1. To compute theinverse modulo m one can use the relation r̄ = rϕ(p)−1 (mod m). Alternatively one may use the extendedEuclidean algorithm, which for two integers a and b finds d = gcd(a, b) and u, v so that a u + b v = d.Feeding a = r, b = m into the algorithm gives u as the inverse: r u + m v ≡ r u ≡ 1 (mod m).While the notion of the Fourier transform as a ‘decomposition into frequencies’ seems to be meaninglessfor NTTs the algorithms are denoted with ‘decimation in time/frequency’ in analogy to those in thecomplex domain.The nice feature of NTTs is that there is no loss of precision in the transform (as there is always with thecomplex FFTs).
Using the analogue of trigonometric recursion (in its most naive form) is mandatory, asthe computation of roots of unity is expensive.4.3.1Radix 2 DIT NTTCode 4.9 (radix 2 DIT NTT) Pseudo code for the radix 2 decimation in time NTT (to be called withldn=log2(n)):procedure mod_fft_dit2(f[], ldn, is)// mod_type f[0..2**ldn-1]{n := 2**ldnrn := element_of_order(n) // (mod_type)if is<0 then rn := rn**(-1)revbin_permute(f[], n)for ldm:=1 to ldn{m := 2**ldmmh := m/2dw := rn**(2**(ldn-ldm))// (mod_type)w := 1// (mod_type)for j:=0 to mh-1{for r:=0 to n-1 step m{t1 := r+jt2 := t1+mhv := f[t2]*w // (mod_type)u := f[t1]// (mod_type)f[t1] := u+vCHAPTER 4.
NUMBER THEORETIC TRANSFORMS (NTTS)}}}f[t2] := u-v}w := w*dwLike in 1.4.2 it is a good idea to extract the ldm==1 stage of the outermost loop:Replacefor ldm:=1 to ldn{byfor r:=0 to n-1 step 2{{f[r], f[r+1]} := {f[r]+f[r+1], f[r]-f[r+1]}}for ldm:=2 to ldn{[FXT: ntt dit2 in mod/nttdit2.cc]4.3.2Radix 2 DIF NTTCode 4.10 (radix 2 DIF NTT) Pseudo code for the radix 2 decimation in frequency NTT:procedure mod_fft_dif2(f[], ldn, is)// mod_type f[0..2**ldn-1]{n := 2**ldndw := element_of_order(n) // (mod_type)if is<0 then dw := rn**(-1)for ldm:=ldn to 1 step -1{m := 2**ldmmh := m/2w := 1 // (mod_type)for j:=0 to mh-1{for r:=0 to n-1 step m{t1 := r+jt2 := t1+mhv := f[t2] // (mod_type)u := f[t1] // (mod_type)f[t1] := u+vf[t2] := (u-v)*w}w := w*dw}dw := dw*dw}revbin_permute(f[], n)}As in section 1.4.3 extract the ldm==1 stage of the outermost loop:Replace the lineforldm:=ldn to 1 step -1forldm:=ldn to 2 step -1by74CHAPTER 4.
NUMBER THEORETIC TRANSFORMS (NTTS)and insertfor r:=0 to n-1 step 2{{f[r], f[r+1]} := {f[r]+f[r+1], f[r]-f[r+1]}}before the call of revbin_permute(f[],n).[FXT: ntt dif2 in mod/nttdif2.cc]4.3.3Radix 4 NTTsC++ code for a radix-4 decimation in time NTT:static const ulong LX = 2;void ntt_dit4(mod *f, ulong ldn, int is)//// radix 4 decimation in time NTT//{const ulong n = (1UL<<ldn);revbin_permute(f, n);// n is not a power of 4, need a radix 2 step:if ( ldn & 1 ){for (ulong i=0; i<n; i+=2) sumdiff(f[i], f[i+1]);}const mod imag = mod::root2pow( is>0 ? 2 : -2 );ulong ldm = LX + (ldn&1);for ( ; ldm<=ldn ; ldm+=LX){const ulong m = (1UL<<ldm);const ulong m4 = (m>>LX);const mod dw = mod::root2pow( is>0 ? ldm : -ldm );mod w = (mod::one);mod w2 = w;mod w3 = w;for (ulong j=0; j<m4; j++){for (ulong r=0, i0=j+r; r<n; r+=m, i0+=m){const ulong i1 = i0 + m4;const ulong i2 = i1 + m4;const ulong i3 = i2 + m4;mod a0 = f[i0];mod a2 = f[i1] * w2;mod a1 = f[i2] * w;mod a3 = f[i3] * w3;mod t02 = a0 + a2;mod t13 = a1 + a3;f[i0] = t02 + t13;f[i2] = t02 - t13;t02 = a0 - a2;t13 = a1 - a3;t13 *= imag;f[i1] = t02 + t13;f[i3] = t02 - t13;}w *= dw;w2 = w * w;w3 = w * w2;}}}75CHAPTER 4.
NUMBER THEORETIC TRANSFORMS (NTTS)The function mod::root2pow(x) returns a primitive root of order 2x .mod/nttdit4.cc].76This is [FXT: ntt dit4 inThe radix-4 DIF variant is [FXT: ntt dif4 in mod/nttdif4.cc].For applications of the NTT see the survey article [35] which also has a good bibliography.4.4Convolution with NTTsThe NTTs are natural candidates for (exact) integer convolutions, as used e.g. in (high precision) multiplications. One must keep in mind that ‘everything is mod p’, the largest value that can be representedis p − 1. As an example consider the multiplication of n-digit radix R numbers3 .
The largest possiblevalue in the convolution is the ‘central’ one, it can be as large as M = n (R − 1)2 (which will occur ifboth numbers consist of ‘nines’ only4 ).One has to choose p > M to get rid of this problem. If p does not fit into a single machine wordthis may slow down the computation unacceptably. The way out is to choose p as the product of severaldistinct primes that are all just below machine word size and use the Chinese Remainder Theorem (CRT)afterwards.If using length-n FFTs for convolution there must be an inverse element for n. This imposes the conditiongcd(n, modulus) = 1, i.e. the modulus must be prime to n.
Usually5 modulus must be an odd number.Integer convolution: Split input mod m1, m2, do 2 FFT convolutions, combine with CRT.4.5The Chinese Remainder Theorem (CRT)The Chinese remainder theorem (CRT):Let m1 , m2 , . . . , mf be pairwise relatively6 prime (i.e. gcd(mi , mj ) = 1, ∀i 6= j)If x ≡ xi (mod mi ) i = 1, 2, . . . , f then x is unique modulo the product m1 · m2 · . . . · mf .For only two moduli m1 , m2 compute x as follows7 :Code 4.11 (CRT for two moduli) pseudo code to find unique x (mod m1 m2 ) with x ≡ x1 (mod m1 )x ≡ x2 (mod m2 ):function crt2(x1,m1,x2,m2){c := m1**(-1) mod m2// inverse of m1 modulo m2s := ((x2-x1)*c) mod m2return x1 + s*m1}For repeated CRT calculations with the same moduli one will use precomputed c.For more more than two moduli use the above algorithm repeatedly.Code 4.12 (CRT) Code to perform the CRT for several moduli:function crt(x[],m[],f){x1 := x[0]m1 := m[0]3 Multiplicationis a convolution of the digits followed by the ‘carry’ operations.radix R ‘nine’ is R − 1, nine in radix 10 is 9.5 for length-2k FFTs6 note that it is not assumed that any of the m is primei7 cf.
[3]4ACHAPTER 4. NUMBER THEORETIC TRANSFORMS (NTTS)77i := 1do{x2 := x[i]m2 := m[i]x1 := crt2(x1,m1,x2,m2)m1 := m1 * m2i := i + 1}while i<freturn x1}To see why these functions really work we have to formulate a more general CRT procedure that specializesto the functions above.DefineTiY:=mk(4.10)k!=iandηi:= Ti−1mod mi(4.11)then forXi:=one hasxi ηi Ti½Ximod mj=xi0(4.12)for j = ielse(4.13)and soXXk= ximod mi(4.14)kFor the special case of two moduli m1 , m2 one hasT1T2η1η2====m2m1m−12m−11mod m1(4.15)(4.16)(4.17)mod m2(4.18)= 1(4.19)which are related by8η1 m2 + η2 m1XXk=x1 η1 T1 + x2 η2 T2(4.20)==x1 η1 m2 + x2 η2 m1x1 (1 − η2 m1 ) + x2 η2 m1(4.21)(4.22)=x1 + (x2 − x1 ) (m−11(4.23)kmod m2 ) m1as given in the code.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
