* fast reduction for generalized mersenne numbers (GM)

#include "os.h"
#include <mp.h>
#include "dat.h"

/*
* fast reduction for generalized mersenne numbers (GM)
* using a series of additions and subtractions.
*/

enum {
MAXDIG = 1024/Dbits,
};

typedef struct GMfield GMfield;
struct GMfield
{
Mfield;

mpint m2[1];

int nadd;
int nsub;
int indx[256];
};

static int
gmreduce(Mfield *m, mpint *a, mpint *r)
{
GMfield *g = (GMfield*)m;
mpdigit d0, t[MAXDIG];
int i, j, d, *x;

if(mpmagcmp(a, g->m2) >= 0)
return -1;

if(a != r)
mpassign(a, r);

d = g->top;
mpbits(r, (d+1)*Dbits*2);
memmove(t+d, r->p+d, d*Dbytes);

r->sign = 1;
r->top = d;
r->p[d] = 0;

if(g->nsub > 0)
mpvecdigmuladd(g->p, d, g->nsub, r->p);

x = g->indx;
for(i=0; i<g->nadd; i++){
t[0] = 0;
d0 = t[*x++];
for(j=1; j<d; j++)
t[j] = t[*x++];
t[0] = d0;

mpvecadd(r->p, d+1, t, d, r->p);
}

for(i=0; i<g->nsub; i++){
t[0] = 0;
d0 = t[*x++];
for(j=1; j<d; j++)
t[j] = t[*x++];
t[0] = d0;

mpvecsub(r->p, d+1, t, d, r->p);
}

mpvecdigmulsub(g->p, d, r->p[d], r->p);
r->p[d] = 0;

mpvecsub(r->p, d+1, g->p, d, r->p+d+1);
d0 = r->p[2*d+1];
for(j=0; j<d; j++)
r->p[j] = (r->p[j] & d0) | (r->p[j+d+1] & ~d0);

mpnorm(r);

return 0;
}

Mfield*
gmfield(mpint *N)
{
int i,j,d, s, *C, *X, *x, *e;
mpint *M, *T;
GMfield *g;

d = N->top;
if(d <= 2 || d > MAXDIG/2 || (mpsignif(N) % Dbits) != 0)
return nil;
g = nil;
T = mpnew(0);
M = mpcopy(N);
C = malloc(sizeof(int)*(d+1));
X = malloc(sizeof(int)*(d*d));
if(C == nil || X == nil)
goto out;

for(i=0; i<=d; i++){
if((M->p[i]>>8) != 0 && (~M->p[i]>>8) != 0)
goto out;
j = M->p[i];
C[d - i] = -j;
itomp(j, T);
mpleft(T, i*Dbits, T);
mpsub(M, T, M);
}
for(j=0; j<d; j++)
X[j] = C[d-j];
for(i=1; i<d; i++){
X[d*i] = X[d*(i-1) + d-1]*C[d];
for(j=1; j<d; j++)
X[d*i + j] = X[d*(i-1) + j-1] + X[d*(i-1) + d-1]*C[d-j];
}
g = mallocz(sizeof(GMfield) + (d+1)*sizeof(mpdigit)*2, 1);
if(g == nil)
goto out;

g->m2->p = (mpdigit*)&g[1];
g->m2->size = d*2+1;
mpmul(N, N, g->m2);
mpassign(N, g);
g->reduce = gmreduce;
g->flags |= MPfield;

s = 0;
x = g->indx;
e = x + nelem(g->indx) - d;
for(g->nadd=0; x <= e; x += d, g->nadd++){
s = 0;
for(i=0; i<d; i++){
for(j=0; j<d; j++){
if(X[d*i+j] > 0 && x[j] == 0){
X[d*i+j]--;
x[j] = d+i;
s = 1;
break;
}
}
}
if(s == 0)
break;
}
for(g->nsub=0; x <= e; x += d, g->nsub++){
s = 0;
for(i=0; i<d; i++){
for(j=0; j<d; j++){
if(X[d*i+j] < 0 && x[j] == 0){
X[d*i+j]++;
x[j] = d+i;
s = 1;
break;
}
}
}
if(s == 0)
break;
}
if(s != 0){
mpfree(g);
g = nil;
}
out:
free(C);
free(X);
mpfree(M);
mpfree(T);
return g;
}