![[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.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.10 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/int.c,v 1.9 2001/06/07 04:54:38 noro 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[] = {
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},
1.8 murao 103:
1.1 noro 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 */
1.8 murao 109:
1.1 noro 110: {"mt_save",Pmt_save,1},
111: {"mt_load",Pmt_load,1},
1.2 noro 112: {"ntoint32",Pntoint32,1},
113: {"int32ton",Pint32ton,1},
1.1 noro 114: {0,0,0},
115: };
116:
117: static int is_prime_small(unsigned int);
118: static unsigned int gcd_small(unsigned int,unsigned int);
119: int TypeT_NB_check(unsigned int, unsigned int);
120: int mpi_mag;
1.2 noro 121:
1.10 ! noro 122: void Pntoint32(NODE arg,USINT *rp)
1.2 noro 123: {
124: Q q;
1.3 noro 125: unsigned int t;
1.2 noro 126:
127: asir_assert(ARG0(arg),O_N,"ntoint32");
128: q = (Q)ARG0(arg);
129: if ( !q ) {
130: MKUSINT(*rp,0);
131: return;
132: }
1.3 noro 133: if ( NID(q)!=N_Q || !INT(q) || PL(NM(q))>1 )
1.2 noro 134: error("ntoint32 : invalid argument");
1.3 noro 135: t = BD(NM(q))[0];
136: if ( SGN(q) < 0 )
1.10 ! noro 137: t = -(int)t;
1.3 noro 138: MKUSINT(*rp,t);
1.2 noro 139: }
140:
1.10 ! noro 141: void Pint32ton(NODE arg,Q *rp)
1.2 noro 142: {
1.3 noro 143: int t;
1.2 noro 144:
145: asir_assert(ARG0(arg),O_USINT,"int32ton");
1.3 noro 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: {
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);
157: }
158:
1.10 ! noro 159: void Psmall_jacobi(NODE arg,Q *rp)
1.1 noro 160: {
161: Q a,m;
162: int a0,m0,s;
163:
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: }
178: }
179:
1.10 ! noro 180: int small_jacobi(int a,int m)
1.1 noro 181: {
182: int m4,m8,a4,j1,i,s;
183:
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: }
204: }
205:
1.10 ! noro 206: void Ptype_t_NB(NODE arg,Q *rp)
1.1 noro 207: {
208: if ( TypeT_NB_check(QTOS((Q)ARG0(arg)),QTOS((Q)ARG1(arg))))
209: *rp = ONE;
210: else
211: *rp = 0;
212: }
213:
214: int TypeT_NB_check(unsigned int m, unsigned int t)
215: {
216: unsigned int p,k,u,h,d;
217:
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;
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: {
241: unsigned int b,t,i;
242:
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;
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: {
264: unsigned int t;
265:
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: }
275: }
276:
1.10 ! noro 277: void Pmt_save(NODE arg,Q *rp)
1.1 noro 278: {
279: int ret;
280:
281: ret = mt_save(BDY((STRING)ARG0(arg)));
282: STOQ(ret,*rp);
283: }
284:
1.10 ! noro 285: void Pmt_load(NODE arg,Q *rp)
1.1 noro 286: {
287: int ret;
288:
289: ret = mt_load(BDY((STRING)ARG0(arg)));
290: STOQ(ret,*rp);
291: }
292:
1.10 ! noro 293: void Pisqrt(NODE arg,Q *rp)
1.1 noro 294: {
295: Q a;
296: N r;
297:
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: }
308: }
309:
1.10 ! noro 310: void Pidiv(NODE arg,Obj *rp)
1.1 noro 311: {
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: }
327: }
328:
1.10 ! noro 329: void Pirem(NODE arg,Obj *rp)
1.1 noro 330: {
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: }
346: }
347:
1.10 ! noro 348: void Piqr(NODE arg,LIST *rp)
1.1 noro 349: {
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);
370: }
371:
1.10 ! noro 372: void Pinttorat(NODE arg,LIST *rp)
1.1 noro 373: {
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: }
409: }
410:
1.10 ! noro 411: void Pigcd(NODE arg,Q *rp)
1.1 noro 412: {
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: }
431: }
432:
1.10 ! noro 433: int comp_n(N *a,N *b)
1.1 noro 434: {
435: return cmpn(*a,*b);
436: }
437:
1.10 ! noro 438: void iqrv(VECT a,Q dvr,LIST *rp)
1.1 noro 439: {
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: {
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: {
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);
511: for ( i = 2; i < n; i++ ) {
512: divn(c[i],g,&q,&r);
513: gcdn(g,r,>);
514: if ( !cmpn(g,gt) ) {
515: for ( j = i+1, nz = 0; j < n; j++ ) {
516: divn(c[j],g,&q,&r); c[j] = r;
517: if ( r )
518: nz = 1;
519: }
520: } else
521: g = gt;
522: }
523: NTOQ(g,1,*rp);
524: }
525:
1.10 ! noro 526: void igcdv_estimate(VECT a,Q *rp)
1.1 noro 527: {
528: int n,i,m;
529: N s,t,u,g;
530: Q *q;
531:
532: n = a->len; q = (Q *)a->body;
533: if ( n == 1 ) {
534: NTOQ(NM(q[0]),1,*rp); return;
535: }
536:
537: m = n/2;
538: for ( i = 0 , s = 0; i < m; i++ ) {
539: addn(s,NM(q[i]),&u); s = u;
540: }
541: for ( t = 0; i < n; i++ ) {
542: addn(t,NM(q[i]),&u); t = u;
543: }
544: gcdn(s,t,&g); NTOQ(g,1,*rp);
545: }
546:
1.10 ! noro 547: void Pilcm(NODE arg,Obj *rp)
1.1 noro 548: {
549: N g,q,r,l;
550: Q n1,n2,a;
551:
552: n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
553: asir_assert(n1,O_N,"ilcm");
554: asir_assert(n2,O_N,"ilcm");
555: if ( !n1 || !n2 )
556: *rp = 0;
557: else {
558: gcdn(NM(n1),NM(n2),&g); divn(NM(n1),g,&q,&r);
559: muln(q,NM(n2),&l); NTOQ(l,1,a); *rp = (Obj)a;
560: }
561: }
562:
1.10 ! noro 563: void Piand(NODE arg,Q *rp)
1.1 noro 564: {
565: N g;
566: Q n1,n2,a;
567:
568: n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
569: asir_assert(n1,O_N,"iand");
570: asir_assert(n2,O_N,"iand");
571: if ( !n1 || !n2 )
572: *rp = 0;
573: else {
574: iand(NM(n1),NM(n2),&g);
575: NTOQ(g,1,a); *rp = a;
576: }
577: }
578:
1.10 ! noro 579: void Pior(NODE arg,Q *rp)
1.1 noro 580: {
581: N g;
582: Q n1,n2,a;
583:
584: n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
585: asir_assert(n1,O_N,"ior");
586: asir_assert(n2,O_N,"ior");
587: if ( !n1 )
588: *rp = n2;
589: else if ( !n2 )
590: *rp = n1;
591: else {
592: ior(NM(n1),NM(n2),&g);
593: NTOQ(g,1,a); *rp = a;
594: }
595: }
596:
1.10 ! noro 597: void Pixor(NODE arg,Q *rp)
1.1 noro 598: {
599: N g;
600: Q n1,n2,a;
601:
602: n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
603: asir_assert(n1,O_N,"ixor");
604: asir_assert(n2,O_N,"ixor");
605: if ( !n1 )
606: *rp = n2;
607: else if ( !n2 )
608: *rp = n1;
609: else {
610: ixor(NM(n1),NM(n2),&g);
611: NTOQ(g,1,a); *rp = a;
612: }
613: }
614:
1.10 ! noro 615: void Pishift(NODE arg,Q *rp)
1.1 noro 616: {
617: N g;
618: Q n1,s,a;
619:
620: n1 = (Q)ARG0(arg); s = (Q)ARG1(arg);
621: asir_assert(n1,O_N,"ixor");
622: asir_assert(s,O_N,"ixor");
623: if ( !n1 )
624: *rp = 0;
625: else if ( !s )
626: *rp = n1;
627: else {
628: bshiftn(NM(n1),QTOS(s),&g);
629: NTOQ(g,1,a); *rp = a;
630: }
631: }
632:
1.10 ! noro 633: void isqrt(N a,N *r)
1.1 noro 634: {
635: int k;
636: N x,t,x2,xh,quo,rem;
637:
638: if ( !a )
639: *r = 0;
640: else if ( UNIN(a) )
641: *r = ONEN;
642: else {
643: k = n_bits(a); /* a <= 2^k-1 */
644: bshiftn(ONEN,-((k>>1)+(k&1)),&x); /* a <= x^2 */
645: while ( 1 ) {
646: pwrn(x,2,&t);
647: if ( cmpn(t,a) <= 0 ) {
648: *r = x; return;
649: } else {
650: if ( BD(x)[0] & 1 )
651: addn(x,a,&t);
652: else
653: t = a;
654: bshiftn(x,-1,&x2); divn(t,x2,&quo,&rem);
655: bshiftn(x,1,&xh); addn(quo,xh,&x);
656: }
657: }
658: }
659: }
660:
1.10 ! noro 661: void iand(N n1,N n2,N *r)
1.1 noro 662: {
663: int d1,d2,d,i;
664: N nr;
665: int *p1,*p2,*pr;
666:
667: d1 = PL(n1); d2 = PL(n2);
668: d = MIN(d1,d2);
669: nr = NALLOC(d);
670: for ( i = 0, p1 = BD(n1), p2 = BD(n2), pr = BD(nr); i < d; i++ )
671: pr[i] = p1[i] & p2[i];
672: for ( i = d-1; i >= 0 && !pr[i]; i-- );
673: if ( i < 0 )
674: *r = 0;
675: else {
676: PL(nr) = i+1; *r = nr;
677: }
678: }
679:
1.10 ! noro 680: void ior(N n1,N n2,N *r)
1.1 noro 681: {
682: int d1,d2,i;
683: N nr;
684: int *p1,*p2,*pr;
685:
686: if ( PL(n1) < PL(n2) ) {
687: nr = n1; n1 = n2; n2 = nr;
688: }
689: d1 = PL(n1); d2 = PL(n2);
690: *r = nr = NALLOC(d1);
691: for ( i = 0, p1 = BD(n1), p2 = BD(n2), pr = BD(nr); i < d2; i++ )
692: pr[i] = p1[i] | p2[i];
693: for ( ; i < d1; i++ )
694: pr[i] = p1[i];
695: for ( i = d1-1; i >= 0 && !pr[i]; i-- );
696: if ( i < 0 )
697: *r = 0;
698: else {
699: PL(nr) = i+1; *r = nr;
700: }
701: }
702:
1.10 ! noro 703: void ixor(N n1,N n2,N *r)
1.1 noro 704: {
705: int d1,d2,i;
706: N nr;
707: int *p1,*p2,*pr;
708:
709: if ( PL(n1) < PL(n2) ) {
710: nr = n1; n1 = n2; n2 = nr;
711: }
712: d1 = PL(n1); d2 = PL(n2);
713: *r = nr = NALLOC(d1);
714: for ( i = 0, p1 = BD(n1), p2 = BD(n2), pr = BD(nr); i < d2; i++ )
715: pr[i] = p1[i] ^ p2[i];
716: for ( ; i < d1; i++ )
717: pr[i] = p1[i];
718: for ( i = d1-1; i >= 0 && !pr[i]; i-- );
719: if ( i < 0 )
720: *r = 0;
721: else {
722: PL(nr) = i+1; *r = nr;
723: }
724: }
725:
1.10 ! noro 726: void Pup2_init_eg(Obj *rp)
1.1 noro 727: {
728: up2_init_eg();
729: *rp = 0;
730: }
731:
1.10 ! noro 732: void Pup2_show_eg(Obj *rp)
1.1 noro 733: {
734: up2_show_eg();
735: *rp = 0;
736: }
737:
1.10 ! noro 738: void Pgf2nton(NODE arg,Q *rp)
1.1 noro 739: {
740: if ( !ARG0(arg) )
741: *rp = 0;
742: else
743: up2ton(((GF2N)ARG0(arg))->body,rp);
744: }
745:
1.10 ! noro 746: void Pntogf2n(NODE arg,GF2N *rp)
1.1 noro 747: {
748: UP2 t;
749:
750: ntoup2(ARG0(arg),&t);
751: MKGF2N(t,*rp);
752: }
753:
1.10 ! noro 754: void Pup2_inv(NODE arg,P *rp)
1.1 noro 755: {
756: UP2 a,b,t;
757:
758: ptoup2(ARG0(arg),&a);
759: ptoup2(ARG1(arg),&b);
760: invup2(a,b,&t);
761: up2top(t,rp);
762: }
763:
1.10 ! noro 764: void Pinv(NODE arg,Num *rp)
1.1 noro 765: {
766: Num n;
767: Q mod;
768: MQ r;
769: int inv;
770:
771: n = (Num)ARG0(arg); mod = (Q)ARG1(arg);
772: asir_assert(n,O_N,"inv");
773: asir_assert(mod,O_N,"inv");
774: if ( !n || !mod )
775: error("inv: invalid input");
776: else
777: switch ( NID(n) ) {
778: case N_Q:
779: invl((Q)n,mod,(Q *)rp);
780: break;
781: case N_M:
782: inv = invm(CONT((MQ)n),QTOS(mod));
783: STOMQ(inv,r);
784: *rp = (Num)r;
785: break;
786: default:
787: error("inv: invalid input");
788: }
789: }
790:
1.10 ! noro 791: void Pfac(NODE arg,Q *rp)
1.1 noro 792: {
793: asir_assert(ARG0(arg),O_N,"fac");
794: factorial(QTOS((Q)ARG0(arg)),rp);
795: }
796:
1.10 ! noro 797: void Plrandom(NODE arg,Q *rp)
1.1 noro 798: {
799: N r;
800:
801: asir_assert(ARG0(arg),O_N,"lrandom");
802: randomn(QTOS((Q)ARG0(arg)),&r);
803: NTOQ(r,1,*rp);
804: }
805:
1.10 ! noro 806: void Prandom(NODE arg,Q *rp)
1.1 noro 807: {
808: unsigned int r;
809:
810: #if 0
811: #if defined(_PA_RISC1_1)
812: r = mrand48()&BMASK;
813: #else
814: if ( arg )
815: srandom(QTOS((Q)ARG0(arg)));
816: r = random()&BMASK;
817: #endif
818: #endif
819: if ( arg )
820: mt_sgenrand(QTOS((Q)ARG0(arg)));
821: r = mt_genrand();
822: UTOQ(r,*rp);
823: }
824:
1.6 noro 825: #if defined(VISUAL)
1.1 noro 826: void srandom(unsigned int);
827:
828: static unsigned int R_Next;
829:
830: unsigned int random() {
831: if ( !R_Next )
832: R_Next = 1;
833: return R_Next = (R_Next * 1103515245 + 12345);
834: }
835:
1.10 ! noro 836: void srandom(unsigned int s)
1.1 noro 837: {
838: if ( s )
839: R_Next = s;
840: else if ( !R_Next )
841: R_Next = 1;
842: }
843: #endif
844:
1.10 ! noro 845: void Pprime(NODE arg,Q *rp)
1.1 noro 846: {
847: int i;
848:
849: asir_assert(ARG0(arg),O_N,"prime");
850: i = QTOS((Q)ARG0(arg));
851: if ( i < 0 || i >= 1900 )
852: *rp = 0;
853: else
854: STOQ(sprime[i],*rp);
855: }
856:
1.10 ! noro 857: void Plprime(NODE arg,Q *rp)
1.1 noro 858: {
1.9 noro 859: int i,p;
1.1 noro 860:
861: asir_assert(ARG0(arg),O_N,"lprime");
862: i = QTOS((Q)ARG0(arg));
1.9 noro 863: if ( i < 0 )
1.1 noro 864: *rp = 0;
1.9 noro 865: else
866: p = get_lprime(i);
867: STOQ(p,*rp);
1.1 noro 868: }
869:
870: extern int up_kara_mag, up_tkara_mag, up_fft_mag;
871:
1.10 ! noro 872: void Pset_upfft(NODE arg,Q *rp)
1.1 noro 873: {
874: if ( arg ) {
875: asir_assert(ARG0(arg),O_N,"set_upfft");
876: up_fft_mag = QTOS((Q)ARG0(arg));
877: }
878: STOQ(up_fft_mag,*rp);
879: }
880:
1.10 ! noro 881: void Pset_upkara(NODE arg,Q *rp)
1.1 noro 882: {
883: if ( arg ) {
884: asir_assert(ARG0(arg),O_N,"set_upkara");
885: up_kara_mag = QTOS((Q)ARG0(arg));
886: }
887: STOQ(up_kara_mag,*rp);
888: }
889:
1.10 ! noro 890: void Pset_uptkara(NODE arg,Q *rp)
1.1 noro 891: {
892: if ( arg ) {
893: asir_assert(ARG0(arg),O_N,"set_uptkara");
894: up_tkara_mag = QTOS((Q)ARG0(arg));
895: }
896: STOQ(up_tkara_mag,*rp);
897: }
898:
899: extern int up2_kara_mag;
900:
1.10 ! noro 901: void Pset_up2kara(NODE arg,Q *rp)
1.1 noro 902: {
903: if ( arg ) {
904: asir_assert(ARG0(arg),O_N,"set_up2kara");
905: up2_kara_mag = QTOS((Q)ARG0(arg));
906: }
907: STOQ(up2_kara_mag,*rp);
908: }
909:
1.10 ! noro 910: void Pigcdbin(NODE arg,Obj *rp)
1.1 noro 911: {
912: N g;
913: Q n1,n2,a;
914:
915: n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
916: asir_assert(n1,O_N,"igcd");
917: asir_assert(n2,O_N,"igcd");
918: if ( !n1 )
919: *rp = (Obj)n2;
920: else if ( !n2 )
921: *rp = (Obj)n1;
922: else {
923: gcdbinn(NM(n1),NM(n2),&g);
924: NTOQ(g,1,a); *rp = (Obj)a;
925: }
926: }
927:
1.10 ! noro 928: void Pigcdbmod(NODE arg,Obj *rp)
1.1 noro 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,"igcdbmod");
935: asir_assert(n2,O_N,"igcdbmod");
936: if ( !n1 )
937: *rp = (Obj)n2;
938: else if ( !n2 )
939: *rp = (Obj)n1;
940: else {
941: gcdbmodn(NM(n1),NM(n2),&g);
942: NTOQ(g,1,a); *rp = (Obj)a;
943: }
944: }
945:
1.10 ! noro 946: void Pigcdacc(NODE arg,Obj *rp)
1.1 noro 947: {
948: N g;
949: Q n1,n2,a;
950:
951: n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
952: asir_assert(n1,O_N,"igcdacc");
953: asir_assert(n2,O_N,"igcdacc");
954: if ( !n1 )
955: *rp = (Obj)n2;
956: else if ( !n2 )
957: *rp = (Obj)n1;
958: else {
959: gcdaccn(NM(n1),NM(n2),&g);
960: NTOQ(g,1,a); *rp = (Obj)a;
961: }
962: }
963:
1.10 ! noro 964: void PigcdEuc(NODE arg,Obj *rp)
1.1 noro 965: {
966: N g;
967: Q n1,n2,a;
968:
969: n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
970: asir_assert(n1,O_N,"igcdbmod");
971: asir_assert(n2,O_N,"igcdbmod");
972: if ( !n1 )
973: *rp = (Obj)n2;
974: else if ( !n2 )
975: *rp = (Obj)n1;
976: else {
977: gcdEuclidn(NM(n1),NM(n2),&g);
978: NTOQ(g,1,a); *rp = (Obj)a;
979: }
980: }
981:
982: extern int igcd_algorithm;
983: /* == 0 : Euclid,
984: * == 1 : binary,
985: * == 2 : bmod,
986: * >= 3 : (Weber's accelerated)/(Jebelean's generalized binary) algorithm,
987: */
988: extern int igcd_thre_inidiv;
989: /*
990: * In the non-Euclidean algorithms, if the ratio of the lengths (number
991: * of words) of two integers is >= igcd_thre_inidiv, we first perform
992: * remainder calculation.
993: * If == 0, this remainder calculation is not performed.
994: */
995: extern int igcdacc_thre;
996: /*
997: * In the accelerated algorithm, if the bit-lengths of two integers is
998: * > igcdacc_thre, "bmod" reduction is done.
999: */
1000:
1.10 ! noro 1001: void Pigcdcntl(NODE arg,Obj *rp)
1.1 noro 1002: {
1003: Obj p;
1004: Q a;
1005: int k, m;
1006:
1007: if ( arg ) {
1008: p = (Obj)ARG0(arg);
1009: if ( !p ) {
1010: igcd_algorithm = 0;
1011: *rp = p;
1012: return;
1013: } else if ( OID(p) == O_N ) {
1014: k = QTOS((Q)p);
1015: a = (Q)p;
1016: if ( k >= 0 ) igcd_algorithm = k;
1017: else if ( k == -1 ) {
1018: ret_thre:
1019: k = - igcd_thre_inidiv - igcdacc_thre*10000;
1020: STOQ(k,a);
1021: *rp = (Obj)a;
1022: return;
1023: } else {
1024: if ( (m = (-k)%10000) != 0 ) igcd_thre_inidiv = m;
1025: if ( (m = -k/10000) != 0 ) igcdacc_thre = m;
1026: goto ret_thre;
1027: }
1028: } else if ( OID(p) == O_STR ) {
1029: char *n = BDY((STRING) p);
1030:
1031: if ( !strcmp( n, "binary" ) || !strcmp( n, "Binary" )
1032: || !strcmp( n, "bin" ) || !strcmp( n, "Bin" ) )
1033: k = igcd_algorithm = 1;
1034: else if ( !strcmp( n, "bmod" ) || !strcmp( n, "Bmod" ) )
1035: igcd_algorithm = 2;
1036: else if ( !strcmp( n, "euc" ) || !strcmp( n, "Euc" )
1037: || !strcmp( n, "euclid" ) || !strcmp( n, "Euclid" ) )
1038: igcd_algorithm = 0;
1039: else if ( !strcmp( n, "acc" ) || !strcmp( n, "Acc" )
1040: || !strcmp( n, "gen" ) || !strcmp( n, "Gen" )
1041: || !strcmp( n, "genbin" ) || !strcmp( n, "GenBin" ) )
1042: igcd_algorithm = 3;
1043: else error( "igcdhow : invalid algorithm specification" );
1044: }
1045: }
1046: STOQ(igcd_algorithm,a);
1047: *rp = (Obj)a;
1048: return;
1049: }