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

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.5     ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2018/builtin/int.c,v 1.4 2018/09/28 08:20:27 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();
1.5     ! noro       72: void Plprime64();
1.1       noro       73:
                     74: struct ftab int_tab[] = {
                     75:   {"dp_set_mpi",Pdp_set_mpi,-1},
                     76:   {"isqrt",Pisqrt,1},
                     77:   {"idiv",Pidiv,2},
                     78:   {"irem",Pirem,2},
                     79:   {"iqr",Piqr,2},
                     80:   {"igcd",Pigcd,-2},
                     81:   {"ilcm",Pilcm,2},
                     82:   {"up2_inv",Pup2_inv,2},
                     83:   {"up2_init_eg",Pup2_init_eg,0},
                     84:   {"up2_show_eg",Pup2_show_eg,0},
                     85:   {"type_t_NB",Ptype_t_NB,2},
                     86:   {"gf2nton",Pgf2nton,1},
                     87:   {"ntogf2n",Pntogf2n,1},
                     88:   {"set_upkara",Pset_upkara,-1},
                     89:   {"set_uptkara",Pset_uptkara,-1},
                     90:   {"set_up2kara",Pset_up2kara,-1},
                     91:   {"set_upfft",Pset_upfft,-1},
                     92:   {"inv",Pinv,2},
                     93:   {"inttorat",Pinttorat,3},
                     94:   {"fac",Pfac,1},
                     95:   {"prime",Pprime,1},
                     96:   {"lprime",Plprime,1},
1.5     ! noro       97: #if SIZEOF_LONG==8
        !            98:   {"lprime64",Plprime64,1},
        !            99: #endif
1.1       noro      100:   {"random",Prandom,-1},
                    101:   {"lrandom",Plrandom,1},
                    102:   {"iand",Piand,2},
                    103:   {"ior",Pior,2},
                    104:   {"ixor",Pixor,2},
                    105:   {"ishift",Pishift,2},
                    106:   {"small_jacobi",Psmall_jacobi,2},
                    107:
                    108:   {"igcdbin",Pigcdbin,2},    /* HM@CCUT extension */
                    109:   {"igcdbmod",Pigcdbmod,2},  /* HM@CCUT extension */
                    110:   {"igcdeuc",PigcdEuc,2},    /* HM@CCUT extension */
                    111:   {"igcdacc",Pigcdacc,2},    /* HM@CCUT extension */
                    112:   {"igcdcntl",Pigcdcntl,-1},  /* HM@CCUT extension */
                    113:
                    114:   {"mt_save",Pmt_save,1},
                    115:   {"mt_load",Pmt_load,1},
                    116:   {"ntoint32",Pntoint32,1},
                    117:   {"int32ton",Pint32ton,1},
                    118:   {0,0,0},
                    119: };
                    120:
                    121: static int is_prime_small(unsigned int);
                    122: static unsigned int gcd_small(unsigned int,unsigned int);
                    123: int TypeT_NB_check(unsigned int, unsigned int);
                    124: int mpi_mag;
                    125:
                    126: void Pntoint32(NODE arg,USINT *rp)
                    127: {
                    128:   Z q,z;
                    129:   unsigned int t;
                    130:
                    131:   asir_assert(ARG0(arg),O_N,"ntoint32");
                    132:   q = (Z)ARG0(arg);
                    133:   if ( !q ) {
                    134:     MKUSINT(*rp,0);
                    135:     return;
                    136:   }
                    137:   if ( !INT(q) || !smallz(q) )
                    138:     error("ntoint32 : invalid argument");
                    139:   absz(q,&z);
1.4       noro      140:   t = ZTOS(z);
1.1       noro      141:   if ( sgnz(q) < 0 )
                    142:     t = -(int)t;
                    143:   MKUSINT(*rp,t);
                    144: }
                    145:
                    146: void Pint32ton(NODE arg,Z *rp)
                    147: {
                    148:   int t;
                    149:
                    150:   asir_assert(ARG0(arg),O_USINT,"int32ton");
                    151:   t = (int)BDY((USINT)ARG0(arg));
1.4       noro      152:   STOZ(t,*rp);
1.1       noro      153: }
                    154:
                    155: void Pdp_set_mpi(NODE arg,Z *rp)
                    156: {
                    157:   if ( arg ) {
                    158:     asir_assert(ARG0(arg),O_N,"dp_set_mpi");
1.4       noro      159:     mpi_mag = ZTOS((Q)ARG0(arg));
1.1       noro      160:   }
1.4       noro      161:   STOZ(mpi_mag,*rp);
1.1       noro      162: }
                    163:
                    164: void Psmall_jacobi(NODE arg,Z *rp)
                    165: {
                    166:   Z a,m;
                    167:   int a0,m0,s;
                    168:
                    169:   a = (Z)ARG0(arg);
                    170:   m = (Z)ARG1(arg);
                    171:   asir_assert(a,O_N,"small_jacobi");
                    172:   asir_assert(m,O_N,"small_jacobi");
                    173:   if ( !a )
                    174:      *rp = ONE;
                    175:   else if ( !m || !INT(m) || !INT(a)
                    176:     || !smallz(m) || !smallz(a) || sgnz(m) < 0 || evenz(m) )
                    177:     error("small_jacobi : invalid input");
                    178:   else {
1.4       noro      179:     a0 = ZTOS(a); m0 = ZTOS(m);
1.1       noro      180:     s = small_jacobi(a0,m0);
1.4       noro      181:     STOZ(s,*rp);
1.1       noro      182:   }
                    183: }
                    184:
                    185: int small_jacobi(int a,int m)
                    186: {
                    187:   int m4,m8,a4,j1,i,s;
                    188:
                    189:   a %= m;
                    190:   if ( a == 0 || a == 1 )
                    191:     return 1;
                    192:   else if ( a < 0 ) {
                    193:     j1 = small_jacobi(-a,m);
                    194:     m4 = m%4;
                    195:     return m4==1?j1:-j1;
                    196:   } else {
                    197:     for ( i = 0; a && !(a&1); i++, a >>= 1 );
                    198:     if ( i&1 ) {
                    199:       m8 = m%8;
                    200:       s = (m8==1||m8==7) ? 1 : -1;
                    201:     } else
                    202:       s = 1;
                    203:     /* a, m are odd */
                    204:     j1 = small_jacobi(m%a,a);
                    205:     m4 = m%4; a4 = a%4;
                    206:     s *= (m4==1||a4==1) ? 1 : -1;
                    207:     return j1*s;
                    208:   }
                    209: }
                    210:
                    211: void Ptype_t_NB(NODE arg,Z *rp)
                    212: {
1.4       noro      213:   if ( TypeT_NB_check(ZTOS((Q)ARG0(arg)),ZTOS((Q)ARG1(arg))))
1.1       noro      214:     *rp = ONE;
                    215:   else
                    216:     *rp = 0;
                    217: }
                    218:
                    219: int TypeT_NB_check(unsigned int m, unsigned int t)
                    220: {
                    221:   unsigned int p,k,u,h,d;
                    222:
                    223:   if ( !(m%8) )
                    224:     return 0;
                    225:   p = t*m+1;
                    226:   if ( !is_prime_small(p) )
                    227:     return 0;
                    228:   for ( k = 1, u = 2%p; ; k++ )
                    229:     if ( u == 1 )
                    230:       break;
                    231:     else
                    232:       u = (2*u)%p;
                    233:   h = t*m/k;
                    234:   d = gcd_small(h,m);
                    235:   return d == 1 ? 1 : 0;
                    236: }
                    237:
                    238: /*
                    239:  * a simple prime checker
                    240:  * return value: 1  ---  prime number
                    241:  *               0  ---  composite number
                    242:  */
                    243:
                    244: static int is_prime_small(unsigned int a)
                    245: {
                    246:   unsigned int b,t,i;
                    247:
                    248:   if ( !(a%2) ) return 0;
                    249:   for ( t = a, i = 0; t; i++, t >>= 1 );
                    250:   /* b >= sqrt(a) */
                    251:   b = 1<<((i+1)/2);
                    252:
                    253:   /* divisibility test by all odd numbers <= b */
                    254:   for ( i = 3; i <= b; i += 2 )
                    255:     if ( !(a%i) )
                    256:       return 0;
                    257:   return 1;
                    258: }
                    259:
                    260: /*
                    261:  * gcd for unsigned int as integers
                    262:  * return value: GCD(a,b)
                    263:  *
                    264:  */
                    265:
                    266:
                    267: static unsigned int gcd_small(unsigned int a,unsigned int b)
                    268: {
                    269:   unsigned int t;
                    270:
                    271:   if ( b > a ) {
                    272:     t = a; a = b; b = t;
                    273:   }
                    274:   /* Euclid's algorithm */
                    275:   while ( 1 )
                    276:     if ( !(t = a%b) ) return b;
                    277:     else {
                    278:       a = b; b = t;
                    279:     }
                    280: }
                    281:
                    282: void Pmt_save(NODE arg,Z *rp)
                    283: {
                    284:   int ret;
                    285:
                    286:   ret = mt_save(BDY((STRING)ARG0(arg)));
1.4       noro      287:   STOZ(ret,*rp);
1.1       noro      288: }
                    289:
                    290: void Pmt_load(NODE arg,Z *rp)
                    291: {
                    292:   int ret;
                    293:
                    294:   ret = mt_load(BDY((STRING)ARG0(arg)));
1.4       noro      295:   STOZ(ret,*rp);
1.1       noro      296: }
                    297:
                    298: void isqrt(Z a,Z *r);
                    299:
                    300: void Pisqrt(NODE arg,Z *rp)
                    301: {
                    302:   Z a;
                    303:   Z r;
                    304:
                    305:   a = (Z)ARG0(arg);
                    306:   asir_assert(a,O_N,"isqrt");
                    307:   if ( !a )
                    308:     *rp = 0;
                    309:   else if ( sgnz(a) < 0 )
                    310:     error("isqrt : negative argument");
                    311:   else {
                    312:     isqrt(a,rp);
                    313:   }
                    314: }
                    315:
                    316: void Pidiv(NODE arg,Z *rp)
                    317: {
                    318:   Z r;
                    319:   Z dnd,dvr;
                    320:
                    321:   dnd = (Z)ARG0(arg); dvr = (Z)ARG1(arg);
                    322:   asir_assert(dnd,O_N,"idiv");
                    323:   asir_assert(dvr,O_N,"idiv");
                    324:   if ( !dvr )
                    325:     error("idiv: division by 0");
                    326:   else if ( !dnd )
                    327:     *rp = 0;
                    328:   else
                    329:     divqrz(dnd,dvr,rp,&r);
                    330: }
                    331:
                    332: void Pirem(NODE arg,Z *rp)
                    333: {
                    334:   Z q,dnd,dvr;
                    335:
                    336:   dnd = (Z)ARG0(arg); dvr = (Z)ARG1(arg);
                    337:   asir_assert(dnd,O_N,"irem");
                    338:   asir_assert(dvr,O_N,"irem");
                    339:   if ( !dvr )
                    340:     error("irem: division by 0");
                    341:   else if ( !dnd )
                    342:     *rp = 0;
                    343:   else
                    344:     divqrz(dnd,dvr,&q,rp);
                    345: }
                    346:
                    347: void iqrv(VECT a,Z dvr,LIST *rp);
                    348:
                    349: void Piqr(NODE arg,LIST *rp)
                    350: {
                    351:   Z dnd,dvr,a,b;
                    352:   NODE node;
                    353:
                    354:   dnd = (Z)ARG0(arg); dvr = (Z)ARG1(arg);
                    355:   if ( !dvr )
                    356:     error("iqr: division by 0");
                    357:   else if ( !dnd )
                    358:     a = b = 0;
                    359:   else if ( OID(dnd) == O_VECT ) {
                    360:     iqrv((VECT)dnd,dvr,rp); return;
                    361:   } else {
                    362:     asir_assert(dnd,O_N,"iqr");
                    363:     asir_assert(dvr,O_N,"iqr");
                    364:     divqrz(dnd,dvr,&a,&b);
                    365:   }
                    366:   node = mknode(2,a,b); MKLIST(*rp,node);
                    367: }
                    368:
                    369: void Pinttorat(NODE arg,LIST *rp)
                    370: {
1.2       noro      371:   int ret;
1.3       noro      372:   Z c,m,t,b,nm,dn;
1.1       noro      373:   NODE node;
                    374:
                    375:   asir_assert(ARG0(arg),O_N,"inttorat");
                    376:   asir_assert(ARG1(arg),O_N,"inttorat");
                    377:   asir_assert(ARG2(arg),O_N,"inttorat");
                    378:   c = (Z)ARG0(arg); m = (Z)ARG1(arg); b = (Z)ARG2(arg);
1.3       noro      379:   remz(c,m,&t); c = t;
1.2       noro      380:   ret = inttorat(c,m,b,&nm,&dn);
                    381:   if ( !ret )
                    382:     *rp = 0;
                    383:   else {
                    384:     node = mknode(2,nm,dn); MKLIST(*rp,node);
                    385:   }
1.1       noro      386: }
                    387:
                    388: void Pigcd(NODE arg,Z *rp)
                    389: {
                    390:   Z n1,n2;
                    391:
                    392:   if ( argc(arg) == 1 ) {
                    393:     gcdvz((VECT)ARG0(arg),rp);
                    394:     return;
                    395:   }
                    396:   n1 = (Z)ARG0(arg); n2 = (Z)ARG1(arg);
                    397:   asir_assert(n1,O_N,"igcd");
                    398:   asir_assert(n2,O_N,"igcd");
                    399:   gcdz(n1,n2,rp);
                    400: }
                    401:
                    402: void iqrv(VECT a,Z dvr,LIST *rp)
                    403: {
                    404:   int i,n;
                    405:   VECT q,r;
                    406:   Z *b;
                    407:   NODE n0;
                    408:
                    409:   if ( !dvr )
                    410:     error("iqrv: division by 0");
                    411:   n = a->len; b = (Z *)BDY(a);
                    412:   MKVECT(q,n); MKVECT(r,n);
                    413:   for ( i = 0; i < n; i++ )
                    414:     divqrz(b[i],dvr,(Z *)&BDY(q)[i],(Z *)&BDY(r)[i]);
                    415:   n0 = mknode(2,q,r); MKLIST(*rp,n0);
                    416: }
                    417:
                    418: /*
                    419:  * gcd = GCD(a,b), ca = a/g, cb = b/g
                    420:  */
                    421:
                    422: void igcd_cofactor(Z a,Z b,Z *gcd,Z *ca,Z *cb)
                    423: {
                    424:   Z g;
                    425:
                    426:   if ( !a ) {
                    427:     if ( !b )
                    428:       error("igcd_cofactor : invalid input");
                    429:     else {
                    430:       *ca = 0;
                    431:       *cb = ONE;
                    432:       *gcd = b;
                    433:     }
                    434:   } else if ( !b ) {
                    435:     *ca = ONE;
                    436:     *cb = 0;
                    437:     *gcd = a;
                    438:   } else {
                    439:     gcdz(a,b,&g);
                    440:     divsz(a,g,ca);
                    441:     divsz(b,g,cb);
                    442:     *gcd = g;
                    443:   }
                    444: }
                    445:
                    446: void Pilcm(NODE arg,Z *rp)
                    447: {
                    448:   Z n1,n2,g,q,l;
                    449:
                    450:   n1 = (Z)ARG0(arg); n2 = (Z)ARG1(arg);
                    451:   asir_assert(n1,O_N,"ilcm");
                    452:   asir_assert(n2,O_N,"ilcm");
                    453:   if ( !n1 || !n2 )
                    454:     *rp = 0;
                    455:   else {
                    456:     gcdz(n1,n2,&g); divsz(n1,g,&q);
                    457:     mulz(q,n2,&l); absz(l,rp);
                    458:   }
                    459: }
                    460:
                    461: void Piand(NODE arg,Z *rp)
                    462: {
                    463:   mpz_t t;
                    464:   Z n1,n2;
                    465:
                    466:   n1 = (Z)ARG0(arg); n2 = (Z)ARG1(arg);
                    467:   asir_assert(n1,O_N,"iand");
                    468:   asir_assert(n2,O_N,"iand");
                    469:   if ( !n1 || !n2 ) *rp = 0;
                    470:   else {
                    471:     mpz_init(t);
                    472:     mpz_and(t,BDY(n1),BDY(n2));
                    473:     MPZTOZ(t,*rp);
                    474:   }
                    475: }
                    476:
                    477: void Pior(NODE arg,Z *rp)
                    478: {
                    479:   Z n1,n2;
                    480:   mpz_t t;
                    481:
                    482:   n1 = (Z)ARG0(arg); n2 = (Z)ARG1(arg);
                    483:   asir_assert(n1,O_N,"ior");
                    484:   asir_assert(n2,O_N,"ior");
                    485:   if ( !n1 ) *rp = n2;
                    486:   else if ( !n2 ) *rp = n1;
                    487:   else {
                    488:     mpz_init(t);
                    489:     mpz_ior(t,BDY(n1),BDY(n2));
                    490:     MPZTOZ(t,*rp);
                    491:   }
                    492: }
                    493:
                    494: void Pixor(NODE arg,Z *rp)
                    495: {
                    496:   Z n1,n2;
                    497:   mpz_t t;
                    498:
                    499:   n1 = (Z)ARG0(arg); n2 = (Z)ARG1(arg);
                    500:   asir_assert(n1,O_N,"ixor");
                    501:   asir_assert(n2,O_N,"ixor");
                    502:   if ( !n1 ) *rp = n2;
                    503:   else if ( !n2 ) *rp = n1;
                    504:   else {
                    505:     mpz_init(t);
                    506:     mpz_xor(t,BDY(n1),BDY(n2));
                    507:     MPZTOZ(t,*rp);
                    508:   }
                    509: }
                    510:
                    511: void Pishift(NODE arg,Z *rp)
                    512: {
                    513:   int i;
                    514:   Z n1,s;
                    515:   mpz_t t;
                    516:
                    517:   n1 = (Z)ARG0(arg);
                    518:   s = (Z)ARG1(arg);
                    519:   asir_assert(n1,O_N,"ixor");
                    520:   asir_assert(s,O_N,"ixor");
1.4       noro      521:   bshiftz(n1,ZTOS(s),rp);
1.1       noro      522: }
                    523:
                    524: void isqrt(Z a,Z *r)
                    525: {
                    526:   int k;
                    527:   Z x,t,x2,xh,quo,rem;
                    528:
                    529:   if ( !a )
                    530:     *r = 0;
                    531:   else if ( UNIZ(a) )
                    532:     *r = ONE;
                    533:   else {
                    534:     k = z_bits((Q)a); /* a <= 2^k-1 */
                    535:     bshiftz(ONE,-((k>>1)+(k&1)),&x); /* a <= x^2 */
                    536:     while ( 1 ) {
                    537:       mulz(x,x,&t);
                    538:       if ( cmpz(t,a) <= 0 ) {
                    539:         *r = x; return;
                    540:       } else {
                    541:         if ( tstbitz(x,0) )
                    542:           addz(x,a,&t);
                    543:         else
                    544:           t = a;
                    545:         bshiftz(x,-1,&x2); divqrz(t,x2,&quo,&rem);
                    546:         bshiftz(x,1,&xh); addz(quo,xh,&x);
                    547:       }
                    548:     }
                    549:   }
                    550: }
                    551:
                    552: void Pup2_init_eg(Obj *rp)
                    553: {
                    554:   up2_init_eg();
                    555:   *rp = 0;
                    556: }
                    557:
                    558: void Pup2_show_eg(Obj *rp)
                    559: {
                    560:   up2_show_eg();
                    561:   *rp = 0;
                    562: }
                    563:
                    564: void Pgf2nton(NODE arg,Z *rp)
                    565: {
                    566:   if ( !ARG0(arg) )
                    567:     *rp = 0;
                    568:   else
                    569:     up2toz(((GF2N)ARG0(arg))->body,rp);
                    570: }
                    571:
                    572: void Pntogf2n(NODE arg,GF2N *rp)
                    573: {
                    574:   UP2 t;
                    575:
                    576:   ztoup2((Z)ARG0(arg),&t);
                    577:   MKGF2N(t,*rp);
                    578: }
                    579:
                    580: void Pup2_inv(NODE arg,P *rp)
                    581: {
                    582:   UP2 a,b,t;
                    583:
                    584:   ptoup2(ARG0(arg),&a);
                    585:   ptoup2(ARG1(arg),&b);
                    586:   invup2(a,b,&t);
                    587:   up2top(t,rp);
                    588: }
                    589:
                    590: void Pinv(NODE arg,Num *rp)
                    591: {
                    592:   Num n;
                    593:   Z mod;
                    594:   MQ r;
                    595:   int inv;
                    596:
                    597:   n = (Num)ARG0(arg); mod = (Z)ARG1(arg);
                    598:   asir_assert(n,O_N,"inv");
                    599:   asir_assert(mod,O_N,"inv");
                    600:   if ( !n || !mod )
                    601:     error("inv: invalid input");
                    602:   else
                    603:     switch ( NID(n) ) {
                    604:       case N_Q:
                    605:         invz((Z)n,mod,(Z *)rp);
                    606:         break;
                    607:       case N_M:
1.4       noro      608:         inv = invm(CONT((MQ)n),ZTOS(mod));
1.1       noro      609:         STOMQ(inv,r);
                    610:         *rp = (Num)r;
                    611:         break;
                    612:       default:
                    613:         error("inv: invalid input");
                    614:     }
                    615: }
                    616:
                    617: void Pfac(NODE arg,Z *rp)
                    618: {
                    619:   asir_assert(ARG0(arg),O_N,"fac");
1.4       noro      620:   factorialz(ZTOS((Q)ARG0(arg)),rp);
1.1       noro      621: }
                    622:
                    623: void Plrandom(NODE arg,Z *rp)
                    624: {
                    625:   asir_assert(ARG0(arg),O_N,"lrandom");
1.4       noro      626:   randomz(ZTOS((Q)ARG0(arg)),rp);
1.1       noro      627: }
                    628:
                    629: void Prandom(NODE arg,Z *rp)
                    630: {
                    631:   unsigned int r;
                    632:
                    633: #if 0
                    634: #if defined(_PA_RISC1_1)
                    635:   r = mrand48()&BMASK;
                    636: #else
                    637:   if ( arg )
1.4       noro      638:     srandom(ZTOS((Q)ARG0(arg)));
1.1       noro      639:   r = random()&BMASK;
                    640: #endif
                    641: #endif
                    642:   if ( arg )
1.4       noro      643:     mt_sgenrand(ZTOS((Q)ARG0(arg)));
1.1       noro      644:   r = mt_genrand();
1.4       noro      645:   UTOZ(r,*rp);
1.1       noro      646: }
                    647:
                    648: #if defined(VISUAL) || defined(__MINGW32__)
                    649: void srandom(unsigned int);
                    650:
                    651: static unsigned int R_Next;
                    652:
                    653: unsigned int random() {
                    654:         if ( !R_Next )
                    655:                 R_Next = 1;
                    656:         return R_Next = (R_Next * 1103515245 + 12345);
                    657: }
                    658:
                    659: void srandom(unsigned int s)
                    660: {
                    661:     if ( s )
                    662:       R_Next = s;
                    663:         else if ( !R_Next )
                    664:             R_Next = 1;
                    665: }
                    666: #endif
                    667:
                    668: void Pprime(NODE arg,Z *rp)
                    669: {
                    670:   int i;
                    671:
                    672:   asir_assert(ARG0(arg),O_N,"prime");
1.4       noro      673:   i = ZTOS((Q)ARG0(arg));
1.1       noro      674:   if ( i < 0 || i >= 1900 )
                    675:     *rp = 0;
                    676:   else
1.4       noro      677:     STOZ(sprime[i],*rp);
1.1       noro      678: }
                    679:
                    680: void Plprime(NODE arg,Z *rp)
                    681: {
                    682:   int i,p;
                    683:
                    684:   asir_assert(ARG0(arg),O_N,"lprime");
1.4       noro      685:   i = ZTOS((Q)ARG0(arg));
1.1       noro      686:   if ( i < 0 )
                    687:     *rp = 0;
                    688:   else
                    689:     p = get_lprime(i);
1.4       noro      690:   STOZ(p,*rp);
1.1       noro      691: }
                    692:
