[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.28

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

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