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