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

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

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