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