version 1.4, 2000/08/22 05:04:06 |
version 1.5, 2001/10/09 01:36:13 |
|
|
* DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, |
* DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, |
* PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. |
* PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. |
* |
* |
* $OpenXM: OpenXM_contrib2/asir2000/engine/up.c,v 1.3 2000/08/21 08:31:28 noro Exp $ |
* $OpenXM: OpenXM_contrib2/asir2000/engine/up.c,v 1.4 2000/08/22 05:04:06 noro Exp $ |
*/ |
*/ |
#include "ca.h" |
#include "ca.h" |
#include <math.h> |
#include <math.h> |
Line 76 extern int lm_lazy; |
|
Line 76 extern int lm_lazy; |
|
extern int current_ff; |
extern int current_ff; |
extern int GC_dont_gc; |
extern int GC_dont_gc; |
|
|
void monicup(a,b) |
void monicup(UP a,UP *b) |
UP a; |
|
UP *b; |
|
{ |
{ |
UP w; |
UP w; |
|
|
|
|
} |
} |
} |
} |
|
|
void simpup(a,b) |
void simpup(UP a,UP *b) |
UP a; |
|
UP *b; |
|
{ |
{ |
int i,d; |
int i,d; |
UP c; |
UP c; |
|
|
} |
} |
} |
} |
|
|
void simpnum(a,b) |
void simpnum(Num a,Num *b) |
Num a; |
|
Num *b; |
|
{ |
{ |
LM lm; |
LM lm; |
GF2N gf; |
GF2N gf; |
|
|
} |
} |
} |
} |
|
|
void uremp(p1,p2,rp) |
void uremp(P p1,P p2,P *rp) |
P p1,p2; |
|
P *rp; |
|
{ |
{ |
UP n1,n2,r; |
UP n1,n2,r; |
|
|
|
|
} |
} |
} |
} |
|
|
void ugcdp(p1,p2,rp) |
void ugcdp(P p1,P p2,P *rp) |
P p1,p2; |
|
P *rp; |
|
{ |
{ |
UP n1,n2,r; |
UP n1,n2,r; |
|
|
|
|
uptop(r,rp); |
uptop(r,rp); |
} |
} |
|
|
void reversep(p1,d,rp) |
void reversep(P p1,Q d,P *rp) |
P p1; |
|
Q d; |
|
P *rp; |
|
{ |
{ |
UP n1,r; |
UP n1,r; |
|
|
|
|
} |
} |
} |
} |
|
|
void invmodp(p1,d,rp) |
void invmodp(P p1,Q d,P *rp) |
P p1; |
|
Q d; |
|
P *rp; |
|
{ |
{ |
UP n1,r; |
UP n1,r; |
|
|
|
|
} |
} |
} |
} |
|
|
void addup(n1,n2,nr) |
void addup(UP n1,UP n2,UP *nr) |
UP n1,n2; |
|
UP *nr; |
|
{ |
{ |
UP r,t; |
UP r,t; |
int i,d1,d2; |
int i,d1,d2; |
|
|
} |
} |
} |
} |
|
|
void subup(n1,n2,nr) |
void subup(UP n1,UP n2,UP *nr) |
UP n1,n2; |
|
UP *nr; |
|
{ |
{ |
UP r; |
UP r; |
int i,d1,d2,d; |
int i,d1,d2,d; |
|
|
} |
} |
} |
} |
|
|
void chsgnup(n1,nr) |
void chsgnup(UP n1,UP *nr) |
UP n1; |
|
UP *nr; |
|
{ |
{ |
UP r; |
UP r; |
int d1,i; |
int d1,i; |
|
|
} |
} |
} |
} |
|
|
void hybrid_mulup(ff,n1,n2,nr) |
void hybrid_mulup(int ff,UP n1,UP n2,UP *nr) |
int ff; |
|
UP n1,n2; |
|
UP *nr; |
|
{ |
{ |
if ( !n1 || !n2 ) |
if ( !n1 || !n2 ) |
*nr = 0; |
*nr = 0; |
|
|
} |
} |
} |
} |
|
|
void hybrid_squareup(ff,n1,nr) |
void hybrid_squareup(int ff,UP n1,UP *nr) |
int ff; |
|
UP n1; |
|
UP *nr; |
|
{ |
{ |
if ( !n1 ) |
if ( !n1 ) |
*nr = 0; |
*nr = 0; |
|
|
} |
} |
} |
} |
|
|
void hybrid_tmulup(ff,n1,n2,d,nr) |
void hybrid_tmulup(int ff,UP n1,UP n2,int d,UP *nr) |
int ff; |
|
UP n1,n2; |
|
int d; |
|
UP *nr; |
|
{ |
{ |
if ( !n1 || !n2 ) |
if ( !n1 || !n2 ) |
*nr = 0; |
*nr = 0; |
|
|
} |
} |
} |
} |
|
|
void mulup(n1,n2,nr) |
void mulup(UP n1,UP n2,UP *nr) |
UP n1,n2; |
|
UP *nr; |
|
{ |
{ |
UP r; |
UP r; |
Num *pc1,*pc,*c1,*c2,*c; |
Num *pc1,*pc,*c1,*c2,*c; |
|
|
|
|
/* nr = c*n1 */ |
/* nr = c*n1 */ |
|
|
void mulcup(c,n1,nr) |
void mulcup(Num c,UP n1,UP *nr) |
Num c; |
|
UP n1; |
|
UP *nr; |
|
{ |
{ |
int d; |
int d; |
UP r; |
UP r; |
|
|
} |
} |
} |
} |
|
|
void tmulup(n1,n2,d,nr) |
void tmulup(UP n1,UP n2,int d,UP *nr) |
UP n1,n2; |
|
int d; |
|
UP *nr; |
|
{ |
{ |
UP r; |
UP r; |
Num *pc1,*pc,*c1,*c2,*c; |
Num *pc1,*pc,*c1,*c2,*c; |
|
|
} |
} |
} |
} |
|
|
void squareup(n1,nr) |
void squareup(UP n1,UP *nr) |
UP n1; |
|
UP *nr; |
|
{ |
{ |
UP r; |
UP r; |
Num *c1,*c; |
Num *c1,*c; |
|
|
} |
} |
} |
} |
|
|
void remup(n1,n2,nr) |
void remup(UP n1,UP n2,UP *nr) |
UP n1,n2; |
|
UP *nr; |
|
{ |
{ |
UP w,r; |
UP w,r; |
|
|
|
|
} |
} |
} |
} |
|
|
void remup_destructive(n1,n2) |
void remup_destructive(UP n1,UP n2) |
UP n1,n2; |
|
{ |
{ |
Num *c1,*c2; |
Num *c1,*c2; |
int d1,d2,i,j; |
int d1,d2,i,j; |
|
|
n1->d = i; |
n1->d = i; |
} |
} |
|
|
void qrup(n1,n2,nq,nr) |
void qrup(UP n1,UP n2,UP *nq,UP *nr) |
UP n1,n2; |
|
UP *nq,*nr; |
|
{ |
{ |
UP w,r,q; |
UP w,r,q; |
struct oUP t; |
struct oUP t; |
|
|
} |
} |
} |
} |
|
|
void qrup_destructive(n1,n2) |
void qrup_destructive(UP n1,UP n2) |
UP n1,n2; |
|
{ |
{ |
Num *c1,*c2; |
Num *c1,*c2; |
int d1,d2,i,j; |
int d1,d2,i,j; |
|
|
n1->d = i; |
n1->d = i; |
} |
} |
|
|
void gcdup(n1,n2,nr) |
void gcdup(UP n1,UP n2,UP *nr) |
UP n1,n2; |
|
UP *nr; |
|
{ |
{ |
UP w1,w2,t,r; |
UP w1,w2,t,r; |
int d1,d2; |
int d1,d2; |
|
|
|
|
/* compute r s.t. a*r = 1 mod m */ |
/* compute r s.t. a*r = 1 mod m */ |
|
|
void extended_gcdup(a,m,r) |
void extended_gcdup(UP a,UP m,UP *r) |
UP a,m; |
|
UP *r; |
|
{ |
{ |
UP one,g1,g2,a1,a2,a3,b1,b2,b3,inv,quo,rem,t; |
UP one,g1,g2,a1,a2,a3,b1,b2,b3,inv,quo,rem,t; |
Num i; |
Num i; |
|
|
mulup(b2,inv,r); |
mulup(b2,inv,r); |
} |
} |
|
|
void reverseup(n1,d,nr) |
void reverseup(UP n1,int d,UP *nr) |
UP n1; |
|
int d; |
|
UP *nr; |
|
{ |
{ |
Num *c1,*c; |
Num *c1,*c; |
int i,d1; |
int i,d1; |
|
|
} |
} |
} |
} |
|
|
void invmodup(n1,d,nr) |
void invmodup(UP n1,int d,UP *nr) |
UP n1; |
|
int d; |
|
UP *nr; |
|
{ |
{ |
UP r; |
UP r; |
Num s,t,u,hinv; |
Num s,t,u,hinv; |
|
|
} |
} |
} |
} |
|
|
void pwrup(n,e,nr) |
void pwrup(UP n,Q e,UP *nr) |
UP n; |
|
Q e; |
|
UP *nr; |
|
{ |
{ |
UP y,z,t; |
UP y,z,t; |
N en,qn; |
N en,qn; |
|
|
} |
} |
} |
} |
|
|
int compup(n1,n2) |
int compup(UP n1,UP n2) |
UP n1,n2; |
|
{ |
{ |
int i,r; |
int i,r; |
|
|
|
|
} |
} |
} |
} |
|
|
void kmulp(vl,n1,n2,nr) |
void kmulp(VL vl,P n1,P n2,P *nr) |
VL vl; |
|
P n1,n2; |
|
P *nr; |
|
{ |
{ |
UP b1,b2,br; |
UP b1,b2,br; |
|
|
|
|
} |
} |
} |
} |
|
|
void ksquarep(vl,n1,nr) |
void ksquarep(VL vl,P n1,P *nr) |
VL vl; |
|
P n1; |
|
P *nr; |
|
{ |
{ |
UP b1,br; |
UP b1,br; |
|
|
|
|
} |
} |
} |
} |
|
|
void kmulup(n1,n2,nr) |
void kmulup(UP n1,UP n2,UP *nr) |
UP n1,n2,*nr; |
|
{ |
{ |
UP n,t,s,m,carry; |
int d1,d2; |
int d,d1,d2,len,i,l; |
|
pointer *r,*r0; |
|
|
|
if ( !n1 || !n2 ) { |
if ( !n1 || !n2 ) { |
*nr = 0; return; |
*nr = 0; return; |
|
|
kmulupmain(n1,n2,nr); |
kmulupmain(n1,n2,nr); |
} |
} |
|
|
void ksquareup(n1,nr) |
void ksquareup(UP n1,UP *nr) |
UP n1,*nr; |
|
{ |
{ |
int d1; |
int d1; |
extern Q TWO; |
extern Q TWO; |
|
|
ksquareupmain(n1,nr); |
ksquareupmain(n1,nr); |
} |
} |
|
|
void copyup(n1,n2) |
void copyup(UP n1,UP n2) |
UP n1,n2; |
|
{ |
{ |
n2->d = n1->d; |
n2->d = n1->d; |
bcopy((char *)n1->c,(char *)n2->c,(n1->d+1)*sizeof(Q)); |
bcopy((char *)n1->c,(char *)n2->c,(n1->d+1)*sizeof(Q)); |
} |
} |
|
|
void c_copyup(n,len,p) |
void c_copyup(UP n,int len,pointer *p) |
UP n; |
|
int len; |
|
pointer *p; |
|
{ |
{ |
if ( n ) |
if ( n ) |
bcopy((char *)COEF(n),(char *)p,MIN((DEG(n)+1),len)*sizeof(pointer)); |
bcopy((char *)COEF(n),(char *)p,MIN((DEG(n)+1),len)*sizeof(pointer)); |
} |
} |
|
|
void kmulupmain(n1,n2,nr) |
void kmulupmain(UP n1,UP n2,UP *nr) |
UP n1,n2,*nr; |
|
{ |
{ |
int d1,d2,h,len; |
int d1,d2,h; |
UP n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2; |
UP n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2; |
|
|
d1 = DEG(n1)+1; d2 = DEG(n2)+1; h = (MAX(d1,d2)+1)/2; |
d1 = DEG(n1)+1; d2 = DEG(n2)+1; h = (MAX(d1,d2)+1)/2; |
|
|
addup(t1,t2,nr); |
addup(t1,t2,nr); |
} |
} |
|
|
void ksquareupmain(n1,nr) |
void ksquareupmain(UP n1,UP *nr) |
UP n1,*nr; |
|
{ |
{ |
int d1,h,len; |
int d1,h; |
UP n1lo,n1hi,hi,lo,mid1,mid2,mid,s1,t1,t2; |
UP n1lo,n1hi,hi,lo,mid1,mid2,mid,s1,t1,t2; |
|
|
d1 = DEG(n1)+1; h = (d1+1)/2; |
d1 = DEG(n1)+1; h = (d1+1)/2; |
|
|
addup(t1,t2,nr); |
addup(t1,t2,nr); |
} |
} |
|
|
void rembymulup(n1,n2,nr) |
void rembymulup(UP n1,UP n2,UP *nr) |
UP n1,n2; |
|
UP *nr; |
|
{ |
{ |
int d1,d2,d; |
int d1,d2,d; |
UP r1,r2,inv2,t,s,q; |
UP r1,r2,inv2,t,s,q; |
|
|
inv2 = inverse of reversep(n2) mod x^d |
inv2 = inverse of reversep(n2) mod x^d |
*/ |
*/ |
|
|
void hybrid_rembymulup_special(ff,n1,n2,inv2,nr) |
void hybrid_rembymulup_special(int ff,UP n1,UP n2,UP inv2,UP *nr) |
int ff; |
|
UP n1,n2,inv2; |
|
UP *nr; |
|
{ |
{ |
int d1,d2,d; |
int d1,d2,d; |
UP r1,t,s,q,u; |
UP r1,t,s,q; |
|
|
if ( !n2 ) |
if ( !n2 ) |
error("hybrid_rembymulup : division by 0"); |
error("hybrid_rembymulup : division by 0"); |
|
|
} |
} |
} |
} |
|
|
void rembymulup_special(n1,n2,inv2,nr) |
void rembymulup_special(UP n1,UP n2,UP inv2,UP *nr) |
UP n1,n2,inv2; |
|
UP *nr; |
|
{ |
{ |
int d1,d2,d; |
int d1,d2,d; |
UP r1,t,s,q,u; |
UP r1,t,s,q; |
|
|
if ( !n2 ) |
if ( !n2 ) |
error("rembymulup : division by 0"); |
error("rembymulup : division by 0"); |
|
|
|
|
/* *nr = n1*n2 mod x^d */ |
/* *nr = n1*n2 mod x^d */ |
|
|
void tkmulup(n1,n2,d,nr) |
void tkmulup(UP n1,UP n2,int d,UP *nr) |
UP n1,n2; |
|
int d; |
|
UP *nr; |
|
{ |
{ |
int m; |
int m; |
UP n1l,n1h,n2l,n2h,l,h,t,s,u,afo; |
UP n1l,n1h,n2l,n2h,l,h,t,s,u; |
|
|
if ( d < 0 ) |
if ( d < 0 ) |
error("tkmulup : invalid argument"); |
error("tkmulup : invalid argument"); |
|
|
|
|
/* n->n*x^d */ |
/* n->n*x^d */ |
|
|
void shiftup(n,d,nr) |
void shiftup(UP n,int d,UP *nr) |
UP n; |
|
int d; |
|
UP *nr; |
|
{ |
{ |
int dr; |
int dr; |
UP r; |
UP r; |
|
|
} |
} |
} |
} |
|
|
void fft_rembymulup_special(n1,n2,inv2,nr) |
void fft_rembymulup_special(UP n1,UP n2,UP inv2,UP *nr) |
UP n1,n2,inv2; |
|
UP *nr; |
|
{ |
{ |
int d1,d2,d; |
int d1,d2,d; |
UP r1,t,s,q,u; |
UP r1,t,s,q,u; |
|
|
} |
} |
} |
} |
|
|
void set_degreeup(n,d) |
void set_degreeup(UP n,int d) |
UP n; |
|
int d; |
|
{ |
{ |
int i; |
int i; |
|
|
|
|
|
|
/* n -> n0 + x^d n1 */ |
/* n -> n0 + x^d n1 */ |
|
|
void decompup(n,d,n0,n1) |
void decompup(UP n,int d,UP *n0,UP *n1) |
UP n; |
|
int d; |
|
UP *n0,*n1; |
|
{ |
{ |
int dn; |
int dn; |
UP r0,r1; |
UP r0,r1; |
|
|
|
|
/* n -> n mod x^d */ |
/* n -> n mod x^d */ |
|
|
void truncup(n1,d,nr) |
void truncup(UP n1,int d,UP *nr) |
UP n1; |
|
int d; |
|
UP *nr; |
|
{ |
{ |
int i; |
int i; |
UP r; |
UP r; |
|
|
} |
} |
} |
} |
|
|
int int_bits(t) |
int int_bits(int t) |
int t; |
|
{ |
{ |
int k; |
int k; |
|
|
|
|
|
|
/* n is assumed to be LM or integer coefficient */ |
/* n is assumed to be LM or integer coefficient */ |
|
|
int maxblenup(n) |
int maxblenup(UP n) |
UP n; |
|
{ |
{ |
int m,r,i,d; |
int m,r,i,d; |
Num *c; |
Num *c; |
|
|
return r; |
return r; |
} |
} |
|
|
void uptofmarray(mod,n,f) |
void uptofmarray(int mod,UP n,ModNum *f) |
int mod; |
|
UP n; |
|
ModNum *f; |
|
{ |
{ |
int d,i; |
int d,i; |
unsigned int r; |
unsigned int r; |
|
|
} |
} |
} |
} |
|
|
void fmarraytoup(f,d,nr) |
void fmarraytoup(ModNum *f,int d,UP *nr) |
ModNum *f; |
|
int d; |
|
UP *nr; |
|
{ |
{ |
int i; |
int i; |
Q *c; |
Q *c; |
|
|
|
|
/* f[i]: an array of length n */ |
/* f[i]: an array of length n */ |
|
|
void uiarraytoup(f,n,d,nr) |
void uiarraytoup(unsigned int **f,int n,int d,UP *nr) |
unsigned int **f; |
|
int n,d; |
|
UP *nr; |
|
{ |
{ |
int i,j; |
int i,j; |
unsigned int *fi; |
unsigned int *fi; |
|
|
} |
} |
} |
} |
|
|
void adj_coefup(n,m,m2,nr) |
void adj_coefup(UP n,N m,N m2,UP *nr) |
UP n; |
|
N m,m2; |
|
UP *nr; |
|
{ |
{ |
int d; |
int d; |
Q *c,*cr; |
Q *c,*cr; |
|
|
|
|
/* n is assumed to have positive integer coefficients. */ |
/* n is assumed to have positive integer coefficients. */ |
|
|
void remcup(n,mod,nr) |
void remcup(UP n,N mod,UP *nr) |
UP n; |
|
N mod; |
|
UP *nr; |
|
{ |
{ |
int i,d; |
int i,d; |
Q *c,*cr; |
Q *c,*cr; |
|
|
|
|
void fft_mulup_main(UP,UP,int,UP *); |
void fft_mulup_main(UP,UP,int,UP *); |
|
|
void fft_mulup(n1,n2,nr) |
void fft_mulup(UP n1,UP n2,UP *nr) |
UP n1,n2; |
|
UP *nr; |
|
{ |
{ |
int d1,d2,d,b1,b2,h; |
int d1,d2,d,b1,b2,h; |
UP n1lo,n1hi,n2lo,n2hi,lo,hi,t1,t2,mid1,mid2,mid,s1,s2; |
UP n1lo,n1hi,n2lo,n2hi,lo,hi,t1,t2,mid1,mid2,mid,s1,s2; |
|
|
} |
} |
} |
} |
|
|
void trunc_fft_mulup(n1,n2,dbd,nr) |
void trunc_fft_mulup(UP n1,UP n2,int dbd,UP *nr) |
UP n1,n2; |
|
int dbd; |
|
UP *nr; |
|
{ |
{ |
int d1,d2,b1,b2,m; |
int d1,d2,b1,b2,m; |
UP n1l,n1h,n2l,n2h,l,h,t,s,u; |
UP n1l,n1h,n2l,n2h,l,h,t,s,u; |
|
|
} |
} |
} |
} |
|
|
void fft_squareup(n1,nr) |
void fft_squareup(UP n1,UP *nr) |
UP n1; |
|
UP *nr; |
|
{ |
{ |
int d1,d,h,len,b1; |
int d1,d,h,b1; |
UP n1lo,n1hi,hi,lo,mid1,mid2,mid,s1,t1,t2; |
UP n1lo,n1hi,hi,lo,mid1,mid2,mid,s1,t1,t2; |
|
|
if ( !n1 ) |
if ( !n1 ) |
|
|
* n1 == n2 => squaring |
* n1 == n2 => squaring |
*/ |
*/ |
|
|
void fft_mulup_main(n1,n2,dbd,nr) |
void fft_mulup_main(UP n1,UP n2,int dbd,UP *nr) |
UP n1,n2; |
|
UP *nr; |
|
{ |
{ |
ModNum *f1,*f2,*w,*fr; |
ModNum *f1,*f2,*w,*fr; |
ModNum **frarray,**fa; |
ModNum **frarray,**fa; |
|
|
#if 0 |
#if 0 |
/* inefficient version */ |
/* inefficient version */ |
|
|
void crup(f,d,mod,index,m,r) |
void crup(ModNum **f,int d,int *mod,int index,N m,UP *r) |
ModNum **f; |
|
int d; |
|
int *mod; |
|
int index; |
|
N m; |
|
UP *r; |
|
{ |
{ |
N *cof,*c; |
N *cof,*c; |
int *inv; |
int *inv; |
|
|
#else |
#else |
/* improved version */ |
/* improved version */ |
|
|
void crup(f,d,mod,index,m,r) |
void crup(ModNum **f,int d,int *mod,int index,N m,UP *r) |
ModNum **f; |
|
int d; |
|
int *mod; |
|
int index; |
|
N m; |
|
UP *r; |
|
{ |
{ |
N cof,c,t,w,w1; |
N cof,c,t,w,w1; |
struct oN fc; |
struct oN fc; |
|
|
* return: n1*n2 mod Primes[modind[0]]*.prime...*Primes[modind[nmod-1]].prime |
* return: n1*n2 mod Primes[modind[0]]*.prime...*Primes[modind[nmod-1]].prime |
*/ |
*/ |
|
|
void fft_mulup_specialmod_main(n1,n2,dbd,modind,nmod,nr) |
void fft_mulup_specialmod_main(UP n1,UP n2,int dbd,int *modind,int nmod,UP *nr) |
UP n1,n2; |
|
int dbd; |
|
int *modind; |
|
int nmod; |
|
UP *nr; |
|
{ |
{ |
ModNum *f1,*f2,*w,*fr; |
ModNum *f1,*f2,*w,*fr; |
ModNum **frarray,**fa; |
ModNum **frarray; |
N m,m1,m2; |
N m,m1,m2; |
unsigned int *modarray; |
unsigned int *modarray; |
int d1,d2,dmin,i,mod,root,d,cond,bound; |
int d1,d2,dmin,i,root,d,cond,bound; |
UP r; |
|
|
|
if ( !n1 || !n2 ) { |
if ( !n1 || !n2 ) { |
*nr = 0; return; |
*nr = 0; return; |