===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/poly.c,v
retrieving revision 1.10
retrieving revision 1.26
diff -u -p -r1.10 -r1.26
--- OpenXM_contrib2/asir2000/builtin/poly.c	2001/05/09 01:41:41	1.10
+++ OpenXM_contrib2/asir2000/builtin/poly.c	2017/02/28 07:06:28	1.26
@@ -45,7 +45,7 @@
  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
  *
- * $OpenXM: OpenXM_contrib2/asir2000/builtin/poly.c,v 1.9 2001/05/02 09:03:52 noro Exp $ 
+ * $OpenXM: OpenXM_contrib2/asir2000/builtin/poly.c,v 1.25 2016/03/31 05:30:32 noro Exp $ 
 */
 #include "ca.h"
 #include "parse.h"
@@ -53,9 +53,11 @@
 
 void Pranp();
 
+void Pheadsgn();
+void Pmul_trunc(),Pquo_trunc();
 void Pumul(),Pumul_ff(),Pusquare(),Pusquare_ff(),Putmul(),Putmul_ff();
 void Pkmul(),Pksquare(),Pktmul();
-void Pord(), Pcoef0(), Pcoef(), Pdeg(), Pmindeg(), Psetmod();
+void Pord(), Premove_vars(), Pcoef0(), Pcoef(), Pdeg(), Pmindeg(), Psetmod();
 void Pcoef_gf2n();
 void getcoef(), getdeglist(), mergedeglist(), change_mvar(), restore_mvar();
 
@@ -63,7 +65,8 @@ void Pp_mag(),Pmaxblen();
 void Pmergelist(), Pch_mv(), Pre_mv(), Pdeglist();
 void Pptomp(),Pmptop();
 void Pptolmp(),Plmptop();
-void Psfptop();
+void Pptosfp(),Psfptop(),Psf_galois_action(),Psf_embed(),Psf_find_root();
+void Psf_minipoly(),Psf_log(),Psfptopsfp();
 void Pptogf2n(),Pgf2ntop(),Pgf2ntovect();
 void Pptogfpn(),Pgfpntop();
 void Pfind_root_gf2n();
@@ -107,13 +110,12 @@ void simp_ff(Obj,Obj *);
 void ranp(int,UP *);
 void field_order_ff(N *);
 
