===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/gfs.c,v
retrieving revision 1.8
retrieving revision 1.14
diff -u -p -r1.8 -r1.14
--- OpenXM_contrib2/asir2000/engine/gfs.c	2001/06/29 09:08:53	1.8
+++ OpenXM_contrib2/asir2000/engine/gfs.c	2002/12/18 06:15:40	1.14
@@ -45,9 +45,10 @@
  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
  *
- * $OpenXM: OpenXM_contrib2/asir2000/engine/gfs.c,v 1.7 2001/06/21 07:47:02 noro Exp $
+ * $OpenXM: OpenXM_contrib2/asir2000/engine/gfs.c,v 1.13 2002/09/30 06:13:07 noro Exp $
 */
 #include "ca.h"
+#include "inline.h"
 
 /* q = p^n */
 
@@ -59,9 +60,6 @@ int *current_gfs_plus1;
 int *current_gfs_ntoi;
 int *current_gfs_iton;
 
-void chsgngfs();
-void generate_defpoly_um();
-
 struct prim_root_info {
 	int p;
 	int extdeg;
@@ -69,6 +67,10 @@ struct prim_root_info {
 	int prim_root;
 };
 
+/* if p >= SF_THRESHOLD, usual representation is used */
+
+#define SF_THRESHOLD 16384
+
 struct prim_root_info prim_root_info_tab[] = {
 {2,1,0,0}, {2,2,7,2}, {2,3,11,2}, {2,4,19,2}, {2,5,37,2}, {2,6,67,2},
 {2,7,131,2}, {2,8,283,3}, {2,9,515,7}, {2,10,1033,2}, {2,11,2053,2},
@@ -473,9 +475,7 @@ struct prim_root_info prim_root_info_tab[] = {
 {16349,1,0,2}, {16361,1,0,3}, {16363,1,0,2}, {16369,1,0,7}, {16381,1,0,2},
 };
 
-void dec_um(p,a,u)
-int p,a;
-UM u;
+void dec_um(int p,int a,UM u)
 {
 	int i;
 
@@ -491,12 +491,24 @@ UM u;
  * current_gfs_iton[p-1] = 0
  */
 
-void setmod_sf(p,n)
-int p,n;
+void setmod_sf(int p,int n)
 {
 	int r,i,q1,q,t,t1;
 	UM dp;
 
+	if ( p >= SF_THRESHOLD ) {
+		if ( n > 1 )
+			error("setmod_ff : p^n is too large");
+		current_gfs_p = p;
+		current_gfs_q = p;
+		current_gfs_q1 = p-1;
+		current_gfs_ext = 0;
+		current_gfs_ntoi = 0;
+		current_gfs_iton = 0;
+		current_gfs_plus1 = 0;
+		return;
+	}
+
 	for ( i = 0, q = 1; i < n; i++ )
 		q *= p;
 	dp = UMALLOC(n);
@@ -504,6 +516,8 @@ int p,n;
 	if ( !r ) {
 		generate_defpoly_um(p,n,dp);
 		r = generate_primitive_root_enc(p,n,dp);
+		if ( !r )
+			error("setmod_sf : primitive root not found");
 	}
 	current_gfs_p = p;
 	current_gfs_q = q;
@@ -531,9 +545,7 @@ int p,n;
 	}
 }
 
-int search_defpoly_and_primitive_root(p,n,dp)
-int p,n;
-UM dp;
+int search_defpoly_and_primitive_root(int p,int n,UM dp)
 {
 	int l,min,max,mid,p1,i,ind,t;
 
@@ -570,11 +582,9 @@ UM dp;
 	return prim_root_info_tab[i].prim_root;
 }
 
-void generate_defpoly_um(p,n,dp)
-int p,n;
-UM dp;
+void generate_defpoly_um(int p,int n,UM dp)
 {
-	int i,j,a,c,q;
+	int i,j,a,q;
 	UM wf,wdf,wgcd;
 
 	wf = W_UMALLOC(n);
@@ -600,9 +610,7 @@ UM dp;
 	}
 }
 
-int generate_primitive_root_enc(p,n,dp)
-int p,n;
-UM dp;
+int generate_primitive_root_enc(int p,int n,UM dp)
 {
 	int i,r,rj,j,q;
 
@@ -618,14 +626,13 @@ UM dp;
 		if ( j == q-1 )
 			return r;
 	}
+	/* not found */
+	return 0;
 }
 
 /* [a(p)]*[b(p)] in GF(p^n) -> [a(x)*b(x) mod dp(x)]_{x->p} */
 
