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