NTT相关问题一则，有兴趣的参与

ssikkiss

这是我用来找NTT参数的代码,由于我的MATHEMATICA刚学不久,有什么不对的地方大家指正!
这里我没有找到具有较好的n次跟的NTT,水平有限啊!

mulinv[a_,b_]:=Module[{p,f},p=
  ExtendedGCD[a,b];p=p[[2]];f=Mod[p[[1]]*a,a*b];Return[f]]
nextprime[m_,step_]:=Module[{p},
    p=m+step;
    While[PrimeQ[p]\[Equal]False,p+=step];
    Return[p];
    ]
nroot[n_,p_]:=Module[{k,a,b,n2=n/2,p2=p-1,i},
    k=2;
    i=0;
    Do[If[PowerMod[k,n,p]\[Equal]1 && PowerMod[k,n2,p]\[Equal]p2,Break[]];
    i++;k*=2,{16}];
    If[i<16,Return[k]];
    k=2;
    i=0;
    Do[If[PowerMod[k+1,n,p]\[Equal]1 && PowerMod[k+1,n2,p]\[Equal]p2,Break[]];
      i++;k*=2,{16}];
    If[i<16,Return[k+1]];
    k=4;
    i=0;
    Do[If[PowerMod[k-1,n,p]\[Equal]1 && PowerMod[k-1,n2,p]\[Equal]p2,
    Break[]];i++;k*=2,{16}];
    If[i<16,Return[k-1]];
    k=1;
    While[(PowerMod[k,n2,p]!=p2)  ||  (PowerMod[k,n,p]!=1),k++];
    Return[k];
    ]
findntt[e_,w_]:=Module[{n=2^e,m,a,ia,invn,i},
    m=n*w*w+1-n;
    m=nextprime[m,n];
    a=nroot[n,m];
    ia=PowerMod[a,-1,m];
    invn=m-(m-1)/n;
    Return[{n,m,a,ia,invn}];
    ]
nttlist[maxe_,w_]:=Module[{},Return[Table[findntt[i,w],{i,1,maxe,1}]]]
Put["ntt paramtarer,for base=10","ntt1.txt" ] ;
PutAppend[nttlist[32,10],"ntt1.txt"];

ssikkiss

这个也是我用来找NTT参数的,不过它同时找2组参数,并且增加对在NTT里使用CRT(中国余数定理)的支持

mulinv[a_,b_]:=Module[{p,f},p=
  ExtendedGCD[a,b];p=p[[2]];f=Mod[p[[1]]*a,a*b];Return[f]]
findntt[e_,w_,fn_]:=Module[{n=2^e,i,j,z,
  m,m1,m2,a1,a2,b1,b2,invm1,invm2,k1,k2,tmp,t1,t2},t1=n/2;
    t2=n-1;
    z=w*Sqrt[n];
    z=Floor[z];
    m=z;
    z=Mod[z,n];
    m-=z;
    m+=1;
    While[PrimeQ[m]≠True,m+=n];
    m1=m;
    invm1=(m*n-m+1)/n;
    a1=findroot[m1,n];
    b1=PowerMod[a1,t2,m];
    m+=n;
    While[PrimeQ[m]≠True,m+=n];
    m2=m;
    invm2=(m*n-m+1)/n;
    a2=findroot[m2,n];
    b2=PowerMod[a2,t2,m];
    m=m1*m2;
    k1=mulinv[m1,m2];
    k2=mulinv[m2,m1];
    PutAppend["___________________________________",fn];
    PutAppend[{"  N = ",n},fn];
    PutAppend["___________________________________",fn];
    PutAppend[{"M1=",m1," a1=",a1,"  1/a1=",b1,"  invm1=",invm1},fn];
    PutAppend[{"M2=",m2,"  a2=",a2,"  1/a2=",b2,"  invm2=",invm2},fn];
    PutAppend[{"M1*M2=",m,"  k1=",k1,"  k2=",k2},fn];
    PutAppend[{"M1M2/NWW=",N[m/(n*w*w),4]},fn];
    PutAppend["___________________________________",fn];]
