===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/poly.c,v
retrieving revision 1.1.1.1
retrieving revision 1.18
diff -u -p -r1.1.1.1 -r1.18
--- OpenXM_contrib2/asir2000/builtin/poly.c	1999/12/03 07:39:07	1.1.1.1
+++ OpenXM_contrib2/asir2000/builtin/poly.c	2003/01/16 00:33:27	1.18
@@ -1,4 +1,52 @@
-/* $OpenXM: OpenXM/src/asir99/builtin/poly.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
+/*
+ * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED 
+ * All rights reserved.
+ * 
+ * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
+ * non-exclusive and royalty-free license to use, copy, modify and
+ * redistribute, solely for non-commercial and non-profit purposes, the
+ * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
+ * conditions of this Agreement. For the avoidance of doubt, you acquire
+ * only a limited right to use the SOFTWARE hereunder, and FLL or any
+ * third party developer retains all rights, including but not limited to
+ * copyrights, in and to the SOFTWARE.
+ * 
+ * (1) FLL does not grant you a license in any way for commercial
+ * purposes. You may use the SOFTWARE only for non-commercial and
+ * non-profit purposes only, such as academic, research and internal
+ * business use.
+ * (2) The SOFTWARE is protected by the Copyright Law of Japan and
+ * international copyright treaties. If you make copies of the SOFTWARE,
+ * with or without modification, as permitted hereunder, you shall affix
+ * to all such copies of the SOFTWARE the above copyright notice.
+ * (3) An explicit reference to this SOFTWARE and its copyright owner
+ * shall be made on your publication or presentation in any form of the
+ * results obtained by use of the SOFTWARE.
+ * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
+ * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
+ * for such modification or the source code of the modified part of the
+ * SOFTWARE.
+ * 
+ * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
+ * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
+ * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
+ * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
+ * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
+ * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
+ * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
+ * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
+ * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
+ * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
+ * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
+ * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
+ * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
+ * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
+ * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
+ * 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.17 2002/09/27 04:42:59 noro Exp $ 
+*/
 #include "ca.h"
 #include "parse.h"
 #include "base.h"
@@ -11,13 +59,16 @@ void Pord(), Pcoef0(), Pcoef(), Pdeg(), Pmindeg(), Pse
 void Pcoef_gf2n();
 void getcoef(), getdeglist(), mergedeglist(), change_mvar(), restore_mvar();
 
-void Pp_mag();
+void Pp_mag(),Pmaxblen();
 void Pmergelist(), Pch_mv(), Pre_mv(), Pdeglist();
 void Pptomp(),Pmptop();
 void Pptolmp(),Plmptop();
+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();
+void Pumul_specialmod(),Pusquare_specialmod(),Putmul_specialmod();
 
 void Pureverse(),Putrunc(),Pudecomp(),Purembymul(),Purembymul_precomp();
 void Puinvmod(),Purevinvmod();
@@ -48,20 +99,25 @@ void Pbininv_gf2n();
 void Prinvtest_gf2n();
 void Pis_irred_gf2();
 void Pis_irred_ddd_gf2();
