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