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

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

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