-
+void Pget_next_fft_prime();
+void Puadj_coef();
+void Preorder();
+void Phomogeneous_part();
+void Phomogeneous_deg();
 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[] = {
+	{"homogeneous_deg",Phomogeneous_deg,-2},
+	{"homogeneous_part",Phomogeneous_part,-3},
+	{"reorder",Preorder,3},
+	{"uadj_coef",Puadj_coef,3},
 	{"ranp",Pranp,2},
 	{"p_mag",Pp_mag,1},
+	{"maxblen",Pmaxblen,1},
 	{"ord",Pord,-1},
 	{"coef0",Pcoef0,-3},
 	{"coef",Pcoef,-3},
@@ -72,7 +128,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},
@@ -91,6 +147,15 @@ struct ftab poly_tab[] = {
 	{"ptolmp",Pptolmp,1},
 	{"lmptop",Plmptop,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},
@@ -107,6 +172,10 @@ struct ftab poly_tab[] = {
 	{"ureverse_inv_as_power_series",Purevinvmod,2},
 	{"uinv_as_power_series",Puinvmod,2},
 
+	{"umul_specialmod",Pumul_specialmod,3},
+	{"usquare_specialmod",Pusquare_specialmod,2},
+	{"utmul_specialmod",Putmul_specialmod,4},
+
 	{"utmul",Putmul,3},
 	{"umul_ff",Pumul_ff,2},
 	{"usquare_ff",Pusquare_ff,1},
@@ -118,7 +187,7 @@ struct ftab poly_tab[] = {
 
 	{"utrunc",Putrunc,2},
 	{"udecomp",Pudecomp,2},
-	{"ureverse",Pureverse,1},
+	{"ureverse",Pureverse,-2},
 	{"urembymul",Purembymul,2},
 	{"urembymul_precomp",Purembymul_precomp,3},
 
@@ -150,15 +219,112 @@ struct ftab poly_tab[] = {
 	{"bininv_gf2n",Pbininv_gf2n,2},
 	{"invtest_gf2n",Pinvtest_gf2n,1},
 	{"rinvtest_gf2n",Prinvtest_gf2n,0},
+	{"get_next_fft_prime",Pget_next_fft_prime,2},
 	{0,0,0},
 };
 
-extern V up_var;
+void Phomogeneous_part(NODE arg,P *rp)
+{
+	if ( argc(arg) == 2 )
+		exthp(CO,(P)ARG0(arg),QTOS((Q)ARG1(arg)),rp);
+	else
+		exthpc_generic(CO,(P)ARG0(arg),QTOS((Q)ARG2(arg)),
+			VR((P)ARG1(arg)),rp);
+}
 
-void Pranp(arg,rp)
-NODE arg;
-P *rp;
+void Phomogeneous_deg(NODE arg,Q *rp)
 {
+	int d;
+
+	if ( argc(arg) == 1 )
+		d = homdeg((P)ARG0(arg));
+	else
+		d = getchomdeg(VR((P)ARG1(arg)),(P)ARG0(arg));
+	STOQ(d,*rp);
+}
+
+/* 
+	p1 = reorder(p,ovl,nvl) => p1 is 'sorted accoding to nvl.
+*/
+
+void Preorder(NODE arg,P *rp)
+{
+	VL ovl,nvl,tvl;
+	NODE n;
+
+	for ( ovl = 0, n = BDY((LIST)ARG1(arg)); n; n = NEXT(n) ) {
+		if ( !ovl ) {
+			NEWVL(ovl); tvl = ovl;
+		} else {
+			NEWVL(NEXT(tvl)); tvl = NEXT(tvl);
+		}
+			VR(tvl) = VR((P)BDY(n));
+	}
+	for ( nvl = 0, n = BDY((LIST)ARG2(arg)); n; n = NEXT(n) ) {
+		if ( !nvl ) {
+			NEWVL(nvl); tvl = nvl;
+		} else {
+			NEWVL(NEXT(tvl)); tvl = NEXT(tvl);
+		}
+			VR(tvl) = VR((P)BDY(n));
+	}
+	reorderp(nvl,ovl,(P)ARG0(arg),rp);
+}
+
+/*
+	uadj_coef(F,M,M2)
+	if ( F is a non-negative integer )
+		return F > M2 ? F-M : M2;
+	else
+		F = CN*V^N+...+C0
+		return uadj_coef(CN,M,M2)*V^N+...+uadj_coef(C0,M,M2);
+*/
+
+void Puadj_coef(NODE arg,P *rp)
+{
+	UP f,r;
+	N m,m2;
+
+	ptoup((P)ARG0(arg),&f);
+	m = NM((Q)ARG1(arg));
+	m2 = NM((Q)ARG2(arg));
+	adj_coefup(f,m,m2,&r);
+	uptop(r,rp);
+}
+
+/*
+	get_next_fft_prime(StartIndex,Bits)
+	tries to find smallest Index >= StartIndex s.t.
+		2^(Bits-1)|FFTprime[Index]-1
+	return [Index,Mod] or 0 (not exist)
+*/
+
+void Pget_next_fft_prime(NODE arg,LIST *rp)
+{
+	unsigned int mod,d;
+	int start,bits,i;
+	NODE n;
+	Q q,ind;
+
+	start = QTOS((Q)ARG0(arg));
+	bits = QTOS((Q)ARG1(arg));
+	for ( i = start; ; i++ ) {
+		get_fft_prime(i,&mod,&d);
+		if ( !mod ) {
+			*rp = 0; return;
+		}
+		if ( bits <= (int)d ) {
+			UTOQ(mod,q);
+			UTOQ(i,ind);
+			n = mknode(2,ind,q);
+			MKLIST(*rp,n);
+			return;
+		}
+	}
+}
+
+void Pranp(NODE arg,P *rp)
+{
 	int n;
 	UP c;
 
@@ -171,9 +337,7 @@ P *rp;
 		*rp = 0;
 }
 
-void ranp(n,nr)
-int n;
-UP *nr;
+void ranp(int n,UP *nr)
 {
 	int i;
 	unsigned int r;
@@ -193,18 +357,21 @@ UP *nr;
 		*nr = 0;
 }
 
-void Pp_mag(arg,rp)
-NODE arg;
-Q *rp;
+void Pmaxblen(NODE arg,Q *rp)
 {
 	int l;
+	l = maxblenp(ARG0(arg));
+	STOQ(l,*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;
 	LIST l;
@@ -254,9 +421,7 @@ LIST *listp;
 	NEXT(tn) = 0; MKLIST(l,n); *listp = l;
 }
 
-void Pcoef0(arg,rp)
-NODE arg;
-Obj *rp;
+void Pcoef0(NODE arg,Obj *rp)
 {
 	Obj t,n;
 	P s;
@@ -297,9 +462,7 @@ Obj *rp;
 	}
 }
 
-void Pcoef(arg,rp)
-NODE arg;
-Obj *rp;
+void Pcoef(NODE arg,Obj *rp)
 {
 	Obj t,n;
 	P s;
@@ -339,9 +502,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;
@@ -362,9 +523,7 @@ Obj *rp;
 		*rp = 0;
 }
 
-void Pdeg(arg,rp)
-NODE arg;
-Q *rp;
+void Pdeg(NODE arg,Q *rp)
 {
 	Obj t,v;
 	int d;
@@ -390,9 +549,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;
 
@@ -402,9 +559,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");
@@ -413,9 +568,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;
 
@@ -428,9 +581,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;
@@ -441,9 +592,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;
@@ -454,9 +603,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;
@@ -467,9 +614,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;
@@ -480,8 +625,7 @@ GF2N *rp;
 	MKGF2N(inv,*rp);
 }
 
-void Prinvtest_gf2n(rp)
-Real *rp;
+void Prinvtest_gf2n(Real *rp)
 {
 	GF2N *a;
 	GF2N c;
@@ -501,9 +645,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
@@ -519,9 +661,7 @@ GF2N *rp;
 #endif
 }
 
-void Pis_irred_gf2(arg,rp)
-NODE arg;
-Q *rp;
+void Pis_irred_gf2(NODE arg,Q *rp)
 {
 	UP2 t;
 
@@ -529,9 +669,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;
@@ -541,17 +679,17 @@ Q *rp;
 	STOQ(ret,*rp);
 }
 
-void Psetmod_ff(arg,rp)
-NODE arg;
-Obj *rp;
+void Psetmod_ff(NODE arg,Obj *rp)
 {
 	int ac;
+	int d;
 	Obj mod,defpoly;
 	N n;
 	UP up;
 	UP2 up2;
-	Q q;
-	P p;
+	UM dp;
+	Q q,r;
+	P p,p1,y;
 	NODE n0,n1;
 	LIST list;
 
@@ -563,7 +701,8 @@ Obj *rp;
 		switch ( OID(mod) ) {
 			case O_N:
 				current_ff = FF_GFP;
-				setmod_lm(NM((Q)mod)); break;
+				setmod_lm(NM((Q)mod));
+				break;
 			case O_P:
 				current_ff = FF_GF2N;
 				setmod_gf2n((P)mod); break;
@@ -571,13 +710,27 @@ Obj *rp;
 				error("setmod_ff : invalid argument");
 		}
 	} else if ( ac == 2 ) {
-		current_ff = FF_GFPN;
-		defpoly = (Obj)ARG0(arg);
-		mod = (Obj)ARG1(arg);
-		if ( !mod || !defpoly )
-			error("setmod_ff : invalid argument");
-		setmod_lm(NM((Q)mod));
-		setmod_gfpn((P)defpoly);
+		if ( OID(ARG0(arg)) == O_N ) {
+			/* small finite field; primitive root representation */
+			current_ff = FF_GFS;
+			setmod_sf(QTOS((Q)ARG0(arg)),QTOS((Q)ARG1(arg)));
+		} else {
+			mod = (Obj)ARG1(arg);
+			current_ff = FF_GFPN;
+			defpoly = (Obj)ARG0(arg);
+			if ( !mod || !defpoly )
+				error("setmod_ff : invalid argument");
+			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:
@@ -590,17 +743,45 @@ Obj *rp;
 			MKNODE(n1,q,0); MKNODE(n0,p,n1);
 			MKLIST(list,n0);
 			*rp = (Obj)list; break;
+		case FF_GFS:
+		case FF_GFSN:
+			STOQ(current_gfs_p,q);
+			if ( current_gfs_ext )
+				enc_to_p(current_gfs_p,current_gfs_iton[1],
+					VR(current_gfs_ext),&p);
+			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:
 			*rp = 0; break;
 	}
 }
 
-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:
@@ -609,13 +790,22 @@ Q *rp;
 			getmod_gf2n(&up2); d = degup2(up2); STOQ(d,*rp); break;
 		case FF_GFPN:
 			getmod_gfpn(&up); STOQ(up->d,*rp); break;
+		case FF_GFS:
+			if ( !current_gfs_ext )
+				*rp = ONE;
+			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;
 
@@ -625,19 +815,20 @@ Q *rp;
 			getmod_lm(&lm); NTOQ(lm,1,*rp); break;
 		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;
 
@@ -645,39 +836,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;
-		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:
@@ -686,75 +851,22 @@ Obj *rp;
 			randomgf2n(&g); *rp = (Obj)g; break;
 		case FF_GFPN:
 			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;
-	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;
-			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;
@@ -779,9 +891,7 @@ P *r;
 	}
 }
 
-void Pdeglist(arg,rp)
-NODE arg;
-LIST *rp;
+void Pdeglist(NODE arg,LIST *rp)
 {
 	NODE d;
 
@@ -789,25 +899,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;
@@ -825,11 +927,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;
@@ -846,10 +944,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;
@@ -868,9 +963,8 @@ NODE *d;
 		*d = d0;
 	}
 }
