Polynomial Multiplication with NTT

[Milkman]

Hello everyone, I'm new to this forum, and I've registered to seek advice on this topic. What I'm attempting to do is multiply two very large numbers together, represented as polynomials, by the use of the number theory transform (ntt). I've done much research on this topic and I've had some success, but I have run into problems as well.

So far, I've coded a (hopefully) working ntt algorithm in c++ which is based off of some of the code found in the apfloat package (http://www.apfloat.org/apfloat/) . With this algorithm, I can successfully ntt a polynomial (represented as an array of numbers), and then inverse ntt the polynomial back to its original state.

Seeing as I could do that, I proceded to the next step in my project. I ntt'd a polynomial, squared each value, and then inverse ntt'd the new polynomial. This should yield a polynomial that represents the square of the original number, as this method works with a fast Fourier transform. Here's where my trouble lies. For some reason, sometimes I do get the correct polynomial as a result, and sometimes I don't. There seems to be no real pattern in the results that I get. If I use a different base for the polynomial, sometimes certain numbers will work when they didn't in another base, but no one base always yields correct results. I would like to know if the error lies in my ntt algorithm, or elsewhere in my method. I was attempting this with a modulus of 5, a length of 4, and a primitive root of 2 for the ntt, but i get similar results for a modulus of 17, a length of 16, and a primitive root of 3 as well.

As I said before, I know that with the fast Fourier transform, you can fft a polynomial, square each value, and inverse fft the new polynomial to get the square of the original number. My whole project relies on the premise that this also works for the number theory transform, which I believe it should, shouldn't it? Below is the code for my number theory transform. Any help with figuring out what is wrong would be much appreciated.
Thank you,
Mike

(tested and compiles with dev c++)

#include <iostream>
using namespace std;



// evaluates exponents, with a modulus 
long power(long base,long power,long mod)
{
       long final=1;
       for (int i=0; i<power; i++)
       {
           final = ((final%mod) * (base%mod)) % mod;
       }
       return final;
} 



// prints an array
void print(long x[],int N)
{
     for (int i=0; i<N; i++)
     {
         cout<<x[i]<<endl;
     }
     cout<<endl;
}



// number theory transform algorithm
void ntt(long data[],long N,long mod,long pri)
{
    long k;
    long j;
    long w;
    long wj=1;
    long wk=1;
    long temp[N];

    w = power(pri,((mod - 1)/N) % mod,mod) % mod;

    for (j=0; j<N; j++)
    {
        temp[j]=0;
        wk=1;
        for (k=0; k<N; k++)
        {
            temp[j] = ((temp[j]%mod) + ((wk%mod)*(data[k]%mod))%mod) % mod;
            wk = ((wk%mod)*(wj%mod)) % mod;
        }
        wj = ((wj%mod)*(w%mod)) % mod;
    }

    for (k=0; k<N; k++)
    {
        data[k] = temp[k] % mod;
    }

}



// inverse number theory transform algorithm
void intt(long data[],long N,long mod,long pri)
{
     long k;
     long j;
     long w;
     long wk=1;
     long wj=1;
     long temp[N];
     long rN;
     
     w = power(pri,(mod-1-((mod-1)/N)) % mod,mod) % mod;
     rN = power(N,(mod-2)%mod,mod) % mod;
     
     for (j=0; j<N; j++)
     {
         temp[j]=0;
         wk=1;
         for (k=0; k<N; k++)
         {
            temp[j] = ((temp[j]%mod) + ((wk%mod)*(data[k]%mod))%mod) % mod;
            wk = ((wk%mod)*(wj%mod)) % mod;
         }
         wj = ((wj%mod)*(w%mod)) % mod;
     }
     
     for (k=0; k<N; k++)
     {
         data[k]=((temp[k]%mod)*(rN%mod)) % mod;
     }
     
}



// main
int main()
{
    long N=16;
    long mod=17;
    long pri=3;
    long data[N];

    data[0]=3;
    data[1]=1;
    data[2]=0;
    data[3]=0;
    for (int i=4; i<N; i++)
    {
        data[i]=0;
    }
    
    print(data,N);
    ntt(data,N,mod,pri);
    print(data,N);
    for (int i=0; i<N; i++)
    {
        data[i]=data[i]*data[i];
    }
    print(data,N);
    intt(data,N,mod,pri);
    print(data,N);
    
    cout<<endl<<endl<<endl;
    system("PAUSE");
    return 0; 
}

 

[Hurkyl]

Can you give examples of times you get the wrong answer?

 

[Milkman]

Yes, here is a list of results that i get in base 16 and 17, for an ntt with a modulus of 17, length of 16, and primitive root of 3.

               base16  base17
