diff --git a/NTT_fromR.nb b/NTT_fromR.nb new file mode 100644 index 0000000..0189b7c --- /dev/null +++ b/NTT_fromR.nb @@ -0,0 +1,89 @@ +(* Content-type: application/vnd.wolfram.mathematica *) + +(*** Wolfram Notebook File ***) +(* http://www.wolfram.com/nb *) + +(* CreatedBy='Mathematica 11.2' *) + +(*CacheID: 234*) +(* Internal cache information: +NotebookFileLineBreakTest +NotebookFileLineBreakTest +NotebookDataPosition[ 158, 7] +NotebookDataLength[ 2708, 81] +NotebookOptionsPosition[ 2399, 67] +NotebookOutlinePosition[ 2755, 83] +CellTagsIndexPosition[ 2712, 80] +WindowFrame->Normal*) + +(* Beginning of Notebook Content *) +Notebook[{ +Cell[BoxData[{ + RowBox[{"SetDirectory", "[", + "\"\\"", + "]"}], "\[IndentingNewLine]", + RowBox[{"<<", "NumberTheoryFunctions.m"}], "\[IndentingNewLine]", + RowBox[{"<<", " ", "numberth.m"}], "\[IndentingNewLine]", + RowBox[{"Print", "[", "1", "]"}], "\[IndentingNewLine]", + RowBox[{"init", "[", "]"}], "\n", + RowBox[{ + RowBox[{"U", " ", "=", " ", + RowBox[{"ReadList", "[", "\n", " ", + RowBox[{ + "\"\\"", + ",", " ", "\n", " ", "Number", ",", " ", + RowBox[{"RecordLists", " ", "->", " ", "True"}]}], "]"}]}], + ";"}], "\n", + RowBox[{ + RowBox[{"A", " ", "=", " ", + RowBox[{"U", "[", + RowBox[{"[", "1", "]"}], "]"}]}], ";"}], "\n", + RowBox[{ + RowBox[{"B", " ", "=", " ", + RowBox[{"U", "[", + RowBox[{"[", "2", "]"}], "]"}]}], ";"}], "\n", + RowBox[{ + RowBox[{"t1", " ", "=", " ", + RowBox[{ + RowBox[{"fnt", "[", "A", "]"}], "*", + RowBox[{"fnt", "[", "B", "]"}]}]}], ";"}], "\n", + RowBox[{ + RowBox[{"c1", " ", "=", " ", + RowBox[{"ifnt", "[", "t1", "]"}]}], ";"}], "\n", "c1", "\n", + RowBox[{"str", " ", "=", " ", + RowBox[{ + "OpenWrite", "[", "\n", " ", + "\"\\"", + "]"}]}], "\n", + RowBox[{"Write", "[", + RowBox[{"str", ",", " ", "c1"}], "]"}], "\n", + RowBox[{"Close", "[", "str", "]"}]}], "Input", + CellChangeTimes->{{3.735253411577821*^9, 3.735253456196847*^9}, { + 3.735253572224071*^9, 3.73525357352667*^9}, 3.7352543189796953`*^9, { + 3.735267261221972*^9, 3.735267263104089*^9}, {3.735268133707547*^9, + 3.735268166321697*^9}, 3.735268365769558*^9, 3.735268476400447*^9, + 3.735268511237224*^9}, + CellLabel->"In[32]:=",ExpressionUUID->"f446b4b9-745a-4d9f-8637-212fe444c6d8"] +}, +WindowSize->{808, 713}, +WindowMargins->{{Automatic, 297}, {-46, Automatic}}, +FrontEndVersion->"11.3 for Mac OS X x86 (32-bit, 64-bit Kernel) (March 5, \ +2018)", +StyleDefinitions->"Default.nb" +] +(* End of Notebook Content *) + +(* Internal cache information *) +(*CellTagsOutline +CellTagsIndex->{} +*) +(*CellTagsIndex +CellTagsIndex->{} +*) +(*NotebookFileOutline +Notebook[{ +Cell[558, 20, 1837, 45, 367, "Input",ExpressionUUID->"f446b4b9-745a-4d9f-8637-212fe444c6d8"] +} +] +*) + diff --git a/NTT_fromR.wls b/NTT_fromR.wls new file mode 100644 index 0000000..d9fc77f --- /dev/null +++ b/NTT_fromR.wls @@ -0,0 +1,17 @@ +#!/usr/bin/env wolframscript +(* ::Package:: *) + +SetDirectory["/Users/maryam/Google Drive/Research/RF_improvment"] +< True]; +A = U[[1]]; +B = U[[2]]; +t1 = fnt[A]*fnt[B]; +c1 = ifnt[t1]; +c1; +Export["/Users/maryam/Google Drive/Research/RF_improvment/outNTT.txt",c1] diff --git a/RF_FFT.R b/RF_FFT.R new file mode 100644 index 0000000..a44ae3a --- /dev/null +++ b/RF_FFT.R @@ -0,0 +1,126 @@ +library("ape") +library("phangorn") +library("readtext") + +#Beta function +Beta=function(m){ + if(m<0){return(1)} + return(dfactorial(2*m+1)) +} + +# the function for computing the number of internal edges + +internaledges <- function(tree,ntip){ + intedges=array(0,c(1,ntip-1)) + edges=tree$edge + for (i in (2*ntip-1):(ntip+1)) { + children=which(edges[,1]==i) + child1=edges[children[1],2] + child2=edges[children[2],2] + if((child1 <= ntip)&(child2 <= ntip)){intedges[i-ntip]=0} + else if((child1<= ntip) & (child2 > ntip)) {intedges[i-ntip]=intedges[child2-ntip]+1} + else if((child2<= ntip )& (child1 > ntip)){intedges[i-ntip]=intedges[child1-ntip]+1} + else {intedges[i-ntip]=intedges[child2-ntip]+intedges[child1-ntip]+2} + } + return(intedges) +} + + +internalchildren <- function(tree,v,ntip){ + edges=tree$edge + children=which(edges[,1]==v) + child1=edges[children[1],2] + child2=edges[children[2],2] + if((child1 > ntip) & (child2 > ntip)){result=c(2,child1,child2)} + else if((child1 > ntip) & (child2 <= ntip)){result=c(1,child1)} + else if((child2 > ntip) & (child1 <= ntip)){result=c(1,child2)} + else {result=0} + return(result) +} + + +RF_Convolve=function(tree,n){ + ntip=n-1 + N=tree$Nnode + R=rep(list(matrix(0,(ntip-1),(ntip-1))),N) + edges=internaledges(tree,ntip) + B=c() + for (k in 0:(n-2)) { + B[k+1]=Beta(k) + } + for (v in N:1) { + intchild=internalchildren(tree,v+ntip,ntip) + intedges=edges[v] + if(intchild[1]==0){ + R[[v]][1,1]=1 + } + else if(intchild[1]==1){ + Rchild=R[[intchild[2]-ntip]] + R[[v]][1,intedges+1]=1 + R[[v]][2:(ntip-1),1]=rowSums(t(t(Rchild[1:(ntip-2),])*B[1:(ntip-1)])) + R[[v]][2:(ntip-1),2:(ntip-1)]=Rchild[2:(ntip-1),1:((ntip-2))] + } + else { + Rchild1=R[[intchild[2]-ntip]] + Rchild2=R[[intchild[3]-ntip]] + R[[v]][1,intedges+1]=1 + R[[v]][3,1]=sum(t(t(Rchild1[1,])*B[1:(ntip-1)]))*sum(t(t(Rchild2[1,])*B[1:(ntip-1)])) + for (s in 4:(ntip-1)) { + R[[v]][s,1]=sum(rowSums(t(t(Rchild1[1:(s-2),])*B[1:(ntip-1)]))*rowSums(t(t(Rchild2[(s-2):1,])*B[1:(ntip-1)]))) + } + sum1=matrix(0,(ntip-2),(ntip-2)) + sum1[1,1:(ntip-2)]=sum(t(t(Rchild1[1,])*B[1:(ntip-1)]))*Rchild2[1,1:(ntip-2)] + for (s in 3:(ntip-1)) { + temp=colSums(rowSums(t(t(Rchild1[1:(s-1),])*B[1:(ntip-1)]))*Rchild2[(s-1):1,1:(ntip-2)]) + sum1[s-1,1:(ntip-2)]=temp + } + sum2=matrix(0,(ntip-2),(ntip-2)) + sum2[1,1:(ntip-2)]=sum(t(t(Rchild2[1,])*B[1:(ntip-1)]))*Rchild1[1,1:(ntip-2)] + for (s in 3:(ntip-1)) { + temp=colSums(rowSums(t(t(Rchild2[1:(s-1),])*B[1:(ntip-1)]))*Rchild1[(s-1):1,1:(ntip-2)]) + sum2[s-1,1:(ntip-2)]=temp + } + + R1=Rchild1[1:(ntip-1),1:(ntip-3)] + R2=Rchild2[1:(ntip-1),1:(ntip-3)] + R1aug=cbind(R1,matrix(0,nrow(R1),ncol(R1))) + R2aug=cbind(R2,matrix(0,nrow(R2),ncol(R2))) + U=round(convolve(t(R1aug), rev(t(R2aug)), type="open"),0) + sum3=t(matrix(c(U,0),2*ncol(R1),2*nrow(R1))[1:ncol(R1),1:nrow(R1)]) + sum3=cbind(array(0, dim=c(nrow(R1)-1,1)),sum3[2:nrow(R1),]) + R[[v]][2:(ntip-1),2:(ntip-1)]=sum1+sum2+sum3 + } + } + return(R) +} + + + +#========================================== +RsT=function(R,n,s){ + B=c() + for (k in 0:(n-2)) { + B[k+1]=Beta(k) + } + rst =sum(t(t(R[[1]][s+1,1:(n-2-s)])*B[1:(n-2-s)])) + return(rst) +} + +#Compute the value of q_m(T) +qmT=function(R,n,m){ + qmt=0 + for (s in m:(n-3)) { + rst=RsT(R,n,s) + qmt=qmt+(factorial(s)/(factorial(m)*factorial(s-m)))*rst*(-1)^(s-m) + } + return(qmt) +} + +polynomial=function(tree,n){ + R=RF_Convolve(tree,n) + for (i in seq(0,2*(n-3),2)) { + print(qmT(R,n,n-3-(i/2))) + } +} + + diff --git a/RF_NTT.R b/RF_NTT.R new file mode 100644 index 0000000..4c705dd --- /dev/null +++ b/RF_NTT.R @@ -0,0 +1,142 @@ +library("ape") +library("phangorn") +library("readtext") + +#Beta function +Beta=function(m){ + if(m<0){return(1)} + return(dfactorial(2*m+1)) +} + +# the function for computing the number of internal edges +internaledges <- function(tree,ntip){ + intedges=array(0,c(1,ntip-1)) + edges=tree$edge + for (i in (2*ntip-1):(ntip+1)) { + children=which(edges[,1]==i) + child1=edges[children[1],2] + child2=edges[children[2],2] + if(child1 <= ntip&child2 <= ntip){intedges[i-ntip]=0} + else if(child1<= ntip & child2 > ntip){intedges[i-ntip]=intedges[child2-ntip]+1} + else if(child2<= ntip & child1 > ntip){intedges[i-ntip]=intedges[child1-ntip]+1} + else {intedges[i-ntip]=intedges[child2-ntip]+intedges[child1-ntip]+2} + } + return(intedges) +} + +# the function for computing the number of internal children +internalchildren <- function(tree,v,ntip){ + edges=tree$edge + children=which(edges[,1]==v) + child1=edges[children[1],2] + child2=edges[children[2],2] + if(child1 > ntip & child2 > ntip){result=c(2,child1,child2)} + else if(child1 > ntip & child2 <= ntip){result=c(1,child1)} + else if(child2 > ntip & child1 <= ntip){result=c(1,child2)} + else {result=0} + return(result) +} + +RF_Convolve=function(tree,n){ + ntip=n-1 + N=tree$Nnode + R=rep(list(matrix(0,(ntip-1),(ntip-1))),N) + edges=internaledges(tree,ntip) + B=c() + for (k in 0:(n-2)) { + B[k+1]=Beta(k) + } + for (v in N:1) { + intchild=internalchildren(tree,v+ntip,ntip) + intedges=edges[v] + if(intchild[1]==0){ + R[[v]][1,1]=1 + } + else if(intchild[1]==1){ + Rchild=R[[intchild[2]-ntip]] + R[[v]][1,intedges+1]=1 + R[[v]][2:(ntip-1),1]=rowSums(t(t(Rchild[1:(ntip-2),])*B[1:(ntip-1)])) + R[[v]][2:(ntip-1),2:(ntip-1)]=Rchild[2:(ntip-1),1:((ntip-2))] + } + else { + Rchild1=R[[intchild[2]-ntip]] + Rchild2=R[[intchild[3]-ntip]] + R[[v]][1,intedges+1]=1 + R[[v]][3,1]=sum(t(t(Rchild1[1,])*B[1:(ntip-1)]))*sum(t(t(Rchild2[1,])*B[1:(ntip-1)])) + for (s in 4:(ntip-1)) { + R[[v]][s,1]=sum(rowSums(t(t(Rchild1[1:(s-2),])*B[1:(ntip-1)]))*rowSums(t(t(Rchild2[(s-2):1,])*B[1:(ntip-1)]))) + } + sum1=matrix(0,(ntip-2),(ntip-2)) + sum1[1,1:(ntip-2)]=sum(t(t(Rchild1[1,])*B[1:(ntip-1)]))*Rchild2[1,1:(ntip-2)] + for (s in 3:(ntip-1)) { + temp=colSums(rowSums(t(t(Rchild1[1:(s-1),])*B[1:(ntip-1)]))*Rchild2[(s-1):1,1:(ntip-2)]) + sum1[s-1,1:(ntip-2)]=temp + } + sum2=matrix(0,(ntip-2),(ntip-2)) + sum2[1,1:(ntip-2)]=sum(t(t(Rchild2[1,])*B[1:(ntip-1)]))*Rchild1[1,1:(ntip-2)] + for (s in 3:(ntip-1)) { + temp=colSums(rowSums(t(t(Rchild2[1:(s-1),])*B[1:(ntip-1)]))*Rchild1[(s-1):1,1:(ntip-2)]) + sum2[s-1,1:(ntip-2)]=temp + } + + R1=Rchild1[1:(ntip-1),1:(ntip-3)] + #R1aug=numeric(nrow(R1)*(2*ncol(R1)-1)) + R1aug=numeric(2048) + t=1 + for(i in 1:ncol(R1)){ + R1aug[t:(t+nrow(R1)-1)]=R1[,i] + t=t+3*nrow(R1) + } + R2=Rchild2[1:(ntip-1),1:(ntip-3)] + #R2aug=numeric(nrow(R2)*(2*ncol(R2)-1)) + R2aug=numeric(2048) + t=1 + for(i in 1:ncol(R2)){ + R2aug[t:(t+nrow(R2)-1)]=R2[,i] + t=t+3*nrow(R2) + } + + write(R1aug,"testNTT.txt",ncolumns=2048,append = TRUE) + write(R2aug,"testNTT.txt",ncolumns=2048,append = TRUE) + #read the output of mathematica + system("/Applications/Mathematica.app/Contents/MacOS/MathKernel -script /Users/maryam/Desktop/NTT_fromR.wls") + U=as.matrix(read.csv("/Users/maryam/Google Drive/Research/RF_improvment/outNTT.txt",header = FALSE, sep="'", quote="")) + Matc=ceiling(length(R1aug)/(3*nrow(R1))) + sum3=matrix(c(U,numeric(Matc*3*nrow(R1)-length(R1aug))),nrow=3*nrow(R1))[1:nrow(R1),1:ncol(R1)] + sum3=cbind(array(0, dim=c(nrow(R1)-1,1)),sum3[2:nrow(R1),]) + R[[v]][2:(ntip-1),2:(ntip-1)]=sum1+sum2+sum3 + file.remove("testNTT.txt") + } + } + + return(R) +} + + + +#========================================== +RsT=function(R,n,s){ + B=c() + for (k in 0:(n-2)) { + B[k+1]=Beta(k) + } + rst =sum(t(t(R[[1]][s+1,1:(n-2-s)])*B[1:(n-2-s)])) + return(rst) +} + +#Compute the value of q_m(T) +qmT=function(R,n,m){ + qmt=0 + for (s in m:(n-3)) { + rst=RsT(R,n,s) + qmt=qmt+(factorial(s)/(factorial(m)*factorial(s-m)))*rst*(-1)^(s-m) + } + return(qmt) +} + +polynomial=function(tree,n){ + R=RF_Convolve(tree,n) + for (i in seq(0,2*(n-3),2)) { + print(qmT(R,n,n-3-(i/2))) + } +} diff --git a/RF_OLD.R b/RF_OLD.R new file mode 100644 index 0000000..c0ecc26 --- /dev/null +++ b/RF_OLD.R @@ -0,0 +1,129 @@ +library("ape") +library("phangorn") +library("readtext") + + +#Beta function +Beta=function(m){ + if(m<0){return(1)} + return(dfactorial(2*m+1)) +} + +# the function for computing the number of internal edges + +internaledges <- function(tree,ntip){ + intedges=array(0,c(1,ntip-1)) + edges=tree$edge + for (i in (2*ntip-1):(ntip+1)) { + children=which(edges[,1]==i) + child1=edges[children[1],2] + child2=edges[children[2],2] + if(child1 <= ntip&child2 <= ntip){intedges[i-ntip]=0} + else if(child1<= ntip & child2 > ntip){intedges[i-ntip]=intedges[child2-ntip]+1} + else if(child2<= ntip & child1 > ntip){intedges[i-ntip]=intedges[child1-ntip]+1} + else {intedges[i-ntip]=intedges[child2-ntip]+intedges[child1-ntip]+2} + } + return(intedges) +} + + +internalchildren <- function(tree,v,ntip){ + edges=tree$edge + children=which(edges[,1]==v) + child1=edges[children[1],2] + child2=edges[children[2],2] + if(child1 > ntip & child2 > ntip){result=c(2,child1,child2)} + else if(child1 > ntip & child2 <= ntip){result=c(1,child1)} + else if(child2 > ntip & child1 <= ntip){result=c(1,child2)} + else {result=0} + return(result) +} + +#============================================ +RF_old=function(tree,n){ + ntip=n-1 + N=tree$Nnode + R=rep(list(matrix(0,(ntip-1),(ntip-1))),N) + edges=internaledges(tree,ntip) + B=c() + for (k in 0:(ntip-2)) { + B[k+1]=Beta(k) + } + for (v in N:1) { + intchild=internalchildren(tree,v+ntip,ntip) + intedges=edges[v] + if(intchild[1]==0){ + R[[v]][1,1]=1 + } + else if(intchild[1]==1){ + Rchild=R[[intchild[2]-ntip]] + R[[v]][1,intedges+1]=Beta(intedges) + R[[v]][2:(ntip-1),1]=rowSums(Rchild[1:(ntip-2),]) + R[[v]][2:(ntip-1),2:(ntip-1)]=t(t(Rchild[2:(ntip-1),1:((ntip-2))])*seq(3,(2*ntip-3),2)) + } + else { + Rchild1=R[[intchild[2]-ntip]] + Rchild2=R[[intchild[3]-ntip]] + R[[v]][1,intedges+1]=Beta(intedges) + R[[v]][3,1]=sum(Rchild1[1,])*sum(Rchild2[1,]) + for (s in 4:(ntip-1)) { + R[[v]][s,1]=sum(rowSums(Rchild1[1:(s-2),])*rowSums(Rchild2[(s-2):1,])) + } + sum1=matrix(0,(ntip-2),(ntip-2)) + sum1[1,1:(ntip-2)]=t(t(rowSums(t(Rchild1[1,]))*Rchild2[1,1:(ntip-2)])*seq(3,(2*ntip-3),2)) + for (s in 3:(ntip-1)) { + temp=colSums(t(t(rowSums(Rchild1[1:(s-1),])*Rchild2[(s-1):1,1:(ntip-2)])*seq(3,(2*ntip-3),2))) + sum1[s-1,1:(ntip-2)]=temp + } + sum2=matrix(0,(ntip-2),(ntip-2)) + sum2[1,1:(ntip-2)]=t(t(rowSums(t(Rchild2[1,]))*Rchild1[1,1:(ntip-2)])*seq(3,(2*ntip-3),2)) + for (s in 3:(ntip-1)) { + temp=colSums(t(t(rowSums(Rchild2[1:(s-1),])*Rchild1[(s-1):1,1:(ntip-2)])*seq(3,(2*ntip-3),2))) + sum2[s-1,1:(ntip-2)]=temp + } + sum3=matrix(0,(ntip-2),(ntip-2)) + for (s in 1:(ntip-2)){ + for (k in 2:(ntip-2)) { + total3=0 + for (s1 in 0:(s)) { + for (k1 in 0:(k-2)) { + # if((s1+1) > 0 & (k1+1) > 0 &(s-s1+1) > 0 & (k-2-k1+1) > 0){ + total3=total3+Rchild1[s1+1,k1+1]*Rchild2[s-s1+1,k-2-k1+1]*Beta(k)/(Beta(k1)*Beta(k-2-k1)) + # } + } + sum3[s,k]=total3 + } + + } + } + + R[[v]][2:(ntip-1),2:(ntip-1)]=sum1+sum2+sum3 + } + } + return(R) +} + + + +#========================================== +RsT=function(R,n,s){ + rst =sum(R[[1]][s+1,1:(n-2-s)]) + return(rst) +} + +#Compute the value of q_m(T) +qmT=function(R,n,m){ + qmt=0 + for (s in m:(n-3)) { + rst=RsT(R,n,s) + qmt=qmt+(factorial(s)/(factorial(m)*factorial(s-m)))*rst*(-1)^(s-m) + } + return(qmt) +} + +polynomial=function(tree,n){ + R=RF_old(tree,n) + for (i in seq(0,2*(n-3),2)) { + print(qmT(R,n,n-3-(i/2))) + } +} diff --git a/numberth.m b/numberth.m new file mode 100644 index 0000000..40ae1e2 --- /dev/null +++ b/numberth.m @@ -0,0 +1,101 @@ +(* ::Package:: *) + +(* numberth.m *) +(* Last updated August 8th 1995 *) + +n = 2048; (* transform length *) +l = 2^256; (* computer's word length *) + +inv[n_Integer] := PowerMod[n, -1, m] (* inverse function *) + +init[] := (* init variables *) +Block[{}, + k = Floor[l/n]-1; + While[PrimeQ[k*n+1] == False, k--]; + m= k*n+1; (* the prime modulo *) + r=PrimitiveRoot[m]; (* primitive root of m *) + w=PowerMod[r, k, m]; (* root of order n, is r^((m-1)/n) *) + w1=inv[w]; (* inverse of w *) + ]; + +initf[] := (* init transform matrixes *) +Block[{}, + f=Table[PowerMod[w, i j, m], {i, 0, n-1}, {j, 0, n-1}]; (* transform *) + f1=Table[PowerMod[w1, i j, m], {i, 0, n-1}, {j, 0, n-1}]; (* inverse *) + f1 = Mod[inv[n] * f1, m]; + ]; + +carry[l_List] := (* convert list to number *) +Block[{s = 0}, + For[k=Length[l], k > 0, k--, s=s*10+l[[k]]]; + Return[s]]; + +n2l[a_Integer] := (* convert number to list *) + Join[Reverse[IntegerDigits[a]], Table[0, {n-1-Floor[Log[10, a]]}]]; + +mulf[a_Integer, b_Integer] := (* multiply with slow algorithm *) +carry[Mod[f1.(Mod[f.n2l[a], m]*Mod[f.n2l[b], m]), m]]; + +mul[a_Integer, b_Integer] := (* multiply with fast algorithm *) +carry[ifnt[Mod[fnt[n2l[a]] * fnt[n2l[b]], m]]]; + +permute[a_Integer] := (* bit reversal *) +Block[{t = 2 * n, o = 2 * a, p = 0}, + While[(t /= 2) > 1, p += p + Mod[(o = Floor[o / 2]), 2]]; + Return[p]]; + +scramble[l_List] := (* bit reverse list order *) +Table[l[[permute[k]+1]], {k, 0, n-1}]; + +fnt[l_List] := (* fast number theoretic transform *) +Block[{f = l, i1, i2, i3, i4, loop, loop1, loop2, a, b, y, z}, + y = w1; + + i1 = n / 2; + i2 = 1; + + For[loop = i1, loop > 0, loop = Floor[loop / 2], + Block[{}, + i3 = 0; + i4 = i1; + For[loop1 = 0, loop1 < i2, loop1++, + Block[{}, + z = 1; + For[loop2 = i3, loop2 < i4, loop2++, + Block[{}, + a = f[[loop2 + 1]]; + b = f[[loop2 + i1+ 1]]; + f[[loop2 + 1]] = Mod[a + b, m]; + f[[loop2 + i1 + 1]] = Mod[z * (a - b), m]; + z = Mod[z * y, m]]]; + i3 += 2 * i1; + i4 += 2 * i1]]; + y = Mod[y * y, m]; + i1 = Floor[i1 / 2]; + i2 *= 2]]; + Return[f]] + +ifnt[l_List] := (* inverse fnt *) +Block[{f = l, i, j, k, istep, mmax, wt, wr, wtemp}, + mmax = 1; + While[n > mmax, + Block[{}, + istep = 2 * mmax; + wt = PowerMod[w, n / istep, m]; + wr = 1; + For[k = 0, k < mmax, k++, + Block[{}, + For[i = k, i < n, i += istep, + Block[{}, + j = i + mmax; + wtemp = Mod[wr * f[[j + 1]], m]; + f[[j + 1]] = Mod[f[[i + 1]] - wtemp, m]; + f[[i + 1]] = Mod[f[[i + 1]] + wtemp, m]]]; + wr = Mod[wr * wt, m]]]; + mmax = istep]]; + Return[Mod[inv[n] * f, m]]] + +fntraw[l_List] := scramble[fnt[l]]; (* for debug only *) +ifntraw[l_List] := Mod[ifnt[scramble[l]] n, m]; (* for debug only *) + +(* end numberth.m *)