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

Annotation of OpenXM_contrib2/asir2000/builtin/poly.c, Revision 1.31

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

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