-void Pmergelist(arg,rp)
-NODE arg;
-LIST *rp;
+
+void Pmergelist(NODE arg,LIST *rp)
 {
     NODE n;
 
@@ -880,8 +974,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;
@@ -913,91 +1006,150 @@ 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 Pptogf2n(arg,rp)
-NODE arg;
-GF2N *rp;
+void Psf_galois_action(NODE arg,P *rp)
 {
-	ptogf2n((Obj)ARG0(arg),rp);
+	sf_galois_action(ARG0(arg),ARG1(arg),rp);
 }
 
-void Pgf2ntop(arg,rp)
-NODE arg;
-P *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)
 {
-	extern V up2_var;
+	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 Psfptopsfp(NODE arg,P *rp)
+{
+	sfptopsfp((P)ARG0(arg),VR((P)ARG1(arg)),rp);
+}
+
+void Pptogf2n(NODE arg,GF2N *rp)
+{
+	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;
 
 	ptoup((P)ARG0(arg),&p);
-	reverseup(p,p->d,&r);
+	if ( argc(arg) == 1 )
+		reverseup(p,p->d,&r);
+	else 
+		reverseup(p,QTOS((Q)ARG1(arg)),&r);
 	uptop(r,rp);
 }
 
-void Putrunc(arg,rp)
-NODE arg;
-P *rp;
+void Putrunc(NODE arg,P *rp)
 {
 	UP p,r;
 
@@ -1006,9 +1158,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;
@@ -1022,9 +1172,7 @@ LIST *rp;
 	MKLIST(*rp,n0);
 }
 
-void Purembymul(arg,rp)
-NODE arg;
-P *rp;
+void Purembymul(NODE arg,P *rp)
 {
 	UP p1,p2,r;
 
@@ -1044,9 +1192,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;
 
@@ -1066,9 +1212,7 @@ P *rp;
 	}
 }
 