-extern int current_mod;
-extern GEN_UP2 current_mod_gf2n;
-extern int lm_lazy;
-
 int current_ff;
 
 struct ftab poly_tab[] = {
+	{"headsgn",Pheadsgn,1},
+	{"quo_trunc",Pquo_trunc,2},
+	{"mul_trunc",Pmul_trunc,3},
 	{"homogeneous_deg",Phomogeneous_deg,-2},
 	{"homogeneous_part",Phomogeneous_part,-3},
 	{"reorder",Preorder,3},
@@ -122,6 +124,8 @@ struct ftab poly_tab[] = {
 	{"p_mag",Pp_mag,1},
 	{"maxblen",Pmaxblen,1},
 	{"ord",Pord,-1},
+	{"remove_vars",Premove_vars,1},
+	{"delete_vars",Premove_vars,1},
 	{"coef0",Pcoef0,-3},
 	{"coef",Pcoef,-3},
 	{"coef_gf2n",Pcoef_gf2n,2},
@@ -131,7 +135,7 @@ struct ftab poly_tab[] = {
 
 	{"sparsemod_gf2n",Psparsemod_gf2n,-1},
 
-	{"setmod_ff",Psetmod_ff,-2},
+	{"setmod_ff",Psetmod_ff,-3},
 	{"simp_ff",Psimp_ff,1},
 	{"extdeg_ff",Pextdeg_ff,0},
 	{"characteristic_ff",Pcharacteristic_ff,0},
@@ -150,8 +154,15 @@ struct ftab poly_tab[] = {
 	{"ptolmp",Pptolmp,1},
 	{"lmptop",Plmptop,1},
 
-	{"sfptop",Psfptop,1},
+	{"sf_galois_action",Psf_galois_action,2},
+	{"sf_find_root",Psf_find_root,1},
+	{"sf_minipoly",Psf_minipoly,2},
+	{"sf_embed",Psf_embed,3},
+	{"sf_log",Psf_log,1},
 
+	{"ptosfp",Pptosfp,1},
+	{"sfptop",Psfptop,1},
+	{"sfptopsfp",Psfptopsfp,2},
 	{"ptogf2n",Pptogf2n,1},
 	{"gf2ntop",Pgf2ntop,-2},
 	{"gf2ntovect",Pgf2ntovect,1},
@@ -219,12 +230,74 @@ struct ftab poly_tab[] = {
 	{0,0,0},
 };
 
-extern V up_var;
+void Pheadsgn(NODE arg,Q *rp)
+{
+	int s;
 
-void Phomogeneous_part(arg,rp)
-NODE arg;
-P *rp;
+	s = headsgn((P)ARG0(arg));
+	STOQ(s,*rp);
+}
+
+void Pmul_trunc(NODE arg,P *rp)
 {
+	P p1,p2,p,h;
+	VL vl0,vl1,vl2,tvl,vl;
+	VN vn;
+	int i,n;
+
+	p1 = (P)ARG0(arg);
+	p2 = (P)ARG1(arg);
+	get_vars((Obj)p1,&vl1); get_vars((Obj)p2,&vl2); mergev(CO,vl1,vl2,&tvl);
+	p = (P)ARG2(arg);
+	get_vars((Obj)p,&vl0); mergev(CO,tvl,vl0,&vl);
+	for ( tvl = vl, n = 0; tvl; tvl = NEXT(tvl), n++ );
+	vn = (VN) ALLOCA((n+1)*sizeof(struct oVN));
+	for ( i = 0, tvl = vl; i < n; tvl = NEXT(tvl), i++ ) {
+		vn[i].v = tvl->v;
+		vn[i].n = 0;
+	}
+	vn[i].v = 0;
+	vn[i].n = 0;
+	for ( h = p, i = 0; OID(h) == O_P; h = COEF(DC(h)) ) {
+		for ( ; vn[i].v != VR(h); i++ );
+		vn[i].n = QTOS(DEG(DC(h)));
+	}
+	mulp_trunc(vl,p1,p2,vn,rp);
+}
+
+void Pquo_trunc(NODE arg,P *rp)
+{
+	P p1,p2,p,h;
+	VL vl0,vl1,vl2,tvl,vl;
+	VN vn;
+	int i,n;
+
+	p1 = (P)ARG0(arg);
+	p2 = (P)ARG1(arg);
+	if ( !p1 )
+		*rp = 0;
+	else if ( NUM(p2) )
+		divsp(CO,p1,p2,rp);
+	else {
+		get_vars((Obj)p1,&vl1); get_vars((Obj)p2,&vl2); mergev(CO,vl1,vl2,&vl);
+		for ( tvl = vl, n = 0; tvl; tvl = NEXT(tvl), n++ );
+		vn = (VN) ALLOCA((n+1)*sizeof(struct oVN));
+		for ( i = 0, tvl = vl; i < n; tvl = NEXT(tvl), i++ ) {
+			vn[i].v = tvl->v;
+			vn[i].n = 0;
+		}
+		vn[i].v = 0;
+		vn[i].n = 0;
+		for ( h = p2, i = 0; OID(h) == O_P; h = COEF(DC(h)) ) {
+			for ( ; vn[i].v != VR(h); i++ );
+			vn[i].n = QTOS(DEG(DC(h)));
+		}
+		quop_trunc(vl,p1,p2,vn,rp);
+	}
+}
+
+void Phomogeneous_part(NODE arg,P *rp)
+{
 	if ( argc(arg) == 2 )
 		exthp(CO,(P)ARG0(arg),QTOS((Q)ARG1(arg)),rp);
 	else
@@ -232,9 +305,7 @@ P *rp;
 			VR((P)ARG1(arg)),rp);
 }
 
-void Phomogeneous_deg(arg,rp)
-NODE arg;
-Q *rp;
+void Phomogeneous_deg(NODE arg,Q *rp)
 {
 	int d;
 
@@ -249,9 +320,7 @@ Q *rp;
 	p1 = reorder(p,ovl,nvl) => p1 is 'sorted accoding to nvl.
 */
 
-void Preorder(arg,rp)
-NODE arg;
-P *rp;
+void Preorder(NODE arg,P *rp)
 {
 	VL ovl,nvl,tvl;
 	NODE n;
@@ -284,9 +353,7 @@ P *rp;
 		return uadj_coef(CN,M,M2)*V^N+...+uadj_coef(C0,M,M2);
 */
 
-void Puadj_coef(arg,rp)
-NODE arg;
-P *rp;
+void Puadj_coef(NODE arg,P *rp)
 {
 	UP f,r;
 	N m,m2;
@@ -305,9 +372,7 @@ P *rp;
 	return [Index,Mod] or 0 (not exist)
 */
 
-void Pget_next_fft_prime(arg,rp)
-NODE arg;
-LIST *rp;
+void Pget_next_fft_prime(NODE arg,LIST *rp)
 {
 	unsigned int mod,d;
 	int start,bits,i;
@@ -321,7 +386,7 @@ LIST *rp;
 		if ( !mod ) {
 			*rp = 0; return;
 		}
-		if ( bits <= d ) {
+		if ( bits <= (int)d ) {
 			UTOQ(mod,q);
 			UTOQ(i,ind);
 			n = mknode(2,ind,q);
@@ -331,9 +396,7 @@ LIST *rp;
 	}
 }
 
-void Pranp(arg,rp)
-NODE arg;
-P *rp;
+void Pranp(NODE arg,P *rp)
 {
 	int n;
 	UP c;
@@ -347,9 +410,7 @@ P *rp;
 		*rp = 0;
 }
 
-void ranp(n,nr)
-int n;
-UP *nr;
+void ranp(int n,UP *nr)
 {
 	int i;
 	unsigned int r;
@@ -369,29 +430,26 @@ UP *nr;
 		*nr = 0;
 }
 
-void Pmaxblen(arg,rp)
-NODE arg;
-Q *rp;
+void Pmaxblen(NODE arg,Q *rp)
 {
 	int l;
 	l = maxblenp(ARG0(arg));
 	STOQ(l,*rp);
 }
 
-void Pp_mag(arg,rp)
-NODE arg;
-Q *rp;
+void Pp_mag(NODE arg,Q *rp)
 {
 	int l;
 	l = p_mag(ARG0(arg));
 	STOQ(l,*rp);
 }
 
-void Pord(arg,listp)
-NODE arg;
-LIST *listp;
+void Pord(NODE arg,LIST *listp)
 {
-	NODE n,tn;
+	NODE n,tn,p,opt;
+	char *key;
+	Obj value;
+	int overwrite=0;
 	LIST l;
 	VL vl,tvl,svl;
 	P t;
@@ -399,6 +457,18 @@ LIST *listp;
 	V *va;
 	V v;
 
+	if ( current_option ) {
+		for ( opt = current_option; opt; opt = NEXT(opt) ) {
+			p = BDY((LIST)BDY(opt));
+			key = BDY((STRING)BDY(p));
+			value = (Obj)BDY(NEXT(p));
+			if ( !strcmp(key,"overwrite") && value ) {
+				overwrite = value ? 1 : 0;
+				break;
+			}
+		}
+	}
+
 	if ( argc(arg) ) {
 		asir_assert(ARG0(arg),O_LIST,"ord");
 		for ( vl = 0, i = 0, n = BDY((LIST)ARG0(arg));
@@ -412,22 +482,29 @@ LIST *listp;
 				error("ord : invalid argument");
 			VR(tvl) = VR(t);
 		}
-		va = (V *)ALLOCA(i*sizeof(V));
-		for ( j = 0, svl = vl; j < i; j++, svl = NEXT(svl) )
-			va[j] = VR(svl);
-		for ( svl = CO; svl; svl = NEXT(svl) ) {
-			v = VR(svl);
-			for ( j = 0; j < i; j++ )
-				if ( v == va[j] )
-					break;
-			if ( j == i ) {
-				if ( !vl ) {
-					NEWVL(vl); tvl = vl;
-				} else {
-					NEWVL(NEXT(tvl)); tvl = NEXT(tvl);
+		if ( !overwrite ) {
+			va = (V *)ALLOCA(i*sizeof(V));
+			for ( j = 0, svl = vl; j < i; j++, svl = NEXT(svl) )
+				va[j] = VR(svl);
+			for ( svl = CO; svl; svl = NEXT(svl) ) {
+				v = VR(svl);
+				for ( j = 0; j < i; j++ )
+					if ( v == va[j] )
+						break;
+				if ( j == i ) {
+					if ( !vl ) {
+						NEWVL(vl); tvl = vl;
+					} else {
+						NEWVL(NEXT(tvl)); tvl = NEXT(tvl);
+					}
+					VR(tvl) = v;
 				}
-				VR(tvl) = v;
 			}
+		} else {
+			for ( svl = vl; svl; svl = NEXT(svl) ) {
+				if ( svl->v->attr == (pointer)V_PF )
+					((PFINS)svl->v->priv)->pf->ins = 0;
+			}
 		}
 		if ( vl )
 			NEXT(tvl) = 0;
@@ -439,10 +516,42 @@ LIST *listp;
 	NEXT(tn) = 0; MKLIST(l,n); *listp = l;
 }
 
-void Pcoef0(arg,rp)
-NODE arg;
-Obj *rp;
+void Premove_vars(NODE arg,LIST *listp)
 {
+	NODE l,nd,tnd;
+	V *v,*va;
+	int n,na,i,j;
+	VL vl,vl1;
+	P t;
+	LIST list;
+
+	asir_assert(ARG0(arg),O_LIST,"remove_vars");
+	l = BDY((LIST)ARG0(arg)); n = length(l);
+	v = (V *)ALLOCA(n*sizeof(V));
+	for ( i = 0; i < n; i++, l = NEXT(l) ) 
+		if ( !(t = (P)BDY(l)) || (OID(t) != O_P) )
+			error("ord : invalid argument");
+		else v[i] = VR(t);
+
+	for ( na = 0, vl = CO; vl; vl = NEXT(vl), na++ );
+	va = (V *)ALLOCA(na*sizeof(V));
+	for ( i = 0, vl = CO; i < na; i++, vl = NEXT(vl) ) va[i] = VR(vl);
+	for ( i = 0; i < na; i++ )
+		for ( j = 0; j < n; j++ ) if ( va[i] == v[j] ) va[i] = 0;
+	for ( vl = 0, i = na-1; i >= 0; i-- )
+		if ( va[i] ) {
+			NEWVL(vl1); VR(vl1) = va[i]; NEXT(vl1) = vl; vl = vl1;
+		}
+	CO = vl;
+	for ( nd = 0, vl = CO; vl; vl = NEXT(vl) ) {
+		NEXTNODE(nd,tnd); MKV(VR(vl),t); BDY(tnd) = (pointer)t;
+	}
+	if ( nd ) NEXT(tnd) = 0;
+	MKLIST(list,nd); *listp = list;
+}
+
+void Pcoef0(NODE arg,Obj *rp)
+{
 	Obj t,n;
 	P s;
 	DCP dc;
@@ -482,9 +591,7 @@ Obj *rp;
 	}
 }
 
-void Pcoef(arg,rp)
-NODE arg;
-Obj *rp;
+void Pcoef(NODE arg,Obj *rp)
 {
 	Obj t,n;
 	P s;
@@ -524,9 +631,7 @@ Obj *rp;
 	}
 }
 
-void Pcoef_gf2n(arg,rp)
-NODE arg;
-Obj *rp;
+void Pcoef_gf2n(NODE arg,Obj *rp)
 {
 	Obj t,n;
 	int id,d;
@@ -547,9 +652,7 @@ Obj *rp;
 		*rp = 0;
 }
 
-void Pdeg(arg,rp)
-NODE arg;
-Q *rp;
+void Pdeg(NODE arg,Q *rp)
 {
 	Obj t,v;
 	int d;
@@ -575,9 +678,7 @@ Q *rp;
 		degp(VR((P)ARG1(arg)),(P)ARG0(arg),rp);
 }
 
-void Pmindeg(arg,rp)
-NODE arg;
-Q *rp;
+void Pmindeg(NODE arg,Q *rp)
 {
 	Obj t;
 
@@ -587,9 +688,7 @@ Q *rp;
 		getmindeg(VR((P)ARG1(arg)),(P)ARG0(arg),rp);
 }
 
-void Psetmod(arg,rp)
-NODE arg;
-Q *rp;
+void Psetmod(NODE arg,Q *rp)
 {
 	if ( arg ) {
 		asir_assert(ARG0(arg),O_N,"setmod");
@@ -598,9 +697,7 @@ Q *rp;
 	STOQ(current_mod,*rp);
 }
 
-void Psparsemod_gf2n(arg,rp)
-NODE arg;
-Q *rp;
+void Psparsemod_gf2n(NODE arg,Q *rp)
 {
 	int id;
 
@@ -613,9 +710,7 @@ Q *rp;
 	STOQ(id,*rp);
 }
 
-void Pmultest_gf2n(arg,rp)
-NODE arg;
-GF2N *rp;
+void Pmultest_gf2n(NODE arg,GF2N *rp)
 {
 	GF2N a,b,c;
 	int i;
@@ -626,9 +721,7 @@ GF2N *rp;
 	*rp = c;
 }
 
-void Psquaretest_gf2n(arg,rp)
-NODE arg;
-GF2N *rp;
+void Psquaretest_gf2n(NODE arg,GF2N *rp)
 {
 	GF2N a,c;
 	int i;
@@ -639,9 +732,7 @@ GF2N *rp;
 	*rp = c;
 }
 
-void Pinvtest_gf2n(arg,rp)
-NODE arg;
-GF2N *rp;
+void Pinvtest_gf2n(NODE arg,GF2N *rp)
 {
 	GF2N a,c;
 	int i;
@@ -652,9 +743,7 @@ GF2N *rp;
 	*rp = c;
 }
 
-void Pbininv_gf2n(arg,rp)
-NODE arg;
-GF2N *rp;
+void Pbininv_gf2n(NODE arg,GF2N *rp)
 {
 	UP2 a,inv;
 	int n;
@@ -665,8 +754,7 @@ GF2N *rp;
 	MKGF2N(inv,*rp);
 }
 
-void Prinvtest_gf2n(rp)
-Real *rp;
+void Prinvtest_gf2n(Real *rp)
 {
 	GF2N *a;
 	GF2N c;
@@ -686,9 +774,7 @@ Real *rp;
 	MKReal(r,*rp);
 }
 
-void Pfind_root_gf2n(arg,rp)
-NODE arg;
-GF2N *rp;
+void Pfind_root_gf2n(NODE arg,GF2N *rp)
 {
 
 #if 0
@@ -704,9 +790,7 @@ GF2N *rp;
 #endif
 }
 
-void Pis_irred_gf2(arg,rp)
-NODE arg;
-Q *rp;
+void Pis_irred_gf2(NODE arg,Q *rp)
 {
 	UP2 t;
 
@@ -714,9 +798,7 @@ Q *rp;
 	*rp = irredcheckup2(t) ? ONE : 0;
 }
 
-void Pis_irred_ddd_gf2(arg,rp)
-NODE arg;
-Q *rp;
+void Pis_irred_ddd_gf2(NODE arg,Q *rp)
 {
 	UP2 t;
 	int ret;
@@ -726,17 +808,7 @@ Q *rp;
 	STOQ(ret,*rp);
 }
 
-extern P current_gfs_ext;
-extern int current_gfs_p;
-extern int current_gfs_q;
-extern int current_gfs_q1;
-extern int *current_gfs_plus1;
-extern int *current_gfs_ntoi;
-extern int *current_gfs_iton;
-
-void Psetmod_ff(arg,rp)
-NODE arg;
-Obj *rp;
+void Psetmod_ff(NODE arg,Obj *rp)
 {
 	int ac;
 	int d;
@@ -744,8 +816,9 @@ Obj *rp;
 	N n;
 	UP up;
 	UP2 up2;
+	UM dp;
 	Q q,r;
-	P p;
+	P p,p1,y;
 	NODE n0,n1;
 	LIST list;
 
@@ -753,8 +826,9 @@ Obj *rp;
 	if ( ac == 1 ) {
 		mod = (Obj)ARG0(arg);
 		if ( !mod )
-			error("setmod_ff : invalid argument");
-		switch ( OID(mod) ) {
+            current_ff = FF_NOT_SET;
+        else {
+		  switch ( OID(mod) ) {
 			case O_N:
 				current_ff = FF_GFP;
 				setmod_lm(NM((Q)mod));
@@ -764,7 +838,8 @@ Obj *rp;
 				setmod_gf2n((P)mod); break;
 			default:
 				error("setmod_ff : invalid argument");
-		}
+		  }
+      }
 	} else if ( ac == 2 ) {
 		if ( OID(ARG0(arg)) == O_N ) {
 			/* small finite field; primitive root representation */
@@ -779,6 +854,14 @@ Obj *rp;
 			setmod_lm(NM((Q)mod));
 			setmod_gfpn((P)defpoly);
 		}
+	} else if ( ac == 3 ) {
+		/* finite extension of a small finite field */
+		current_ff = FF_GFS;
+		setmod_sf(QTOS((Q)ARG0(arg)),QTOS((Q)ARG1(arg)));
+		d = QTOS((Q)ARG2(arg));
+		generate_defpoly_sfum(d,&dp);
+		setmod_gfsn(dp);
+		current_ff = FF_GFSN;
 	}
 	switch ( current_ff ) {
 		case FF_GFP:
@@ -792,15 +875,31 @@ Obj *rp;
 			MKLIST(list,n0);
 			*rp = (Obj)list; break;
 		case FF_GFS:
+		case FF_GFSN:
 			STOQ(current_gfs_p,q);
-			if ( current_gfs_ext ) {
+			if ( current_gfs_ext )
 				enc_to_p(current_gfs_p,current_gfs_iton[1],
 					VR(current_gfs_ext),&p);
-				n0 = mknode(3,q,current_gfs_ext,p);
-			} else {
-				STOQ(current_gfs_iton[1],r);
-				n0 = mknode(3,q,current_gfs_ext,r);
+			else {
+				if ( current_gfs_p == 2 )
+					r = ONE;
+				else if ( !current_gfs_ntoi )
+					r = 0;
+				else
+					STOQ(current_gfs_iton[1],r);
+				p = (P)r;
 			}
+			switch ( current_ff ) {
+				case FF_GFS:
+					n0 = mknode(3,q,current_gfs_ext,p);
+					break;
+				case FF_GFSN:
+					getmod_gfsn(&dp);
+					makevar("y",&y);
+					sfumtop(VR(y),dp,&p1);
+					n0 = mknode(4,q,current_gfs_ext,p,p1);
+					break;
+			}
 			MKLIST(list,n0);
 			*rp = (Obj)list; break;
 		default:
@@ -808,12 +907,12 @@ Obj *rp;
 	}
 }
 
-void Pextdeg_ff(rp)
-Q *rp;
+void Pextdeg_ff(Q *rp)
 {
 	int d;
 	UP2 up2;
 	UP up;
+	UM dp;
 
 	switch ( current_ff ) {
 		case FF_GFP:
@@ -828,13 +927,16 @@ Q *rp;
 			else
 				*rp = DEG(DC(current_gfs_ext));
 			break;
+		case FF_GFSN:
+			getmod_gfsn(&dp);
+			STOQ(DEG(dp),*rp);
+			break;
 		default:
 			error("extdeg_ff : current_ff is not set");
 	}
 }
 
-void Pcharacteristic_ff(rp)
-Q *rp;
+void Pcharacteristic_ff(Q *rp)
 {
 	N lm;
 
@@ -845,20 +947,19 @@ Q *rp;
 		case FF_GF2N:
 			STOQ(2,*rp); break;
 		case FF_GFS:
+		case FF_GFSN:
 			STOQ(current_gfs_p,*rp); break;
 		default:
 			error("characteristic_ff : current_ff is not set");
 	}
 }
 
-void Pfield_type_ff(rp)
-Q *rp;
+void Pfield_type_ff(Q *rp)
 {
 	STOQ(current_ff,*rp);
 }
 
-void Pfield_order_ff(rp)
-Q *rp;
+void Pfield_order_ff(Q *rp)
 {
 	N n;
 
@@ -866,42 +967,13 @@ Q *rp;
 	NTOQ(n,1,*rp);
 }
 
-void field_order_ff(order)
-N *order;
+void Prandom_ff(Obj *rp)
 {
-	UP2 up2;
-	UP up;
-	N m;
-	int d,w;
-
-	switch ( current_ff ) {
-		case FF_GFP:
-			getmod_lm(order); break;
-		case FF_GF2N:
-			getmod_gf2n(&up2); d = degup2(up2);
-			w = (d>>5)+1;
-			*order = m = NALLOC(w);
-			PL(m)=w;
-			bzero((char *)BD(m),w*sizeof(unsigned int));
-			BD(m)[d>>5] |= 1<<(d&31);
-			break;
-		case FF_GFPN:
-			getmod_lm(&m);
-			getmod_gfpn(&up); pwrn(m,up->d,order); break;
-		case FF_GFS:
-			STON(current_gfs_q,*order); break;
-		default:
-			error("field_order_ff : current_ff is not set");
-	}
-}
-
-void Prandom_ff(rp)
-Obj *rp;
-{
 	LM l;
 	GF2N g;
 	GFPN p;
 	GFS s;
+	GFSN spn;
 
 	switch ( current_ff ) {
 		case FF_GFP:
@@ -912,84 +984,20 @@ Obj *rp;
 			randomgfpn(&p); *rp = (Obj)p; break;
 		case FF_GFS:
 			randomgfs(&s); *rp = (Obj)s; break;
+		case FF_GFSN:
+			randomgfsn(&spn); *rp = (Obj)spn; break;
 		default:
 			error("random_ff : current_ff is not set");
 	}
 }
 
-void Psimp_ff(arg,rp)
-NODE arg;
-Obj *rp;
+void Psimp_ff(NODE arg,Obj *rp)
 {
-	LM r;
-	GF2N rg;
-	extern lm_lazy;
-
 	simp_ff((Obj)ARG0(arg),rp);
 }
 
-void simp_ff(p,rp)
-Obj p;
-Obj *rp;
+void getcoef(VL vl,P p,V v,Q d,P *r)
 {
-	Num n;
-	LM r,s;
-	DCP dc,dcr0,dcr;
-	GF2N rg,sg;
-	GFPN rpn,spn;
-	GFS rs;
-	P t;
-	Obj obj;
-
-	lm_lazy = 0;
-	if ( !p )
-		*rp = 0;
-	else if ( OID(p) == O_N ) {
-		switch ( current_ff ) {
-			case FF_GFP:
-				ptolmp((P)p,&t); simplm((LM)t,&s); *rp = (Obj)s;
-				break;
-			case FF_GF2N:
-				ptogf2n((Obj)p,&rg); simpgf2n((GF2N)rg,&sg); *rp = (Obj)sg;
-				break;
-			case FF_GFPN:
-				ntogfpn((Obj)p,&rpn); simpgfpn((GFPN)rpn,&spn); *rp = (Obj)spn;
-				break;
-			case FF_GFS:
-				if ( NID((Num)p) == N_GFS )
-					*rp = p;
-				else {
-					ptomp(current_gfs_p,(P)p,&t); mqtogfs(t,&rs);
-					*rp = (Obj)rs;
-				} 
-				break;
-			default:
-				*rp = (Obj)p;
-				break;
-		}
-	} else if ( OID(p) == O_P ) {
-		for ( dc = DC((P)p), dcr0 = 0; dc; dc = NEXT(dc) ) {
-			simp_ff((Obj)COEF(dc),&obj);
-			if ( obj ) {
-				NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = (P)obj;
-			}
-		}
-		if ( !dcr0 ) 
-			*rp = 0;
-		else {
-			NEXT(dcr) = 0; MKP(VR((P)p),dcr0,t); *rp = (Obj)t;
-		}
-	} else
-		error("simp_ff : not implemented yet");
-}
-
-void getcoef(vl,p,v,d,r)
-VL vl;
-P p;
-V v;
-Q d;
-P *r;
-{
 	P s,t,u,a,b,x;
 	DCP dc;
 	V w;
@@ -1014,9 +1022,7 @@ P *r;
 	}
 }
 
-void Pdeglist(arg,rp)
-NODE arg;
-LIST *rp;
+void Pdeglist(NODE arg,LIST *rp)
 {
 	NODE d;
 
@@ -1024,25 +1030,17 @@ LIST *rp;
 	MKLIST(*rp,d);
 }
 
-void Pch_mv(arg,rp)
-NODE arg;
-P *rp;
+void Pch_mv(NODE arg,P *rp)
 {
 	change_mvar(CO,(P)ARG0(arg),VR((P)ARG1(arg)),rp);
 }
 
-void Pre_mv(arg,rp)
-NODE arg;
-P *rp;
+void Pre_mv(NODE arg,P *rp)
 {
 	restore_mvar(CO,(P)ARG0(arg),VR((P)ARG1(arg)),rp);
 }
 
-void change_mvar(vl,p,v,r)
-VL vl;
-P p;
-V v;
-P *r;
+void change_mvar(VL vl,P p,V v,P *r)
 {
 	Q d;
 	DCP dc,dc0;
@@ -1060,11 +1058,7 @@ P *r;
 	}
 }
 
-void restore_mvar(vl,p,v,r)
-VL vl;
-P p;
-V v;
-P *r;
+void restore_mvar(VL vl,P p,V v,P *r)
 {
 	P s,u,a,b,x;
 	DCP dc;
@@ -1081,10 +1075,7 @@ P *r;
 	}
 }
 
-void getdeglist(p,v,d)
-P p;
-V v;
-NODE *d;
+void getdeglist(P p,V v,NODE *d)
 {
 	NODE n,n0,d0,d1,d2;
 	DCP dc;
@@ -1103,9 +1094,8 @@ NODE *d;
 		*d = d0;
 	}
 }
-void Pmergelist(arg,rp)
-NODE arg;
-LIST *rp;
+
+void Pmergelist(NODE arg,LIST *rp)
 {
     NODE n;
 
@@ -1115,8 +1105,7 @@ LIST *rp;
 	MKLIST(*rp,n);
 }
 
-void mergedeglist(d0,d1,dr)
-NODE d0,d1,*dr;
+void mergedeglist(NODE d0,NODE d1,NODE *dr)
 {
 	NODE t0,t,dt;
 	Q d;
@@ -1148,87 +1137,138 @@ NODE d0,d1,*dr;
 	}
 }
 
-void Pptomp(arg,rp)
-NODE arg;
-P *rp;
+void Pptomp(NODE arg,P *rp)
 {
 	ptomp(QTOS((Q)ARG1(arg)),(P)ARG0(arg),rp);
 }
 
-void Pmptop(arg,rp)
-NODE arg;
-P *rp;
+void Pmptop(NODE arg,P *rp)
 {
 	mptop((P)ARG0(arg),rp);
 }
 
-void Pptolmp(arg,rp)
-NODE arg;
-P *rp;
+void Pptolmp(NODE arg,P *rp)
 {
 	ptolmp((P)ARG0(arg),rp);
 }
 
-void Plmptop(arg,rp)
-NODE arg;
-P *rp;
+void Plmptop(NODE arg,P *rp)
 {
 	lmptop((P)ARG0(arg),rp);
 }
 
-void Psfptop(arg,rp)
-NODE arg;
-P *rp;
+void Psf_galois_action(NODE arg,P *rp)
 {
+	sf_galois_action(ARG0(arg),ARG1(arg),rp);
+}
+
+/*
+  sf_embed(F,B,PM)
+  F : an element of GF(pn)
+  B : the image of the primitive root of GF(pn)
+  PM : order of GF(pm)
+*/
+
+void Psf_embed(NODE arg,P *rp)
+{
+	int k,pm;
+
+	/* GF(pn)={0,1,a,a^2,...}->GF(pm)={0,1,b,b^2,...}; a->b^k */
+	k = CONT((GFS)ARG1(arg));
+	pm = QTOS((Q)ARG2(arg));
+	sf_embed((P)ARG0(arg),k,pm,rp);
+}
+
+void Psf_log(NODE arg,Q *rp)
+{
+	int k;
+
+	if ( !ARG0(arg) )
+		error("sf_log : invalid armument");
+	k = CONT((GFS)ARG0(arg));
+	STOQ(k,*rp);
+}
+
+void Psf_find_root(NODE arg,GFS *rp)
+{
+	P p;
+	Obj t;
+	int d;
+	UM u;
+	int *root;
+
+	p = (P)ARG0(arg);
+	simp_ff((Obj)p,&t); p = (P)t;
+	d = getdeg(VR(p),p);
+	u = W_UMALLOC(d);
+	ptosfum(p,u);
+	root = (int *)ALLOCA(d*sizeof(int));
+	find_rootsf(u,root);
+	MKGFS(IFTOF(root[0]),*rp);
+}
+
+void Psf_minipoly(NODE arg,P *rp)
+{
+	Obj t;
+	P p1,p2;
+	int d1,d2;
+	UM up1,up2,m;
+
+	p1 = (P)ARG0(arg); simp_ff((Obj)p1,&t); p1 = (P)t;
+	p2 = (P)ARG1(arg); simp_ff((Obj)p2,&t); p2 = (P)t;
+	d1 = getdeg(VR(p1),p1); up1 = W_UMALLOC(d1); ptosfum(p1,up1);
+	d2 = getdeg(VR(p2),p2); up2 = W_UMALLOC(d2); ptosfum(p2,up2);
+	m = W_UMALLOC(d2);
+	minipolysf(up1,up2,m);
+	sfumtop(VR(p2),m,&p1);
+	sfptop(p1,rp);
+}
+
+void Pptosfp(NODE arg,P *rp)
+{
+	ptosfp(ARG0(arg),rp);
+}
+
+void Psfptop(NODE arg,P *rp)
+{
 	sfptop((P)ARG0(arg),rp);
 }
 
-void Pptogf2n(arg,rp)
-NODE arg;
-GF2N *rp;
+void Psfptopsfp(NODE arg,P *rp)
 {
-	ptogf2n((Obj)ARG0(arg),rp);
+	sfptopsfp((P)ARG0(arg),VR((P)ARG1(arg)),rp);
 }
 
-void Pgf2ntop(arg,rp)
-NODE arg;
-P *rp;
+void Pptogf2n(NODE arg,GF2N *rp)
 {
-	extern V up2_var;
+	ptogf2n((Obj)ARG0(arg),rp);
+}
 
+void Pgf2ntop(NODE arg,P *rp)
+{
 	if ( argc(arg) == 2 )
 		up2_var = VR((P)ARG1(arg));
 	gf2ntop((GF2N)ARG0(arg),rp);
 }
 
-void Pgf2ntovect(arg,rp)
-NODE arg;
-VECT *rp;
+void Pgf2ntovect(NODE arg,VECT *rp)
 {
 	gf2ntovect((GF2N)ARG0(arg),rp);
 }
 
-void Pptogfpn(arg,rp)
-NODE arg;
-GF2N *rp;
+void Pptogfpn(NODE arg,GFPN *rp)
 {
 	ptogfpn((Obj)ARG0(arg),rp);
 }
 
-void Pgfpntop(arg,rp)
-NODE arg;
-P *rp;
+void Pgfpntop(NODE arg,P *rp)
 {
-	extern V up_var;
-
 	if ( argc(arg) == 2 )
 		up_var = VR((P)ARG1(arg));
 	gfpntop((GFPN)ARG0(arg),rp);
 }
 
-void Pureverse(arg,rp)
-NODE arg;
-P *rp;
+void Pureverse(NODE arg,P *rp)
 {
 	UP p,r;
 
@@ -1240,9 +1280,7 @@ P *rp;
 	uptop(r,rp);
 }
 
-void Putrunc(arg,rp)
-NODE arg;
-P *rp;
+void Putrunc(NODE arg,P *rp)
 {
 	UP p,r;
 
@@ -1251,9 +1289,7 @@ P *rp;
 	uptop(r,rp);
 }
 
-void Pudecomp(arg,rp)
-NODE arg;
-LIST *rp;
+void Pudecomp(NODE arg,LIST *rp)
 {
 	P u,l;
 	UP p,up,low;
@@ -1267,9 +1303,7 @@ LIST *rp;
 	MKLIST(*rp,n0);
 }
 
-void Purembymul(arg,rp)
-NODE arg;
-P *rp;
+void Purembymul(NODE arg,P *rp)
 {
 	UP p1,p2,r;
 
@@ -1289,9 +1323,7 @@ P *rp;
  * p2*inv = 1 mod x^d2
  */
 
-void Purembymul_precomp(arg,rp)
-NODE arg;
-P *rp;
+void Purembymul_precomp(NODE arg,P *rp)
 {
 	UP p1,p2,inv,r;
 
@@ -1311,9 +1343,7 @@ P *rp;
 	}
 }
 
-void Puinvmod(arg,rp)
-NODE arg;
-P *rp;
+void Puinvmod(NODE arg,P *rp)
 {
 	UP p,r;
 
@@ -1322,9 +1352,7 @@ P *rp;
 	uptop(r,rp);
 }
 
-void Purevinvmod(arg,rp)
-NODE arg;
-P *rp;
+void Purevinvmod(NODE arg,P *rp)
 {
 	UP p,pr,r;
 
@@ -1334,9 +1362,7 @@ P *rp;
 	uptop(r,rp);
 }
 
-void Ppwrmod_ff(arg,rp)
-NODE arg;
-P *rp;
+void Ppwrmod_ff(NODE arg,P *rp)
 {
 	UP p1,p2;
 
@@ -1348,6 +1374,7 @@ P *rp;
 			powermodup_gf2n(p1,&p2); break;
 		case FF_GFPN:
 		case FF_GFS:
+		case FF_GFSN:
 			powermodup(p1,&p2); break;
 		default:
 			error("pwrmod_ff : current_ff is not set");
@@ -1355,9 +1382,7 @@ P *rp;
 	uptop(p2,rp);
 }
 
-void Pgeneric_pwrmod_ff(arg,rp)
-NODE arg;
-P *rp;
+void Pgeneric_pwrmod_ff(NODE arg,P *rp)
 {
 	UP g,f,r;
 
@@ -1370,6 +1395,7 @@ P *rp;
 			generic_powermodup_gf2n(g,f,(Q)ARG2(arg),&r); break;
 		case FF_GFPN:
 		case FF_GFS:
+		case FF_GFSN:
 			generic_powermodup(g,f,(Q)ARG2(arg),&r); break;
 		default:
 			error("generic_pwrmod_ff : current_ff is not set");
@@ -1377,9 +1403,7 @@ P *rp;
 	uptop(r,rp);
 }
 
-void Ppwrtab_ff(arg,rp)
-NODE arg;
-VECT *rp;
+void Ppwrtab_ff(NODE arg,VECT *rp)
 {
 	UP f,xp;
 	UP *tab;
@@ -1398,6 +1422,7 @@ VECT *rp;
 			powertabup_gf2n(f,xp,tab); break;
 		case FF_GFPN:
 		case FF_GFS:
+		case FF_GFSN:
 			powertabup(f,xp,tab); break;
 		default:
 			error("pwrtab_ff : current_ff is not set");
@@ -1407,9 +1432,7 @@ VECT *rp;
 		uptop(tab[i],(P *)&BDY(r)[i]);
 }
 
-void Pkpwrmod_lm(arg,rp)
-NODE arg;
-P *rp;
+void Pkpwrmod_lm(NODE arg,P *rp)
 {
 	UP p1,p2;
 
@@ -1418,9 +1441,7 @@ P *rp;
 	uptop(p2,rp);
 }
 
-void Pkgeneric_pwrmod_lm(arg,rp)
-NODE arg;
-P *rp;
+void Pkgeneric_pwrmod_lm(NODE arg,P *rp)
 {
 	UP g,f,r;
 
@@ -1430,9 +1451,7 @@ P *rp;
 	uptop(r,rp);
 }
 
-void Pkpwrtab_lm(arg,rp)
-NODE arg;
-VECT *rp;
+void Pkpwrtab_lm(NODE arg,VECT *rp)
 {
 	UP f,xp;
 	UP *tab;
@@ -1450,17 +1469,13 @@ VECT *rp;
 		uptop(tab[i],(P *)&BDY(r)[i]);
 }
 
-void Plazy_lm(arg,rp)
-NODE arg;
-Q *rp;
+void Plazy_lm(NODE arg,Q *rp)
 {
 	lm_lazy = QTOS((Q)ARG0(arg));
 	*rp = 0;
 }
 
-void Pkmul(arg,rp)
-NODE arg;
-P *rp;
+void Pkmul(NODE arg,P *rp)
 {
 	P n1,n2;
 	
@@ -1470,9 +1485,7 @@ P *rp;
 	kmulp(CO,n1,n2,rp);
 }
 
-void Pksquare(arg,rp)
-NODE arg;
-P *rp;
+void Pksquare(NODE arg,P *rp)
 {
 	P n1;
 	
@@ -1481,9 +1494,7 @@ P *rp;
 	ksquarep(CO,n1,rp);
 }
 
-void Pktmul(arg,rp)
-NODE arg;
-P *rp;
+void Pktmul(NODE arg,P *rp)
 {
 	UP p1,p2,r;
 
@@ -1493,9 +1504,7 @@ P *rp;
 	uptop(r,rp);
 }
 
-void Pumul(arg,rp)
-NODE arg;
-P *rp;
+void Pumul(NODE arg,P *rp)
 {
 	P a1,a2;
 	UP p1,p2,r;
@@ -1504,7 +1513,7 @@ P *rp;
 	if ( !a1 || !a2 || NUM(a1) || NUM(a2) )
 		mulp(CO,a1,a2,rp);
 	else {
-		if ( !uzpcheck(a1) || !uzpcheck(a2) || VR(a1) != VR(a2) )
+		if ( !uzpcheck((Obj)a1) || !uzpcheck((Obj)a2) || VR(a1) != VR(a2) )
 			error("umul : invalid argument");
 		ptoup(a1,&p1);
 		ptoup(a2,&p2);
@@ -1513,20 +1522,16 @@ P *rp;
 	}
 }
 
-void Pusquare(arg,rp)
-NODE arg;
-P *rp;
+void Pusquare(NODE arg,P *rp)
 {
-	UP p1,p2,r;
+	UP p1,r;
 
 	ptoup((P)ARG0(arg),&p1);
 	hybrid_squareup(0,p1,&r);
 	uptop(r,rp);
 }
 
-void Putmul(arg,rp)
-NODE arg;
-P *rp;
+void Putmul(NODE arg,P *rp)
 {
 	UP p1,p2,r;
 
@@ -1536,9 +1541,7 @@ P *rp;
 	uptop(r,rp);
 }
 
-void Pumul_ff(arg,rp)
-NODE arg;
-Obj *rp;
+void Pumul_ff(NODE arg,Obj *rp)
 {
 	P a1,a2;
 	UP p1,p2,r;
@@ -1552,11 +1555,9 @@ Obj *rp;
 	simp_ff((Obj)p,rp);
 }
 
-void Pusquare_ff(arg,rp)
-NODE arg;
-Obj *rp;
+void Pusquare_ff(NODE arg,Obj *rp)
 {
-	UP p1,p2,r;
+	UP p1,r;
 	P p;
 
 	ptoup((P)ARG0(arg),&p1);
@@ -1565,9 +1566,7 @@ Obj *rp;
 	simp_ff((Obj)p,rp);
 }
 
-void Putmul_ff(arg,rp)
-NODE arg;
-Obj *rp;
+void Putmul_ff(NODE arg,Obj *rp)
 {
 	UP p1,p2,r;
 	P p;
@@ -1579,9 +1578,7 @@ Obj *rp;
 	simp_ff((Obj)p,rp);
 }
 
-void Phfmul_lm(arg,rp)
-NODE arg;
-P *rp;
+void Phfmul_lm(NODE arg,P *rp)
 {
 	UP p1,p2,r;
 	UP hi,lo,mid,t,s,p10,p11,p20,p21;
@@ -1622,9 +1619,7 @@ P *rp;
 	uptop(r,rp);
 }
 
-void Pfmultest(arg,rp)
-NODE arg;
-LIST *rp;
+void Pfmultest(NODE arg,LIST *rp)
 {
 	P p1,p2,r;
 	int d1,d2;
@@ -1669,9 +1664,7 @@ LIST *rp;
 	}
 }
 
-void Pkmulum(arg,rp)
-NODE arg;
-P *rp;
+void Pkmulum(NODE arg,P *rp)
 {
 	P p1,p2;
 	int d1,d2,mod;
@@ -1686,9 +1679,7 @@ P *rp;
 	umtop(VR(p1),wr,rp);
 }
 
-void Pksquareum(arg,rp)
-NODE arg;
-P *rp;
+void Pksquareum(NODE arg,P *rp)
 {
 	P p1;
 	int d1,mod;
@@ -1703,9 +1694,7 @@ P *rp;
 	umtop(VR(p1),wr,rp);
 }
 
-void Ptracemod_gf2n(arg,rp)
-NODE arg;
-P *rp;
+void Ptracemod_gf2n(NODE arg,P *rp)
 {
 	UP g,f,r;
 
@@ -1715,9 +1704,7 @@ P *rp;
 	uptop(r,rp);
 }
 
-void Pumul_specialmod(arg,rp)
-NODE arg;
-P *rp;
+void Pumul_specialmod(NODE arg,P *rp)
 {
 	P a1,a2;
 	UP p1,p2,r;
@@ -1729,7 +1716,7 @@ P *rp;
 	if ( !a1 || !a2 || NUM(a1) || NUM(a2) )
 		mulp(CO,a1,a2,rp);
 	else {
-		if ( !uzpcheck(a1) || !uzpcheck(a2) || VR(a1) != VR(a2) )
+		if ( !uzpcheck((Obj)a1) || !uzpcheck((Obj)a2) || VR(a1) != VR(a2) )
 			error("umul_specialmod : invalid argument");
 		ptoup(a1,&p1);
 		ptoup(a2,&p2);
@@ -1743,9 +1730,7 @@ P *rp;
 	}
 }
 
-void Pusquare_specialmod(arg,rp)
-NODE arg;
-P *rp;
+void Pusquare_specialmod(NODE arg,P *rp)
 {
 	P a1;
 	UP p1,r;
@@ -1757,7 +1742,7 @@ P *rp;
 	if ( !a1 || NUM(a1) )
 		mulp(CO,a1,a1,rp);
 	else {
-		if ( !uzpcheck(a1) )
+		if ( !uzpcheck((Obj)a1) )
 			error("usquare_specialmod : invalid argument");
 		ptoup(a1,&p1);
 		n = BDY((LIST)ARG1(arg));
@@ -1770,9 +1755,7 @@ P *rp;
 	}
 }
 
-void Putmul_specialmod(arg,rp)
-NODE arg;
-P *rp;
+void Putmul_specialmod(NODE arg,P *rp)
 {
 	P a1,a2;
 	UP p1,p2,r;
@@ -1784,7 +1767,7 @@ P *rp;
 	if ( !a1 || !a2 || NUM(a1) || NUM(a2) )
 		mulp(CO,a1,a2,rp);
 	else {
-		if ( !uzpcheck(a1) || !uzpcheck(a2) || VR(a1) != VR(a2) )
+		if ( !uzpcheck((Obj)a1) || !uzpcheck((Obj)a2) || VR(a1) != VR(a2) )
 			error("utmul_specialmod : invalid argument");
 		ptoup(a1,&p1);
 		ptoup(a2,&p2);