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

Annotation of OpenXM_contrib2/asir2018/builtin/int.c, Revision 1.2

1.1       noro        1: /*
                      2:  * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
                      3:  * All rights reserved.
                      4:  *
                      5:  * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
                      6:  * non-exclusive and royalty-free license to use, copy, modify and
                      7:  * redistribute, solely for non-commercial and non-profit purposes, the
                      8:  * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
                      9:  * conditions of this Agreement. For the avoidance of doubt, you acquire
                     10:  * only a limited right to use the SOFTWARE hereunder, and FLL or any
                     11:  * third party developer retains all rights, including but not limited to
                     12:  * copyrights, in and to the SOFTWARE.
                     13:  *
                     14:  * (1) FLL does not grant you a license in any way for commercial
                     15:  * purposes. You may use the SOFTWARE only for non-commercial and
                     16:  * non-profit purposes only, such as academic, research and internal
                     17:  * business use.
                     18:  * (2) The SOFTWARE is protected by the Copyright Law of Japan and
                     19:  * international copyright treaties. If you make copies of the SOFTWARE,
                     20:  * with or without modification, as permitted hereunder, you shall affix
                     21:  * to all such copies of the SOFTWARE the above copyright notice.
                     22:  * (3) An explicit reference to this SOFTWARE and its copyright owner
                     23:  * shall be made on your publication or presentation in any form of the
                     24:  * results obtained by use of the SOFTWARE.
                     25:  * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
                     26:  * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
                     27:  * for such modification or the source code of the modified part of the
                     28:  * SOFTWARE.
                     29:  *
                     30:  * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
                     31:  * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
                     32:  * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
                     33:  * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
                     34:  * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
                     35:  * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
                     36:  * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
                     37:  * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
                     38:  * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
                     39:  * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
                     40:  * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
                     41:  * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
                     42:  * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
                     43:  * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
                     44:  * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
                     45:  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
                     46:  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
                     47:  *
1.2     ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2018/builtin/int.c,v 1.1 2018/09/19 05:45:06 noro Exp $
1.1       noro       49: */
                     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 Plrandom();
                     61: void Pset_upkara(), Pset_uptkara(), Pset_up2kara(), Pset_upfft();
                     62: void Pmt_save(), Pmt_load();
                     63: void Psmall_jacobi();
                     64: void Pdp_set_mpi();
                     65: void Pntoint32(),Pint32ton();
                     66:
                     67: void Pigcdbin(), Pigcdbmod(), PigcdEuc(), Pigcdacc(), Pigcdcntl();
                     68:
                     69: void Pihex();
                     70: void Pimaxrsh(), Pilen();
                     71: void Ptype_t_NB();
                     72:
                     73: struct ftab int_tab[] = {
                     74:   {"dp_set_mpi",Pdp_set_mpi,-1},
                     75:   {"isqrt",Pisqrt,1},
                     76:   {"idiv",Pidiv,2},
                     77:   {"irem",Pirem,2},
                     78:   {"iqr",Piqr,2},
                     79:   {"igcd",Pigcd,-2},
                     80:   {"ilcm",Pilcm,2},
                     81:   {"up2_inv",Pup2_inv,2},
                     82:   {"up2_init_eg",Pup2_init_eg,0},
                     83:   {"up2_show_eg",Pup2_show_eg,0},
                     84:   {"type_t_NB",Ptype_t_NB,2},
                     85:   {"gf2nton",Pgf2nton,1},
                     86:   {"ntogf2n",Pntogf2n,1},
                     87:   {"set_upkara",Pset_upkara,-1},
                     88:   {"set_uptkara",Pset_uptkara,-1},
                     89:   {"set_up2kara",Pset_up2kara,-1},
                     90:   {"set_upfft",Pset_upfft,-1},
                     91:   {"inv",Pinv,2},
                     92:   {"inttorat",Pinttorat,3},
                     93:   {"fac",Pfac,1},
                     94:   {"prime",Pprime,1},
                     95:   {"lprime",Plprime,1},
                     96:   {"random",Prandom,-1},
                     97:   {"lrandom",Plrandom,1},
                     98:   {"iand",Piand,2},
                     99:   {"ior",Pior,2},
                    100:   {"ixor",Pixor,2},
                    101:   {"ishift",Pishift,2},
                    102:   {"small_jacobi",Psmall_jacobi,2},
                    103:
                    104:   {"igcdbin",Pigcdbin,2},    /* HM@CCUT extension */
                    105:   {"igcdbmod",Pigcdbmod,2},  /* HM@CCUT extension */
                    106:   {"igcdeuc",PigcdEuc,2},    /* HM@CCUT extension */
                    107:   {"igcdacc",Pigcdacc,2},    /* HM@CCUT extension */
                    108:   {"igcdcntl",Pigcdcntl,-1},  /* HM@CCUT extension */
                    109:
                    110:   {"mt_save",Pmt_save,1},
                    111:   {"mt_load",Pmt_load,1},
                    112:   {"ntoint32",Pntoint32,1},
                    113:   {"int32ton",Pint32ton,1},
                    114:   {0,0,0},
                    115: };
                    116:
                    117: static int is_prime_small(unsigned int);
                    118: static unsigned int gcd_small(unsigned int,unsigned int);
                    119: int TypeT_NB_check(unsigned int, unsigned int);
                    120: int mpi_mag;
                    121:
                    122: void Pntoint32(NODE arg,USINT *rp)
                    123: {
                    124:   Z q,z;
                    125:   unsigned int t;
                    126:
                    127:   asir_assert(ARG0(arg),O_N,"ntoint32");
                    128:   q = (Z)ARG0(arg);
                    129:   if ( !q ) {
                    130:     MKUSINT(*rp,0);
                    131:     return;
                    132:   }
                    133:   if ( !INT(q) || !smallz(q) )
                    134:     error("ntoint32 : invalid argument");
                    135:   absz(q,&z);
                    136:   t = QTOS(z);
                    137:   if ( sgnz(q) < 0 )
                    138:     t = -(int)t;
                    139:   MKUSINT(*rp,t);
                    140: }
                    141:
                    142: void Pint32ton(NODE arg,Z *rp)
                    143: {
                    144:   int t;
                    145:
                    146:   asir_assert(ARG0(arg),O_USINT,"int32ton");
                    147:   t = (int)BDY((USINT)ARG0(arg));
                    148:   STOQ(t,*rp);
                    149: }
                    150:
                    151: void Pdp_set_mpi(NODE arg,Z *rp)
                    152: {
                    153:   if ( arg ) {
                    154:     asir_assert(ARG0(arg),O_N,"dp_set_mpi");
                    155:     mpi_mag = QTOS((Q)ARG0(arg));
                    156:   }
                    157:   STOQ(mpi_mag,*rp);
                    158: }
                    159:
                    160: void Psmall_jacobi(NODE arg,Z *rp)
                    161: {
                    162:   Z a,m;
                    163:   int a0,m0,s;
                    164:
                    165:   a = (Z)ARG0(arg);
                    166:   m = (Z)ARG1(arg);
                    167:   asir_assert(a,O_N,"small_jacobi");
                    168:   asir_assert(m,O_N,"small_jacobi");
                    169:   if ( !a )
                    170:      *rp = ONE;
                    171:   else if ( !m || !INT(m) || !INT(a)
                    172:     || !smallz(m) || !smallz(a) || sgnz(m) < 0 || evenz(m) )
                    173:     error("small_jacobi : invalid input");
                    174:   else {
                    175:     a0 = QTOS(a); m0 = QTOS(m);
                    176:     s = small_jacobi(a0,m0);
                    177:     STOQ(s,*rp);
                    178:   }
                    179: }
                    180:
                    181: int small_jacobi(int a,int m)
                    182: {
                    183:   int m4,m8,a4,j1,i,s;
                    184:
                    185:   a %= m;
                    186:   if ( a == 0 || a == 1 )
                    187:     return 1;
                    188:   else if ( a < 0 ) {
                    189:     j1 = small_jacobi(-a,m);
                    190:     m4 = m%4;
                    191:     return m4==1?j1:-j1;
                    192:   } else {
                    193:     for ( i = 0; a && !(a&1); i++, a >>= 1 );
                    194:     if ( i&1 ) {
                    195:       m8 = m%8;
                    196:       s = (m8==1||m8==7) ? 1 : -1;
                    197:     } else
                    198:       s = 1;
                    199:     /* a, m are odd */
                    200:     j1 = small_jacobi(m%a,a);
                    201:     m4 = m%4; a4 = a%4;
                    202:     s *= (m4==1||a4==1) ? 1 : -1;
                    203:     return j1*s;
                    204:   }
                    205: }
                    206:
                    207: void Ptype_t_NB(NODE arg,Z *rp)
                    208: {
                    209:   if ( TypeT_NB_check(QTOS((Q)ARG0(arg)),QTOS((Q)ARG1(arg))))
                    210:     *rp = ONE;
                    211:   else
                    212:     *rp = 0;
                    213: }
                    214:
                    215: int TypeT_NB_check(unsigned int m, unsigned int t)
                    216: {
                    217:   unsigned int p,k,u,h,d;
                    218:
                    219:   if ( !(m%8) )
                    220:     return 0;
                    221:   p = t*m+1;
                    222:   if ( !is_prime_small(p) )
                    223:     return 0;
                    224:   for ( k = 1, u = 2%p; ; k++ )
                    225:     if ( u == 1 )
                    226:       break;
                    227:     else
                    228:       u = (2*u)%p;
                    229:   h = t*m/k;
                    230:   d = gcd_small(h,m);
                    231:   return d == 1 ? 1 : 0;
                    232: }
                    233:
                    234: /*
                    235:  * a simple prime checker
                    236:  * return value: 1  ---  prime number
                    237:  *               0  ---  composite number
                    238:  */
                    239:
                    240: static int is_prime_small(unsigned int a)
                    241: {
                    242:   unsigned int b,t,i;
                    243:
                    244:   if ( !(a%2) ) return 0;
                    245:   for ( t = a, i = 0; t; i++, t >>= 1 );
                    246:   /* b >= sqrt(a) */
                    247:   b = 1<<((i+1)/2);
                    248:
                    249:   /* divisibility test by all odd numbers <= b */
                    250:   for ( i = 3; i <= b; i += 2 )
                    251:     if ( !(a%i) )
                    252:       return 0;
                    253:   return 1;
                    254: }
                    255:
                    256: /*
                    257:  * gcd for unsigned int as integers
                    258:  * return value: GCD(a,b)
                    259:  *
                    260:  */
                    261:
                    262:
                    263: static unsigned int gcd_small(unsigned int a,unsigned int b)
                    264: {
                    265:   unsigned int t;
                    266:
                    267:   if ( b > a ) {
                    268:     t = a; a = b; b = t;
                    269:   }
                    270:   /* Euclid's algorithm */
                    271:   while ( 1 )
                    272:     if ( !(t = a%b) ) return b;
                    273:     else {
                    274:       a = b; b = t;
                    275:     }
                    276: }
                    277:
                    278: void Pmt_save(NODE arg,Z *rp)
                    279: {
                    280:   int ret;
                    281:
                    282:   ret = mt_save(BDY((STRING)ARG0(arg)));
                    283:   STOQ(ret,*rp);
                    284: }
                    285:
                    286: void Pmt_load(NODE arg,Z *rp)
                    287: {
                    288:   int ret;
                    289:
                    290:   ret = mt_load(BDY((STRING)ARG0(arg)));
                    291:   STOQ(ret,*rp);
                    292: }
                    293:
                    294: void isqrt(Z a,Z *r);
                    295:
                    296: void Pisqrt(NODE arg,Z *rp)
                    297: {
                    298:   Z a;
                    299:   Z r;
                    300:
                    301:   a = (Z)ARG0(arg);
                    302:   asir_assert(a,O_N,"isqrt");
                    303:   if ( !a )
                    304:     *rp = 0;
                    305:   else if ( sgnz(a) < 0 )
                    306:     error("isqrt : negative argument");
                    307:   else {
                    308:     isqrt(a,rp);
                    309:   }
                    310: }
                    311:
                    312: void Pidiv(NODE arg,Z *rp)
                    313: {
                    314:   Z r;
                    315:   Z dnd,dvr;
                    316:
                    317:   dnd = (Z)ARG0(arg); dvr = (Z)ARG1(arg);
                    318:   asir_assert(dnd,O_N,"idiv");
                    319:   asir_assert(dvr,O_N,"idiv");
                    320:   if ( !dvr )
                    321:     error("idiv: division by 0");
                    322:   else if ( !dnd )
                    323:     *rp = 0;
                    324:   else
                    325:     divqrz(dnd,dvr,rp,&r);
                    326: }
                    327:
                    328: void Pirem(NODE arg,Z *rp)
                    329: {
                    330:   Z q,dnd,dvr;
                    331:
                    332:   dnd = (Z)ARG0(arg); dvr = (Z)ARG1(arg);
                    333:   asir_assert(dnd,O_N,"irem");
                    334:   asir_assert(dvr,O_N,"irem");
                    335:   if ( !dvr )
                    336:     error("irem: division by 0");
                    337:   else if ( !dnd )
                    338:     *rp = 0;
                    339:   else
                    340:     divqrz(dnd,dvr,&q,rp);
                    341: }
                    342:
                    343: void iqrv(VECT a,Z dvr,LIST *rp);
                    344:
                    345: void Piqr(NODE arg,LIST *rp)
                    346: {
                    347:   Z dnd,dvr,a,b;
                    348:   NODE node;
                    349:
                    350:   dnd = (Z)ARG0(arg); dvr = (Z)ARG1(arg);
                    351:   if ( !dvr )
                    352:     error("iqr: division by 0");
                    353:   else if ( !dnd )
                    354:     a = b = 0;
                    355:   else if ( OID(dnd) == O_VECT ) {
                    356:     iqrv((VECT)dnd,dvr,rp); return;
                    357:   } else {
                    358:     asir_assert(dnd,O_N,"iqr");
                    359:     asir_assert(dvr,O_N,"iqr");
                    360:     divqrz(dnd,dvr,&a,&b);
                    361:   }
                    362:   node = mknode(2,a,b); MKLIST(*rp,node);
                    363: }
                    364:
                    365: void Pinttorat(NODE arg,LIST *rp)
                    366: {
1.2     ! noro      367:   int ret;
1.1       noro      368:   Z c,m,b,nm,dn;
                    369:   NODE node;
                    370:
                    371:   asir_assert(ARG0(arg),O_N,"inttorat");
                    372:   asir_assert(ARG1(arg),O_N,"inttorat");
                    373:   asir_assert(ARG2(arg),O_N,"inttorat");
                    374:   c = (Z)ARG0(arg); m = (Z)ARG1(arg); b = (Z)ARG2(arg);
1.2     ! noro      375:   ret = inttorat(c,m,b,&nm,&dn);
        !           376:   if ( !ret )
        !           377:     *rp = 0;
        !           378:   else {
        !           379:     node = mknode(2,nm,dn); MKLIST(*rp,node);
        !           380:   }
1.1       noro      381: }
                    382:
                    383: void Pigcd(NODE arg,Z *rp)
                    384: {
                    385:   Z n1,n2;
                    386:
                    387:   if ( argc(arg) == 1 ) {
                    388:     gcdvz((VECT)ARG0(arg),rp);
                    389:     return;
                    390:   }
                    391:   n1 = (Z)ARG0(arg); n2 = (Z)ARG1(arg);
                    392:   asir_assert(n1,O_N,"igcd");
                    393:   asir_assert(n2,O_N,"igcd");
                    394:   gcdz(n1,n2,rp);
                    395: }
                    396:
                    397: void iqrv(VECT a,Z dvr,LIST *rp)
                    398: {
                    399:   int i,n;
                    400:   VECT q,r;
                    401:   Z *b;
                    402:   NODE n0;
                    403:
                    404:   if ( !dvr )
                    405:     error("iqrv: division by 0");
                    406:   n = a->len; b = (Z *)BDY(a);
                    407:   MKVECT(q,n); MKVECT(r,n);
                    408:   for ( i = 0; i < n; i++ )
                    409:     divqrz(b[i],dvr,(Z *)&BDY(q)[i],(Z *)&BDY(r)[i]);
                    410:   n0 = mknode(2,q,r); MKLIST(*rp,n0);
                    411: }
                    412:
                    413: /*
                    414:  * gcd = GCD(a,b), ca = a/g, cb = b/g
                    415:  */
                    416:
                    417: void igcd_cofactor(Z a,Z b,Z *gcd,Z *ca,Z *cb)
                    418: {
                    419:   Z g;
                    420:
                    421:   if ( !a ) {
                    422:     if ( !b )
                    423:       error("igcd_cofactor : invalid input");
                    424:     else {
                    425:       *ca = 0;
                    426:       *cb = ONE;
                    427:       *gcd = b;
                    428:     }
                    429:   } else if ( !b ) {
                    430:     *ca = ONE;
                    431:     *cb = 0;
                    432:     *gcd = a;
                    433:   } else {
                    434:     gcdz(a,b,&g);
                    435:     divsz(a,g,ca);
                    436:     divsz(b,g,cb);
                    437:     *gcd = g;
                    438:   }
                    439: }
                    440:
                    441: void Pilcm(NODE arg,Z *rp)
                    442: {
                    443:   Z n1,n2,g,q,l;
                    444:
                    445:   n1 = (Z)ARG0(arg); n2 = (Z)ARG1(arg);
                    446:   asir_assert(n1,O_N,"ilcm");
                    447:   asir_assert(n2,O_N,"ilcm");
                    448:   if ( !n1 || !n2 )
                    449:     *rp = 0;
                    450:   else {
                    451:     gcdz(n1,n2,&g); divsz(n1,g,&q);
                    452:     mulz(q,n2,&l); absz(l,rp);
                    453:   }
                    454: }
                    455:
                    456: void Piand(NODE arg,Z *rp)
                    457: {
                    458:   mpz_t t;
                    459:   Z n1,n2;
                    460:
                    461:   n1 = (Z)ARG0(arg); n2 = (Z)ARG1(arg);
                    462:   asir_assert(n1,O_N,"iand");
                    463:   asir_assert(n2,O_N,"iand");
                    464:   if ( !n1 || !n2 ) *rp = 0;
                    465:   else {
                    466:     mpz_init(t);
                    467:     mpz_and(t,BDY(n1),BDY(n2));
                    468:     MPZTOZ(t,*rp);
                    469:   }
                    470: }
                    471:
                    472: void Pior(NODE arg,Z *rp)
                    473: {
                    474:   Z n1,n2;
                    475:   mpz_t t;
                    476:
                    477:   n1 = (Z)ARG0(arg); n2 = (Z)ARG1(arg);
                    478:   asir_assert(n1,O_N,"ior");
                    479:   asir_assert(n2,O_N,"ior");
                    480:   if ( !n1 ) *rp = n2;
                    481:   else if ( !n2 ) *rp = n1;
                    482:   else {
                    483:     mpz_init(t);
                    484:     mpz_ior(t,BDY(n1),BDY(n2));
                    485:     MPZTOZ(t,*rp);
                    486:   }
                    487: }
                    488:
                    489: void Pixor(NODE arg,Z *rp)
                    490: {
                    491:   Z n1,n2;
                    492:   mpz_t t;
                    493:
                    494:   n1 = (Z)ARG0(arg); n2 = (Z)ARG1(arg);
                    495:   asir_assert(n1,O_N,"ixor");
                    496:   asir_assert(n2,O_N,"ixor");
                    497:   if ( !n1 ) *rp = n2;
                    498:   else if ( !n2 ) *rp = n1;
                    499:   else {
                    500:     mpz_init(t);
                    501:     mpz_xor(t,BDY(n1),BDY(n2));
                    502:     MPZTOZ(t,*rp);
                    503:   }
                    504: }
                    505:
                    506: void Pishift(NODE arg,Z *rp)
                    507: {
                    508:   int i;
                    509:   Z n1,s;
                    510:   mpz_t t;
                    511:
                    512:   n1 = (Z)ARG0(arg);
                    513:   s = (Z)ARG1(arg);
                    514:   asir_assert(n1,O_N,"ixor");
                    515:   asir_assert(s,O_N,"ixor");
                    516:   bshiftz(n1,QTOS(s),rp);
                    517: }
                    518:
                    519: void isqrt(Z a,Z *r)
                    520: {
                    521:   int k;
                    522:   Z x,t,x2,xh,quo,rem;
                    523:
                    524:   if ( !a )
                    525:     *r = 0;
                    526:   else if ( UNIZ(a) )
                    527:     *r = ONE;
                    528:   else {
                    529:     k = z_bits((Q)a); /* a <= 2^k-1 */
                    530:     bshiftz(ONE,-((k>>1)+(k&1)),&x); /* a <= x^2 */
                    531:     while ( 1 ) {
                    532:       mulz(x,x,&t);
                    533:       if ( cmpz(t,a) <= 0 ) {
                    534:         *r = x; return;
                    535:       } else {
                    536:         if ( tstbitz(x,0) )
                    537:           addz(x,a,&t);
                    538:         else
                    539:           t = a;
                    540:         bshiftz(x,-1,&x2); divqrz(t,x2,&quo,&rem);
                    541:         bshiftz(x,1,&xh); addz(quo,xh,&x);
                    542:       }
                    543:     }
                    544:   }
                    545: }
                    546:
                    547: void Pup2_init_eg(Obj *rp)
                    548: {
                    549:   up2_init_eg();
                    550:   *rp = 0;
                    551: }
                    552:
                    553: void Pup2_show_eg(Obj *rp)
                    554: {
                    555:   up2_show_eg();
                    556:   *rp = 0;
                    557: }
                    558:
                    559: void Pgf2nton(NODE arg,Z *rp)
                    560: {
                    561:   if ( !ARG0(arg) )
                    562:     *rp = 0;
                    563:   else
                    564:     up2toz(((GF2N)ARG0(arg))->body,rp);
                    565: }
                    566:
                    567: void Pntogf2n(NODE arg,GF2N *rp)
                    568: {
                    569:   UP2 t;
                    570:
                    571:   ztoup2((Z)ARG0(arg),&t);
                    572:   MKGF2N(t,*rp);
                    573: }
                    574:
                    575: void Pup2_inv(NODE arg,P *rp)
                    576: {
                    577:   UP2 a,b,t;
                    578:
                    579:   ptoup2(ARG0(arg),&a);
                    580:   ptoup2(ARG1(arg),&b);
                    581:   invup2(a,b,&t);
                    582:   up2top(t,rp);
                    583: }
                    584:
                    585: void Pinv(NODE arg,Num *rp)
                    586: {
                    587:   Num n;
                    588:   Z mod;
                    589:   MQ r;
                    590:   int inv;
                    591:
                    592:   n = (Num)ARG0(arg); mod = (Z)ARG1(arg);
                    593:   asir_assert(n,O_N,"inv");
                    594:   asir_assert(mod,O_N,"inv");
                    595:   if ( !n || !mod )
                    596:     error("inv: invalid input");
                    597:   else
                    598:     switch ( NID(n) ) {
                    599:       case N_Q:
                    600:         invz((Z)n,mod,(Z *)rp);
                    601:         break;
                    602:       case N_M:
                    603:         inv = invm(CONT((MQ)n),QTOS(mod));
                    604:         STOMQ(inv,r);
                    605:         *rp = (Num)r;
                    606:         break;
                    607:       default:
                    608:         error("inv: invalid input");
                    609:     }
                    610: }
                    611:
                    612: void Pfac(NODE arg,Z *rp)
                    613: {
                    614:   asir_assert(ARG0(arg),O_N,"fac");
                    615:   factorialz(QTOS((Q)ARG0(arg)),rp);
                    616: }
                    617:
                    618: void Plrandom(NODE arg,Z *rp)
                    619: {
                    620:   asir_assert(ARG0(arg),O_N,"lrandom");
                    621:   randomz(QTOS((Q)ARG0(arg)),rp);
                    622: }
                    623:
                    624: void Prandom(NODE arg,Z *rp)
                    625: {
                    626:   unsigned int r;
                    627:
                    628: #if 0
                    629: #if defined(_PA_RISC1_1)
                    630:   r = mrand48()&BMASK;
                    631: #else
                    632:   if ( arg )
                    633:     srandom(QTOS((Q)ARG0(arg)));
                    634:   r = random()&BMASK;
                    635: #endif
                    636: #endif
                    637:   if ( arg )
                    638:     mt_sgenrand(QTOS((Q)ARG0(arg)));
                    639:   r = mt_genrand();
                    640:   UTOQ(r,*rp);
                    641: }
                    642:
                    643: #if defined(VISUAL) || defined(__MINGW32__)
                    644: void srandom(unsigned int);
                    645:
                    646: static unsigned int R_Next;
                    647:
                    648: unsigned int random() {
                    649:         if ( !R_Next )
                    650:                 R_Next = 1;
                    651:         return R_Next = (R_Next * 1103515245 + 12345);
                    652: }
                    653:
                    654: void srandom(unsigned int s)
                    655: {
                    656:     if ( s )
                    657:       R_Next = s;
                    658:         else if ( !R_Next )
                    659:             R_Next = 1;
                    660: }
                    661: #endif
                    662:
                    663: void Pprime(NODE arg,Z *rp)
                    664: {
                    665:   int i;
                    666:
                    667:   asir_assert(ARG0(arg),O_N,"prime");
                    668:   i = QTOS((Q)ARG0(arg));
                    669:   if ( i < 0 || i >= 1900 )
                    670:     *rp = 0;
                    671:   else
                    672:     STOQ(sprime[i],*rp);
                    673: }
                    674:
                    675: void Plprime(NODE arg,Z *rp)
                    676: {
                    677:   int i,p;
                    678:
                    679:   asir_assert(ARG0(arg),O_N,"lprime");
                    680:   i = QTOS((Q)ARG0(arg));
                    681:   if ( i < 0 )
                    682:     *rp = 0;
                    683:   else
                    684:     p = get_lprime(i);
                    685:   STOQ(p,*rp);
                    686: }
                    687:
                    688: extern int up_kara_mag, up_tkara_mag, up_fft_mag;
                    689:
                    690: void Pset_upfft(NODE arg,Z *rp)
                    691: {
                    692:   if ( arg ) {
                    693:     asir_assert(ARG0(arg),O_N,"set_upfft");
                    694:     up_fft_mag = QTOS((Q)ARG0(arg));
                    695:   }
                    696:   STOQ(up_fft_mag,*rp);
                    697: }
                    698:
                    699: void Pset_upkara(NODE arg,Z *rp)
                    700: {
                    701:   if ( arg ) {
                    702:     asir_assert(ARG0(arg),O_N,"set_upkara");
                    703:     up_kara_mag = QTOS((Q)ARG0(arg));
                    704:   }
                    705:   STOQ(up_kara_mag,*rp);
                    706: }
                    707:
                    708: void Pset_uptkara(NODE arg,Z *rp)
                    709: {
                    710:   if ( arg ) {
                    711:     asir_assert(ARG0(arg),O_N,"set_uptkara");
                    712:     up_tkara_mag = QTOS((Q)ARG0(arg));
                    713:   }
                    714:   STOQ(up_tkara_mag,*rp);
                    715: }
                    716:
                    717: extern int up2_kara_mag;
                    718:
                    719: void Pset_up2kara(NODE arg,Z *rp)
                    720: {
                    721:   if ( arg ) {
                    722:     asir_assert(ARG0(arg),O_N,"set_up2kara");
                    723:     up2_kara_mag = QTOS((Q)ARG0(arg));
                    724:   }
                    725:   STOQ(up2_kara_mag,*rp);
                    726: }
                    727:
                    728: void Pigcdbin(NODE arg,Z *rp)
                    729: {
                    730:   Pigcd(arg,rp);
                    731: }
                    732:
                    733: void Pigcdbmod(NODE arg,Z *rp)
                    734: {
                    735:   Pigcd(arg,rp);
                    736: }
                    737:
                    738: void Pigcdacc(NODE arg,Z *rp)
                    739: {
                    740:   Pigcd(arg,rp);
                    741: }
                    742:
                    743: void PigcdEuc(NODE arg,Z *rp)
                    744: {
                    745:   Pigcd(arg,rp);
                    746: }
                    747:
                    748: void Pigcdcntl(NODE arg,Z *rp)
                    749: {
                    750:   *rp = ONE;
                    751: }

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