![[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 */