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