-void Puinvmod(arg,rp)
-NODE arg;
-P *rp;
+void Puinvmod(NODE arg,P *rp)
 {
 	UP p,r;
 
@@ -1077,9 +1221,7 @@ P *rp;
 	uptop(r,rp);
 }
 
-void Purevinvmod(arg,rp)
-NODE arg;
-P *rp;
+void Purevinvmod(NODE arg,P *rp)
 {
 	UP p,pr,r;
 
@@ -1089,9 +1231,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;
 
@@ -1102,6 +1242,8 @@ P *rp;
 		case FF_GF2N:
 			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");
@@ -1109,9 +1251,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;
 
@@ -1123,6 +1263,8 @@ P *rp;
 		case FF_GF2N:
 			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");
@@ -1130,9 +1272,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;
@@ -1150,6 +1290,8 @@ VECT *rp;
 		case FF_GF2N:
 			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");
@@ -1159,9 +1301,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;
 
@@ -1170,9 +1310,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;
 
@@ -1182,9 +1320,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;
@@ -1202,17 +1338,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;
 	
@@ -1222,9 +1354,7 @@ P *rp;
 	kmulp(CO,n1,n2,rp);
 }
 
-void Pksquare(arg,rp)
-NODE arg;
-P *rp;
+void Pksquare(NODE arg,P *rp)
 {
 	P n1;
 	
@@ -1233,9 +1363,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;
 
@@ -1245,9 +1373,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;
@@ -1256,7 +1382,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);
@@ -1265,20 +1391,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;
 