1.5     ! noro      693: #if SIZEOF_LONG==8
        !           694: void Plprime64(NODE arg,Z *rp)
        !           695: {
        !           696:   int i;
        !           697:   mp_limb_t p;
        !           698:
        !           699:   asir_assert(ARG0(arg),O_N,"lprime64");
        !           700:   i = ZTOS((Q)ARG0(arg));
        !           701:   if ( i < 0 )
        !           702:     *rp = 0;
        !           703:   else
        !           704:     p = get_lprime64(i);
        !           705:   STOZ(p,*rp);
        !           706: }
        !           707: #endif
        !           708:
1.1       noro      709: extern int up_kara_mag, up_tkara_mag, up_fft_mag;
                    710:
                    711: void Pset_upfft(NODE arg,Z *rp)
                    712: {
                    713:   if ( arg ) {
                    714:     asir_assert(ARG0(arg),O_N,"set_upfft");
1.4       noro      715:     up_fft_mag = ZTOS((Q)ARG0(arg));
1.1       noro      716:   }
1.4       noro      717:   STOZ(up_fft_mag,*rp);
1.1       noro      718: }
                    719:
                    720: void Pset_upkara(NODE arg,Z *rp)
                    721: {
                    722:   if ( arg ) {
                    723:     asir_assert(ARG0(arg),O_N,"set_upkara");
1.4       noro      724:     up_kara_mag = ZTOS((Q)ARG0(arg));
1.1       noro      725:   }
1.4       noro      726:   STOZ(up_kara_mag,*rp);
1.1       noro      727: }
                    728:
                    729: void Pset_uptkara(NODE arg,Z *rp)
                    730: {
                    731:   if ( arg ) {
                    732:     asir_assert(ARG0(arg),O_N,"set_uptkara");
1.4       noro      733:     up_tkara_mag = ZTOS((Q)ARG0(arg));
1.1       noro      734:   }
1.4       noro      735:   STOZ(up_tkara_mag,*rp);
1.1       noro      736: }
                    737:
                    738: extern int up2_kara_mag;
                    739:
                    740: void Pset_up2kara(NODE arg,Z *rp)
                    741: {
                    742:   if ( arg ) {
                    743:     asir_assert(ARG0(arg),O_N,"set_up2kara");
1.4       noro      744:     up2_kara_mag = ZTOS((Q)ARG0(arg));
1.1       noro      745:   }
1.4       noro      746:   STOZ(up2_kara_mag,*rp);
1.1       noro      747: }
                    748:
                    749: void Pigcdbin(NODE arg,Z *rp)
                    750: {
                    751:   Pigcd(arg,rp);
                    752: }
                    753:
                    754: void Pigcdbmod(NODE arg,Z *rp)
                    755: {
                    756:   Pigcd(arg,rp);
                    757: }
                    758:
                    759: void Pigcdacc(NODE arg,Z *rp)
                    760: {
                    761:   Pigcd(arg,rp);
                    762: }
                    763:
                    764: void PigcdEuc(NODE arg,Z *rp)
                    765: {
                    766:   Pigcd(arg,rp);
                    767: }
                    768:
                    769: void Pigcdcntl(NODE arg,Z *rp)
                    770: {
                    771:   *rp = ONE;
                    772: }

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