===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/poly.c,v
retrieving revision 1.1.1.1
retrieving revision 1.27
diff -u -p -r1.1.1.1 -r1.27
--- OpenXM_contrib2/asir2000/builtin/poly.c	1999/12/03 07:39:07	1.1.1.1
+++ OpenXM_contrib2/asir2000/builtin/poly.c	2017/09/06 06:25:26	1.27
@@ -1,23 +1,76 @@
-/* $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.26 2017/02/28 07:06:28 noro Exp $ 
+*/
 #include "ca.h"
 #include "parse.h"
 #include "base.h"
 
 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();
 
-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,21 +101,31 @@ 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[] = {
+	{"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},
+	{"uadj_coef",Puadj_coef,3},
 	{"ranp",Pranp,2},
 	{"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},
@@ -72,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},
@@ -85,12 +148,21 @@ struct ftab poly_tab[] = {
 	{"ch_mv",Pch_mv,2},
 	{"re_mv",Pre_mv,2},
 
-	{"ptomp",Pptomp,2},
+	{"ptomp",Pptomp,-2},
 	{"mptop",Pmptop,1},
 
 	{"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 +179,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 +194,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 +226,178 @@ 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 Pheadsgn(NODE arg,Q *rp)
+{
+	int s;
 
-void Pranp(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
+		exthpc_generic(CO,(P)ARG0(arg),QTOS((Q)ARG2(arg)),
+			VR((P)ARG1(arg)),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 +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;
@@ -193,20 +430,26 @@ 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;
+	NODE n,tn,p,opt;
+	char *key;
+	Obj value;
+	int overwrite=0;
 	LIST l;
 	VL vl,tvl,svl;
 	P t;
@@ -214,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));
@@ -227,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;
@@ -254,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;
@@ -297,9 +591,7 @@ Obj *rp;
 	}
 }
 
-void Pcoef(arg,rp)
-NODE arg;
-Obj *rp;
+void Pcoef(NODE arg,Obj *rp)
 {
 	Obj t,n;
 	P s;
@@ -339,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;
@@ -362,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;
@@ -390,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;
 
@@ -402,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");
@@ -413,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;
 
@@ -428,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;
@@ -441,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;
@@ -454,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;
@@ -467,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;
@@ -480,8 +754,7 @@ GF2N *rp;
 	MKGF2N(inv,*rp);
 }
 
-void Prinvtest_gf2n(rp)
-Real *rp;
+void Prinvtest_gf2n(Real *rp)
 {
 	GF2N *a;
 	GF2N c;
@@ -501,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
@@ -519,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;
 
@@ -529,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;
@@ -541,17 +808,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;
 
@@ -559,25 +826,42 @@ 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)); break;
+				setmod_lm(NM((Q)mod));
+				break;
 			case O_P:
 				current_ff = FF_GF2N;
 				setmod_gf2n((P)mod); break;
 			default:
 				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 +874,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 +921,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 +946,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 +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;
-		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 +982,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 +1022,7 @@ P *r;
 	}
 }
 
-void Pdeglist(arg,rp)
-NODE arg;
-LIST *rp;
+void Pdeglist(NODE arg,LIST *rp)
 {
 	NODE d;
 
@@ -789,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;
@@ -825,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;
@@ -846,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;
@@ -868,9 +1094,8 @@ NODE *d;
 		*d = d0;
 	}
 }
-void Pmergelist(arg,rp)
-NODE arg;
-LIST *rp;
+
+void Pmergelist(NODE arg,LIST *rp)
 {
     NODE n;
 
@@ -880,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;
@@ -913,91 +1137,159 @@ 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);
+  int mod;
+
+  if ( argc(arg) == 1 ) {
+    if ( !current_mod )
+      error("ptomp : current_mod is not set");
+    else
+      mod = current_mod;
+  } else
+    mod = QTOS((Q)ARG1(arg));
+	ptomp(mod,(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 +1298,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 +1312,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 +1332,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 +1352,7 @@ P *rp;
 	}
 }
 
-void Puinvmod(arg,rp)
-NODE arg;
-P *rp;
+void Puinvmod(NODE arg,P *rp)
 {
 	UP p,r;
 
@@ -1077,9 +1361,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 +1371,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 +1382,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 +1391,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 +1403,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 +1412,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 +1430,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 +1441,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 +1450,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 +1460,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 +1478,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 +1494,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 +1503,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 +1513,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 +1522,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 +1531,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 +1550,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 +1564,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 +1575,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 +1587,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 +1628,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 +1673,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 +1688,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 +1703,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 +1711,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);
+	}
 }