Arndt - Algorithms for Programmers (523138), страница 13
Текст из файла (страница 13)
THE HARTLEY TRANSFORM (HT)66Code 3.13 (DST via DCT) Pseudo code for the computation of the DST via DCT:procedure dst(x[],ldn)// real x[0..n-1] input,result{n := 2**nnh := n/2for k:=1 to n-1 step 2{x[k] := -x[k]}dct(x,ldn)for k:=0 to nh-1{swap(x[k],x[n-1-k])}}[FXT: dsth in dctdst/dsth.cc]Code 3.14 (IDST via IDCT) Pseudo code for the computation of the inverse sine transform (IDST)using the inverse cosine transform (IDCT):procedure idst(x[],ldn)// real x[0..n-1] input,result{n := 2**nnh := n/2for k:=0 to nh-1{swap(x[k],x[n-1-k])}idct(x,ldn)for k:=1 to n-1 step 2{x[k] := -x[k]}}[FXT: idsth in dctdst/dsth.cc]3.8Convolution via FHTThe convolution property of the Hartley transform isH [a ~ b] =´1³H [a] H [b] − H [a] H [b] + H [a] H [b] + H [a] H [b]2(3.24)or, written element-wise:H [a ~ b]k==¢1¡ck dk − ck dk + ck dk + ck dk2¢1¡ck (dk + dk ) + ck (dk − dk )where c = H [a],2d = H [b](3.25)Code 3.15 (cyclic convolution via FHT) Pseudo code for the cyclic convolution of two real valuedsequences x[] and y[], n must be even, result is found in y[]:procedure fht_cyclic_convolution(x[], y[], n)// real x[0..n-1] input, modifiedCHAPTER 3.
THE HARTLEY TRANSFORM (HT)67// real y[0..n-1] result{// transform data:fht(x[], n)fht(y[], n)// convolution in transformed domain:j := n-1for i:=1 to n/2-1{xi := x[i]xj := x[j]yp := y[i] + y[j]// =y[j] + y[i]ym := y[i] - y[j]// = -(y[j] - y[i])y[i] := (xi*yp + xj*ym)/2y[j] := (xj*yp - xi*ym)/2j := j-1}y[0] := y[0]*y[0]if n>1 then y[n/2] := y[n/2]*y[n/2]// transform back:fht(y[], n)// normalise:for i:=0 to n-1{y[i] := y[i] / n}}It is assumed that the procedure fht() does no normalization.fht/fhtcnvl.cc]Cf.
[FXT: fht convolution inEquation 3.25 (slightly optimized) for the auto convolution isH [a ~ a]k==1(ck (ck + ck ) + ck (ck − ck ))2¢1¡ 2ck ck +c − ck 2where c = H [a]2 k(3.26)Code 3.16 (cyclic auto convolution via FHT) Pseudo code for an auto convolution that uses a fastHartley transform, n must be even:procedure cyclic_self_convolution(x[], n)// real x[0..n-1] input, result{// transform data:fht(x[], n)// convolution in transformed domain:j := n-1for i:=1 to n/2-1{ci := x[i]cj := x[j]t1 := ci*cj// = cj*cit2 := 1/2*(ci*ci-cj*cj)// = -1/2*(cj*cj-ci*ci)x[i] := t1 + t2x[j] := t1 - t2j := j-1}x[0]:= x[0]*x[0]if n>1 then x[n/2] := x[n/2]*x[n/2]// transform back:fht(x[], n)// normalise:for i:=0 to n-1{CHAPTER 3.
THE HARTLEY TRANSFORM (HT)}}68x[i] := x[i] / nFor odd n replace the linefor i:=1 to n/2-1byfor i:=1 to (n-1)/2and omit the lineif n>1 thenx[n/2] := x[n/2]*x[n/2]in both procedures above. Cf. [FXT: fht auto convolution in fht/fhtcnvla.cc]3.9Negacyclic convolution via FHTCode 3.17 (negacyclic auto convolution via FHT) Code for the computation of the negacyclic(auto-) convolution:procedure negacyclic_self_convolution(x[], n)// real x[0..n-1] input, result{// preprocessing:hartley_shift_05(x, n)// transform data:fht(x, n)// convolution in transformed domain:j := n-1for i:=0 to n/2-1 // here i starts from zero{a := x[i]b := x[j]x[i] := a*b+(a*a-b*b)/2x[j] := a*b-(a*a-b*b)/2j := j-1}// transform back:fht(x, n)// postprocessing:hartley_shift_05(x, n)}(The code for hartley_shift_05() was given on page 56.)Cf. [FXT: fht negacyclic auto convolution in fht/fhtnegacnvla.cc]Code for the negacyclic convolution (without the ’self’):[FXT: fht negacyclic convolution in fht/fhtnegacnvl.cc]The underlying idea can be derived by closely looking at the convolution of real sequences by the radix-2FHT.The FHT-based negacyclic convolution turns out to be extremely useful for the computation of weightedtransforms, e.g.
in the MFA-based convolution for real input.Chapter 4Number theoretic transforms(NTTs)How to make a number theoretic transform out of your FFT:‘Replace exp(± 2 π i/n) by a primitive n-th root of unity, done.’We want to do FFTs in Z/mZ (the ring of integers modulo some integer m) instead of C, the (field ofthe) complex numbers.
These FFTs are called number theoretic transforms (NTTs), mod m FFTs or (ifm is a prime) prime modulus transforms.There is a restriction for the choice of m: For a length n NTT we need a primitive n-th root of unity.A number r is called an n-th root of unity if rn = 1. It is called a primitive n-th root if rk 6= 1 ∀ k < n.In C matters are simple: e± 2 π i/n is a primitive n-th root of unity for arbitrary n. e2 π i/21 is a 21-th rootof unity.
r = e2 π i/3 is also 21-th root of unity but not a primitive root, because r3 = 1. A primitive n-throot of 1 in Z/mZ is also called an element of order n. The ‘cyclic’ property of the elements r of ordern lies in the heart of all FFT algorithms: rn+k = rk .In Z/mZ things are not that simple since primitive roots of unity do not exist for arbitrary n, they existfor some maximal order R only. Roots of unity of an order different from R are available only for thedivisors di of R: rR/di is a di -th root of unity because (rR/di )di = rR = 1.Therefore n must divide R, the first condition for NTTs:n\R⇐⇒ ∃√n1(4.1)The operations needed in FFTs are addition, subtraction and multiplication.
Division is not needed,except for division by n for the final normalization after transform and back-transform. Division by n ismultiplication by the inverse of n. Hence n must be invertible in Z/mZ: n must be coprime1 to m, thesecond condition for NTTs:n⊥m⇐⇒ ∃ n−1 in Z/mZ(4.2)Cf. [1], [3], [19] or [2] and books on number theory.4.1Prime modulus: Z/pZ = FpIf the modulus is a prime p then Z/pZ is the field Fp : All elements except 0 have inverses and ‘division ispossible’ in Z/pZ.
Thereby the second condition is trivially fulfilled for all FFT lengths n < p: a primep is coprime to all integers n < p.1ncoprime to m ⇐⇒ gcd(n, m) = 169CHAPTER 4. NUMBER THEORETIC TRANSFORMS (NTTS)70Roots of unity are available for the maximal order R = p−1 and its divisors: Therefore the first conditionon n for a length-n mod p FFT being possible is that n divides p − 1. This restricts the choice for p toprimes of the form p = v n + 1: For length-n = 2k FFTs one will use primes like p = 3 · 5 · 227 + 1 (31bits), p = 13 · 228 + 1 (32 bits), p = 3 · 29 · 256 + 1 (63 bits) or p = 27 · 259 + 1 (64 bits)2 .
The elementsof maximal order in Z/pZ are called primitive elements, generators or primitive roots modulo p. If r is agenerator, then every element in Fp different from 0 is equal to some power re (1 ≤ e < p) of r and itsorder is R/e. To test whether r is a primitive n-th root of unity in Fp one does not need to check rk 6= 1for all k < n. It suffices to do the check for exponents k that are prime factors of n.
This is because theorder of any element divides the maximal order. To find a primitive root in Fp proceed as indicated bythe following pseudo code:Code 4.1 (Primitive root modulo p) Return a primitive root in Fpfunction primroot(p){if p==2 then return 1f[] := distinct_prime_factors(p-1)for r:=2 to p-1{x := TRUEforeach q in f[]{if r**((p-1)/q)==1 then x:=FALSE}if x==TRUE then return r}error("no primitive root found") // p cannot be prime !}An element of order n is returned by this function:Code 4.2 (Find element of order n) Return an element of order n in Fp :function element_of_order(n,p){R := p-1 // maxorderif (R/n)*n != R then error("order n must divide maxorder p-1")r := primroot(p)x := r**(R/n)return x}4.2Composite modulus: Z/mZIn what follows we will need the function ϕ(), the so-called ‘totient’ function.
ϕ(m) counts the numberof integers prime to and less than m. For m = p prime ϕ(p) = p − 1. For m composite ϕ(m) is alwaysless than m − 1. For m = pk a prime powerϕ(pk )= pk − pk−1(4.3)e.g. ϕ(2k ) = 2k−1 . ϕ(1) = 1. For coprime p1 , p2 (p1 , p2 not necessarily primes) ϕ(p1 p2 ) = ϕ(p1 ) ϕ(p2 ),ϕ() is a so-called multiplicative function.For the computation of ϕ(m) for m a prime power one can use this simple piece of codeCode 4.3 (Compute phi(m) for m a prime power) Return ϕ(px )2 Primesof that form are not ‘exceptional’, cf.
Lipson [3]CHAPTER 4. NUMBER THEORETIC TRANSFORMS (NTTS)71function phi_pp(p,x){if x==1 then return p - 1elsereturn p**x - p**(x-1)}Pseudo code to compute ϕ(m) for general m:Code 4.4 (Compute phi(m)) Return ϕ(m)function phi(m){{n, p[], x[]} := factorization(m)ph := 1for i:=0 to n-1{ph := ph * phi_pp(p[i],x[i])}}// m==product(i=0..n-1,p[i]**x[i])Further we need the notion of Z/mZ∗ , the ring of units in Z/mZ. Z/mZ∗ contains all invertible elements(‘units’) of Z/mZ, i.e. those which are coprime to m. Evidently the total number of units is given byϕ(m):|Z/mZ∗ | =ϕ(m)(4.4)kIf m factorizes as m = 2k0 · pk11 · .
. . · pq q then|Z/mZ∗ |= ϕ(2k0 ) · ϕ(pk11 ) · . . . · ϕ(pkq q )(4.5)It turns out that the maximal order R of an element can be equal to or less than |Z/mZ∗ |, the ringZ/mZ∗ is then called cyclic or noncyclic, respectively. For m a power of an odd prime p the maximalorder R in Z/mZ∗ (and also in Z/mZ) isR(pk ) =ϕ(pk )while for m a power of two a tiny irregularity enters: 12R(2k ) = k−22for k = 1for k = 2for k ≥ 3(4.6)(4.7)i.e. for powers of two greater than 4 the maximal order deviates from ϕ(2k ) = 2k−1 by a factor of 2.
Forkthe general modulus m = 2k0 · pk11 · . . . · pq q the maximal order isR(m) =lcm(R(2k0 ), R(pk11 ), . . . , R(pkq q ))where lcm() denotes the least common multiple.Pseudo code to compute R(m):Code 4.5 (Maximal order modulo m) Return R(m), the maximal order in Z/mZfunction maxorder(m){{n, p[], k[]} := factorization(m) // m==product(i=0..n-1,p[i]**k[i])R := 1for i:=0 to n-1{t := phi_pp(p[i],k[i])if p[i]==2 AND k[i]>=3 then t := t / 2R := lcm(R,t)}return R}(4.8)CHAPTER 4. NUMBER THEORETIC TRANSFORMS (NTTS)72Now we can see for which m the ring Z/mZ∗ will be cyclic:Z/mZ∗cyclic form = 2, 4, pk , 2 · pk(4.9)where p is an odd prime.If m contains two different odd primes pa , pb then R(m) =lcm(. . .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