findroot[p_,n_]:=Module[{i,j,l,q,b,v,f,d},
    j=2;
    i=p-1;
    i=i/n;
    q=0;
    l=PowerMod[j,i,p];
    While[( l<2) || (ispow[l]\[Equal]2),j++;q++;l=PowerMod[j,i,p];If[q>
    200000,Break[]]];
    If[q<200000,Return[l]];
    q=0;
    l=2;
    b=n/2;
    v=p-1;
    f=PowerMod[tt,n,p];
    d=PowerMod[tt,b,p];
    While[(!(f== 1 )) || (!(d == v))  ,l++;q++;f=PowerMod[l,n,p];d=PowerMod[l,
    b,p];If[q>1024,Break[]];];
    If[(f\[Equal]1)&&(d\[Equal]v),Return[l]];
    q=0;
    l=PowerMod[j,i,p];
    While[( l<2) || (ispow[l]\[Equal]2),j++;
    q++;l=PowerMod[j,i,p];If[q>200000,Break[]]];
    If[q<200000,Return[l]];
    l=PowerMod[j,i,p];
    While[( l<2) ,j++;q++;l=PowerMod[j,i,p];];
    Return[l];
    ]

ispow[n_]:=Module[{xx=
  n},If[xx\[Equal]1,Return[2]];While[Mod[xx,2]\[Equal]0,xx/=
      2];If[xx\[Equal]1,Return[1],Return[2];];]
startfindntt[wordSize_,filename_]:=Module[{},
    Put["Search For NTT Parameter",filename];
    PutAppend[{"It's used for ",wordSize," word width,"},filename];
    Do[findntt[i,wordSize,filename],{i,1,32}];]
startfindntt[10,"ntt_1.txt"];
startfindntt[100,"ntt_2.txt"];
startfindntt[1000,"ntt_3.txt"];
startfindntt[10000,"ntt_4.txt"];
startfindntt[100000,"ntt_5.txt"];
startfindntt[1000000,"ntt_6.txt"];
startfindntt[10000000,"ntt_7.txt"];
startfindntt[100000000,"ntt_8.txt"];
startfindntt[1000000000,"ntt_9.txt"];

无心人

程序有点难懂

呵呵

能写出程序运行结果么？

另外有其他形式的表示

稍后会开贴讨论
不过日期可能无限向后推迟

ssikkiss

传上个附件,里面是我用MATHEMATICA的计算结果
现在对这个附件里的文件的解释:
"Search For NTT Parameter"   每个文件都以这句打头
{"It's used for ", 10, " word width,"}  表示这个文件里的参数,例如,这里是 用于 10^1  进制
"___________________________________"  
{"  N = ", 2}                             表示紧接着的下面的参数,这里是 可用于长n=2的变换
"___________________________________"
{"M1=", 17, " a1=", 16, "  1/a1=", 16, "  invm1=", 9}  表示 其中一个模M=17,根为16,根的倒数也是16,1/2是9
{"M2=", 19, "  a2=", 18, "  1/a2=", 18, "  invm2=", 10} 第2个模M=19,........................
{"M1*M2=", 323, "  k1=", 171, "  k2=", 153}  M1*M2是这两个模的积,k1和k2是用于CRT的2个参数
{"M1M2/NWW=", 1.615\`4.000000000000001} 这个的意义,如果这个值小于1,那么此组参数将无法进行NTT变换
"___________________________________"

ssikkiss

再传各附件,这是用前面第一个程序得到的,里面仅计算一组NTT参数,没有对CRT的支持
对此附件里的内容的解释:
{2, 211, 210, 210, 106},  表示一个长度n=2的变换,其模M=211,根是210,根的倒数也是210,1/2是106
{4, 401, 20, 381, 301},
{8, 809, 44, 570, 708}, 表示一个长n=8的变换,模M=809,根是44,根的倒数是570,1/8是708
..
..
..

无心人

多谢分享
不过，怎么说呢

有实际用途的参数不多的

无心人

你做下模是2^b + 1，可能b不是2^t形式，但是s*2^t形式，s为奇数
模不是素数
变换长度是2的方幂，根也最好是2的方幂的NTT的搜索

这个实用性强些

ssikkiss

你说的是伪费马变换吧!我上面的参数全是s*2^t形式的,但我注意过s是不是奇数.
目前为止我还没找到有效的好的NTT,那种根是2或2得幂的,正探索中....
顺带问个小问题,我做NTT变换时,是用刚才的程序先生成一堆(其实就32组)NTT参数,作ntt变换的时候根据数据的长短调用相应的参数
不知楼主是怎么做的?不可能对不同长度的数据应用相同德一组NTT参数吧?

无心人

s必须是奇数，否则也可转换成这种形式

看NTT自身的限制

在我的介绍NTT的帖子里有限制条件

如果用途就是卷积，则就是了

如果是大数运算，则还需要修改下

无心人

可以对不同长度的数据应用相同的NTT参数

关键是补零

在我帖子里开头有证明