Karatsuba multiplication algorithm with modular reduction

Question

I am trying to estimate the time it takes to multiply two large bit numbers modulo some large prime. My computational power is limited to adding, multiplying, and storing 32-bit numbers.

Due to this I wanted to use the Karatsuba algorithm since the multiplication is reduced to single-digit operations. However, the algorithm does not incorporate modular reduction.

My other thought would be to use Montgomery reduction but I'm not sure how this would perform given my computational limitations.

Which algorithm would you suggest I use?

Answer 1

This is not an answer but this would be unreadable inside comment

• it is just illustration how Karatsuba code looks like

This is mine implementation of Karatsuba multiplication for arbitrary floating point numbers

//---------------------------------------------------------------------------
void arbnum::_mul_karatsuba(DWORD *z,DWORD *x,DWORD *y,int n)
    {
    // recursion for karatsuba
    // z[2n]=x[n]*y[n]; n=2^m
    int i;    for (i=0;i<n;i++) if (x[i]) { i=-1; break; }      // x==0 ?
    if (i< 0) for (i=0;i<n;i++) if (y[i]) { i=-1; break; }      // y==0 ?
    if (i>=0){for (i=0;i<n+n;i++) z[i]=0; return; }             // 0.? = 0
    if (n==1) { alu.mul(z[0],z[1],x[0],y[0]); return; }
    if (n< 1) return;
    int n2=n>>1;
    _mul_karatsuba(z+n,x+n2,y+n2,n2);                           // z0 = x0.y0
    _mul_karatsuba(z  ,x   ,y   ,n2);                           // z2 = x1.y1
    DWORD *q=new DWORD[n<<1],*q0,*q1,*qq; BYTE cx,cy;
    if (q==NULL) { _error(_arbnum_error_NotEnoughMemory); return; }
    #define _add { alu.add(qq[i],q0[i],q1[i]); for (i--;i>=0;i--) alu.adc(qq[i],q0[i],q1[i]); } // qq=q0+q1 ...[i..0]
    #define _sub { alu.sub(qq[i],q0[i],q1[i]); for (i--;i>=0;i--) alu.sbc(qq[i],q0[i],q1[i]); } // qq=q0-q1 ...[i..0]
    qq=q;    q0=x+n2; q1=x;   i=n2-1; _add; cx=alu.cy;          // =x0+x1
    qq=q+n2; q0=y+n2; q1=y;   i=n2-1; _add; cy=alu.cy;          // =y0+y1
    _mul_karatsuba(q+n,q+n2,q,n2);                              // =(x0+x1)(y0+y1) mod ((2^N)-1)
    if (cx) { qq=q+n; q0=qq; q1=q+n2; i=n2-1; _add; cx=alu.cy; }// +=cx*(y0+y1)<<n2
    if (cy) { qq=q+n; q0=qq; q1=q   ; i=n2-1; _add; cy=alu.cy; }// +=cy*(x0+x1)<<n2
    qq=q+n;  q0=qq;   q1=z+n; i=n-1;  _sub;                     // -=z0
    qq=q+n;  q0=qq;   q1=z;   i=n-1;  _sub;                     // -=z2
    qq=z+n2; q0=qq;   q1=q+n; i=n-1;  _add;                     // z1=(x0+x1)(y0+y1)-z0-z2
    for (i=n2-1;i>=0;i--) if (alu.cy) alu.inc(z[i]); else break;                // adc(1)
    alu.cy=cx|cy; for (i=n2-1;i>=0;i--) if (alu.cy) alu.inc(z[i]); else break;  // adc(1) if +=cx*(y0+y1)<<n2 or +=cy*(x0+x1)<<n2 overflowed
    delete[] q;
    #undef _add
    #undef _sub
    }
//---------------------------------------------------------------------------
void arbnum::mul_karatsuba(const arbnum &x,const arbnum &y)
    {
    // O(3*(N)^log2(3)) ~ O(3*(N^1.585))
    // Karatsuba multiplication
    int s=x.sig*y.sig;
    arbnum a,b; a=x; b=y; a.sig=+1; b.sig=+1;
    int i,n;
    for (n=1;(n<a.siz)||(n<b.siz);n<<=1);
    a._realloc(n);
    b._realloc(n);
    _alloc(n+n); for (i=0;i<siz;i++) dat[i]=0;
    _mul_karatsuba(dat,a.dat,b.dat,n);
    bits=siz<<5;
    sig=s;
    exp=a.exp+b.exp+((siz-a.siz-b.siz)<<5)+1;
//  _normalize();
    }
//---------------------------------------------------------------------------

• no modulo is present
• as you can see you will need +,- operations with variable bit width both are encoded as #define so no need for another function calls ...
• it uses 32bit ALU in x86 C++ as a block to build variable bit width arithmetics
• result is in this object
• _mul_karatsuba is just internal recursive function
• mul_karatsuba is the main api ...

Number representation:

// dat is MSDW first ... LSDW last
DWORD *dat; int siz,exp,sig,bits;

• dat[siz] is the mantisa LSDW means least significant DWORD
• exp is exponent of MSB of dat[0]

• first nonzero bit is present in mantisa !!!

  // |-----|---------------------------|---------------|------|
  // | sig | MSB      mantisa      LSB |   exponent    | bits |
  // |-----|---------------------------|---------------|------|
  // | +1  | 0.(0      ...          0) | 2^0           |   0  | +zero
  // | -1  | 0.(0      ...          0) | 2^0           |   0  | -zero
  // |-----|---------------------------|---------------|------|
  // | +1  | 1.(dat[0] ... dat[siz-1]) | 2^exp         |   n  | +number
  // | -1  | 1.(dat[0] ... dat[siz-1]) | 2^exp         |   n  | -number
  // |-----|---------------------------|---------------|------|
  // | +1  | 1.0                       | 2^+0x7FFFFFFE |   1  | +infinity
  // | -1  | 1.0                       | 2^+0x7FFFFFFE |   1  | -infinity
  // |-----|---------------------------|---------------|------|