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