| 13 | } |
| 14 | |
| 15 | void FFT_Data::init(const Ring& Rg,const Zp_Data& PrD) |
| 16 | { |
| 17 | R=Rg; |
| 18 | prData=PrD; |
| 19 | root.resize(2); |
| 20 | |
| 21 | // Find which case we are in |
| 22 | int hwt=Hwt(Rg.m()); |
| 23 | bigint prm1=PrD.pr-1; |
| 24 | twop=1<<(numBits(Rg.m())+1); |
| 25 | if (hwt==1) |
| 26 | { twop=0; } |
| 27 | else if ((prm1%twop)!=0) |
| 28 | { twop=-twop; } |
| 29 | |
| 30 | if (twop==0) |
| 31 | { int nb=numBits(Rg.m()); |
| 32 | nb=1<<(nb-1); |
| 33 | if (nb==Rg.m()) |
| 34 | { //cout << Rg.m() << " " << PrD.pr << endl; |
| 35 | root[0]=Find_Primitive_Root_2power(Rg.m(),PrD); |
| 36 | Inv(root[1],root[0],PrD); |
| 37 | to_modp(iphi,Rg.phi_m(),PrD); |
| 38 | Inv(iphi,iphi,PrD); |
| 39 | compute_roots(Rg.m()); |
| 40 | } |
| 41 | } |
| 42 | else |
| 43 | { bigint pr=PrD.pr; |
| 44 | if ((pr-1)%(2*Rg.m())!=0) |
| 45 | { throw invalid_params(); } |
| 46 | root[0]=Find_Primitive_Root_2m(Rg.m(),Rg.Phi(),PrD); |
| 47 | Inv(root[1],root[0],PrD); |
| 48 | compute_roots(2 * Rg.m()); |
| 49 | |
| 50 | int ptwop=twop; if (twop<0) { ptwop=-twop; } |
| 51 | |
| 52 | powers.resize(2,vector<modp>(Rg.m())); |
| 53 | powers_i.resize(2,vector<modp>(Rg.m())); |
| 54 | b.resize(2,vector<modp>(ptwop)); |
| 55 | |
| 56 | modp rInv,bi; |
| 57 | bigint ee=ptwop; ee=ee*Rg.m(); |
| 58 | for (int r=0; r<2; r++) |
| 59 | { assignOne(powers[r][0],PrD); |
| 60 | if (r==0) |
| 61 | { to_modp(powers_i[0][0],ptwop,PrD); } |
| 62 | else |
| 63 | { to_modp(powers_i[1][0],ee,PrD); } |
| 64 | Inv(powers_i[r][0],powers_i[r][0],PrD); |
| 65 | |
| 66 | Inv(rInv,root[r],PrD); |
| 67 | assignOne(b[r][Rg.m()-1],PrD); |
| 68 | for (long i=1; i<Rg.m(); i++) |
| 69 | { long iSqr=(i*i)%(2*Rg.m()); |
| 70 | Power(powers[r][i],root[r],iSqr,PrD); |
| 71 | Mul(powers_i[r][i],powers[r][i],powers_i[r][0],PrD); |
| 72 | Power(bi,rInv,iSqr,PrD); |
no test coverage detected