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

Annotation of OpenXM_contrib2/asir2000/builtin/int.c, Revision 1.5

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.5     ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2000/builtin/int.c,v 1.4 2000/08/21 08:31:20 noro Exp $
1.4       noro       49: */
1.1       noro       50: #include "ca.h"
                     51: #include "parse.h"
                     52: #include "base.h"
                     53:
                     54: void Pidiv(), Pirem(), Pigcd(), Pilcm(), Pfac(), Prandom(), Pinv();
                     55: void Pup2_inv(),Pgf2nton(), Pntogf2n();
                     56: void Pup2_init_eg(), Pup2_show_eg();
                     57: void Piqr(), Pprime(), Plprime(), Pinttorat();
                     58: void Piand(), Pior(), Pixor(), Pishift();
                     59: void Pisqrt();
                     60: void iand(), ior(), ixor();
                     61: void isqrt();
                     62: void Plrandom();
                     63: void Pset_upkara(), Pset_uptkara(), Pset_up2kara(), Pset_upfft();
                     64: void Pmt_save(), Pmt_load();
                     65: void Psmall_jacobi();
                     66: void Pdp_set_mpi();
1.2       noro       67: void Pntoint32(),Pint32ton();
1.1       noro       68:
                     69: #ifdef HMEXT
                     70: void Pigcdbin(), Pigcdbmod(), PigcdEuc(), Pigcdacc(), Pigcdcntl();
                     71:
                     72: void Pihex();
                     73: void Pimaxrsh(), Pilen();
                     74: void Ptype_t_NB();
                     75:
                     76: #endif /* HMEXT */
                     77:
                     78: struct ftab int_tab[] = {
                     79:        {"dp_set_mpi",Pdp_set_mpi,-1},
                     80:        {"isqrt",Pisqrt,1},
                     81:        {"idiv",Pidiv,2},
                     82:        {"irem",Pirem,2},
                     83:        {"iqr",Piqr,2},
                     84:        {"igcd",Pigcd,-2},
                     85:        {"ilcm",Pilcm,2},
                     86:        {"up2_inv",Pup2_inv,2},
                     87:        {"up2_init_eg",Pup2_init_eg,0},
                     88:        {"up2_show_eg",Pup2_show_eg,0},
                     89:        {"type_t_NB",Ptype_t_NB,2},
                     90:        {"gf2nton",Pgf2nton,1},
                     91:        {"ntogf2n",Pntogf2n,1},
                     92:        {"set_upkara",Pset_upkara,-1},
                     93:        {"set_uptkara",Pset_uptkara,-1},
                     94:        {"set_up2kara",Pset_up2kara,-1},
                     95:        {"set_upfft",Pset_upfft,-1},
                     96:        {"inv",Pinv,2},
                     97:        {"inttorat",Pinttorat,3},
                     98:        {"fac",Pfac,1},
                     99:        {"prime",Pprime,1},
                    100:        {"lprime",Plprime,1},
                    101:        {"random",Prandom,-1},
                    102:        {"lrandom",Plrandom,1},
                    103:        {"iand",Piand,2},
                    104:        {"ior",Pior,2},
                    105:        {"ixor",Pixor,2},
                    106:        {"ishift",Pishift,2},
                    107:        {"small_jacobi",Psmall_jacobi,2},
                    108: #ifdef HMEXT
                    109:        {"igcdbin",Pigcdbin,2},         /* HM@CCUT extension */
                    110:        {"igcdbmod",Pigcdbmod,2},       /* HM@CCUT extension */
                    111:        {"igcdeuc",PigcdEuc,2},         /* HM@CCUT extension */
                    112:        {"igcdacc",Pigcdacc,2},         /* HM@CCUT extension */
                    113:        {"igcdcntl",Pigcdcntl,-1},      /* HM@CCUT extension */
                    114:        {"ihex",Pihex,1},       /* HM@CCUT extension */
                    115:        {"imaxrsh",Pimaxrsh,1}, /* HM@CCUT extension */
                    116:        {"ilen",Pilen,1},       /* HM@CCUT extension */
                    117: #endif /* HMEXT */
                    118:        {"mt_save",Pmt_save,1},
                    119:        {"mt_load",Pmt_load,1},
1.2       noro      120:        {"ntoint32",Pntoint32,1},
                    121:        {"int32ton",Pint32ton,1},
1.1       noro      122:        {0,0,0},
                    123: };
                    124:
                    125: static int is_prime_small(unsigned int);
                    126: static unsigned int gcd_small(unsigned int,unsigned int);
                    127: int TypeT_NB_check(unsigned int, unsigned int);
                    128: int mpi_mag;
