=================================================================== 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); + } }