@@ -1288,9 +1410,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;
@@ -1304,11 +1424,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);
@@ -1317,9 +1435,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;
@@ -1331,9 +1447,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;
@@ -1374,9 +1488,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;
@@ -1421,9 +1533,7 @@ LIST *rp;
 	}
 }
 
-void Pkmulum(arg,rp)
-NODE arg;
-P *rp;
+void Pkmulum(NODE arg,P *rp)
 {
 	P p1,p2;
 	int d1,d2,mod;
@@ -1438,9 +1548,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;
@@ -1455,9 +1563,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;
 
@@ -1465,4 +1571,81 @@ P *rp;
 	ptoup((P)ARG1(arg),&f);
 	tracemodup_gf2n(g,f,(Q)ARG2(arg),&r);
 	uptop(r,rp);
+}
+
+void Pumul_specialmod(NODE arg,P *rp)
+{
+	P a1,a2;
+	UP p1,p2,r;
+	int i,nmod;
+	int *modind;
+	NODE t,n;
+
+	a1 = (P)ARG0(arg); a2 = (P)ARG1(arg);
+	if ( !a1 || !a2 || NUM(a1) || NUM(a2) )
+		mulp(CO,a1,a2,rp);
+	else {
+		if ( !uzpcheck((Obj)a1) || !uzpcheck((Obj)a2) || VR(a1) != VR(a2) )
+			error("umul_specialmod : invalid argument");
+		ptoup(a1,&p1);
+		ptoup(a2,&p2);
+		n = BDY((LIST)ARG2(arg));
+		nmod = length(n);
+		modind = (int *)MALLOC_ATOMIC(nmod*sizeof(int));
+		for ( i = 0, t = n; i < nmod; i++, t = NEXT(t) )
+			modind[i] = QTOS((Q)BDY(t));
+		fft_mulup_specialmod_main(p1,p2,0,modind,nmod,&r);
+		uptop(r,rp);
+	}
+}
+
+void Pusquare_specialmod(NODE arg,P *rp)
+{
+	P a1;
+	UP p1,r;
+	int i,nmod;
+	int *modind;
+	NODE t,n;
+
+	a1 = (P)ARG0(arg);
+	if ( !a1 || NUM(a1) )
+		mulp(CO,a1,a1,rp);
+	else {
+		if ( !uzpcheck((Obj)a1) )
+			error("usquare_specialmod : invalid argument");
+		ptoup(a1,&p1);
+		n = BDY((LIST)ARG1(arg));
+		nmod = length(n);
+		modind = (int *)MALLOC_ATOMIC(nmod*sizeof(int));
+		for ( i = 0, t = n; i < nmod; i++, t = NEXT(t) )
+			modind[i] = QTOS((Q)BDY(t));
+		fft_mulup_specialmod_main(p1,p1,0,modind,nmod,&r);
+		uptop(r,rp);
+	}
+}
+
+void Putmul_specialmod(NODE arg,P *rp)
+{
+	P a1,a2;
+	UP p1,p2,r;
+	int i,nmod;
+	int *modind;
+	NODE t,n;
+
+	a1 = (P)ARG0(arg); a2 = (P)ARG1(arg);
+	if ( !a1 || !a2 || NUM(a1) || NUM(a2) )
+		mulp(CO,a1,a2,rp);
+	else {
+		if ( !uzpcheck((Obj)a1) || !uzpcheck((Obj)a2) || VR(a1) != VR(a2) )
+			error("utmul_specialmod : invalid argument");
+		ptoup(a1,&p1);
+		ptoup(a2,&p2);
+		n = BDY((LIST)ARG3(arg));
+		nmod = length(n);
+		modind = (int *)MALLOC_ATOMIC(nmod*sizeof(int));
+		for ( i = 0, t = n; i < nmod; i++, t = NEXT(t) )
+			modind[i] = QTOS((Q)BDY(t));
+		fft_mulup_specialmod_main(p1,p2,QTOS((Q)ARG2(arg))+1,modind,nmod,&r);
+		uptop(r,rp);
+	}
 }