1.2       noro      129:
                    130: void Pntoint32(arg,rp)
                    131: NODE arg;
                    132: USINT *rp;
                    133: {
                    134:        Q q;
1.3       noro      135:        unsigned int t;
1.2       noro      136:
                    137:        asir_assert(ARG0(arg),O_N,"ntoint32");
                    138:        q = (Q)ARG0(arg);
                    139:        if ( !q ) {
                    140:                MKUSINT(*rp,0);
                    141:                return;
                    142:        }
1.3       noro      143:        if ( NID(q)!=N_Q || !INT(q) || PL(NM(q))>1 )
1.2       noro      144:                error("ntoint32 : invalid argument");
1.3       noro      145:        t = BD(NM(q))[0];
                    146:        if ( SGN(q) < 0 )
                    147:                t = -t;
                    148:        MKUSINT(*rp,t);
1.2       noro      149: }
                    150:
                    151: void Pint32ton(arg,rp)
                    152: NODE arg;
                    153: Q *rp;
                    154: {
1.3       noro      155:        int t;
1.2       noro      156:
                    157:        asir_assert(ARG0(arg),O_USINT,"int32ton");
1.3       noro      158:        t = (int)BDY((USINT)ARG0(arg));
                    159:        STOQ(t,*rp);
1.2       noro      160: }
1.1       noro      161:
                    162: void Pdp_set_mpi(arg,rp)
                    163: NODE arg;
                    164: Q *rp;
                    165: {
                    166:        if ( arg ) {
                    167:                asir_assert(ARG0(arg),O_N,"dp_set_mpi");
                    168:                mpi_mag = QTOS((Q)ARG0(arg));
                    169:        }
                    170:        STOQ(mpi_mag,*rp);
                    171: }
                    172:
                    173: void Psmall_jacobi(arg,rp)
                    174: NODE arg;
                    175: Q *rp;
                    176: {
                    177:        Q a,m;
                    178:        int a0,m0,s;
                    179:
                    180:        a = (Q)ARG0(arg);
                    181:        m = (Q)ARG1(arg);
                    182:        asir_assert(a,O_N,"small_jacobi");
                    183:        asir_assert(m,O_N,"small_jacobi");
                    184:        if ( !a )
                    185:                 *rp = ONE;
                    186:        else if ( !m || !INT(m) || !INT(a)
                    187:                || PL(NM(m))>1 || PL(NM(a))>1 || SGN(m) < 0 || EVENN(NM(m)) )
                    188:                error("small_jacobi : invalid input");
                    189:        else {
                    190:                a0 = QTOS(a); m0 = QTOS(m);
                    191:                s = small_jacobi(a0,m0);
                    192:                STOQ(s,*rp);
                    193:        }
                    194: }
                    195:
                    196: int small_jacobi(a,m)
                    197: int a,m;
                    198: {
                    199:        int m4,m8,a4,j1,i,s;
                    200:
                    201:        a %= m;
                    202:        if ( a == 0 || a == 1 )
                    203:                return 1;
                    204:        else if ( a < 0 ) {
                    205:                j1 = small_jacobi(-a,m);
                    206:                m4 = m%4;
                    207:                return m4==1?j1:-j1;
                    208:        } else {
                    209:                for ( i = 0; a && !(a&1); i++, a >>= 1 );
                    210:                if ( i&1 ) {
                    211:                        m8 = m%8;
                    212:                        s = (m8==1||m8==7) ? 1 : -1;
                    213:                } else
                    214:                        s = 1;
                    215:                /* a, m are odd */
                    216:                j1 = small_jacobi(m%a,a);
                    217:                m4 = m%4; a4 = a%4;
                    218:                s *= (m4==1||a4==1) ? 1 : -1;
                    219:                return j1*s;
                    220:        }
                    221: }
                    222:
                    223: void Ptype_t_NB(arg,rp)
                    224: NODE arg;
                    225: Q *rp;
                    226: {
                    227:        if ( TypeT_NB_check(QTOS((Q)ARG0(arg)),QTOS((Q)ARG1(arg))))
                    228:                *rp = ONE;
                    229:        else
                    230:                *rp = 0;
                    231: }
                    232:
                    233: int TypeT_NB_check(unsigned int m, unsigned int t)
                    234: {
                    235:        unsigned int p,k,u,h,d;
                    236:
                    237:        if ( !(m%8) )
                    238:                return 0;
                    239:        p = t*m+1;
                    240:        if ( !is_prime_small(p) )
                    241:                return 0;
                    242:        for ( k = 1, u = 2%p; ; k++ )
                    243:                if ( u == 1 )
                    244:                        break;
                    245:                else
                    246:                        u = (2*u)%p;
                    247:        h = t*m/k;
                    248:        d = gcd_small(h,m);
                    249:        return d == 1 ? 1 : 0;
                    250: }
                    251:
                    252: /*
                    253:  * a simple prime checker
                    254:  * return value: 1  ---  prime number
                    255:  *               0  ---  composite number
                    256:  */
                    257:
                    258: static int is_prime_small(unsigned int a)
                    259: {
                    260:        unsigned int b,t,i;
                    261:
                    262:        if ( !(a%2) ) return 0;
                    263:        for ( t = a, i = 0; t; i++, t >>= 1 );
                    264:        /* b >= sqrt(a) */
                    265:        b = 1<<((i+1)/2);
                    266:
                    267:        /* divisibility test by all odd numbers <= b */
                    268:        for ( i = 3; i <= b; i += 2 )
                    269:                if ( !(a%i) )
                    270:                        return 0;
                    271:        return 1;
                    272: }
                    273:
                    274: /*
                    275:  * gcd for unsigned int as integers
                    276:  * return value: GCD(a,b)
                    277:  *
                    278:  */
                    279:
                    280:
                    281: static unsigned int gcd_small(unsigned int a,unsigned int b)
                    282: {
                    283:        unsigned int t;
                    284:
                    285:        if ( b > a ) {
                    286:                t = a; a = b; b = t;
                    287:        }
                    288:        /* Euclid's algorithm */
                    289:        while ( 1 )
                    290:                if ( !(t = a%b) ) return b;
                    291:                else {
                    292:                        a = b; b = t;
                    293:                }
                    294: }
                    295:
                    296: void Pmt_save(arg,rp)
                    297: NODE arg;
                    298: Q *rp;
                    299: {
                    300:        int ret;
                    301:
                    302:        ret = mt_save(BDY((STRING)ARG0(arg)));
                    303:        STOQ(ret,*rp);
                    304: }
                    305:
                    306: void Pmt_load(arg,rp)
                    307: NODE arg;
                    308: Q *rp;
                    309: {
                    310:        int ret;
                    311:
                    312:        ret = mt_load(BDY((STRING)ARG0(arg)));
                    313:        STOQ(ret,*rp);
                    314: }
                    315:
                    316: void Pisqrt(arg,rp)
                    317: NODE arg;
                    318: Q *rp;
                    319: {
                    320:        Q a;
                    321:        N r;
                    322:
                    323:        a = (Q)ARG0(arg);
                    324:        asir_assert(a,O_N,"isqrt");
                    325:        if ( !a )
                    326:                *rp = 0;
                    327:        else if ( SGN(a) < 0 )
                    328:                error("isqrt : negative argument");
                    329:        else {
                    330:                isqrt(NM(a),&r);
                    331:                NTOQ(r,1,*rp);
                    332:        }
                    333: }
                    334:
                    335: void Pidiv(arg,rp)
                    336: NODE arg;
                    337: Obj *rp;
                    338: {
                    339:        N q,r;
                    340:        Q a;
                    341:        Q dnd,dvr;
                    342:
                    343:        dnd = (Q)ARG0(arg); dvr = (Q)ARG1(arg);
                    344:        asir_assert(dnd,O_N,"idiv");
                    345:        asir_assert(dvr,O_N,"idiv");
                    346:        if ( !dvr )
                    347:                error("idiv: division by 0");
                    348:        else if ( !dnd )
                    349:                *rp = 0;
                    350:        else {
                    351:                divn(NM(dnd),NM(dvr),&q,&r);
                    352:                NTOQ(q,SGN(dnd)*SGN(dvr),a); *rp = (Obj)a;
                    353:        }
                    354: }
                    355:
                    356: void Pirem(arg,rp)
                    357: NODE arg;
                    358: Obj *rp;
                    359: {
                    360:        N q,r;
                    361:        Q a;
                    362:        Q dnd,dvr;
                    363:
                    364:        dnd = (Q)ARG0(arg); dvr = (Q)ARG1(arg);
                    365:        asir_assert(dnd,O_N,"irem");
                    366:        asir_assert(dvr,O_N,"irem");
                    367:        if ( !dvr )
                    368:                error("irem: division by 0");
                    369:        else if ( !dnd )
                    370:                *rp = 0;
                    371:        else {
                    372:                divn(NM(dnd),NM(dvr),&q,&r);
                    373:                NTOQ(r,SGN(dnd),a); *rp = (Obj)a;
                    374:        }
                    375: }
                    376:
                    377: void Piqr(arg,rp)
                    378: NODE arg;
                    379: LIST *rp;
                    380: {
                    381:        N q,r;
                    382:        Q a,b;
                    383:        Q dnd,dvr;
                    384:        NODE node1,node2;
                    385:
                    386:        dnd = (Q)ARG0(arg); dvr = (Q)ARG1(arg);
                    387:        if ( !dvr )
                    388:                error("iqr: division by 0");
                    389:        else if ( !dnd )
                    390:                a = b = 0;
                    391:        else if ( OID(dnd) == O_VECT ) {
                    392:                iqrv((VECT)dnd,dvr,rp); return;
                    393:        } else {
                    394:                asir_assert(dnd,O_N,"iqr");
                    395:                asir_assert(dvr,O_N,"iqr");
                    396:                divn(NM(dnd),NM(dvr),&q,&r);
                    397:                NTOQ(q,SGN(dnd)*SGN(dvr),a);
                    398:                NTOQ(r,SGN(dnd),b);
                    399:        }
                    400:        MKNODE(node2,b,0); MKNODE(node1,a,node2); MKLIST(*rp,node1);
                    401: }
                    402:
                    403: void Pinttorat(arg,rp)
                    404: NODE arg;
                    405: LIST *rp;
                    406: {
                    407:        Q cq,qq,t,u1,v1,r1,nm;
                    408:        N m,b,q,r,c,u2,v2,r2;
                    409:        NODE node1,node2;
                    410:        P p;
                    411:
                    412:        asir_assert(ARG0(arg),O_N,"inttorat");
                    413:        asir_assert(ARG1(arg),O_N,"inttorat");
                    414:        asir_assert(ARG2(arg),O_N,"inttorat");
                    415:        cq = (Q)ARG0(arg); m = NM((Q)ARG1(arg)); b = NM((Q)ARG2(arg));
                    416:        if ( !cq ) {
                    417:                MKNODE(node2,(pointer)ONE,0); MKNODE(node1,0,node2); MKLIST(*rp,node1);
                    418:                return;
                    419:        }
                    420:        divn(NM(cq),m,&q,&r);
                    421:        if ( !r ) {
                    422:                MKNODE(node2,(pointer)ONE,0); MKNODE(node1,0,node2); MKLIST(*rp,node1);
                    423:                return;
                    424:        } else if ( SGN(cq) < 0 ) {
                    425:                subn(m,r,&c);
                    426:        } else
                    427:                c = r;
                    428:        u1 = 0; v1 = ONE; u2 = m; v2 = c;
                    429:        while ( cmpn(v2,b) >= 0 ) {
                    430:                divn(u2,v2,&q,&r2); u2 = v2; v2 = r2;
                    431:                NTOQ(q,1,qq); mulq(qq,v1,&t); subq(u1,t,&r1); u1 = v1; v1 = r1;
                    432:        }
                    433:        if ( cmpn(NM(v1),b) >= 0 )
                    434:                *rp = 0;
                    435:        else {
                    436:                if ( SGN(v1) < 0 ) {
                    437:                        chsgnp((P)v1,&p); v1 = (Q)p; NTOQ(v2,-1,nm);
                    438:                } else
                    439:                        NTOQ(v2,1,nm);
                    440:                MKNODE(node2,v1,0); MKNODE(node1,nm,node2); MKLIST(*rp,node1);
                    441:        }
                    442: }
                    443:
                    444: void Pigcd(arg,rp)
                    445: NODE arg;
                    446: Q *rp;
                    447: {
                    448:        N g;
                    449:        Q n1,n2,a;
                    450:
                    451:        if ( argc(arg) == 1 ) {
                    452:                igcdv((VECT)ARG0(arg),rp);
                    453:                return;
                    454:        }
                    455:        n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
                    456:        asir_assert(n1,O_N,"igcd");
                    457:        asir_assert(n2,O_N,"igcd");
                    458:        if ( !n1 )
                    459:                *rp = n2;
                    460:        else if ( !n2 )
                    461:                *rp = n1;
                    462:        else {
                    463:                gcdn(NM(n1),NM(n2),&g);
                    464:                NTOQ(g,1,a); *rp = a;
                    465:        }
                    466: }
                    467:
                    468: int comp_n(a,b)
                    469: N *a,*b;
                    470: {
                    471:        return cmpn(*a,*b);
                    472: }
                    473:
                    474: void iqrv(a,dvr,rp)
                    475: VECT a;
                    476: Q dvr;
                    477: LIST *rp;
                    478: {
                    479:        int i,n;
                    480:        VECT q,r;
                    481:        Q dnd,qi,ri;
                    482:        Q *b;
                    483:        N qn,rn;
                    484:        NODE n0,n1;
                    485:
                    486:        if ( !dvr )
                    487:                error("iqrv: division by 0");
                    488:        n = a->len; b = (Q *)BDY(a);
                    489:        MKVECT(q,n); MKVECT(r,n);
                    490:        for ( i = 0; i < n; i++ ) {
                    491:                dnd = b[i];
                    492:                if ( !dnd ) {
                    493:                        qi = ri = 0;
                    494:                } else {
                    495:                        divn(NM(dnd),NM(dvr),&qn,&rn);
                    496:                        NTOQ(qn,SGN(dnd)*SGN(dvr),qi);
                    497:                        NTOQ(rn,SGN(dnd),ri);
                    498:                }
                    499:                BDY(q)[i] = (pointer)qi; BDY(r)[i] = (pointer)ri;
                    500:        }
                    501:        MKNODE(n1,r,0); MKNODE(n0,q,n1); MKLIST(*rp,n0);
                    502: }
                    503:
                    504: void igcdv(a,rp)
                    505: VECT a;
                    506: Q *rp;
                    507: {
                    508:        int i,j,n,nz;
                    509:        N g,gt,q,r;
                    510:        N *c;
                    511:
                    512:        n = a->len;
                    513:        c = (N *)ALLOCA(n*sizeof(N));
                    514:        for ( i = 0; i < n; i++ )
                    515:                c[i] = BDY(a)[i]?NM((Q)(BDY(a)[i])):0;
                    516:        qsort(c,n,sizeof(N),(int (*) (const void *,const void *))comp_n);
                    517:        for ( ; n && ! *c; n--, c++ );
                    518:
                    519:        if ( !n ) {
                    520:                *rp = 0; return;
                    521:        } else if ( n == 1 ) {
                    522:                NTOQ(*c,1,*rp); return;
                    523:        }
                    524:        gcdn(c[0],c[1],&g);
                    525:        for ( i = 2; i < n; i++ ) {
                    526:                divn(c[i],g,&q,&r);
                    527:                gcdn(g,r,&gt);
                    528:                if ( !cmpn(g,gt) ) {
                    529:                        for ( j = i+1, nz = 0; j < n; j++ ) {
                    530:                                divn(c[j],g,&q,&r); c[j] = r;
                    531:                                if ( r )
                    532:                                        nz = 1;
                    533:                        }
                    534:                } else
                    535:                        g = gt;
                    536:        }
                    537:        NTOQ(g,1,*rp);
                    538: }
                    539:
                    540: void igcdv_estimate(a,rp)
                    541: VECT a;
                    542: Q *rp;
                    543: {
                    544:        int n,i,m;
                    545:        N s,t,u,g;
                    546:        Q *q;
                    547:
                    548:        n = a->len; q = (Q *)a->body;
                    549:        if ( n == 1 ) {
                    550:                NTOQ(NM(q[0]),1,*rp); return;
                    551:        }
                    552:
                    553:        m = n/2;
                    554:        for ( i = 0 , s = 0; i < m; i++ ) {
                    555:                addn(s,NM(q[i]),&u); s = u;
                    556:        }
                    557:        for ( t = 0; i < n; i++ ) {
                    558:                addn(t,NM(q[i]),&u); t = u;
                    559:        }
                    560:        gcdn(s,t,&g); NTOQ(g,1,*rp);
                    561: }
                    562:
                    563: void Pilcm(arg,rp)
                    564: NODE arg;
                    565: Obj *rp;
                    566: {
                    567:        N g,q,r,l;
                    568:        Q n1,n2,a;
                    569:
                    570:        n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
                    571:        asir_assert(n1,O_N,"ilcm");
                    572:        asir_assert(n2,O_N,"ilcm");
                    573:        if ( !n1 || !n2 )
                    574:                *rp = 0;
                    575:        else {
                    576:                gcdn(NM(n1),NM(n2),&g); divn(NM(n1),g,&q,&r);
                    577:                muln(q,NM(n2),&l); NTOQ(l,1,a); *rp = (Obj)a;
                    578:        }
                    579: }
                    580:
                    581: void Piand(arg,rp)
                    582: NODE arg;
                    583: Q *rp;
                    584: {
                    585:        N g;
                    586:        Q n1,n2,a;
                    587:
                    588:        n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
                    589:        asir_assert(n1,O_N,"iand");
                    590:        asir_assert(n2,O_N,"iand");
                    591:        if ( !n1 || !n2 )
                    592:                *rp = 0;
                    593:        else {
                    594:                iand(NM(n1),NM(n2),&g);
                    595:                NTOQ(g,1,a); *rp = a;
                    596:        }
                    597: }
                    598:
                    599: void Pior(arg,rp)
                    600: NODE arg;
                    601: Q *rp;
                    602: {
                    603:        N g;
                    604:        Q n1,n2,a;
                    605:
                    606:        n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
                    607:        asir_assert(n1,O_N,"ior");
                    608:        asir_assert(n2,O_N,"ior");
                    609:        if ( !n1 )
                    610:                *rp = n2;
                    611:        else if ( !n2 )
                    612:                *rp = n1;
                    613:        else {
                    614:                ior(NM(n1),NM(n2),&g);
                    615:                NTOQ(g,1,a); *rp = a;
                    616:        }
                    617: }
                    618:
                    619: void Pixor(arg,rp)
                    620: NODE arg;
                    621: Q *rp;
                    622: {
                    623:        N g;
                    624:        Q n1,n2,a;
                    625:
                    626:        n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
                    627:        asir_assert(n1,O_N,"ixor");
                    628:        asir_assert(n2,O_N,"ixor");
                    629:        if ( !n1 )
                    630:                *rp = n2;
                    631:        else if ( !n2 )
                    632:                *rp = n1;
                    633:        else {
                    634:                ixor(NM(n1),NM(n2),&g);
                    635:                NTOQ(g,1,a); *rp = a;
                    636:        }
                    637: }
                    638:
                    639: void Pishift(arg,rp)
                    640: NODE arg;
                    641: Q *rp;
                    642: {
                    643:        N g;
                    644:        Q n1,s,a;
                    645:
                    646:        n1 = (Q)ARG0(arg); s = (Q)ARG1(arg);
                    647:        asir_assert(n1,O_N,"ixor");
                    648:        asir_assert(s,O_N,"ixor");
                    649:        if ( !n1 )
                    650:                *rp = 0;
                    651:        else if ( !s )
                    652:                *rp = n1;
                    653:        else {
                    654:                bshiftn(NM(n1),QTOS(s),&g);
                    655:                NTOQ(g,1,a); *rp = a;
                    656:        }
                    657: }
                    658:
                    659: void isqrt(a,r)
                    660: N a,*r;
                    661: {
                    662:        int k;
                    663:        N x,t,x2,xh,quo,rem;
                    664:
                    665:        if ( !a )
                    666:                *r = 0;
                    667:        else if ( UNIN(a) )
                    668:                *r = ONEN;
                    669:        else {
                    670:                k = n_bits(a); /* a <= 2^k-1 */
                    671:                bshiftn(ONEN,-((k>>1)+(k&1)),&x); /* a <= x^2 */
                    672:                while ( 1 ) {
                    673:                        pwrn(x,2,&t);
                    674:                        if ( cmpn(t,a) <= 0 ) {
                    675:                                *r = x; return;
                    676:                        } else {
                    677:                                if ( BD(x)[0] & 1 )
                    678:                                        addn(x,a,&t);
                    679:                                else
                    680:                                        t = a;
                    681:                                bshiftn(x,-1,&x2); divn(t,x2,&quo,&rem);
                    682:                                bshiftn(x,1,&xh); addn(quo,xh,&x);
                    683:                        }
                    684:                }
                    685:        }
                    686: }
                    687:
                    688: void iand(n1,n2,r)
                    689: N n1,n2,*r;
                    690: {
                    691:        int d1,d2,d,i;
                    692:        N nr;
                    693:        int *p1,*p2,*pr;
                    694:
                    695:        d1 = PL(n1); d2 = PL(n2);
                    696:        d = MIN(d1,d2);
                    697:        nr = NALLOC(d);
                    698:        for ( i = 0, p1 = BD(n1), p2 = BD(n2), pr = BD(nr); i < d; i++ )
                    699:                pr[i] = p1[i] & p2[i];
                    700:        for ( i = d-1; i >= 0 && !pr[i]; i-- );
                    701:        if ( i < 0 )
                    702:                *r = 0;
                    703:        else {
                    704:                PL(nr) = i+1; *r = nr;
                    705:        }
                    706: }
                    707:
                    708: void ior(n1,n2,r)
                    709: N n1,n2,*r;
                    710: {
                    711:        int d1,d2,i;
                    712:        N nr;
                    713:        int *p1,*p2,*pr;
                    714:
                    715:        if ( PL(n1) < PL(n2) ) {
                    716:                nr = n1; n1 = n2; n2 = nr;
                    717:        }
                    718:        d1 = PL(n1); d2 = PL(n2);
                    719:        *r = nr = NALLOC(d1);
                    720:        for ( i = 0, p1 = BD(n1), p2 = BD(n2), pr = BD(nr); i < d2; i++ )
                    721:                pr[i] = p1[i] | p2[i];
                    722:        for ( ; i < d1; i++ )
                    723:                pr[i] = p1[i];
                    724:        for ( i = d1-1; i >= 0 && !pr[i]; i-- );
                    725:        if ( i < 0 )
                    726:                *r = 0;
                    727:        else {
                    728:                PL(nr) = i+1; *r = nr;
                    729:        }
                    730: }
                    731:
                    732: void ixor(n1,n2,r)
                    733: N n1,n2,*r;
                    734: {
                    735:        int d1,d2,i;
                    736:        N nr;
                    737:        int *p1,*p2,*pr;
                    738:
                    739:        if ( PL(n1) < PL(n2) ) {
                    740:                nr = n1; n1 = n2; n2 = nr;
                    741:        }
                    742:        d1 = PL(n1); d2 = PL(n2);
                    743:        *r = nr = NALLOC(d1);
                    744:        for ( i = 0, p1 = BD(n1), p2 = BD(n2), pr = BD(nr); i < d2; i++ )
                    745:                pr[i] = p1[i] ^ p2[i];
                    746:        for ( ; i < d1; i++ )
                    747:                pr[i] = p1[i];
                    748:        for ( i = d1-1; i >= 0 && !pr[i]; i-- );
                    749:        if ( i < 0 )
                    750:                *r = 0;
                    751:        else {
                    752:                PL(nr) = i+1; *r = nr;
                    753:        }
                    754: }
                    755:
                    756: void Pup2_init_eg(rp)
                    757: Obj *rp;
                    758: {
                    759:        up2_init_eg();
                    760:        *rp = 0;
                    761: }
                    762:
                    763: void Pup2_show_eg(rp)
                    764: Obj *rp;
                    765: {
                    766:        up2_show_eg();
                    767:        *rp = 0;
                    768: }
                    769:
                    770: void Pgf2nton(arg,rp)
                    771: NODE arg;
                    772: Q *rp;
                    773: {
                    774:        if ( !ARG0(arg) )
                    775:                *rp = 0;
                    776:        else
                    777:                up2ton(((GF2N)ARG0(arg))->body,rp);
                    778: }
                    779:
                    780: void Pntogf2n(arg,rp)
                    781: NODE arg;
                    782: GF2N *rp;
                    783: {
                    784:        UP2 t;
                    785:
                    786:        ntoup2(ARG0(arg),&t);
                    787:        MKGF2N(t,*rp);
                    788: }
                    789:
                    790: void Pup2_inv(arg,rp)
                    791: NODE arg;
                    792: P *rp;
                    793: {
                    794:        UP2 a,b,t;
                    795:
                    796:        ptoup2(ARG0(arg),&a);
                    797:        ptoup2(ARG1(arg),&b);
                    798:        invup2(a,b,&t);
                    799:        up2top(t,rp);
                    800: }
                    801:
                    802: void Pinv(arg,rp)
                    803: NODE arg;
                    804: Num *rp;
                    805: {
                    806:        Num n;
                    807:        Q mod;
                    808:        MQ r;
                    809:        int inv;
                    810:
                    811:        n = (Num)ARG0(arg); mod = (Q)ARG1(arg);
                    812:        asir_assert(n,O_N,"inv");
                    813:        asir_assert(mod,O_N,"inv");
                    814:        if ( !n || !mod )
                    815:                error("inv: invalid input");
                    816:        else
                    817:                switch ( NID(n) ) {
                    818:                        case N_Q:
                    819:                                invl((Q)n,mod,(Q *)rp);
                    820:                                break;
                    821:                        case N_M:
                    822:                                inv = invm(CONT((MQ)n),QTOS(mod));
                    823:                                STOMQ(inv,r);
                    824:                                *rp = (Num)r;
                    825:                                break;
                    826:                        default:
                    827:                                error("inv: invalid input");
                    828:                }
                    829: }
                    830:
                    831: void Pfac(arg,rp)
                    832: NODE arg;
                    833: Q *rp;
                    834: {
                    835:        asir_assert(ARG0(arg),O_N,"fac");
                    836:        factorial(QTOS((Q)ARG0(arg)),rp);
                    837: }
                    838:
                    839: void Plrandom(arg,rp)
                    840: NODE arg;
                    841: Q *rp;
                    842: {
                    843:        N r;
                    844:
                    845:        asir_assert(ARG0(arg),O_N,"lrandom");
                    846:        randomn(QTOS((Q)ARG0(arg)),&r);
                    847:        NTOQ(r,1,*rp);
                    848: }
                    849:
                    850: void Prandom(arg,rp)
                    851: NODE arg;
                    852: Q *rp;
                    853: {
                    854:        unsigned int r;
                    855:
                    856: #if 0
                    857: #if defined(_PA_RISC1_1)
                    858:        r = mrand48()&BMASK;
                    859: #else
                    860:        if ( arg )
                    861:                srandom(QTOS((Q)ARG0(arg)));
                    862:        r = random()&BMASK;
                    863: #endif
                    864: #endif
                    865:        if ( arg )
                    866:                mt_sgenrand(QTOS((Q)ARG0(arg)));
                    867:        r = mt_genrand();
                    868:        UTOQ(r,*rp);
                    869: }
                    870:
                    871: #if defined(VISUAL) || defined(THINK_C)
                    872: void srandom(unsigned int);
                    873:
                    874: static unsigned int R_Next;
                    875:
                    876: unsigned int random() {
                    877:         if ( !R_Next )
                    878:                 R_Next = 1;
                    879:         return R_Next = (R_Next * 1103515245 + 12345);
                    880: }
                    881:
                    882: void srandom(s)
                    883: unsigned int s;
                    884: {
                    885:                if ( s )
                    886:                        R_Next = s;
                    887:         else if ( !R_Next )
                    888:             R_Next = 1;
                    889: }
                    890: #endif
                    891:
                    892: void Pprime(arg,rp)
                    893: NODE arg;
                    894: Q *rp;
                    895: {
                    896:        int i;
                    897:
                    898:        asir_assert(ARG0(arg),O_N,"prime");
                    899:        i = QTOS((Q)ARG0(arg));
                    900:        if ( i < 0 || i >= 1900 )
                    901:                *rp = 0;
                    902:        else
                    903:                STOQ(sprime[i],*rp);
                    904: }
                    905:
                    906: void Plprime(arg,rp)
                    907: NODE arg;
                    908: Q *rp;
                    909: {
                    910:        int i;
                    911:
                    912:        asir_assert(ARG0(arg),O_N,"lprime");
                    913:        i = QTOS((Q)ARG0(arg));
                    914:        if ( i < 0 || i >= 1000 )
                    915:                *rp = 0;
                    916:        else
                    917:                STOQ(lprime[i],*rp);
                    918: }
                    919:
                    920: extern int up_kara_mag, up_tkara_mag, up_fft_mag;
                    921:
                    922: void Pset_upfft(arg,rp)
                    923: NODE arg;
                    924: Q *rp;
                    925: {
                    926:        if ( arg ) {
                    927:                asir_assert(ARG0(arg),O_N,"set_upfft");
                    928:                up_fft_mag = QTOS((Q)ARG0(arg));
                    929:        }
                    930:        STOQ(up_fft_mag,*rp);
                    931: }
                    932:
                    933: void Pset_upkara(arg,rp)
                    934: NODE arg;
                    935: Q *rp;
                    936: {
                    937:        if ( arg ) {
                    938:                asir_assert(ARG0(arg),O_N,"set_upkara");
                    939:                up_kara_mag = QTOS((Q)ARG0(arg));
                    940:        }
                    941:        STOQ(up_kara_mag,*rp);
                    942: }
                    943:
                    944: void Pset_uptkara(arg,rp)
                    945: NODE arg;
                    946: Q *rp;
                    947: {
                    948:        if ( arg ) {
                    949:                asir_assert(ARG0(arg),O_N,"set_uptkara");
                    950:                up_tkara_mag = QTOS((Q)ARG0(arg));
                    951:        }
                    952:        STOQ(up_tkara_mag,*rp);
                    953: }
                    954:
                    955: extern int up2_kara_mag;
                    956:
                    957: void Pset_up2kara(arg,rp)
                    958: NODE arg;
                    959: Q *rp;
                    960: {
                    961:        if ( arg ) {
                    962:                asir_assert(ARG0(arg),O_N,"set_up2kara");
                    963:                up2_kara_mag = QTOS((Q)ARG0(arg));
                    964:        }
                    965:        STOQ(up2_kara_mag,*rp);
                    966: }
                    967:
                    968: #ifdef HMEXT
                    969: void Pigcdbin(arg,rp)
                    970: NODE arg;
                    971: Obj *rp;
                    972: {
                    973:        N g;
                    974:        Q n1,n2,a;
                    975:
                    976:        n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
                    977:        asir_assert(n1,O_N,"igcd");
                    978:        asir_assert(n2,O_N,"igcd");
                    979:        if ( !n1 )
                    980:                *rp = (Obj)n2;
                    981:        else if ( !n2 )
                    982:                *rp = (Obj)n1;
                    983:        else {
                    984:                gcdbinn(NM(n1),NM(n2),&g);
                    985:                NTOQ(g,1,a); *rp = (Obj)a;
                    986:        }
                    987: }
                    988:
                    989: void Pigcdbmod(arg,rp)
                    990: NODE arg;
                    991: Obj *rp;
                    992: {
                    993:        N g;
                    994:        Q n1,n2,a;
                    995:
                    996:        n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
                    997:        asir_assert(n1,O_N,"igcdbmod");
                    998:        asir_assert(n2,O_N,"igcdbmod");
                    999:        if ( !n1 )
                   1000:                *rp = (Obj)n2;
                   1001:        else if ( !n2 )
                   1002:                *rp = (Obj)n1;
                   1003:        else {
                   1004:                gcdbmodn(NM(n1),NM(n2),&g);
                   1005:                NTOQ(g,1,a); *rp = (Obj)a;
                   1006:        }
                   1007: }
                   1008:
                   1009: void Pigcdacc(arg,rp)
                   1010: NODE arg;
                   1011: Obj *rp;
                   1012: {
                   1013:        N g;
                   1014:        Q n1,n2,a;
                   1015:
                   1016:        n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
                   1017:        asir_assert(n1,O_N,"igcdacc");
                   1018:        asir_assert(n2,O_N,"igcdacc");
                   1019:        if ( !n1 )
                   1020:                *rp = (Obj)n2;
                   1021:        else if ( !n2 )
                   1022:                *rp = (Obj)n1;
                   1023:        else {
                   1024:                gcdaccn(NM(n1),NM(n2),&g);
                   1025:                NTOQ(g,1,a); *rp = (Obj)a;
                   1026:        }
                   1027: }
                   1028:
                   1029: void PigcdEuc(arg,rp)
                   1030: NODE arg;
                   1031: Obj *rp;
                   1032: {
                   1033:        N g;
                   1034:        Q n1,n2,a;
                   1035:
                   1036:        n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
                   1037:        asir_assert(n1,O_N,"igcdbmod");
                   1038:        asir_assert(n2,O_N,"igcdbmod");
                   1039:        if ( !n1 )
                   1040:                *rp = (Obj)n2;
                   1041:        else if ( !n2 )
                   1042:                *rp = (Obj)n1;
                   1043:        else {
                   1044:                gcdEuclidn(NM(n1),NM(n2),&g);
                   1045:                NTOQ(g,1,a); *rp = (Obj)a;
                   1046:        }
                   1047: }
                   1048:
                   1049: extern int igcd_algorithm;
                   1050:     /* == 0 : Euclid,
                   1051:      * == 1 : binary,
                   1052:      * == 2 : bmod,
                   1053:      * >= 3 : (Weber's accelerated)/(Jebelean's generalized binary) algorithm,
                   1054:      */
                   1055: extern int igcd_thre_inidiv;
                   1056:     /*
                   1057:      *  In the non-Euclidean algorithms, if the ratio of the lengths (number
                   1058:      *  of words) of two integers is >= igcd_thre_inidiv, we first perform
                   1059:      *  remainder calculation.
                   1060:      *  If == 0, this remainder calculation is not performed.
                   1061:      */
                   1062: extern int igcdacc_thre;
                   1063:     /*
                   1064:      *  In the accelerated algorithm, if the bit-lengths of two integers is
                   1065:      *  > igcdacc_thre, "bmod" reduction is done.
                   1066:      */
                   1067:
                   1068: void Pigcdcntl(arg,rp)
                   1069: NODE arg;
                   1070: Obj *rp;
                   1071: {
                   1072:        Obj p;
                   1073:        Q a;
                   1074:        int k, m;
                   1075:
                   1076:        if ( arg ) {
                   1077:                p = (Obj)ARG0(arg);
                   1078:                if ( !p ) {
                   1079:                        igcd_algorithm = 0;
                   1080:                        *rp = p;
                   1081:                        return;
                   1082:                } else if ( OID(p) == O_N ) {
                   1083:                        k = QTOS((Q)p);
                   1084:                        a = (Q)p;
                   1085:                        if ( k >= 0 ) igcd_algorithm = k;
                   1086:                        else if ( k == -1 ) {
                   1087:                        ret_thre:
                   1088:                                k = - igcd_thre_inidiv - igcdacc_thre*10000;
                   1089:                                STOQ(k,a);
                   1090:                                *rp = (Obj)a;
                   1091:                                return;
                   1092:                        } else {
                   1093:                                if ( (m = (-k)%10000) != 0 ) igcd_thre_inidiv = m;
                   1094:                                if ( (m = -k/10000) != 0 ) igcdacc_thre = m;
                   1095:                                goto ret_thre;
                   1096:                        }
                   1097:                } else if ( OID(p) == O_STR ) {
                   1098:                        char *n = BDY((STRING) p);
                   1099:
                   1100:                        if ( !strcmp( n, "binary" ) || !strcmp( n, "Binary" )
                   1101:                             || !strcmp( n, "bin" ) || !strcmp( n, "Bin" ) )
                   1102:                                k = igcd_algorithm = 1;
                   1103:                        else if ( !strcmp( n, "bmod" ) || !strcmp( n, "Bmod" ) )
                   1104:                                igcd_algorithm = 2;
                   1105:                        else if ( !strcmp( n, "euc" ) || !strcmp( n, "Euc" )
                   1106:                                  || !strcmp( n, "euclid" ) || !strcmp( n, "Euclid" ) )
                   1107:                                igcd_algorithm = 0;
                   1108:                        else if ( !strcmp( n, "acc" ) || !strcmp( n, "Acc" )
                   1109:                             || !strcmp( n, "gen" ) || !strcmp( n, "Gen" )
                   1110:                             || !strcmp( n, "genbin" ) || !strcmp( n, "GenBin" ) )
                   1111:                                igcd_algorithm = 3;
                   1112:                        else error( "igcdhow : invalid algorithm specification" );
                   1113:                }
                   1114:        }
                   1115:        STOQ(igcd_algorithm,a);
                   1116:        *rp = (Obj)a;
                   1117:        return;
                   1118: }
                   1119:
                   1120:
                   1121: void Pimaxrsh(arg,rp)
                   1122: NODE arg;
                   1123: LIST *rp;
                   1124: {
                   1125:        N n1, n2;
                   1126:        Q q;
                   1127:        NODE node1, node2;
                   1128:        int bits;
                   1129:        N maxrshn();
                   1130:
                   1131:        q = (Q)ARG0(arg);
                   1132:        asir_assert(q,O_N,"imaxrsh");
                   1133:        if ( !q ) n1 = n2 = 0;
                   1134:        else {
                   1135:                n1 = maxrshn( NM(q), &bits );
                   1136:                STON( bits, n2 );
                   1137:        }
                   1138:        NTOQ( n2, 1, q );
                   1139:        MKNODE( node2, q, 0 );
                   1140:        NTOQ( n1, 1, q );
                   1141:        MKNODE( node1, q, node2 );
                   1142:        MKLIST( *rp, node1 );
                   1143: }
                   1144: void Pilen(arg,rp)
                   1145: NODE arg;
                   1146: Obj *rp;
                   1147: {
                   1148:        Q q;
                   1149:        int l;
                   1150:
                   1151:        q = (Q)ARG0(arg);
                   1152:        asir_assert(q,O_N,"ilenh");
                   1153:        if ( !q ) l = 0;
                   1154:        else l = PL(NM(q));
                   1155:        STOQ(l,q);
                   1156:        *rp = (Obj)q;
                   1157: }
                   1158:
                   1159: void Pihex(arg,rp)
                   1160: NODE arg;
                   1161: Obj *rp;
                   1162: {
                   1163:        Q n1;
                   1164:
                   1165:        n1 = (Q)ARG0(arg);
                   1166:        asir_assert(n1,O_N,"ihex");
                   1167:        if ( n1 ) {
                   1168:                int i, l = PL(NM(n1)), *p = BD(NM(n1));
                   1169:
                   1170:                for ( i = 0; i < l; i++ ) printf( " 0x%08x", p[i] );
                   1171:                printf( "\n" );
                   1172:        } else printf( " 0x%08x\n", 0 );
                   1173:        *rp = (Obj)n1;
                   1174: }
                   1175: #endif /* HMEXT */

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