[BACK]Return to poly.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2018 / builtin

Annotation of OpenXM_contrib2/asir2018/builtin/poly.c, Revision 1.1

1.1     ! noro        1: /*
        !             2:  * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
        !             3:  * All rights reserved.
        !             4:  *
        !             5:  * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
        !             6:  * non-exclusive and royalty-free license to use, copy, modify and
        !             7:  * redistribute, solely for non-commercial and non-profit purposes, the
        !             8:  * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
        !             9:  * conditions of this Agreement. For the avoidance of doubt, you acquire
        !            10:  * only a limited right to use the SOFTWARE hereunder, and FLL or any
        !            11:  * third party developer retains all rights, including but not limited to
        !            12:  * copyrights, in and to the SOFTWARE.
        !            13:  *
        !            14:  * (1) FLL does not grant you a license in any way for commercial
        !            15:  * purposes. You may use the SOFTWARE only for non-commercial and
        !            16:  * non-profit purposes only, such as academic, research and internal
        !            17:  * business use.
        !            18:  * (2) The SOFTWARE is protected by the Copyright Law of Japan and
        !            19:  * international copyright treaties. If you make copies of the SOFTWARE,
        !            20:  * with or without modification, as permitted hereunder, you shall affix
        !            21:  * to all such copies of the SOFTWARE the above copyright notice.
        !            22:  * (3) An explicit reference to this SOFTWARE and its copyright owner
        !            23:  * shall be made on your publication or presentation in any form of the
        !            24:  * results obtained by use of the SOFTWARE.
        !            25:  * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
        !            26:  * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
        !            27:  * for such modification or the source code of the modified part of the
        !            28:  * SOFTWARE.
        !            29:  *
        !            30:  * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
        !            31:  * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
        !            32:  * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
        !            33:  * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
        !            34:  * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
        !            35:  * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
        !            36:  * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
        !            37:  * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
        !            38:  * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
        !            39:  * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
        !            40:  * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
        !            41:  * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
        !            42:  * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
        !            43:  * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
        !            44:  * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
        !            45:  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
        !            46:  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
        !            47:  *
        !            48:  * $OpenXM$
        !            49: */
        !            50: #include "ca.h"
        !            51: #include "parse.h"
        !            52: #include "base.h"
        !            53:
        !            54: void Pranp();
        !            55:
        !            56: void Pheadsgn();
        !            57: void Pmul_trunc(),Pquo_trunc();
        !            58: void Pumul(),Pumul_ff(),Pusquare(),Pusquare_ff(),Putmul(),Putmul_ff();
        !            59: void Pkmul(),Pksquare(),Pktmul();
        !            60: void Pord(), Premove_vars(), Pcoef0(), Pcoef(), Pdeg(), Pmindeg(), Psetmod();
        !            61: void Pcoef_gf2n();
        !            62: void getcoef(), getdeglist(), mergedeglist(), change_mvar(), restore_mvar();
        !            63:
        !            64: void Pp_mag(),Pmaxblen();
        !            65: void Pmergelist(), Pch_mv(), Pre_mv(), Pdeglist();
        !            66: void Pptomp(),Pmptop();
        !            67: void Pptolmp(),Plmptop();
        !            68: void Pptosfp(),Psfptop(),Psf_galois_action(),Psf_embed(),Psf_find_root();
        !            69: void Psf_minipoly(),Psf_log(),Psfptopsfp();
        !            70: void Pptogf2n(),Pgf2ntop(),Pgf2ntovect();
        !            71: void Pptogfpn(),Pgfpntop();
        !            72: void Pfind_root_gf2n();
        !            73: void Pumul_specialmod(),Pusquare_specialmod(),Putmul_specialmod();
        !            74:
        !            75: void Pureverse(),Putrunc(),Pudecomp(),Purembymul(),Purembymul_precomp();
        !            76: void Puinvmod(),Purevinvmod();
        !            77: void Ppwrmod_ff(),Ppwrtab_ff(),Pgeneric_pwrmod_ff();
        !            78: void Pkpwrmod_lm(),Pkpwrtab_lm(),Pkgeneric_pwrmod_lm();
        !            79:
        !            80: void Pkmulum();
        !            81: void Pksquareum();
        !            82:
        !            83: void Pfmultest();
        !            84: void Phfmul_lm();
        !            85: void Plazy_lm();
        !            86:
        !            87: void Psetmod_ff();
        !            88: void Psimp_ff();
        !            89: void Pextdeg_ff();
        !            90: void Pcharacteristic_ff();
        !            91: void Pfield_type_ff();
        !            92: void Pfield_order_ff();
        !            93: void Prandom_ff();
        !            94:
        !            95: void Ptracemod_gf2n();
        !            96: void Psparsemod_gf2n();
        !            97: void Pmultest_gf2n();
        !            98: void Psquaretest_gf2n();
        !            99: void Pinvtest_gf2n();
        !           100: void Pbininv_gf2n();
        !           101: void Prinvtest_gf2n();
        !           102: void Pis_irred_gf2();
        !           103: void Pis_irred_ddd_gf2();
        !           104: void Pget_next_fft_prime();
        !           105: void Puadj_coef();
        !           106: void Preorder();
        !           107: void Phomogeneous_part();
        !           108: void Phomogeneous_deg();
        !           109: void simp_ff(Obj,Obj *);
        !           110: void ranp(int,UP *);
        !           111: void field_order_ff(Z *);
        !           112:
        !           113: int current_ff;
        !           114:
        !           115: struct ftab poly_tab[] = {
        !           116:   {"headsgn",Pheadsgn,1},
        !           117:   {"quo_trunc",Pquo_trunc,2},
        !           118:   {"mul_trunc",Pmul_trunc,3},
        !           119:   {"homogeneous_deg",Phomogeneous_deg,-2},
        !           120:   {"homogeneous_part",Phomogeneous_part,-3},
        !           121:   {"reorder",Preorder,3},
        !           122:   {"uadj_coef",Puadj_coef,3},
        !           123:   {"ranp",Pranp,2},
        !           124:   {"p_mag",Pp_mag,1},
        !           125:   {"maxblen",Pmaxblen,1},
        !           126:   {"ord",Pord,-1},
        !           127:   {"remove_vars",Premove_vars,1},
        !           128:   {"delete_vars",Premove_vars,1},
        !           129:   {"coef0",Pcoef0,-3},
        !           130:   {"coef",Pcoef,-3},
        !           131:   {"coef_gf2n",Pcoef_gf2n,2},
        !           132:   {"deg",Pdeg,2},
        !           133:   {"mindeg",Pmindeg,2},
        !           134:   {"setmod",Psetmod,-1},
        !           135:
        !           136:   {"sparsemod_gf2n",Psparsemod_gf2n,-1},
        !           137:
        !           138:   {"setmod_ff",Psetmod_ff,-3},
        !           139:   {"simp_ff",Psimp_ff,1},
        !           140:   {"extdeg_ff",Pextdeg_ff,0},
        !           141:   {"characteristic_ff",Pcharacteristic_ff,0},
        !           142:   {"field_type_ff",Pfield_type_ff,0},
        !           143:   {"field_order_ff",Pfield_order_ff,0},
        !           144:   {"random_ff",Prandom_ff,0},
        !           145:
        !           146:   {"deglist",Pdeglist,2},
        !           147:   {"mergelist",Pmergelist,2},
        !           148:   {"ch_mv",Pch_mv,2},
        !           149:   {"re_mv",Pre_mv,2},
        !           150:
        !           151:   {"ptomp",Pptomp,-2},
        !           152:   {"mptop",Pmptop,1},
        !           153:
        !           154:   {"ptolmp",Pptolmp,1},
        !           155:   {"lmptop",Plmptop,1},
        !           156:
        !           157:   {"sf_galois_action",Psf_galois_action,2},
        !           158:   {"sf_find_root",Psf_find_root,1},
        !           159:   {"sf_minipoly",Psf_minipoly,2},
        !           160:   {"sf_embed",Psf_embed,3},
        !           161:   {"sf_log",Psf_log,1},
        !           162:
        !           163:   {"ptosfp",Pptosfp,1},
        !           164:   {"sfptop",Psfptop,1},
        !           165:   {"sfptopsfp",Psfptopsfp,2},
        !           166:   {"ptogf2n",Pptogf2n,1},
        !           167:   {"gf2ntop",Pgf2ntop,-2},
        !           168:   {"gf2ntovect",Pgf2ntovect,1},
        !           169:
        !           170:   {"ptogfpn",Pptogfpn,1},
        !           171:   {"gfpntop",Pgfpntop,-2},
        !           172:
        !           173:   {"kmul",Pkmul,2},
        !           174:   {"ksquare",Pksquare,1},
        !           175:   {"ktmul",Pktmul,3},
        !           176:
        !           177:   {"umul",Pumul,2},
        !           178:   {"usquare",Pusquare,1},
        !           179:   {"ureverse_inv_as_power_series",Purevinvmod,2},
        !           180:   {"uinv_as_power_series",Puinvmod,2},
        !           181:
        !           182:   {"umul_specialmod",Pumul_specialmod,3},
        !           183:   {"usquare_specialmod",Pusquare_specialmod,2},
        !           184:   {"utmul_specialmod",Putmul_specialmod,4},
        !           185:
        !           186:   {"utmul",Putmul,3},
        !           187:   {"umul_ff",Pumul_ff,2},
        !           188:   {"usquare_ff",Pusquare_ff,1},
        !           189:   {"utmul_ff",Putmul_ff,3},
        !           190:
        !           191:   /* for historical reason */
        !           192:   {"trunc",Putrunc,2},
        !           193:   {"decomp",Pudecomp,2},
        !           194:
        !           195:   {"utrunc",Putrunc,2},
        !           196:   {"udecomp",Pudecomp,2},
        !           197:   {"ureverse",Pureverse,-2},
        !           198:   {"urembymul",Purembymul,2},
        !           199:   {"urembymul_precomp",Purembymul_precomp,3},
        !           200:
        !           201:   {"lazy_lm",Plazy_lm,1},
        !           202:   {"lazy_ff",Plazy_lm,1},
        !           203:
        !           204:   {"pwrmod_ff",Ppwrmod_ff,1},
        !           205:   {"generic_pwrmod_ff",Pgeneric_pwrmod_ff,3},
        !           206:   {"pwrtab_ff",Ppwrtab_ff,2},
        !           207:
        !           208:   {"tracemod_gf2n",Ptracemod_gf2n,3},
        !           209:   {"b_find_root_gf2n",Pfind_root_gf2n,1},
        !           210:
        !           211:   {"is_irred_gf2",Pis_irred_gf2,1},
        !           212:   {"is_irred_ddd_gf2",Pis_irred_ddd_gf2,1},
        !           213:
        !           214:   {"kpwrmod_lm",Pkpwrmod_lm,1},
        !           215:   {"kgeneric_pwrmod_lm",Pkgeneric_pwrmod_lm,3},
        !           216:   {"kpwrtab_lm",Pkpwrtab_lm,2},
        !           217:
        !           218:   {"kmulum",Pkmulum,3},
        !           219:   {"ksquareum",Pksquareum,2},
        !           220:
        !           221:   {"fmultest",Pfmultest,3},
        !           222:   {"hfmul_lm",Phfmul_lm,2},
        !           223:
        !           224:   {"multest_gf2n",Pmultest_gf2n,2},
        !           225:   {"squaretest_gf2n",Psquaretest_gf2n,1},
        !           226:   {"bininv_gf2n",Pbininv_gf2n,2},
        !           227:   {"invtest_gf2n",Pinvtest_gf2n,1},
        !           228:   {"rinvtest_gf2n",Prinvtest_gf2n,0},
        !           229:   {"get_next_fft_prime",Pget_next_fft_prime,2},
        !           230:   {0,0,0},
        !           231: };
        !           232:
        !           233: void Pheadsgn(NODE arg,Z *rp)
        !           234: {
        !           235:   int s;
        !           236:
        !           237:   s = headsgn((P)ARG0(arg));
        !           238:   STOQ(s,*rp);
        !           239: }
        !           240:
        !           241: void Pmul_trunc(NODE arg,P *rp)
        !           242: {
        !           243:   P p1,p2,p,h;
        !           244:   VL vl0,vl1,vl2,tvl,vl;
        !           245:   VN vn;
        !           246:   int i,n;
        !           247:
        !           248:   p1 = (P)ARG0(arg);
        !           249:   p2 = (P)ARG1(arg);
        !           250:   get_vars((Obj)p1,&vl1); get_vars((Obj)p2,&vl2); mergev(CO,vl1,vl2,&tvl);
        !           251:   p = (P)ARG2(arg);
        !           252:   get_vars((Obj)p,&vl0); mergev(CO,tvl,vl0,&vl);
        !           253:   for ( tvl = vl, n = 0; tvl; tvl = NEXT(tvl), n++ );
        !           254:   vn = (VN) ALLOCA((n+1)*sizeof(struct oVN));
        !           255:   for ( i = 0, tvl = vl; i < n; tvl = NEXT(tvl), i++ ) {
        !           256:     vn[i].v = tvl->v;
        !           257:     vn[i].n = 0;
        !           258:   }
        !           259:   vn[i].v = 0;
        !           260:   vn[i].n = 0;
        !           261:   for ( h = p, i = 0; OID(h) == O_P; h = COEF(DC(h)) ) {
        !           262:     for ( ; vn[i].v != VR(h); i++ );
        !           263:     vn[i].n = QTOS(DEG(DC(h)));
        !           264:   }
        !           265:   mulp_trunc(vl,p1,p2,vn,rp);
        !           266: }
        !           267:
        !           268: void Pquo_trunc(NODE arg,P *rp)
        !           269: {
        !           270:   P p1,p2,p,h;
        !           271:   VL vl0,vl1,vl2,tvl,vl;
        !           272:   VN vn;
        !           273:   int i,n;
        !           274:
        !           275:   p1 = (P)ARG0(arg);
        !           276:   p2 = (P)ARG1(arg);
        !           277:   if ( !p1 )
        !           278:     *rp = 0;
        !           279:   else if ( NUM(p2) )
        !           280:     divsp(CO,p1,p2,rp);
        !           281:   else {
        !           282:     get_vars((Obj)p1,&vl1); get_vars((Obj)p2,&vl2); mergev(CO,vl1,vl2,&vl);
        !           283:     for ( tvl = vl, n = 0; tvl; tvl = NEXT(tvl), n++ );
        !           284:     vn = (VN) ALLOCA((n+1)*sizeof(struct oVN));
        !           285:     for ( i = 0, tvl = vl; i < n; tvl = NEXT(tvl), i++ ) {
        !           286:       vn[i].v = tvl->v;
        !           287:       vn[i].n = 0;
        !           288:     }
        !           289:     vn[i].v = 0;
        !           290:     vn[i].n = 0;
        !           291:     for ( h = p2, i = 0; OID(h) == O_P; h = COEF(DC(h)) ) {
        !           292:       for ( ; vn[i].v != VR(h); i++ );
        !           293:       vn[i].n = QTOS(DEG(DC(h)));
        !           294:     }
        !           295:     quop_trunc(vl,p1,p2,vn,rp);
        !           296:   }
        !           297: }
        !           298:
        !           299: void Phomogeneous_part(NODE arg,P *rp)
        !           300: {
        !           301:   if ( argc(arg) == 2 )
        !           302:     exthp(CO,(P)ARG0(arg),QTOS((Q)ARG1(arg)),rp);
        !           303:   else
        !           304:     exthpc_generic(CO,(P)ARG0(arg),QTOS((Q)ARG2(arg)),
        !           305:       VR((P)ARG1(arg)),rp);
        !           306: }
        !           307:
        !           308: void Phomogeneous_deg(NODE arg,Z *rp)
        !           309: {
        !           310:   int d;
        !           311:
        !           312:   if ( argc(arg) == 1 )
        !           313:     d = homdeg((P)ARG0(arg));
        !           314:   else
        !           315:     d = getchomdeg(VR((P)ARG1(arg)),(P)ARG0(arg));
        !           316:   STOQ(d,*rp);
        !           317: }
        !           318:
        !           319: /*
        !           320:   p1 = reorder(p,ovl,nvl) => p1 is 'sorted accoding to nvl.
        !           321: */
        !           322:
        !           323: void Preorder(NODE arg,P *rp)
        !           324: {
        !           325:   VL ovl,nvl,tvl;
        !           326:   NODE n;
        !           327:
        !           328:   for ( ovl = 0, n = BDY((LIST)ARG1(arg)); n; n = NEXT(n) ) {
        !           329:     if ( !ovl ) {
        !           330:       NEWVL(ovl); tvl = ovl;
        !           331:     } else {
        !           332:       NEWVL(NEXT(tvl)); tvl = NEXT(tvl);
        !           333:     }
        !           334:       VR(tvl) = VR((P)BDY(n));
        !           335:   }
        !           336:   for ( nvl = 0, n = BDY((LIST)ARG2(arg)); n; n = NEXT(n) ) {
        !           337:     if ( !nvl ) {
        !           338:       NEWVL(nvl); tvl = nvl;
        !           339:     } else {
        !           340:       NEWVL(NEXT(tvl)); tvl = NEXT(tvl);
        !           341:     }
        !           342:       VR(tvl) = VR((P)BDY(n));
        !           343:   }
        !           344:   reorderp(nvl,ovl,(P)ARG0(arg),rp);
        !           345: }
        !           346:
        !           347: /*
        !           348:   uadj_coef(F,M,M2)
        !           349:   if ( F is a non-negative integer )
        !           350:     return F > M2 ? F-M : M2;
        !           351:   else
        !           352:     F = CN*V^N+...+C0
        !           353:     return uadj_coef(CN,M,M2)*V^N+...+uadj_coef(C0,M,M2);
        !           354: */
        !           355:
        !           356: void Puadj_coef(NODE arg,P *rp)
        !           357: {
        !           358:   UP f,r;
        !           359:   Z m,m2;
        !           360:
        !           361:   ptoup((P)ARG0(arg),&f);
        !           362:   m = (Z)ARG1(arg);
        !           363:   m2 = (Z)ARG2(arg);
        !           364:   adj_coefup(f,m,m2,&r);
        !           365:   uptop(r,rp);
        !           366: }
        !           367:
        !           368: /*
        !           369:   get_next_fft_prime(StartIndex,Bits)
        !           370:   tries to find smallest Index >= StartIndex s.t.
        !           371:     2^(Bits-1)|FFTprime[Index]-1
        !           372:   return [Index,Mod] or 0 (not exist)
        !           373: */
        !           374:
        !           375: void Pget_next_fft_prime(NODE arg,LIST *rp)
        !           376: {
        !           377:   unsigned int mod,d;
        !           378:   int start,bits,i;
        !           379:   NODE n;
        !           380:   Z q,ind;
        !           381:
        !           382:   start = QTOS((Q)ARG0(arg));
        !           383:   bits = QTOS((Q)ARG1(arg));
        !           384:   for ( i = start; ; i++ ) {
        !           385:     get_fft_prime(i,&mod,&d);
        !           386:     if ( !mod ) {
        !           387:       *rp = 0; return;
        !           388:     }
        !           389:     if ( bits <= (int)d ) {
        !           390:       UTOQ(mod,q);
        !           391:       UTOQ(i,ind);
        !           392:       n = mknode(2,ind,q);
        !           393:       MKLIST(*rp,n);
        !           394:       return;
        !           395:     }
        !           396:   }
        !           397: }
        !           398:
        !           399: void Pranp(NODE arg,P *rp)
        !           400: {
        !           401:   int n;
        !           402:   UP c;
        !           403:
        !           404:   n = QTOS((Q)ARG0(arg));
        !           405:   ranp(n,&c);
        !           406:   if ( c ) {
        !           407:     up_var = VR((P)ARG1(arg));
        !           408:     uptop(c,rp);
        !           409:   } else
        !           410:     *rp = 0;
        !           411: }
        !           412:
        !           413: void ranp(int n,UP *nr)
        !           414: {
        !           415:   int i;
        !           416:   unsigned int r;
        !           417:   Z q;
        !           418:   UP c;
        !           419:
        !           420:   *nr = c = UPALLOC(n);
        !           421:   for ( i = 0; i <= n; i++ ) {
        !           422:     r = random();
        !           423:     UTOQ(r,q);
        !           424:     c->c[i] = (Num)q;
        !           425:   }
        !           426:   for ( i = n; i >= 0 && !c->c[i]; i-- );
        !           427:   if ( i >= 0 )
        !           428:     c->d = i;
        !           429:   else
        !           430:     *nr = 0;
        !           431: }
        !           432:
        !           433: void Pmaxblen(NODE arg,Z *rp)
        !           434: {
        !           435:   int l;
        !           436:   l = maxblenp(ARG0(arg));
        !           437:   STOQ(l,*rp);
        !           438: }
        !           439:
        !           440: void Pp_mag(NODE arg,Z *rp)
        !           441: {
        !           442:   int l;
        !           443:   l = p_mag(ARG0(arg));
        !           444:   STOQ(l,*rp);
        !           445: }
        !           446:
        !           447: void Pord(NODE arg,LIST *listp)
        !           448: {
        !           449:   NODE n,tn,p,opt;
        !           450:   char *key;
        !           451:   Obj value;
        !           452:   int overwrite=0;
        !           453:   LIST l;
        !           454:   VL vl,tvl,svl;
        !           455:   P t;
        !           456:   int i,j;
        !           457:   V *va;
        !           458:   V v;
        !           459:
        !           460: #if 0
        !           461: printf("LASTCO="); printv(CO,LASTCO->v); printf("\n");
        !           462: #endif
        !           463:   if ( current_option ) {
        !           464:     for ( opt = current_option; opt; opt = NEXT(opt) ) {
        !           465:       p = BDY((LIST)BDY(opt));
        !           466:       key = BDY((STRING)BDY(p));
        !           467:       value = (Obj)BDY(NEXT(p));
        !           468:       if ( !strcmp(key,"overwrite") && value ) {
        !           469:         overwrite = value ? 1 : 0;
        !           470:         break;
        !           471:       }
        !           472:     }
        !           473:   }
        !           474:
        !           475:   if ( argc(arg) ) {
        !           476:     asir_assert(ARG0(arg),O_LIST,"ord");
        !           477:     for ( vl = 0, i = 0, n = BDY((LIST)ARG0(arg));
        !           478:       n; n = NEXT(n), i++ ) {
        !           479:       if ( !vl ) {
        !           480:         NEWVL(vl); tvl = vl;
        !           481:       } else {
        !           482:         NEWVL(NEXT(tvl)); tvl = NEXT(tvl);
        !           483:       }
        !           484:       if ( !(t = (P)BDY(n)) || (OID(t) != O_P) )
        !           485:         error("ord : invalid argument");
        !           486:       VR(tvl) = VR(t);
        !           487:     }
        !           488:     if ( !overwrite ) {
        !           489:       va = (V *)ALLOCA(i*sizeof(V));
        !           490:       for ( j = 0, svl = vl; j < i; j++, svl = NEXT(svl) )
        !           491:         va[j] = VR(svl);
        !           492:       for ( svl = CO; svl; svl = NEXT(svl) ) {
        !           493:         v = VR(svl);
        !           494:         for ( j = 0; j < i; j++ )
        !           495:           if ( v == va[j] )
        !           496:             break;
        !           497:         if ( j == i ) {
        !           498:           if ( !vl ) {
        !           499:             NEWVL(vl); tvl = vl;
        !           500:           } else {
        !           501:             NEWVL(NEXT(tvl)); tvl = NEXT(tvl);
        !           502:           }
        !           503:           VR(tvl) = v;
        !           504:         }
        !           505:       }
        !           506:     } else {
        !           507:       for ( svl = vl; svl; svl = NEXT(svl) ) {
        !           508:         if ( svl->v->attr == (pointer)V_PF )
        !           509:           ((PFINS)svl->v->priv)->pf->ins = 0;
        !           510:       }
        !           511:     }
        !           512:     if ( vl )
        !           513:       NEXT(tvl) = 0;
        !           514:     CO = vl;
        !           515:     update_LASTCO();
        !           516:   }
        !           517:   for ( n = 0, vl = CO; vl; vl = NEXT(vl) ) {
        !           518:     NEXTNODE(n,tn); MKV(VR(vl),t); BDY(tn) = (pointer)t;
        !           519:   }
        !           520:   NEXT(tn) = 0; MKLIST(l,n); *listp = l;
        !           521: }
        !           522:
        !           523: void Premove_vars(NODE arg,LIST *listp)
        !           524: {
        !           525:   NODE l,nd,tnd;
        !           526:   V *v,*va;
        !           527:   int n,na,i,j;
        !           528:   VL vl,vl1;
        !           529:   P t;
        !           530:   LIST list;
        !           531:
        !           532:   asir_assert(ARG0(arg),O_LIST,"remove_vars");
        !           533:   l = BDY((LIST)ARG0(arg)); n = length(l);
        !           534:   v = (V *)ALLOCA(n*sizeof(V));
        !           535:   for ( i = 0; i < n; i++, l = NEXT(l) )
        !           536:     if ( !(t = (P)BDY(l)) || (OID(t) != O_P) )
        !           537:       error("ord : invalid argument");
        !           538:     else v[i] = VR(t);
        !           539:
        !           540:   for ( na = 0, vl = CO; vl; vl = NEXT(vl), na++ );
        !           541:   va = (V *)ALLOCA(na*sizeof(V));
        !           542:   for ( i = 0, vl = CO; i < na; i++, vl = NEXT(vl) ) va[i] = VR(vl);
        !           543:   for ( i = 0; i < na; i++ )
        !           544:     for ( j = 0; j < n; j++ ) if ( va[i] == v[j] ) va[i] = 0;
        !           545:   for ( vl = 0, i = na-1; i >= 0; i-- )
        !           546:     if ( va[i] ) {
        !           547:       NEWVL(vl1); VR(vl1) = va[i]; NEXT(vl1) = vl; vl = vl1;
        !           548:     }
        !           549:   CO = vl;
        !           550:   for ( nd = 0, vl = CO; vl; vl = NEXT(vl) ) {
        !           551:     NEXTNODE(nd,tnd); MKV(VR(vl),t); BDY(tnd) = (pointer)t;
        !           552:   }
        !           553:   if ( nd ) NEXT(tnd) = 0;
        !           554:   MKLIST(list,nd); *listp = list;
        !           555: }
        !           556:
        !           557: void Pcoef0(NODE arg,Obj *rp)
        !           558: {
        !           559:   Obj t,n;
        !           560:   P s;
        !           561:   DCP dc;
        !           562:   int id;
        !           563:   V v;
        !           564:   VL vl;
        !           565:
        !           566:   if ( !(t = (Obj)ARG0(arg)) || ((id = OID(ARG0(arg))) > O_P) )
        !           567:     *rp = 0;
        !           568:   else if ( (n = (Obj)ARG1(arg)) && (OID(n) > O_N) )
        !           569:     *rp = 0;
        !           570:   else if ( id == O_N )
        !           571:     if ( !n )
        !           572:       *rp = t;
        !           573:     else
        !           574:       *rp = 0;
        !           575:   else {
        !           576:     if ( argc(arg) == 3 ) {
        !           577:       if ( (v = VR((P)ARG2(arg))) != VR((P)t) ) {
        !           578:         reordvar(CO,v,&vl); reorderp(vl,CO,(P)t,&s);
        !           579:       } else
        !           580:         s = (P)t;
        !           581:       if ( VR(s) != v ) {
        !           582:         if ( n )
        !           583:           *rp = 0;
        !           584:         else
        !           585:           *rp = t;
        !           586:         return;
        !           587:       }
        !           588:     } else
        !           589:       s = (P)t;
        !           590:     for ( dc = DC(s); dc && cmpz(DEG(dc),(Z)n); dc = NEXT(dc) );
        !           591:     if ( dc )
        !           592:       *rp = (Obj)COEF(dc);
        !           593:     else
        !           594:       *rp = 0;
        !           595:   }
        !           596: }
        !           597:
        !           598: void Pcoef(NODE arg,Obj *rp)
        !           599: {
        !           600:   Obj t,n;
        !           601:   P s;
        !           602:   DCP dc;
        !           603:   int id;
        !           604:   V v;
        !           605:
        !           606:   if ( !(t = (Obj)ARG0(arg)) || ((id = OID(ARG0(arg))) > O_P) )
        !           607:     *rp = 0;
        !           608:   else if ( (n = (Obj)ARG1(arg)) && (OID(n) > O_N) )
        !           609:     *rp = 0;
        !           610:   else if ( id == O_N ) {
        !           611:     if ( !n )
        !           612:       *rp = t;
        !           613:     else
        !           614:       *rp = 0;
        !           615:   } else {
        !           616:     if ( argc(arg) == 3 ) {
        !           617:       if ( (v = VR((P)ARG2(arg))) != VR((P)t) ) {
        !           618:         getcoef(CO,(P)t,v,(Z)n,(P *)rp); return;
        !           619:       } else
        !           620:         s = (P)t;
        !           621:       if ( VR(s) != v ) {
        !           622:         if ( n )
        !           623:           *rp = 0;
        !           624:         else
        !           625:           *rp = t;
        !           626:         return;
        !           627:       }
        !           628:     } else
        !           629:       s = (P)t;
        !           630:     for ( dc = DC(s); dc && cmpz(DEG(dc),(Z)n); dc = NEXT(dc) );
        !           631:     if ( dc )
        !           632:       *rp = (Obj)COEF(dc);
        !           633:     else
        !           634:       *rp = 0;
        !           635:   }
        !           636: }
        !           637:
        !           638: void Pcoef_gf2n(NODE arg,Obj *rp)
        !           639: {
        !           640:   Obj t,n;
        !           641:   int id,d;
        !           642:   UP2 up2;
        !           643:
        !           644:   if ( !(t = (Obj)ARG0(arg)) || ((id = OID(ARG0(arg))) > O_P) )
        !           645:     *rp = 0;
        !           646:   else if ( (n = (Obj)ARG1(arg)) && (OID(n) > O_N) )
        !           647:     *rp = 0;
        !           648:   else if ( id == O_N && NID((Num)t) == N_GF2N ) {
        !           649:     d = QTOS((Q)n);
        !           650:     up2 = ((GF2N)t)->body;
        !           651:     if ( d > degup2(up2) )
        !           652:       *rp = 0;
        !           653:     else
        !           654:       *rp = (Obj)(up2->b[d/BSH]&(((unsigned long)1)<<(d%BSH))?ONE:0);
        !           655:   } else
        !           656:     *rp = 0;
        !           657: }
        !           658:
        !           659: void Pdeg(NODE arg,Z *rp)
        !           660: {
        !           661:   Obj t,v;
        !           662:   int d;
        !           663:
        !           664: #if 0
        !           665:   if ( !(t = (Obj)ARG0(arg)) || (OID(t) != O_P) )
        !           666:     *rp = 0;
        !           667:   else if ( !(v = (Obj)ARG1(arg)) || (VR((P)v) != VR((P)t)) )
        !           668:     *rp = 0;
        !           669:   else
        !           670:     *rp = (Obj)DEG(DC((P)t));
        !           671: #endif
        !           672:   if ( !(t = (Obj)ARG0(arg)) )
        !           673:     STOQ(-1,*rp);
        !           674:   else if ( OID(t) != O_P ) {
        !           675:     if ( OID(t) == O_N && NID(t) == N_GF2N
        !           676:       && (v=(Obj)ARG1(arg)) && OID(v)== O_N && NID(v) == N_GF2N ) {
        !           677:       d = degup2(((GF2N)t)->body);
        !           678:       STOQ(d,*rp);
        !           679:     } else
        !           680:       *rp = 0;
        !           681:   } else
        !           682:     degp(VR((P)ARG1(arg)),(P)ARG0(arg),rp);
        !           683: }
        !           684:
        !           685: void Pmindeg(NODE arg,Z *rp)
        !           686: {
        !           687:   Obj t;
        !           688:
        !           689:   if ( !(t = (Obj)ARG0(arg)) || (OID(t) != O_P) )
        !           690:     *rp = 0;
        !           691:   else
        !           692:     getmindeg(VR((P)ARG1(arg)),(P)ARG0(arg),rp);
        !           693: }
        !           694:
        !           695: void Psetmod(NODE arg,Z *rp)
        !           696: {
        !           697:   if ( arg ) {
        !           698:     asir_assert(ARG0(arg),O_N,"setmod");
        !           699:     current_mod = QTOS((Q)ARG0(arg));
        !           700:   }
        !           701:   STOQ(current_mod,*rp);
        !           702: }
        !           703:
        !           704: void Psparsemod_gf2n(NODE arg,Z *rp)
        !           705: {
        !           706:   int id;
        !           707:
        !           708:   if ( arg && current_mod_gf2n )
        !           709:     current_mod_gf2n->id = ARG0(arg)?1:0;
        !           710:   if ( !current_mod_gf2n )
        !           711:     id = -1;
        !           712:   else
        !           713:     id = current_mod_gf2n->id;
        !           714:   STOQ(id,*rp);
        !           715: }
        !           716:
        !           717: void Pmultest_gf2n(NODE arg,GF2N *rp)
        !           718: {
        !           719:   GF2N a,b,c;
        !           720:   int i;
        !           721:
        !           722:   a = (GF2N)ARG0(arg); b = (GF2N)ARG0(arg);
        !           723:   for ( i = 0; i < 10000; i++ )
        !           724:     mulgf2n(a,b,&c);
        !           725:   *rp = c;
        !           726: }
        !           727:
        !           728: void Psquaretest_gf2n(NODE arg,GF2N *rp)
        !           729: {
        !           730:   GF2N a,c;
        !           731:   int i;
        !           732:
        !           733:   a = (GF2N)ARG0(arg);
        !           734:   for ( i = 0; i < 10000; i++ )
        !           735:     squaregf2n(a,&c);
        !           736:   *rp = c;
        !           737: }
        !           738:
        !           739: void Pinvtest_gf2n(NODE arg,GF2N *rp)
        !           740: {
        !           741:   GF2N a,c;
        !           742:   int i;
        !           743:
        !           744:   a = (GF2N)ARG0(arg);
        !           745:   for ( i = 0; i < 10000; i++ )
        !           746:     invgf2n(a,&c);
        !           747:   *rp = c;
        !           748: }
        !           749:
        !           750: void Pbininv_gf2n(NODE arg,GF2N *rp)
        !           751: {
        !           752:   UP2 a,inv;
        !           753:   int n;
        !           754:
        !           755:   a = ((GF2N)ARG0(arg))->body;
        !           756:   n = QTOS((Q)ARG1(arg));
        !           757:   type1_bin_invup2(a,n,&inv);
        !           758:   MKGF2N(inv,*rp);
        !           759: }
        !           760:
        !           761: void Prinvtest_gf2n(Real *rp)
        !           762: {
        !           763:   GF2N *a;
        !           764:   GF2N c;
        !           765:   double t0,t1,r;
        !           766:   int i;
        !           767:   double get_clock();
        !           768:
        !           769:   a = (GF2N *)ALLOCA(1000*sizeof(GF2N));
        !           770:   for ( i = 0; i < 1000; i++ ) {
        !           771:     randomgf2n(&a[i]);
        !           772:   }
        !           773:   t0 = get_clock();
        !           774:   for ( i = 0; i < 1000; i++ )
        !           775:     invgf2n(a[i],&c);
        !           776:   t1 = get_clock();
        !           777:   r = (t1-t0)/1000;
        !           778:   MKReal(r,*rp);
        !           779: }
        !           780:
        !           781: void Pfind_root_gf2n(NODE arg,GF2N *rp)
        !           782: {
        !           783:
        !           784: #if 0
        !           785:   UP p;
        !           786:
        !           787:   ptoup((P)ARG0(arg),&p);
        !           788:   find_root_gf2n(p,rp);
        !           789: #else
        !           790:   UP2 p;
        !           791:
        !           792:   ptoup2((P)ARG0(arg),&p);
        !           793:   find_root_up2(p,rp);
        !           794: #endif
        !           795: }
        !           796:
        !           797: void Pis_irred_gf2(NODE arg,Z *rp)
        !           798: {
        !           799:   UP2 t;
        !           800:
        !           801:   ptoup2(ARG0(arg),&t);
        !           802:   *rp = irredcheckup2(t) ? ONE : 0;
        !           803: }
        !           804:
        !           805: void Pis_irred_ddd_gf2(NODE arg,Z *rp)
        !           806: {
        !           807:   UP2 t;
        !           808:   int ret;
        !           809:
        !           810:   ptoup2(ARG0(arg),&t);
        !           811:   ret = irredcheck_dddup2(t);
        !           812:   STOQ(ret,*rp);
        !           813: }
        !           814:
        !           815: void Psetmod_ff(NODE arg,Obj *rp)
        !           816: {
        !           817:   int ac;
        !           818:   int d;
        !           819:   Obj mod,defpoly;
        !           820:   Z n;
        !           821:   UP up;
        !           822:   UP2 up2;
        !           823:   UM dp;
        !           824:   Z q,r;
        !           825:   P p,p1,y;
        !           826:   NODE n0,n1;
        !           827:   LIST list;
        !           828:
        !           829:   ac = argc(arg);
        !           830:   if ( ac == 1 ) {
        !           831:     mod = (Obj)ARG0(arg);
        !           832:     if ( !mod )
        !           833:             current_ff = FF_NOT_SET;
        !           834:         else {
        !           835:       switch ( OID(mod) ) {
        !           836:       case O_N:
        !           837:         current_ff = FF_GFP;
        !           838:         setmod_lm((Z)mod);
        !           839:         break;
        !           840:       case O_P:
        !           841:         current_ff = FF_GF2N;
        !           842:         setmod_gf2n((P)mod); break;
        !           843:       default:
        !           844:         error("setmod_ff : invalid argument");
        !           845:       }
        !           846:       }
        !           847:   } else if ( ac == 2 ) {
        !           848:     if ( OID(ARG0(arg)) == O_N ) {
        !           849:       /* small finite field; primitive root representation */
        !           850:       current_ff = FF_GFS;
        !           851:       setmod_sf(QTOS((Q)ARG0(arg)),QTOS((Q)ARG1(arg)));
        !           852:     } else {
        !           853:       mod = (Obj)ARG1(arg);
        !           854:       current_ff = FF_GFPN;
        !           855:       defpoly = (Obj)ARG0(arg);
        !           856:       if ( !mod || !defpoly )
        !           857:         error("setmod_ff : invalid argument");
        !           858:       setmod_lm((Z)mod);
        !           859:       setmod_gfpn((P)defpoly);
        !           860:     }
        !           861:   } else if ( ac == 3 ) {
        !           862:     /* finite extension of a small finite field */
        !           863:     current_ff = FF_GFS;
        !           864:     setmod_sf(QTOS((Q)ARG0(arg)),QTOS((Q)ARG1(arg)));
        !           865:     d = QTOS((Q)ARG2(arg));
        !           866:     generate_defpoly_sfum(d,&dp);
        !           867:     setmod_gfsn(dp);
        !           868:     current_ff = FF_GFSN;
        !           869:   }
        !           870:   switch ( current_ff ) {
        !           871:     case FF_GFP:
        !           872:       getmod_lm(&n); *rp = (Obj)n; break;
        !           873:     case FF_GF2N:
        !           874:       getmod_gf2n(&up2); up2top(up2,&p); *rp = (Obj)p; break;
        !           875:     case FF_GFPN:
        !           876:       getmod_lm(&n);
        !           877:       getmod_gfpn(&up); uptop(up,&p);
        !           878:       MKNODE(n1,n,0); MKNODE(n0,p,n1);
        !           879:       MKLIST(list,n0);
        !           880:       *rp = (Obj)list; break;
        !           881:     case FF_GFS:
        !           882:     case FF_GFSN:
        !           883:       STOQ(current_gfs_p,q);
        !           884:       if ( current_gfs_ext )
        !           885:         enc_to_p(current_gfs_p,current_gfs_iton[1],
        !           886:           VR(current_gfs_ext),&p);
        !           887:       else {
        !           888:         if ( current_gfs_p == 2 )
        !           889:           r = ONE;
        !           890:         else if ( !current_gfs_ntoi )
        !           891:           r = 0;
        !           892:         else
        !           893:           STOQ(current_gfs_iton[1],r);
        !           894:         p = (P)r;
        !           895:       }
        !           896:       switch ( current_ff ) {
        !           897:         case FF_GFS:
        !           898:           n0 = mknode(3,q,current_gfs_ext,p);
        !           899:           break;
        !           900:         case FF_GFSN:
        !           901:           getmod_gfsn(&dp);
        !           902:           makevar("y",&y);
        !           903:           sfumtop(VR(y),dp,&p1);
        !           904:           n0 = mknode(4,q,current_gfs_ext,p,p1);
        !           905:           break;
        !           906:       }
        !           907:       MKLIST(list,n0);
        !           908:       *rp = (Obj)list; break;
        !           909:     default:
        !           910:       *rp = 0; break;
        !           911:   }
        !           912: }
        !           913:
        !           914: void Pextdeg_ff(Z *rp)
        !           915: {
        !           916:   int d;
        !           917:   UP2 up2;
        !           918:   UP up;
        !           919:   UM dp;
        !           920:
        !           921:   switch ( current_ff ) {
        !           922:     case FF_GFP:
        !           923:       *rp = ONE; break;
        !           924:     case FF_GF2N:
        !           925:       getmod_gf2n(&up2); d = degup2(up2); STOQ(d,*rp); break;
        !           926:     case FF_GFPN:
        !           927:       getmod_gfpn(&up); STOQ(up->d,*rp); break;
        !           928:     case FF_GFS:
        !           929:       if ( !current_gfs_ext )
        !           930:         *rp = ONE;
        !           931:       else
        !           932:         *rp = DEG(DC(current_gfs_ext));
        !           933:       break;
        !           934:     case FF_GFSN:
        !           935:       getmod_gfsn(&dp);
        !           936:       STOQ(DEG(dp),*rp);
        !           937:       break;
        !           938:     default:
        !           939:       error("extdeg_ff : current_ff is not set");
        !           940:   }
        !           941: }
        !           942:
        !           943: void Pcharacteristic_ff(Z *rp)
        !           944: {
        !           945:   switch ( current_ff ) {
        !           946:     case FF_GFP:
        !           947:     case FF_GFPN:
        !           948:       getmod_lm(rp); break;
        !           949:     case FF_GF2N:
        !           950:       STOQ(2,*rp); break;
        !           951:     case FF_GFS:
        !           952:     case FF_GFSN:
        !           953:       STOQ(current_gfs_p,*rp); break;
        !           954:     default:
        !           955:       error("characteristic_ff : current_ff is not set");
        !           956:   }
        !           957: }
        !           958:
        !           959: void Pfield_type_ff(Z *rp)
        !           960: {
        !           961:   STOQ(current_ff,*rp);
        !           962: }
        !           963:
        !           964: void Pfield_order_ff(Z *rp)
        !           965: {
        !           966:   field_order_ff(rp);
        !           967: }
        !           968:
        !           969: void Prandom_ff(Obj *rp)
        !           970: {
        !           971:   LM l;
        !           972:   GF2N g;
        !           973:   GFPN p;
        !           974:   GFS s;
        !           975:   GFSN spn;
        !           976:
        !           977:   switch ( current_ff ) {
        !           978:     case FF_GFP:
        !           979:       random_lm(&l); *rp = (Obj)l; break;
        !           980:     case FF_GF2N:
        !           981:       randomgf2n(&g); *rp = (Obj)g; break;
        !           982:     case FF_GFPN:
        !           983:       randomgfpn(&p); *rp = (Obj)p; break;
        !           984:     case FF_GFS:
        !           985:       randomgfs(&s); *rp = (Obj)s; break;
        !           986:     case FF_GFSN:
        !           987:       randomgfsn(&spn); *rp = (Obj)spn; break;
        !           988:     default:
        !           989:       error("random_ff : current_ff is not set");
        !           990:   }
        !           991: }
        !           992:
        !           993: void Psimp_ff(NODE arg,Obj *rp)
        !           994: {
        !           995:   simp_ff((Obj)ARG0(arg),rp);
        !           996: }
        !           997:
        !           998: void getcoef(VL vl,P p,V v,Z d,P *r)
        !           999: {
        !          1000:   P s,t,u,a,b,x;
        !          1001:   DCP dc;
        !          1002:   V w;
        !          1003:
        !          1004:   if ( !p )
        !          1005:     *r = 0;
        !          1006:   else if ( NUM(p) )
        !          1007:     *r = d ? 0 : p;
        !          1008:   else if ( (w=VR(p)) == v ) {
        !          1009:     for ( dc = DC(p); dc && cmpz(DEG(dc),d); dc = NEXT(dc) );
        !          1010:     *r = dc ? COEF(dc) : 0;
        !          1011:   } else {
        !          1012:     MKV(w,x);
        !          1013:     for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
        !          1014:       getcoef(vl,COEF(dc),v,d,&t);
        !          1015:       if ( t ) {
        !          1016:         pwrp(vl,x,DEG(dc),&u); mulp(vl,t,u,&a);
        !          1017:         addp(vl,s,a,&b); s = b;
        !          1018:       }
        !          1019:     }
        !          1020:     *r = s;
        !          1021:   }
        !          1022: }
        !          1023:
        !          1024: void Pdeglist(NODE arg,LIST *rp)
        !          1025: {
        !          1026:   NODE d;
        !          1027:
        !          1028:   getdeglist((P)ARG0(arg),VR((P)ARG1(arg)),&d);
        !          1029:   MKLIST(*rp,d);
        !          1030: }
        !          1031:
        !          1032: void Pch_mv(NODE arg,P *rp)
        !          1033: {
        !          1034:   change_mvar(CO,(P)ARG0(arg),VR((P)ARG1(arg)),rp);
        !          1035: }
        !          1036:
        !          1037: void Pre_mv(NODE arg,P *rp)
        !          1038: {
        !          1039:   restore_mvar(CO,(P)ARG0(arg),VR((P)ARG1(arg)),rp);
        !          1040: }
        !          1041:
        !          1042: void change_mvar(VL vl,P p,V v,P *r)
        !          1043: {
        !          1044:   Z d;
        !          1045:   DCP dc,dc0;
        !          1046:   NODE dl;
        !          1047:
        !          1048:   if ( !p || NUM(p) || (VR(p) == v) )
        !          1049:     *r = p;
        !          1050:   else {
        !          1051:     getdeglist(p,v,&dl);
        !          1052:     for ( dc0 = 0; dl; dl = NEXT(dl) ) {
        !          1053:       NEXTDC(dc0,dc); DEG(dc) = d = (Z)BDY(dl);
        !          1054:       getcoef(vl,p,v,d,&COEF(dc));
        !          1055:     }
        !          1056:     NEXT(dc) = 0; MKP(v,dc0,*r);
        !          1057:   }
        !          1058: }
        !          1059:
        !          1060: void restore_mvar(VL vl,P p,V v,P *r)
        !          1061: {
        !          1062:   P s,u,a,b,x;
        !          1063:   DCP dc;
        !          1064:
        !          1065:   if ( !p || NUM(p) || (VR(p) != v) )
        !          1066:     *r = p;
        !          1067:   else {
        !          1068:     MKV(v,x);
        !          1069:     for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
        !          1070:       pwrp(vl,x,DEG(dc),&u); mulp(vl,COEF(dc),u,&a);
        !          1071:       addp(vl,s,a,&b); s = b;
        !          1072:     }
        !          1073:     *r = s;
        !          1074:   }
        !          1075: }
        !          1076:
        !          1077: void getdeglist(P p,V v,NODE *d)
        !          1078: {
        !          1079:   NODE n,n0,d0,d1,d2;
        !          1080:   DCP dc;
        !          1081:
        !          1082:   if ( !p || NUM(p) ) {
        !          1083:     MKNODE(n,0,0); *d = n;
        !          1084:   } else if ( VR(p) == v ) {
        !          1085:     for ( n0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
        !          1086:       NEXTNODE(n0,n); BDY(n) = (pointer)DEG(dc);
        !          1087:     }
        !          1088:     NEXT(n) = 0; *d = n0;
        !          1089:   } else {
        !          1090:     for ( dc = DC(p), d0 = 0; dc; dc = NEXT(dc) ) {
        !          1091:       getdeglist(COEF(dc),v,&d1); mergedeglist(d0,d1,&d2); d0 = d2;
        !          1092:     }
        !          1093:     *d = d0;
        !          1094:   }
        !          1095: }
        !          1096:
        !          1097: void Pmergelist(NODE arg,LIST *rp)
        !          1098: {
        !          1099:     NODE n;
        !          1100:
        !          1101:   asir_assert(ARG0(arg),O_LIST,"mergelist");
        !          1102:   asir_assert(ARG1(arg),O_LIST,"mergelist");
        !          1103:   mergedeglist(BDY((LIST)ARG0(arg)),BDY((LIST)ARG1(arg)),&n);
        !          1104:   MKLIST(*rp,n);
        !          1105: }
        !          1106:
        !          1107: void mergedeglist(NODE d0,NODE d1,NODE *dr)
        !          1108: {
        !          1109:   NODE t0,t,dt;
        !          1110:   Z d;
        !          1111:   int c;
        !          1112:
        !          1113:   if ( !d0 )
        !          1114:     *dr = d1;
        !          1115:   else {
        !          1116:     while ( d1 ) {
        !          1117:       dt = d1; d1 = NEXT(d1); d = (Z)BDY(dt);
        !          1118:       for ( t0 = 0, t = d0; t; t0 = t, t = NEXT(t) ) {
        !          1119:         c = cmpz(d,(Z)BDY(t));
        !          1120:         if ( !c )
        !          1121:           break;
        !          1122:         else if ( c > 0 ) {
        !          1123:           if ( !t0 ) {
        !          1124:             NEXT(dt) = d0; d0 = dt;
        !          1125:           } else {
        !          1126:             NEXT(t0) = dt; NEXT(dt) = t;
        !          1127:           }
        !          1128:           break;
        !          1129:         }
        !          1130:       }
        !          1131:       if ( !t ) {
        !          1132:         NEXT(t0) = dt; *dr = d0; return;
        !          1133:       }
        !          1134:     }
        !          1135:     *dr = d0;
        !          1136:   }
        !          1137: }
        !          1138:
        !          1139: void Pptomp(NODE arg,P *rp)
        !          1140: {
        !          1141:   int mod;
        !          1142:
        !          1143:   if ( argc(arg) == 1 ) {
        !          1144:     if ( !current_mod )
        !          1145:       error("ptomp : current_mod is not set");
        !          1146:     else
        !          1147:       mod = current_mod;
        !          1148:   } else
        !          1149:     mod = QTOS((Q)ARG1(arg));
        !          1150:   ptomp(mod,(P)ARG0(arg),rp);
        !          1151: }
        !          1152:
        !          1153: void Pmptop(NODE arg,P *rp)
        !          1154: {
        !          1155:   mptop((P)ARG0(arg),rp);
        !          1156: }
        !          1157:
        !          1158: void Pptolmp(NODE arg,P *rp)
        !          1159: {
        !          1160:   ptolmp((P)ARG0(arg),rp);
        !          1161: }
        !          1162:
        !          1163: void Plmptop(NODE arg,P *rp)
        !          1164: {
        !          1165:   lmptop((P)ARG0(arg),rp);
        !          1166: }
        !          1167:
        !          1168: void Psf_galois_action(NODE arg,P *rp)
        !          1169: {
        !          1170:   sf_galois_action(ARG0(arg),ARG1(arg),rp);
        !          1171: }
        !          1172:
        !          1173: /*
        !          1174:   sf_embed(F,B,PM)
        !          1175:   F : an element of GF(pn)
        !          1176:   B : the image of the primitive root of GF(pn)
        !          1177:   PM : order of GF(pm)
        !          1178: */
        !          1179:
        !          1180: void Psf_embed(NODE arg,P *rp)
        !          1181: {
        !          1182:   int k,pm;
        !          1183:
        !          1184:   /* GF(pn)={0,1,a,a^2,...}->GF(pm)={0,1,b,b^2,...}; a->b^k */
        !          1185:   k = CONT((GFS)ARG1(arg));
        !          1186:   pm = QTOS((Q)ARG2(arg));
        !          1187:   sf_embed((P)ARG0(arg),k,pm,rp);
        !          1188: }
        !          1189:
        !          1190: void Psf_log(NODE arg,Z *rp)
        !          1191: {
        !          1192:   int k;
        !          1193:
        !          1194:   if ( !ARG0(arg) )
        !          1195:     error("sf_log : invalid armument");
        !          1196:   k = CONT((GFS)ARG0(arg));
        !          1197:   STOQ(k,*rp);
        !          1198: }
        !          1199:
        !          1200: void Psf_find_root(NODE arg,GFS *rp)
        !          1201: {
        !          1202:   P p;
        !          1203:   Obj t;
        !          1204:   int d;
        !          1205:   UM u;
        !          1206:   int *root;
        !          1207:
        !          1208:   p = (P)ARG0(arg);
        !          1209:   simp_ff((Obj)p,&t); p = (P)t;
        !          1210:   d = getdeg(VR(p),p);
        !          1211:   u = W_UMALLOC(d);
        !          1212:   ptosfum(p,u);
        !          1213:   root = (int *)ALLOCA(d*sizeof(int));
        !          1214:   find_rootsf(u,root);
        !          1215:   MKGFS(IFTOF(root[0]),*rp);
        !          1216: }
        !          1217:
        !          1218: void Psf_minipoly(NODE arg,P *rp)
        !          1219: {
        !          1220:   Obj t;
        !          1221:   P p1,p2;
        !          1222:   int d1,d2;
        !          1223:   UM up1,up2,m;
        !          1224:
        !          1225:   p1 = (P)ARG0(arg); simp_ff((Obj)p1,&t); p1 = (P)t;
        !          1226:   p2 = (P)ARG1(arg); simp_ff((Obj)p2,&t); p2 = (P)t;
        !          1227:   d1 = getdeg(VR(p1),p1); up1 = W_UMALLOC(d1); ptosfum(p1,up1);
        !          1228:   d2 = getdeg(VR(p2),p2); up2 = W_UMALLOC(d2); ptosfum(p2,up2);
        !          1229:   m = W_UMALLOC(d2);
        !          1230:   minipolysf(up1,up2,m);
        !          1231:   sfumtop(VR(p2),m,&p1);
        !          1232:   sfptop(p1,rp);
        !          1233: }
        !          1234:
        !          1235: void Pptosfp(NODE arg,P *rp)
        !          1236: {
        !          1237:   ptosfp(ARG0(arg),rp);
        !          1238: }
        !          1239:
        !          1240: void Psfptop(NODE arg,P *rp)
        !          1241: {
        !          1242:   sfptop((P)ARG0(arg),rp);
        !          1243: }
        !          1244:
        !          1245: void Psfptopsfp(NODE arg,P *rp)
        !          1246: {
        !          1247:   sfptopsfp((P)ARG0(arg),VR((P)ARG1(arg)),rp);
        !          1248: }
        !          1249:
        !          1250: void Pptogf2n(NODE arg,GF2N *rp)
        !          1251: {
        !          1252:   ptogf2n((Obj)ARG0(arg),rp);
        !          1253: }
        !          1254:
        !          1255: void Pgf2ntop(NODE arg,P *rp)
        !          1256: {
        !          1257:   if ( argc(arg) == 2 )
        !          1258:     up2_var = VR((P)ARG1(arg));
        !          1259:   gf2ntop((GF2N)ARG0(arg),rp);
        !          1260: }
        !          1261:
        !          1262: void Pgf2ntovect(NODE arg,VECT *rp)
        !          1263: {
        !          1264:   gf2ntovect((GF2N)ARG0(arg),rp);
        !          1265: }
        !          1266:
        !          1267: void Pptogfpn(NODE arg,GFPN *rp)
        !          1268: {
        !          1269:   ptogfpn((Obj)ARG0(arg),rp);
        !          1270: }
        !          1271:
        !          1272: void Pgfpntop(NODE arg,P *rp)
        !          1273: {
        !          1274:   if ( argc(arg) == 2 )
        !          1275:     up_var = VR((P)ARG1(arg));
        !          1276:   gfpntop((GFPN)ARG0(arg),rp);
        !          1277: }
        !          1278:
        !          1279: void Pureverse(NODE arg,P *rp)
        !          1280: {
        !          1281:   UP p,r;
        !          1282:
        !          1283:   ptoup((P)ARG0(arg),&p);
        !          1284:   if ( argc(arg) == 1 )
        !          1285:     reverseup(p,p->d,&r);
        !          1286:   else
        !          1287:     reverseup(p,QTOS((Q)ARG1(arg)),&r);
        !          1288:   uptop(r,rp);
        !          1289: }
        !          1290:
        !          1291: void Putrunc(NODE arg,P *rp)
        !          1292: {
        !          1293:   UP p,r;
        !          1294:
        !          1295:   ptoup((P)ARG0(arg),&p);
        !          1296:   truncup(p,QTOS((Q)ARG1(arg))+1,&r);
        !          1297:   uptop(r,rp);
        !          1298: }
        !          1299:
        !          1300: void Pudecomp(NODE arg,LIST *rp)
        !          1301: {
        !          1302:   P u,l;
        !          1303:   UP p,up,low;
        !          1304:   NODE n0,n1;
        !          1305:
        !          1306:   ptoup((P)ARG0(arg),&p);
        !          1307:   decompup(p,QTOS((Q)ARG1(arg))+1,&low,&up);
        !          1308:   uptop(low,&l);
        !          1309:   uptop(up,&u);
        !          1310:   MKNODE(n1,u,0); MKNODE(n0,l,n1);
        !          1311:   MKLIST(*rp,n0);
        !          1312: }
        !          1313:
        !          1314: void Purembymul(NODE arg,P *rp)
        !          1315: {
        !          1316:   UP p1,p2,r;
        !          1317:
        !          1318:   if ( !ARG0(arg) || !ARG1(arg) )
        !          1319:     *rp = 0;
        !          1320:   else {
        !          1321:     ptoup((P)ARG0(arg),&p1);
        !          1322:     ptoup((P)ARG1(arg),&p2);
        !          1323:     rembymulup(p1,p2,&r);
        !          1324:     uptop(r,rp);
        !          1325:   }
        !          1326: }
        !          1327:
        !          1328: /*
        !          1329:  * d1 = deg(p1), d2 = deg(p2)
        !          1330:  * d1 <= 2*d2-1
        !          1331:  * p2*inv = 1 mod x^d2
        !          1332:  */
        !          1333:
        !          1334: void Purembymul_precomp(NODE arg,P *rp)
        !          1335: {
        !          1336:   UP p1,p2,inv,r;
        !          1337:
        !          1338:   if ( !ARG0(arg) || !ARG1(arg) )
        !          1339:     *rp = 0;
        !          1340:   else {
        !          1341:     ptoup((P)ARG0(arg),&p1);
        !          1342:     ptoup((P)ARG1(arg),&p2);
        !          1343:     ptoup((P)ARG2(arg),&inv);
        !          1344:     if ( p1->d >= 2*p2->d ) {
        !          1345:       error("urembymul_precomp : degree of 1st arg is too large");
        !          1346: /*      fprintf(stderr,"urembymul_precomp : degree of 1st arg is too large"); */
        !          1347:       remup(p1,p2,&r);
        !          1348:     } else
        !          1349:       hybrid_rembymulup_special(current_ff,p1,p2,inv,&r);
        !          1350:     uptop(r,rp);
        !          1351:   }
        !          1352: }
        !          1353:
        !          1354: void Puinvmod(NODE arg,P *rp)
        !          1355: {
        !          1356:   UP p,r;
        !          1357:
        !          1358:   ptoup((P)ARG0(arg),&p);
        !          1359:   invmodup(p,QTOS((Q)ARG1(arg)),&r);
        !          1360:   uptop(r,rp);
        !          1361: }
        !          1362:
        !          1363: void Purevinvmod(NODE arg,P *rp)
        !          1364: {
        !          1365:   UP p,pr,r;
        !          1366:
        !          1367:   ptoup((P)ARG0(arg),&p);
        !          1368:   reverseup(p,p->d,&pr);
        !          1369:   invmodup(pr,QTOS((Q)ARG1(arg)),&r);
        !          1370:   uptop(r,rp);
        !          1371: }
        !          1372:
        !          1373: void Ppwrmod_ff(NODE arg,P *rp)
        !          1374: {
        !          1375:   UP p1,p2;
        !          1376:
        !          1377:   ptoup((P)ARG0(arg),&p1);
        !          1378:   switch ( current_ff ) {
        !          1379:     case FF_GFP:
        !          1380:       hybrid_powermodup(p1,&p2); break;
        !          1381:     case FF_GF2N:
        !          1382:       powermodup_gf2n(p1,&p2); break;
        !          1383:     case FF_GFPN:
        !          1384:     case FF_GFS:
        !          1385:     case FF_GFSN:
        !          1386:       powermodup(p1,&p2); break;
        !          1387:     default:
        !          1388:       error("pwrmod_ff : current_ff is not set");
        !          1389:   }
        !          1390:   uptop(p2,rp);
        !          1391: }
        !          1392:
        !          1393: void Pgeneric_pwrmod_ff(NODE arg,P *rp)
        !          1394: {
        !          1395:   UP g,f,r;
        !          1396:
        !          1397:   ptoup((P)ARG0(arg),&g);
        !          1398:   ptoup((P)ARG1(arg),&f);
        !          1399:   switch ( current_ff ) {
        !          1400:     case FF_GFP:
        !          1401:       hybrid_generic_powermodup(g,f,(Z)ARG2(arg),&r); break;
        !          1402:     case FF_GF2N:
        !          1403:       generic_powermodup_gf2n(g,f,(Z)ARG2(arg),&r); break;
        !          1404:     case FF_GFPN:
        !          1405:     case FF_GFS:
        !          1406:     case FF_GFSN:
        !          1407:       generic_powermodup(g,f,(Z)ARG2(arg),&r); break;
        !          1408:     default:
        !          1409:       error("generic_pwrmod_ff : current_ff is not set");
        !          1410:   }
        !          1411:   uptop(r,rp);
        !          1412: }
        !          1413:
        !          1414: void Ppwrtab_ff(NODE arg,VECT *rp)
        !          1415: {
        !          1416:   UP f,xp;
        !          1417:   UP *tab;
        !          1418:   VECT r;
        !          1419:   int i,d;
        !          1420:
        !          1421:   ptoup((P)ARG0(arg),&f);
        !          1422:   ptoup((P)ARG1(arg),&xp);
        !          1423:   d = f->d;
        !          1424:
        !          1425:   tab = (UP *)ALLOCA(d*sizeof(UP));
        !          1426:   switch ( current_ff ) {
        !          1427:     case FF_GFP:
        !          1428:       hybrid_powertabup(f,xp,tab); break;
        !          1429:     case FF_GF2N:
        !          1430:       powertabup_gf2n(f,xp,tab); break;
        !          1431:     case FF_GFPN:
        !          1432:     case FF_GFS:
        !          1433:     case FF_GFSN:
        !          1434:       powertabup(f,xp,tab); break;
        !          1435:     default:
        !          1436:       error("pwrtab_ff : current_ff is not set");
        !          1437:   }
        !          1438:   MKVECT(r,d); *rp = r;
        !          1439:   for ( i = 0; i < d; i++ )
        !          1440:     uptop(tab[i],(P *)&BDY(r)[i]);
        !          1441: }
        !          1442:
        !          1443: void Pkpwrmod_lm(NODE arg,P *rp)
        !          1444: {
        !          1445:   UP p1,p2;
        !          1446:
        !          1447:   ptoup((P)ARG0(arg),&p1);
        !          1448:   powermodup(p1,&p2);
        !          1449:   uptop(p2,rp);
        !          1450: }
        !          1451:
        !          1452: void Pkgeneric_pwrmod_lm(NODE arg,P *rp)
        !          1453: {
        !          1454:   UP g,f,r;
        !          1455:
        !          1456:   ptoup((P)ARG0(arg),&g);
        !          1457:   ptoup((P)ARG1(arg),&f);
        !          1458:   generic_powermodup(g,f,(Z)ARG2(arg),&r);
        !          1459:   uptop(r,rp);
        !          1460: }
        !          1461:
        !          1462: void Pkpwrtab_lm(NODE arg,VECT *rp)
        !          1463: {
        !          1464:   UP f,xp;
        !          1465:   UP *tab;
        !          1466:   VECT r;
        !          1467:   int i,d;
        !          1468:
        !          1469:   ptoup((P)ARG0(arg),&f);
        !          1470:   ptoup((P)ARG1(arg),&xp);
        !          1471:   d = f->d;
        !          1472:
        !          1473:   tab = (UP *)ALLOCA(d*sizeof(UP));
        !          1474:   powertabup(f,xp,tab);
        !          1475:   MKVECT(r,d); *rp = r;
        !          1476:   for ( i = 0; i < d; i++ )
        !          1477:     uptop(tab[i],(P *)&BDY(r)[i]);
        !          1478: }
        !          1479:
        !          1480: void Plazy_lm(NODE arg,Q *rp)
        !          1481: {
        !          1482:   lm_lazy = QTOS((Q)ARG0(arg));
        !          1483:   *rp = 0;
        !          1484: }
        !          1485:
        !          1486: void Pkmul(NODE arg,P *rp)
        !          1487: {
        !          1488:   P n1,n2;
        !          1489:
        !          1490:   n1 = (P)ARG0(arg); n2 = (P)ARG1(arg);
        !          1491:   asir_assert(n1,O_P,"kmul");
        !          1492:   asir_assert(n2,O_P,"kmul");
        !          1493:   kmulp(CO,n1,n2,rp);
        !          1494: }
        !          1495:
        !          1496: void Pksquare(NODE arg,P *rp)
        !          1497: {
        !          1498:   P n1;
        !          1499:
        !          1500:   n1 = (P)ARG0(arg);
        !          1501:   asir_assert(n1,O_P,"ksquare");
        !          1502:   ksquarep(CO,n1,rp);
        !          1503: }
        !          1504:
        !          1505: void Pktmul(NODE arg,P *rp)
        !          1506: {
        !          1507:   UP p1,p2,r;
        !          1508:
        !          1509:   ptoup((P)ARG0(arg),&p1);
        !          1510:   ptoup((P)ARG1(arg),&p2);
        !          1511:   tkmulup(p1,p2,QTOS((Q)ARG2(arg))+1,&r);
        !          1512:   uptop(r,rp);
        !          1513: }
        !          1514:
        !          1515: void Pumul(NODE arg,P *rp)
        !          1516: {
        !          1517:   P a1,a2;
        !          1518:   UP p1,p2,r;
        !          1519:
        !          1520:   a1 = (P)ARG0(arg); a2 = (P)ARG1(arg);
        !          1521:   if ( !a1 || !a2 || NUM(a1) || NUM(a2) )
        !          1522:     mulp(CO,a1,a2,rp);
        !          1523:   else {
        !          1524:     if ( !uzpcheck((Obj)a1) || !uzpcheck((Obj)a2) || VR(a1) != VR(a2) )
        !          1525:       error("umul : invalid argument");
        !          1526:     ptoup(a1,&p1);
        !          1527:     ptoup(a2,&p2);
        !          1528:     hybrid_mulup(0,p1,p2,&r);
        !          1529:     uptop(r,rp);
        !          1530:   }
        !          1531: }
        !          1532:
        !          1533: void Pusquare(NODE arg,P *rp)
        !          1534: {
        !          1535:   UP p1,r;
        !          1536:
        !          1537:   ptoup((P)ARG0(arg),&p1);
        !          1538:   hybrid_squareup(0,p1,&r);
        !          1539:   uptop(r,rp);
        !          1540: }
        !          1541:
        !          1542: void Putmul(NODE arg,P *rp)
        !          1543: {
        !          1544:   UP p1,p2,r;
        !          1545:
        !          1546:   ptoup((P)ARG0(arg),&p1);
        !          1547:   ptoup((P)ARG1(arg),&p2);
        !          1548:   hybrid_tmulup(0,p1,p2,QTOS((Q)ARG2(arg))+1,&r);
        !          1549:   uptop(r,rp);
        !          1550: }
        !          1551:
        !          1552: void Pumul_ff(NODE arg,Obj *rp)
        !          1553: {
        !          1554:   P a1,a2;
        !          1555:   UP p1,p2,r;
        !          1556:   P p;
        !          1557:
        !          1558:   a1 = (P)ARG0(arg); a2 = (P)ARG1(arg);
        !          1559:   ptoup(a1,&p1);
        !          1560:   ptoup(a2,&p2);
        !          1561:   hybrid_mulup(current_ff,p1,p2,&r);
        !          1562:   uptop(r,&p);
        !          1563:   simp_ff((Obj)p,rp);
        !          1564: }
        !          1565:
        !          1566: void Pusquare_ff(NODE arg,Obj *rp)
        !          1567: {
        !          1568:   UP p1,r;
        !          1569:   P p;
        !          1570:
        !          1571:   ptoup((P)ARG0(arg),&p1);
        !          1572:   hybrid_squareup(current_ff,p1,&r);
        !          1573:   uptop(r,&p);
        !          1574:   simp_ff((Obj)p,rp);
        !          1575: }
        !          1576:
        !          1577: void Putmul_ff(NODE arg,Obj *rp)
        !          1578: {
        !          1579:   UP p1,p2,r;
        !          1580:   P p;
        !          1581:
        !          1582:   ptoup((P)ARG0(arg),&p1);
        !          1583:   ptoup((P)ARG1(arg),&p2);
        !          1584:   hybrid_tmulup(current_ff,p1,p2,QTOS((Q)ARG2(arg))+1,&r);
        !          1585:   uptop(r,&p);
        !          1586:   simp_ff((Obj)p,rp);
        !          1587: }
        !          1588:
        !          1589: void Phfmul_lm(NODE arg,P *rp)
        !          1590: {
        !          1591:   UP p1,p2,r;
        !          1592:   UP hi,lo,mid,t,s,p10,p11,p20,p21;
        !          1593:   unsigned int m,d;
        !          1594:   int i;
        !          1595:
        !          1596:   ptoup((P)ARG0(arg),&p1);
        !          1597:   ptoup((P)ARG1(arg),&p2);
        !          1598:   d = MAX(p1->d,p2->d);
        !          1599:   for ( m = 1; m < d; m <<= 1 );
        !          1600:   if ( m > d )
        !          1601:     m >>= 1;
        !          1602:   if ( d-m < 10000 ) {
        !          1603:     decompup(p1,m,&p10,&p11); /* p1 = p11*x^m+p10 */
        !          1604:     decompup(p2,m,&p20,&p21); /* p2 = p21*x^m+p20 */
        !          1605:     fft_mulup_lm(p10,p20,&lo);
        !          1606:     kmulup(p11,p21,&hi);
        !          1607:     kmulup(p11,p20,&t); kmulup(p10,p21,&s); addup(t,s,&mid);
        !          1608:     r = UPALLOC(2*d);
        !          1609:     bzero((char *)COEF(r),(2*d+1)*sizeof(Num));
        !          1610:     if ( lo )
        !          1611:       bcopy((char *)COEF(lo),(char *)COEF(r),(DEG(lo)+1)*sizeof(Num));
        !          1612:     if ( hi )
        !          1613:       bcopy((char *)COEF(hi),(char *)(COEF(r)+2*m),(DEG(hi)+1)*sizeof(Num));
        !          1614:     for ( i = 2*d; i >= 0 && !COEF(r)[i]; i-- );
        !          1615:     if ( i < 0 )
        !          1616:       r = 0;
        !          1617:     else {
        !          1618:       DEG(r) = i;
        !          1619:       t = UPALLOC(DEG(mid)+m); DEG(t) = DEG(mid)+m;
        !          1620:       bzero((char *)COEF(t),(DEG(mid)+m+1)*sizeof(Num));
        !          1621:       bcopy((char *)COEF(mid),(char *)(COEF(t)+m),(DEG(mid)+1)*sizeof(Num));
        !          1622:       addup(t,r,&s);
        !          1623:       r = s;
        !          1624:     }
        !          1625:   } else
        !          1626:     fft_mulup_lm(p1,p2,&r);
        !          1627:   uptop(r,rp);
        !          1628: }
        !          1629:
        !          1630: void Pfmultest(NODE arg,LIST *rp)
        !          1631: {
        !          1632:   P p1,p2,r;
        !          1633:   int d1,d2;
        !          1634:   UM w1,w2,wr;
        !          1635:   unsigned int *f1,*f2,*fr,*w;
        !          1636:   int index,mod,root,d,maxint,i;
        !          1637:   int cond;
        !          1638:   Z prime;
        !          1639:   NODE n0,n1;
        !          1640:
        !          1641:   p1 = (P)ARG0(arg); p2 = (P)ARG1(arg); index = QTOS((Q)ARG2(arg));
        !          1642:   FFT_primes(index,&mod,&root,&d);
        !          1643:   maxint = 1<<d;
        !          1644:   d1 = UDEG(p1); d2 = UDEG(p2);
        !          1645:   if ( maxint < d1+d2+1 )
        !          1646:     *rp = 0;
        !          1647:   else {
        !          1648:     w1 = W_UMALLOC(d1); w2 = W_UMALLOC(d2);
        !          1649:     wr = W_UMALLOC(d1+d2);
        !          1650:     ptoum(mod,p1,w1); ptoum(mod,p2,w2);
        !          1651:     f1 = (unsigned int *)ALLOCA(maxint*sizeof(unsigned int));
        !          1652:     f2 = (unsigned int *)ALLOCA(maxint*sizeof(unsigned int));
        !          1653:     fr = (unsigned int *)ALLOCA(maxint*sizeof(unsigned int));
        !          1654:     w = (unsigned int *)ALLOCA(12*maxint*sizeof(unsigned int));
        !          1655:
        !          1656:     for ( i = 0; i <= d1; i++ )
        !          1657:       f1[i] = (unsigned int)w1->c[i];
        !          1658:     for ( i = 0; i <= d2; i++ )
        !          1659:       f2[i] = (unsigned int)w2->c[i];
        !          1660:
        !          1661:     cond = FFT_pol_product(d1,f1,d2,f2,fr,index,w);
        !          1662:     if ( cond )
        !          1663:       error("fmultest : ???");
        !          1664:     wr->d = d1+d2;
        !          1665:     for ( i = 0; i <= d1+d2; i++ )
        !          1666:       wr->c[i] = (unsigned int)fr[i];
        !          1667:     umtop(VR(p1),wr,&r);
        !          1668:     STOQ(mod,prime);
        !          1669:     MKNODE(n1,prime,0);
        !          1670:     MKNODE(n0,r,n1);
        !          1671:     MKLIST(*rp,n0);
        !          1672:   }
        !          1673: }
        !          1674:
        !          1675: void Pkmulum(NODE arg,P *rp)
        !          1676: {
        !          1677:   P p1,p2;
        !          1678:   int d1,d2,mod;
        !          1679:   UM w1,w2,wr;
        !          1680:
        !          1681:   p1 = (P)ARG0(arg); p2 = (P)ARG1(arg); mod = QTOS((Q)ARG2(arg));
        !          1682:   d1 = UDEG(p1); d2 = UDEG(p2);
        !          1683:   w1 = W_UMALLOC(d1); w2 = W_UMALLOC(d2);
        !          1684:   wr = W_UMALLOC(d1+d2);
        !          1685:   ptoum(mod,p1,w1); ptoum(mod,p2,w2);
        !          1686:   kmulum(mod,w1,w2,wr);
        !          1687:   umtop(VR(p1),wr,rp);
        !          1688: }
        !          1689:
        !          1690: void Pksquareum(NODE arg,P *rp)
        !          1691: {
        !          1692:   P p1;
        !          1693:   int d1,mod;
        !          1694:   UM w1,wr;
        !          1695:
        !          1696:   p1 = (P)ARG0(arg); mod = QTOS((Q)ARG1(arg));
        !          1697:   d1 = UDEG(p1);
        !          1698:   w1 = W_UMALLOC(d1);
        !          1699:   wr = W_UMALLOC(2*d1);
        !          1700:   ptoum(mod,p1,w1);
        !          1701:   kmulum(mod,w1,w1,wr);
        !          1702:   umtop(VR(p1),wr,rp);
        !          1703: }
        !          1704:
        !          1705: void Ptracemod_gf2n(NODE arg,P *rp)
        !          1706: {
        !          1707:   UP g,f,r;
        !          1708:
        !          1709:   ptoup((P)ARG0(arg),&g);
        !          1710:   ptoup((P)ARG1(arg),&f);
        !          1711:   tracemodup_gf2n(g,f,(Z)ARG2(arg),&r);
        !          1712:   uptop(r,rp);
        !          1713: }
        !          1714:
        !          1715: void Pumul_specialmod(NODE arg,P *rp)
        !          1716: {
        !          1717:   P a1,a2;
        !          1718:   UP p1,p2,r;
        !          1719:   int i,nmod;
        !          1720:   int *modind;
        !          1721:   NODE t,n;
        !          1722:
        !          1723:   a1 = (P)ARG0(arg); a2 = (P)ARG1(arg);
        !          1724:   if ( !a1 || !a2 || NUM(a1) || NUM(a2) )
        !          1725:     mulp(CO,a1,a2,rp);
        !          1726:   else {
        !          1727:     if ( !uzpcheck((Obj)a1) || !uzpcheck((Obj)a2) || VR(a1) != VR(a2) )
        !          1728:       error("umul_specialmod : invalid argument");
        !          1729:     ptoup(a1,&p1);
        !          1730:     ptoup(a2,&p2);
        !          1731:     n = BDY((LIST)ARG2(arg));
        !          1732:     nmod = length(n);
        !          1733:     modind = (int *)MALLOC_ATOMIC(nmod*sizeof(int));
        !          1734:     for ( i = 0, t = n; i < nmod; i++, t = NEXT(t) )
        !          1735:       modind[i] = QTOS((Q)BDY(t));
        !          1736:     fft_mulup_specialmod_main(p1,p2,0,modind,nmod,&r);
        !          1737:     uptop(r,rp);
        !          1738:   }
        !          1739: }
        !          1740:
        !          1741: void Pusquare_specialmod(NODE arg,P *rp)
        !          1742: {
        !          1743:   P a1;
        !          1744:   UP p1,r;
        !          1745:   int i,nmod;
        !          1746:   int *modind;
        !          1747:   NODE t,n;
        !          1748:
        !          1749:   a1 = (P)ARG0(arg);
        !          1750:   if ( !a1 || NUM(a1) )
        !          1751:     mulp(CO,a1,a1,rp);
        !          1752:   else {
        !          1753:     if ( !uzpcheck((Obj)a1) )
        !          1754:       error("usquare_specialmod : invalid argument");
        !          1755:     ptoup(a1,&p1);
        !          1756:     n = BDY((LIST)ARG1(arg));
        !          1757:     nmod = length(n);
        !          1758:     modind = (int *)MALLOC_ATOMIC(nmod*sizeof(int));
        !          1759:     for ( i = 0, t = n; i < nmod; i++, t = NEXT(t) )
        !          1760:       modind[i] = QTOS((Q)BDY(t));
        !          1761:     fft_mulup_specialmod_main(p1,p1,0,modind,nmod,&r);
        !          1762:     uptop(r,rp);
        !          1763:   }
        !          1764: }
        !          1765:
        !          1766: void Putmul_specialmod(NODE arg,P *rp)
        !          1767: {
        !          1768:   P a1,a2;
        !          1769:   UP p1,p2,r;
        !          1770:   int i,nmod;
        !          1771:   int *modind;
        !          1772:   NODE t,n;
        !          1773:
        !          1774:   a1 = (P)ARG0(arg); a2 = (P)ARG1(arg);
        !          1775:   if ( !a1 || !a2 || NUM(a1) || NUM(a2) )
        !          1776:     mulp(CO,a1,a2,rp);
        !          1777:   else {
        !          1778:     if ( !uzpcheck((Obj)a1) || !uzpcheck((Obj)a2) || VR(a1) != VR(a2) )
        !          1779:       error("utmul_specialmod : invalid argument");
        !          1780:     ptoup(a1,&p1);
        !          1781:     ptoup(a2,&p2);
        !          1782:     n = BDY((LIST)ARG3(arg));
        !          1783:     nmod = length(n);
        !          1784:     modind = (int *)MALLOC_ATOMIC(nmod*sizeof(int));
        !          1785:     for ( i = 0, t = n; i < nmod; i++, t = NEXT(t) )
        !          1786:       modind[i] = QTOS((Q)BDY(t));
        !          1787:     fft_mulup_specialmod_main(p1,p2,QTOS((Q)ARG2(arg))+1,modind,nmod,&r);
        !          1788:     uptop(r,rp);
        !          1789:   }
        !          1790: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>