-int mulremum_enc(p,n,dp,a,b)
-int p,n;
-UM dp;
-int a,b;
+int mulremum_enc(int p,int n,UM dp,int a,int b)
 {
 	int i,dr,r;
 	UM wa,wb,wc,wq;
@@ -650,30 +657,26 @@ int a,b;
 	return r;
 }
 
-void gfs_galois_action(a,e,c)
-GFS a;
-Q e;
-GFS *c;
+/* sigma : alpha -> alpha^q */
+
+void gfs_galois_action(GFS a,Q e,GFS *c)
 {
-	Q p;
+	Q q;
 	int i,k;
 	GFS t,s;
 
 	t = a;
 	k = QTOS(e);
-	STOQ(current_gfs_p,p);
+	STOQ(current_gfs_q,q);
 	for ( i = 0; i < k; i++ ) {
-		pwrgfs(t,p,&s); t = s;
+		pwrgfs(t,q,&s); t = s;
 	}
 	*c = t;
 }
 
 /* GF(pn)={0,1,a,a^2,...} -> GF(pm)={0,1,b,b^2,...}; a->b^k */
 
-void gfs_embed(z,k,pm,c)
-GFS z;
-int k,pm;
-GFS *c;
+void gfs_embed(GFS z,int k,int pm,GFS *c)
 {
 	int t;
 
@@ -684,46 +687,70 @@ GFS *c;
 		MKGFS(t,*c);
 	}
 }
-void qtogfs(a,c)
-Q a;
-GFS *c;
+
+/* 0 <= index <= q-1 */
+
+void indextogfs(int index, GFS *c)
 {
-	int s;
+	if ( index == 0 )
+		*c = 0;
+	else if ( index >= current_gfs_q )
+		error("indextogfs : exhausted");
+	else if ( !current_gfs_ntoi ) {
+		MKGFS(index,*c);
+	} else {
+		MKGFS(index-1,*c);
+	}
+}
 
-	s = QTOS(a)%current_gfs_q;
-	if ( s < 0 )
-		s += current_gfs_q;
-	if ( !s )
+void itogfs(int n, GFS *c)
+{
+	n = _itosf(n);
+	if ( !n )
 		*c = 0;
-	else
-		MKGFS(current_gfs_ntoi[s],*c);	
+	else {
+		n = IFTOF(n);
+		MKGFS(n,*c);
+	}
 }
 
-void mqtogfs(a,c)
-MQ a;
-GFS *c;
+void iftogfs(int n, GFS *c)
 {
-	if ( !a )
+	if ( !n )
 		*c = 0;
 	else {
-		MKGFS(current_gfs_ntoi[CONT(a)],*c);	
+		MKGFS(IFTOF(n),*c);
 	}
 }
 