num   num^2    result  result
16     256      256    
17     289      289     289
18     324      324     324
19     361      361     361
20     400      400     400
21     441      424     441
22     484      450     467
23     529      495     495
24     576      525     542
25     625      285     574
26     676      319     319
27     729      338     355
28     784      376     376
29     841      416     416
30     900      441     458
31     961      468     485
32     1024     1024    514
33     1089     1089    545
34     1156     1156    1156

so...
in base 16, the results deviate from num^2 when num=21,22,23,24,25,26,27,28,29,30,31
in base 17, the results deviate from num^2 when num=22,23,24,25,26,27,28,29,30,31,32,33

I guess there is a pattern with the incorrect results, but i don't know what to make of it or how to fix the problem.

 

[Hurkyl]

What is "num" and "result" in that table?


I was looking for a specific example of a polynomial (and the corresponding N, mod, and pri), and the polynomial your code thought was its square.


By the way, did you notice that all of the differences are divisible by 17? e.g. for num = 26: 676 - 319 = 357 = 17 * 21.


One of the things I was suspecting is that your code gives the right answer modulo 17... which is the best it can hope to do because all the arithmetic is done modulo 17.

I had wanted to see an actual polynomial so I could test that hypothesis, but that results column suggests it too.

 

[shmoe]

When you multiply two polynomials using a number theory transform mod 17, you get the product of these two polynomials mod 17 (please correct me if I'm wrong, I can't say I've ever looked at a ntt before today-keep this in mind when reading this post i could be off). For example, when squaring 21 using base 16, you used 21=1*16+5. This corresponds to squaring the polynomial x+5. The answer you want is x^2+10x+25, but mod 17 the square of this polynomial is x^2+10x+8, which explains your off by 17.

Mod 17 is a different story, 21=1*17+4. The square of x+4 is x^2+8x+16, which is also it's square mod 17. Essentially you get problems when the square of your coefficients is 17 or greater. (other problems will occur when your transform isn't long enough to get the highest powers of the polynomial). By looking at the coefficients in your base 16 and 17 representations of "num" in your table, it should be clear why the rows your mod 17 ntt fail or succeed.

There's a few ways around this. First you can represent your numbers by a small enough base with respect to your transform to ensure the square of the coefficients is never too large. Try squaring those numbers again using the same ntt but working base 5 (e.g. 21=4*5+1, you can use your mod 17 ntt to square 4*x+1 without problems).

Second, you could use multiple transforms to different bases and the chinese remainder theorem to find the coeffcients of your squared polynomials exactly. for example, if you are representing your numbers in base 16 you wanted to know (x+5)^2 to find 21^2. You know it's x^2+10*x+8 mod 17 from your mod 17 ntt. Do another ntt using the primitive root 2 mod 5, you'll get (x+5)^2=x^2 mod 5 for example (okie, that is obvious without the transform). By the Chinese remainder theorem you'd then know the coefficients of (x+5)^2 mod 17*5=85, which is sufficient here.edit-had "mod" where I should have had "base" and a 16 that should have been a 17

 

[Milkman]

Thank you shmoe, I understand exactly what's going wrong now. In order to be absolutely sure that a polynomial can be squared using this method without any problems, I think that this inequality must evaluate true:

(d+1)(m-1)^2 < p

where d is the degree of the polynomial, m is the base used, and p is the prime modulus used in the ntt. This is because m-1 is the largest possible diget value for the given base, and d+1 times the square of the largest possible diget is the size of the largest possible coefficient when the polynomial is squared. The largest possible coefficient must be smaller than p for it to work.

According to this, even a base of 5 would not be good enough to safely square all first degree polynomials when p=17 because (1+1)(5-1)^2 = 32 > 17. The largest base that you can use to safely square polynomials when p=17 is 3! (not factorial ) And that's only talking about first degree polynomials, i.e. only two digets in base 3.

Thanks for the help guys, I'll let you know if I run into any more snags along my journey.

 

[Hurkyl]

Oh, whoops! I missed that your goal is to multiply integers, not polynomials.


Hrm, this is interesting: I seem to remember that when squaring an integer n, you usually want to work with a modulus bigger than n^2... but off hand, I don't see anything wrong with your suggestion on a lower bound for the modulus!

 

[Milkman]

Multiplication works perfectly now, but I'm wondering about division. If i wish to divide two sequences of numbers with an ntt, would I just ntt those two sequences, and multiply the first sequence's terms by the modular inverse of the second sequence's terms? Is division even possible to do with an ntt (i know it's possible with an ftt)? Any help would be appreciated, thanks.

 

[Hurkyl]

I can tell you what I know about polynomials. When using one of these evaluate/interpolate algorithms, you can only do exact division: remainders must be zero.

Since you're really working with integers, but converting to polynomials, you have a very difficult problem: not only does the integer division have to be exact, but so does the polynomial division!

Of course, I would try coding it up to see if it works; sometimes the heavens align and good things happen.