-void gfstomq(a,c)
-GFS a;
-MQ *c;
+void qtogfs(Q a,GFS *c)
 {
+	int s;
+
+	s = QTOS(a)%current_gfs_q;
+	itogfs(s,c);
+}
+
+void mqtogfs(MQ a,GFS *c)
+{
 	if ( !a )
 		*c = 0;
-	else {
+	else
+		itogfs(CONT(a),c);
+}
+
+void gfstomq(GFS a,MQ *c)
+{
+	if ( !a )
+		*c = 0;
+	else if ( !current_gfs_ntoi ) {
+		UTOMQ(CONT(a),*c);
+	} else {
 		UTOMQ(current_gfs_iton[CONT(a)],*c);
 	}
 }
 
-void ntogfs(a,b)
-Obj a;
-GFS *b;
+void ntogfs(Obj a,GFS *b)
 {
 	P t;
 
@@ -732,28 +759,36 @@ GFS *b;
 	if ( !a || (OID(a)==O_N && NID(a) == N_GFS) )
 		*b = (GFS)a;
 	else if ( OID(a) == O_N && NID(a) == N_M )
-		mqtogfs(a,b);
+		mqtogfs((MQ)a,b);
 	else if ( OID(a) == O_N && NID(a) == N_Q ) {
-		ptomp(current_gfs_p,(P)a,&t); mqtogfs(t,b);
+		ptomp(current_gfs_p,(P)a,&t); mqtogfs((MQ)t,b);
 	} else
 		error("ntogfs : invalid argument");
 }
 
-void addgfs(a,b,c)
-GFS a,b;
-GFS *c;
+void addgfs(GFS a,GFS b,GFS *c)
 {
 	int ai,bi,ci;
 	GFS z;
 
-	ntogfs(a,&z); a = z;
-	ntogfs(b,&z); b = z;
+	ntogfs((Obj)a,&z); a = z;
+	ntogfs((Obj)b,&z); b = z;
 	if ( !a )
 		*c = b;
 	else if ( !b )
 		*c = a;
-	else {
+	else if ( !current_gfs_ntoi ) {
 		ai = CONT(a); bi = CONT(b);
+		ci = ai+bi-current_gfs_q;
+		if ( ci == 0 )
+			*c = 0;
+		else {
+			if ( ci < 0 )
+				ci += current_gfs_q;	
+			MKGFS(ci,*c);
+		}
+	} else {
+		ai = CONT(a); bi = CONT(b);
 		if ( ai > bi ) {
 			/* tab[ai]+tab[bi] = tab[bi](tab[ai-bi]+1) */
 			ci = current_gfs_plus1[ai-bi];
@@ -780,14 +815,12 @@ GFS *c;
 	}
 }
 
-void subgfs(a,b,c)
-GFS a,b;
-GFS *c;
+void subgfs(GFS a,GFS b,GFS *c)
 {
 	GFS t,z;
 
-	ntogfs(a,&z); a = z;
-	ntogfs(b,&z); b = z;
+	ntogfs((Obj)a,&z); a = z;
+	ntogfs((Obj)b,&z); b = z;
 	if ( !b )
 		*c = a;
 	else {
@@ -796,18 +829,20 @@ GFS *c;
 	}
 }
 
-void mulgfs(a,b,c)
-GFS a,b;
-GFS *c;
+void mulgfs(GFS a,GFS b,GFS *c)
 {
-	int ai;
+	int ai,bi,ci;
 	GFS z;
 
-	ntogfs(a,&z); a = z;
-	ntogfs(b,&z); b = z;
+	ntogfs((Obj)a,&z); a = z;
+	ntogfs((Obj)b,&z); b = z;
 	if ( !a || !b )
 		*c = 0;
-	else {
+	else if ( !current_gfs_ntoi ) {
+		ai = CONT(a); bi = CONT(b);
+		DMAR(ai,bi,0,current_gfs_q,ci);
+		MKGFS(ci,*c);
+	} else {
 		ai = CONT(a) + CONT(b);
 		if ( ai >= current_gfs_q1 )
 			ai -= current_gfs_q1;
@@ -815,20 +850,22 @@ GFS *c;
 	}
 }
 
-void divgfs(a,b,c)
-GFS a,b;
-GFS *c;
+void divgfs(GFS a,GFS b,GFS *c)
 {
-	int ai;
+	int ai,bi,ci;
 	GFS z;
 
-	ntogfs(a,&z); a = z;
-	ntogfs(b,&z); b = z;
+	ntogfs((Obj)a,&z); a = z;
+	ntogfs((Obj)b,&z); b = z;
 	if ( !b )
 		error("divgfs : division by 0");
 	else if ( !a )
 		*c = 0;
-	else {
+	else if ( !current_gfs_ntoi ) {
+		ai = CONT(a); bi = invm(CONT(b),current_gfs_q);
+		DMAR(ai,bi,0,current_gfs_q,ci);
+		MKGFS(ci,*c);
+	} else {
 		ai = CONT(a) - CONT(b);
 		if ( ai < 0 )
 			ai += current_gfs_q1;
@@ -836,16 +873,18 @@ GFS *c;
 	}
 }
 
-void chsgngfs(a,c)
-GFS a,*c;
+void chsgngfs(GFS a,GFS *c)
 {
 	int ai;
 	GFS z;
 
-	ntogfs(a,&z); a = z;
+	ntogfs((Obj)a,&z); a = z;
 	if ( !a )
 		*c = 0;
-	else if ( current_gfs_q1&1 )
+	else if ( !current_gfs_ntoi ) {
+		ai = current_gfs_q - CONT(a);
+		MKGFS(ai,*c);
+	} else if ( current_gfs_q1&1 )
 		*c = a;
 	else {	
 		/* r^((q-1)/2) = -1 */
@@ -856,40 +895,40 @@ GFS a,*c;
 	}
 }
 
-void pwrgfs(a,b,c)
-GFS a;
-Q b;
-GFS *c;
+void pwrgfs(GFS a,Q b,GFS *c)
 {
 	N an,tn,rn;
 	GFS t,s,z;
+	int ai;
 
-	ntogfs(a,&z); a = z;
+	ntogfs((Obj)a,&z); a = z;
 	if ( !b )
-		MKGFS(0,*c);
+		itogfs(1,c);
 	else if ( !a )
 		*c = 0;
-	else {
+	else if ( !current_gfs_ntoi) {
+		ai = pwrm(current_gfs_q,CONT(a),QTOS(b));
+		MKGFS(ai,*c);
+	} else {
 		STON(CONT(a),an); muln(an,NM(b),&tn);
 		STON(current_gfs_q1,an); remn(tn,an,&rn);
 		if ( !rn )
-			MKGFS(0,*c);
+			itogfs(1,c);
 		else if ( SGN(b) > 0 )
 			MKGFS(BD(rn)[0],*c);
 		else {
-			MKGFS(0,t);
+			itogfs(1,&t);
 			MKGFS(BD(rn)[0],s);
 			divgfs(t,s,c);
 		}
 	}
 }
 
-int cmpgfs(a,b)
-GFS a,b;
+int cmpgfs(GFS a,GFS b)
 {
 	GFS z;
 
-	ntogfs(a,&z); a = z;
+	ntogfs((Obj)a,&z); a = z;
 	if ( !a )
 		return !b ? 0 : -1;
 	else
@@ -905,27 +944,34 @@ GFS a,b;
 		}
 }
 
-void randomgfs(r)
-GFS *r;
+void pthrootgfs(GFS a,GFS *b)
 {
+	Q p;
+	int e,i;
+	GFS t,s;
+
+	STOQ(characteristic_sf(),p);
+	e = extdeg_sf()-1;
+	t = a;
+	for ( i = 0; i < e; i++ ) {
+		pwrgfs(t,p,&s); t = s;
+	}
+	*b = t;
+}
+
+void randomgfs(GFS *r)
+{
 	unsigned int t;
 
 	if ( !current_gfs_q1 )
 		error("addgfs : current_gfs_q is not set");
 	t = mt_genrand()%current_gfs_q;
-	if ( !t )
-		*r = 0;
-	else {
-		if ( t == current_gfs_q1 )
-			t = 0;
-		MKGFS(t,*r);
-	}
+	indextogfs(t,r);
 }
 
 /* arithmetic operations for 'immediate values of GFS */
 
-int _addsf(a,b)
-int a,b;
+int _addsf(int a,int b)
 {
 	if ( !a )
 		return b;
@@ -933,6 +979,18 @@ int a,b;
 		return a;
 
 	a = IFTOF(a); b = IFTOF(b);
+
+	if ( !current_gfs_ntoi ) {
+		a = a+b-current_gfs_q;
+		if ( a == 0 )
+			return 0;
+		else { 
+			if ( a < 0 )
+				a += current_gfs_q;	
+			return FTOIF(a);
+		}
+	}
+
 	if ( a > b ) {
 		/* tab[a]+tab[b] = tab[b](tab[a-b]+1) */
 		a = current_gfs_plus1[a-b];
@@ -958,12 +1016,14 @@ int a,b;
 	}
 }
 
-int _chsgnsf(a)
-int a;
+int _chsgnsf(int a)
 {
 	if ( !a )
 		return 0;
-	else if ( current_gfs_q1&1 )
+	else if ( !current_gfs_ntoi ) {
+		a = current_gfs_q-IFTOF(a);
+		return FTOIF(a);
+	} else if ( current_gfs_q1&1 )
 		return a;
 	else {	
 		/* r^((q-1)/2) = -1 */
@@ -975,8 +1035,7 @@ int a;
 	}
 }
 
-int _subsf(a,b)
-int a,b;
+int _subsf(int a,int b)
 {
 	if ( !a )
 		return _chsgnsf(b);
@@ -986,12 +1045,17 @@ int a,b;
 		return _addsf(a,_chsgnsf(b));
 }
 
-int _mulsf(a,b)
-int a,b;
+int _mulsf(int a,int b)
 {
+	int c;
+
 	if ( !a || !b )
 		return 0;
-	else {
+	else if ( !current_gfs_ntoi ) {
+		a = IFTOF(a); b = IFTOF(b);
+		DMAR(a,b,0,current_gfs_q,c);
+		return FTOIF(c);
+	} else {
 		a = IFTOF(a) + IFTOF(b);
 		if ( a >= current_gfs_q1 )
 			a -= current_gfs_q1;
@@ -999,23 +1063,35 @@ int a,b;
 	}
 }
 
-int _invsf(a)
-int a;
+int _invsf(int a)
 {
-	if ( !a )
+	if ( !a ) {
 		error("_invsf : division by 0");
-	else {
+		/* NOTREACHED */
+		return -1;
+	} else if ( !current_gfs_ntoi ) {
+		a = invm(IFTOF(a),current_gfs_q);
+		return FTOIF(a);
+	} else {
 		a = current_gfs_q1 - IFTOF(a);
 		return FTOIF(a);
 	}
 }
 
-int _divsf(a,b)
-int a,b;
+int _divsf(int a,int b)
 {
-	if ( !b )
+	int c;
+
+	if ( !b ) {
 		error("_divsf : division by 0");
-	else if ( !a )
+		/* NOTREACHED */
+		return -1;
+	} else if ( !current_gfs_ntoi ) {
+		b = invm(IFTOF(b),current_gfs_q);
+		a = IFTOF(a);
+		DMAR(a,b,0,current_gfs_q,c);
+		return FTOIF(c);
+	} else if ( !a )
 		return 0;
 	else {
 		a = IFTOF(a) - IFTOF(b);
@@ -1025,8 +1101,7 @@ int a,b;
 	}
 }
 
-int _pwrsf(a,b)
-int a,b;
+int _pwrsf(int a,int b)
 {
 	GFS at,ct;
 	Q bt;
@@ -1036,9 +1111,11 @@ int a,b;
 		return _onesf();
 	else if ( !a )
 		return 0;
-	else {
-		a = IFTOF(a);
-		MKGFS(a,at);
+	if ( !current_gfs_ntoi ) {
+		a = pwrm(current_gfs_q,IFTOF(a),b);
+		return FTOIF(a);
+	} else {
+		iftogfs(a,&at);
 		STOQ(b,bt);
 		pwrgfs(at,bt,&ct);
 		c = CONT(ct);
@@ -1048,39 +1125,42 @@ int a,b;
 
 int _onesf()
 {
-	return FTOIF(0);
+	return !current_gfs_ntoi ? FTOIF(1) : FTOIF(0);
 }
 
-int _itosf(n)
-int n;
+int _itosf(int n)
 {
 	int i;
 
+	/* XXX */
+#if 0
 	n %= current_gfs_p;
+#else
+	n %= current_gfs_q;
+#endif
 	if ( !n )
 		return 0;
-	i = current_gfs_ntoi[n];
+	i = !current_gfs_ntoi ? n : current_gfs_ntoi[n];
 	i = FTOIF(i);
 	if ( n < 0 )
 		i = _chsgnsf(i);
 	return i;
 }
 
-int _isonesf(a)
-int a;
+int _isonesf(int a)
 {
-	return a == FTOIF(0);
+	return a == _onesf();
 }
 
 int _randomsf()
 {
 	int t;
 
-	t = (int) (mt_genrand() % current_gfs_q1);
-	if ( !t )
-		return 0;
+	t = (int) (mt_genrand() % current_gfs_q);
+	if ( !current_gfs_ntoi )
+		return t ? FTOIF(t) : 0;
 	else
-		return FTOIF(t);
+		return t!=current_gfs_q1 ? FTOIF(t) : 0;
 }
 
 int field_order_sf()
@@ -1095,5 +1175,8 @@ int characteristic_sf()
 
 int extdeg_sf()
 {
-	return UDEG(current_gfs_ext);
+	if ( !current_gfs_ext )
+		return 1;
+	else
+		return UDEG(current_gfs_ext);
 }