![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to N.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine|
Return to N.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.2 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.3 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2 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.4 ! murao 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/N.c,v 1.3 2000/08/22 05:04:04 noro Exp $
1.2 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include "base.h"
52:
53: void addn();
54: int subn();
55: void _bshiftn();
56: void dupn();
57:
58: #if defined(VISUAL) || defined(i386)
59: void addn(n1,n2,nr)
60: N n1,n2,*nr;
61: {
62: unsigned int *m1,*m2,*mr;
63: unsigned int c;
64: N r;
65: int i,d1,d2;
66: unsigned int tmp;
67:
68: if ( !n1 )
69: COPY(n2,*nr);
70: else if ( !n2 )
71: COPY(n1,*nr);
72: else {
73: if ( PL(n1) > PL(n2) ) {
74: d1 = PL(n1); d2 = PL(n2); m1 = BD(n1); m2 = BD(n2);
75: } else {
76: d1 = PL(n2); d2 = PL(n1); m1 = BD(n2); m2 = BD(n1);
77: }
78: *nr = r = NALLOC(d1 + 1); INITRC(r); mr = BD(r);
79:
80: #if defined(VISUAL)
81: __asm {
82: push esi
83: push edi
84: mov esi,m1
85: mov edi,m2
86: mov ebx,mr
87: mov ecx,d2
88: xor eax,eax
89: Lstart_addn:
90: mov eax,DWORD PTR [esi]
91: mov edx,DWORD PTR [edi]
92: adc eax,edx
93: mov DWORD PTR [ebx],eax
94: lea esi,DWORD PTR [esi+4]
95: lea edi,DWORD PTR [edi+4]
96: lea ebx,DWORD PTR [ebx+4]
97: dec ecx
98: jnz Lstart_addn
99: pop edi
100: pop esi
101: mov eax,0
102: adc eax,eax
103: mov c,eax
104: }
105: #else
106: asm volatile("\
107: movl %1,%%esi;\
108: movl %2,%%edi;\
109: movl %3,%%ebx;\
110: movl %4,%%ecx;\
111: testl %%eax,%%eax;\
112: Lstart_addn:;\
113: movl (%%esi),%%eax;\
114: movl (%%edi),%%edx;\
115: adcl %%edx,%%eax;\
116: movl %%eax,(%%ebx);\
117: leal 4(%%esi),%%esi;\
118: leal 4(%%edi),%%edi;\
119: leal 4(%%ebx),%%ebx;\
120: decl %%ecx;\
121: jnz Lstart_addn;\
122: movl $0,%%eax;\
123: adcl %%eax,%%eax;\
124: movl %%eax,%0"\
125: :"=m"(c)\
126: :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
127: :"eax","ebx","ecx","edx","esi","edi");
128: #endif
129: for ( i = d2, m1 += d2, mr += d2; (i < d1) && c ; i++ ) {
130: tmp = *m1++ + c;
131: c = tmp < c ? 1 : 0;
132: *mr++ = tmp;
133: }
134: for ( ; i < d1; i++ )
135: *mr++ = *m1++;
136: *mr = c;
137: PL(r) = (c?d1+1:d1);
138: }
139: }
140:
141: int subn(n1,n2,nr)
142: N n1,n2,*nr;
143: {
144: N r;
145: unsigned int *m1,*m2,*mr,br;
146: unsigned int tmp,t;
147: int d1,d2,sgn,i;
148:
149: if ( !n1 ) {
150: if ( n2 ) {
151: COPY(n2,*nr);
152: return -1;
153: } else {
154: *nr = 0;
155: return 0;
156: }
157: } else if ( !n2 ) {
158: COPY(n1,*nr);
159: return 1;
160: } else {
161: d1 = PL(n1); d2 = PL(n2);
162: m1 = BD(n1); m2 = BD(n2);
163: if ( (d1 = PL(n1)) > (d2 = PL(n2)) )
164: sgn = 1;
165: else if ( d1 < d2 ) {
166: d1 = PL(n2); d2 = PL(n1);
167: m1 = BD(n2); m2 = BD(n1);
168: sgn = -1;
169: } else {
170: for ( i = d1-1; i >= 0 && m1[i] == m2[i]; i-- );
171: if ( i < 0 ) {
172: *nr = 0;
173: return 0;
174: }
175: d1 = d2 = i+1;
176: if ( m1[i] > m2[i] )
177: sgn = 1;
178: else {
179: m1 = BD(n2); m2 = BD(n1);
180: sgn = -1;
181: }
182: }
183: *nr = r = NALLOC(d1); INITRC(r); mr = BD(r);
184:
185: #if defined(VISUAL)
186: __asm {
187: push esi
188: push edi
189: mov esi,m1
190: mov edi,m2
191: mov ebx,mr
192: mov ecx,d2
193: xor eax,eax
194: Lstart_subn:
195: mov eax,DWORD PTR [esi]
196: mov edx,DWORD PTR [edi]
197: sbb eax,edx
198: mov DWORD PTR [ebx],eax
199: lea esi,DWORD PTR [esi+4]
200: lea edi,DWORD PTR [edi+4]
201: lea ebx,DWORD PTR [ebx+4]
202: dec ecx
203: jnz Lstart_subn
204: pop edi
205: pop esi
206: mov eax,0
207: adc eax,eax
208: mov br,eax
209: }
210: #else
211: asm volatile("\
212: movl %1,%%esi;\
213: movl %2,%%edi;\
214: movl %3,%%ebx;\
215: movl %4,%%ecx;\
216: testl %%eax,%%eax;\
217: Lstart_subn:;\
218: movl (%%esi),%%eax;\
219: movl (%%edi),%%edx;\
220: sbbl %%edx,%%eax;\
221: movl %%eax,(%%ebx);\
222: leal 4(%%esi),%%esi;\
223: leal 4(%%edi),%%edi;\
224: leal 4(%%ebx),%%ebx;\
225: decl %%ecx;\
226: jnz Lstart_subn;\
227: movl $0,%%eax;\
228: adcl %%eax,%%eax;\
229: movl %%eax,%0"\
230: :"=m"(br)\
231: :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
232: :"eax","ebx","ecx","edx","esi","edi");
233: #endif
234: for ( i = d2, m1 += d2, mr += d2; (i < d1) && br; i++ ) {
235: t = *m1++;
236: tmp = t - br;
237: br = tmp > t ? 1 : 0;
238: *mr++ = tmp;
239: }
240: for ( ; i < d1; i++ )
241: *mr++ = *m1++;
242: for ( i = d1-1, mr--; i >= 0 && !*mr--; i-- );
243: PL(r) = i + 1;
244: return sgn;
245: }
246: }
247:
248: void _addn(n1,n2,nr)
249: N n1,n2,nr;
250: {
251: unsigned int *m1,*m2,*mr;
252: unsigned int c;
253: int i,d1,d2;
254: unsigned int tmp;
255:
256: if ( !n1 || !PL(n1) )
257: dupn(n2,nr);
258: else if ( !n2 || !PL(n2) )
259: dupn(n1,nr);
260: else {
261: if ( PL(n1) > PL(n2) ) {
262: d1 = PL(n1); d2 = PL(n2); m1 = BD(n1); m2 = BD(n2);
263: } else {
264: d1 = PL(n2); d2 = PL(n1); m1 = BD(n2); m2 = BD(n1);
265: }
266: mr = BD(nr);
267:
268: #if defined(VISUAL)
269: __asm {
270: push esi
271: push edi
272: mov esi,m1
273: mov edi,m2
274: mov ebx,mr
275: mov ecx,d2
276: xor eax,eax
277: Lstart__addn:
278: mov eax,DWORD PTR [esi]
279: mov edx,DWORD PTR [edi]
280: adc eax,edx
281: mov DWORD PTR [ebx],eax
282: lea esi,DWORD PTR [esi+4]
283: lea edi,DWORD PTR [edi+4]
284: lea ebx,DWORD PTR [ebx+4]
285: dec ecx
286: jnz Lstart__addn
287: pop edi
288: pop esi
289: mov eax,0
290: adc eax,eax
291: mov c,eax
292: }
293: #else
294: asm volatile("\
295: movl %1,%%esi;\
296: movl %2,%%edi;\
297: movl %3,%%ebx;\
298: movl %4,%%ecx;\
299: testl %%eax,%%eax;\
300: Lstart__addn:;\
301: movl (%%esi),%%eax;\
302: movl (%%edi),%%edx;\
303: adcl %%edx,%%eax;\
304: movl %%eax,(%%ebx);\
305: leal 4(%%esi),%%esi;\
306: leal 4(%%edi),%%edi;\
307: leal 4(%%ebx),%%ebx;\
308: decl %%ecx;\
309: jnz Lstart__addn;\
310: movl $0,%%eax;\
311: adcl %%eax,%%eax;\
312: movl %%eax,%0"\
313: :"=m"(c)\
314: :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
315: :"eax","ebx","ecx","edx","esi","edi");
316: #endif
317: for ( i = d2, m1 += d2, mr += d2; (i < d1) && c ; i++ ) {
318: tmp = *m1++ + c;
319: c = tmp < c ? 1 : 0;
320: *mr++ = tmp;
321: }
322: for ( ; i < d1; i++ )
323: *mr++ = *m1++;
324: *mr = c;
325: PL(nr) = (c?d1+1:d1);
326: }
327: }
328:
329: int _subn(n1,n2,nr)
330: N n1,n2,nr;
331: {
332: unsigned int *m1,*m2,*mr,br;
333: unsigned int tmp,t;
334: int d1,d2,sgn,i;
335:
336: if ( !n1 || !PL(n1) ) {
337: if ( n2 && PL(n2) ) {
338: dupn(n2,nr);
339: return -1;
340: } else {
341: PL(nr) = 0;
342: return 0;
343: }
344: } else if ( !n2 || !PL(n2) ) {
345: dupn(n1,nr);
346: return 1;
347: } else {
348: d1 = PL(n1); d2 = PL(n2);
349: m1 = BD(n1); m2 = BD(n2);
350: if ( (d1 = PL(n1)) > (d2 = PL(n2)) )
351: sgn = 1;
352: else if ( d1 < d2 ) {
353: d1 = PL(n2); d2 = PL(n1);
354: m1 = BD(n2); m2 = BD(n1);
355: sgn = -1;
356: } else {
357: for ( i = d1-1; i >= 0 && m1[i] == m2[i]; i-- );
358: if ( i < 0 ) {
359: PL(nr) = 0;
360: return 0;
361: }
362: d1 = d2 = i+1;
363: if ( m1[i] > m2[i] )
364: sgn = 1;
365: else {
366: m1 = BD(n2); m2 = BD(n1);
367: sgn = -1;
368: }
369: }
370: mr = BD(nr);
371:
372: #if defined(VISUAL)
373: __asm {
374: push esi
375: push edi
376: mov esi,m1
377: mov edi,m2
378: mov ebx,mr
379: mov ecx,d2
380: xor eax,eax
381: Lstart__subn:
382: mov eax,DWORD PTR [esi]
383: mov edx,DWORD PTR [edi]
384: sbb eax,edx
385: mov DWORD PTR [ebx],eax
386: lea esi,DWORD PTR [esi+4]
387: lea edi,DWORD PTR [edi+4]
388: lea ebx,DWORD PTR [ebx+4]
389: dec ecx
390: jnz Lstart__subn
391: pop edi
392: pop esi
393: mov eax,0
394: adc eax,eax
395: mov br,eax
396: }
397: #else
398: asm volatile("\
399: movl %1,%%esi;\
400: movl %2,%%edi;\
401: movl %3,%%ebx;\
402: movl %4,%%ecx;\
403: testl %%eax,%%eax;\
404: Lstart__subn:;\
405: movl (%%esi),%%eax;\
406: movl (%%edi),%%edx;\
407: sbbl %%edx,%%eax;\
408: movl %%eax,(%%ebx);\
409: leal 4(%%esi),%%esi;\
410: leal 4(%%edi),%%edi;\
411: leal 4(%%ebx),%%ebx;\
412: decl %%ecx;\
413: jnz Lstart__subn;\
414: movl $0,%%eax;\
415: adcl %%eax,%%eax;\
416: movl %%eax,%0"\
417: :"=m"(br)\
418: :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
419: :"eax","ebx","ecx","edx","esi","edi");
420: #endif
421: for ( i = d2, m1 += d2, mr += d2; (i < d1) && br; i++ ) {
422: t = *m1++;
423: tmp = t - br;
424: br = tmp > t ? 1 : 0;
425: *mr++ = tmp;
426: }
427: for ( ; i < d1; i++ )
428: *mr++ = *m1++;
429: for ( i = d1-1, mr--; i >= 0 && !*mr--; i-- );
430: PL(nr) = i + 1;
431: return sgn;
432: }
433: }
434: #else
435:
436: void addn(n1,n2,nr)
437: N n1,n2,*nr;
438: {
439: unsigned int *m1,*m2,*mr,i,c;
440: N r;
441: int d1,d2;
442: unsigned int tmp;
443:
444: if ( !n1 )
445: COPY(n2,*nr);
446: else if ( !n2 )
447: COPY(n1,*nr);
448: else {
449: if ( PL(n1) > PL(n2) ) {
450: d1 = PL(n1); d2 = PL(n2); m1 = BD(n1); m2 = BD(n2);
451: } else {
452: d1 = PL(n2); d2 = PL(n1); m1 = BD(n2); m2 = BD(n1);
453: }
454: *nr = r = NALLOC(d1 + 1); INITRC(r);
455: for ( i = 0, c = 0, mr = BD(r); i < d2; i++, m1++, m2++, mr++ ) {
456: tmp = *m1 + *m2;
457: if ( tmp < *m1 ) {
458: tmp += c;
459: c = 1;
460: } else {
461: tmp += c;
462: c = tmp < c ? 1 : 0;
463: }
464: *mr = tmp;
465: }
466: for ( ; (i < d1) && c ; i++, m1++, mr++ ) {
467: tmp = *m1 + c;
468: c = tmp < c ? 1 : 0;
469: *mr = tmp;
470: }
471: for ( ; i < d1; i++ )
472: *mr++ = *m1++;
473: *mr = c;
474: PL(r) = (c?d1+1:d1);
475: }
476: }
477:
478: int subn(n1,n2,nr)
479: N n1,n2,*nr;
480: {
481: N r;
482: unsigned int *m1,*m2,*mr,i,br;
483: L tmp;
484: int d1,d2,nz,sgn;
485:
486: if ( !n1 ) {
487: if ( n2 ) {
488: COPY(n2,*nr);
489: return -1;
490: } else {
491: *nr = 0;
492: return 0;
493: }
494: } else if ( !n2 ) {
495: COPY(n1,*nr);
496: return 1;
497: } else {
498: switch ( cmpn(n1,n2) ) {
499: case 1:
500: d1 = PL(n1); d2 = PL(n2); m1 = BD(n1); m2 = BD(n2);
501: sgn = 1; break;
502: case -1:
503: d1 = PL(n2); d2 = PL(n1); m1 = BD(n2); m2 = BD(n1);
504: sgn = -1; break;
505: case 0:
506: default:
507: *nr = 0; return ( 0 ); break;
508: }
509: *nr = r = NALLOC(d1); INITRC(r);
510: for ( i = 0, br = 0, nz = -1, mr = BD(r);
511: i < d2; i++ ) {
512: if ( (tmp = (L)*m1++ - (L)*m2++ - (L)br) && ( tmp > -LBASE ) )
513: nz = i;
514: if ( tmp < 0 ) {
515: br = 1; *mr++ = (unsigned int)(tmp + LBASE);
516: } else {
517: br = 0; *mr++ = (unsigned int)tmp;
518: }
519: }
520: for ( ; (i < d1) && br; i++ ) {
521: if ( (tmp = (L)*m1++ - (L)br) && ( tmp > -LBASE ) )
522: nz = i;
523: if ( tmp < 0 ) {
524: br = 1; *mr++ = (unsigned int)(tmp + LBASE);
525: } else {
526: br = 0; *mr++ = (unsigned int)tmp;
527: }
528: }
529: for ( ; i < d1; i++ )
530: if ( *mr++ = *m1++ )
531: nz = i;
532: PL(r) = nz + 1;
533: return sgn;
534: }
535: }
536:
537: void _addn(n1,n2,nr)
538: N n1,n2,nr;
539: {
540: unsigned int *m1,*m2,*mr,i,c;
541: int d1,d2;
542: unsigned int tmp;
543:
544: if ( !n1 || !PL(n1) )
545: dupn(n2,nr);
546: else if ( !n2 || !PL(n2) )
547: dupn(n1,nr);
548: else {
549: if ( PL(n1) > PL(n2) ) {
550: d1 = PL(n1); d2 = PL(n2); m1 = BD(n1); m2 = BD(n2);
551: } else {
552: d1 = PL(n2); d2 = PL(n1); m1 = BD(n2); m2 = BD(n1);
553: }
554: for ( i = 0, c = 0, mr = BD(nr); i < d2; i++, m1++, m2++, mr++ ) {
555: tmp = *m1 + *m2;
556: if ( tmp < *m1 ) {
557: tmp += c;
558: c = 1;
559: } else {
560: tmp += c;
561: c = tmp < c ? 1 : 0;
562: }
563: *mr = tmp;
564: }
565: for ( ; (i < d1) && c ; i++, m1++, mr++ ) {
566: tmp = *m1 + c;
567: c = tmp < c ? 1 : 0;
568: *mr = tmp;
569: }
570: for ( ; i < d1; i++ )
571: *mr++ = *m1++;
572: *mr = c;
573: PL(nr) = (c?d1+1:d1);
574: }
575: }
576:
577: int _subn(n1,n2,nr)
578: N n1,n2,nr;
579: {
580: N r;
581: unsigned int *m1,*m2,*mr,i,br;
582: L tmp;
583: int d1,d2,nz,sgn;
584:
585: if ( !n1 || !PL(n1) ) {
586: if ( n2 && PL(n2) ) {
587: dupn(n2,nr);
588: return -1;
589: } else {
590: PL(nr) = 0;
591: return 0;
592: }
593: } else if ( !n2 || !PL(n2) ) {
594: dupn(n1,nr);
595: return 1;
596: } else {
597: switch ( cmpn(n1,n2) ) {
598: case 1:
599: d1 = PL(n1); d2 = PL(n2); m1 = BD(n1); m2 = BD(n2);
600: sgn = 1; break;
601: case -1:
602: d1 = PL(n2); d2 = PL(n1); m1 = BD(n2); m2 = BD(n1);
603: sgn = -1; break;
604: case 0:
605: default:
606: PL(nr) = 0; return ( 0 ); break;
607: }
608: for ( i = 0, br = 0, nz = -1, mr = BD(nr);
609: i < d2; i++ ) {
610: if ( (tmp = (L)*m1++ - (L)*m2++ - (L)br) && ( tmp > -LBASE ) )
611: nz = i;
612: if ( tmp < 0 ) {
613: br = 1; *mr++ = (unsigned int)(tmp + LBASE);
614: } else {
615: br = 0; *mr++ = (unsigned int)tmp;
616: }
617: }
618: for ( ; (i < d1) && br; i++ ) {
619: if ( (tmp = (L)*m1++ - (L)br) && ( tmp > -LBASE ) )
620: nz = i;
621: if ( tmp < 0 ) {
622: br = 1; *mr++ = (unsigned int)(tmp + LBASE);
623: } else {
624: br = 0; *mr++ = (unsigned int)tmp;
625: }
626: }
627: for ( ; i < d1; i++ )
628: if ( *mr++ = *m1++ )
629: nz = i;
630: PL(nr) = nz + 1;
631: return sgn;
632: }
633: }
634: #endif
635:
636: /* a2 += a1; n2 >= n1 */
637:
638: void addarray_to(a1,n1,a2,n2)
639: unsigned int *a1;
640: int n1;
641: unsigned int *a2;
642: int n2;
643: {
644: int i;
645: unsigned int c,tmp;
646:
647: for ( i = 0, c = 0; i < n2; i++, a1++, a2++ ) {
648: tmp = *a1 + *a2;
649: if ( tmp < *a1 ) {
650: tmp += c;
651: c = 1;
652: } else {
653: tmp += c;
654: c = tmp < c ? 1 : 0;
655: }
656: *a2 = tmp;
657: }
658: for ( ; (i < n2) && c ; i++, a2++ ) {
659: tmp = *a2 + c;
660: c = tmp < c ? 1 : 0;
661: *a2 = tmp;
662: }
663: if ( i == n2 && c )
664: *a2 = c;
665: }
666:
667: void pwrn(n,e,nr)
668: N n,*nr;
669: int e;
670: {
671: N nw,nw1;
672:
673: if ( e == 1 ) {
674: COPY(n,*nr);
675: return;
676: }
677: pwrn(n,e / 2,&nw);
678: if ( e % 2 == 0 )
679: kmuln(nw,nw,nr);
680: else {
681: kmuln(nw,nw,&nw1); kmuln(nw1,n,nr); FREEN(nw1);
682: }
683: FREEN(nw);
684: }
685:
1.4 ! murao 686:
! 687:
1.1 noro 688: extern int igcd_algorithm;
689:
1.4 ! murao 690: void gcdEuclidn(), gcdn_HMEXT();
1.1 noro 691:
692: void gcdn(n1,n2,nr)
693: N n1,n2,*nr;
694: {
695: N m1,m2,g;
696:
697: if ( !igcd_algorithm )
698: gcdEuclidn(n1,n2,nr);
699: else {
1.4 ! murao 700: gcdn_HMEXT(n1,n2,nr);
1.1 noro 701: }
702: }
703:
1.4 ! murao 704: #include "Ngcd.c"
! 705:
1.1 noro 706: void gcdEuclidn(n1,n2,nr)
707: N n1,n2,*nr;
708: {
709: N m1,m2,q,r;
710: unsigned int i1,i2,ir;
711:
712: if ( !n1 )
713: COPY(n2,*nr);
714: else if ( !n2 )
715: COPY(n1,*nr);
716: else {
717: if ( PL(n1) > PL(n2) ) {
718: COPY(n1,m1); COPY(n2,m2);
719: } else {
720: COPY(n1,m2); COPY(n2,m1);
721: }
722: while ( PL(m1) > 1 ) {
723: divn(m1,m2,&q,&r); FREEN(m1); FREEN(q);
724: if ( !r ) {
725: *nr = m2;
726: return;
727: } else {
728: m1 = m2; m2 = r;
729: }
730: }
1.4 ! murao 731: for ( i1 = BD(m1)[0], i2 = BD(m2)[0]; ir = i1 % i2; i2 = ir )
! 732: i1 = i2;
1.1 noro 733: if ( i2 == 1 )
734: COPY(ONEN,*nr);
735: else {
736: *nr = r = NALLOC(1); INITRC(r); PL(r) = 1; BD(r)[0] = i2;
737: }
738: }
739: }
740:
741: int cmpn(n1,n2)
742: N n1,n2;
743: {
744: int i;
745: unsigned int *m1,*m2;
746:
747: if ( !n1 )
748: if ( !n2 )
749: return 0;
750: else
751: return -1;
752: else if ( !n2 )
753: return 1;
754: else if ( PL(n1) > PL(n2) )
755: return 1;
756: else if ( PL(n1) < PL(n2) )
757: return -1;
758: else {
759: for ( i = PL(n1)-1, m1 = BD(n1)+i, m2 = BD(n2)+i;
760: i >= 0; i--, m1--, m2-- )
761: if ( *m1 > *m2 )
762: return 1;
763: else if ( *m1 < *m2 )
764: return -1;
765: return 0;
766: }
767: }
768:
769: void bshiftn(n,b,r)
770: N n;
771: int b;
772: N *r;
773: {
774: int w,l,nl,i,j;
775: N z;
776: unsigned int msw;
777: unsigned int *p,*pz;
778:
779: if ( b == 0 ) {
780: COPY(n,*r); return;
781: }
782: if ( b > 0 ) { /* >> */
783: w = b / BSH; l = PL(n)-w;
784: if ( l <= 0 ) {
785: *r = 0; return;
786: }
787: b %= BSH; p = BD(n)+w;
788: if ( !b ) {
789: *r = z = NALLOC(l); INITRC(z); PL(z) = l;
790: bcopy(p,BD(z),l*sizeof(unsigned int));
791: return;
792: }
793: msw = p[l-1];
794: for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
795: if ( b >= i ) {
796: l--;
797: if ( !l ) {
798: *r = 0; return;
799: }
800: *r = z = NALLOC(l); INITRC(z); PL(z) = l; pz = BD(z);
801: for ( j = 0; j < l; j++, p++ )
802: *pz++ = (*(p+1)<<(BSH-b))|(*p>>b);
803: } else {
804: *r = z = NALLOC(l); INITRC(z); PL(z) = l; pz = BD(z);
805: for ( j = 1; j < l; j++, p++ )
806: *pz++ = (*(p+1)<<(BSH-b))|(*p>>b);
807: *pz = *p>>b;
808: }
809: } else { /* << */
810: b = -b;
811: w = b / BSH; b %= BSH; l = PL(n); p = BD(n);
812: if ( !b ) {
813: nl = l+w; *r = z = NALLOC(nl); INITRC(z); PL(z) = nl;
814: bzero((char *)BD(z),w*sizeof(unsigned int));
815: bcopy(p,BD(z)+w,l*sizeof(unsigned int));
816: return;
817: }
818: msw = p[l-1];
819: for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
820: if ( b + i > BSH ) {
821: nl = l+w+1;
822: *r = z = NALLOC(nl); INITRC(z); PL(z) = nl; pz = BD(z)+w;
823: bzero((char *)BD(z),w*sizeof(unsigned int));
824: *pz++ = *p++<<b;
825: for ( j = 1; j < l; j++, p++ )
826: *pz++ = (*p<<b)|(*(p-1)>>(BSH-b));
827: *pz = *(p-1)>>(BSH-b);
828: } else {
829: nl = l+w;
830: *r = z = NALLOC(nl); INITRC(z); PL(z) = nl; pz = BD(z)+w;
831: bzero((char *)BD(z),w*sizeof(unsigned int));
832: *pz++ = *p++<<b;
833: for ( j = 1; j < l; j++, p++ )
834: *pz++ = (*p<<b)|(*(p-1)>>(BSH-b));
835: }
836: }
837: }
838:
839: #if 0
840: void _bshiftn(n,b,z)
841: N n;
842: int b;
843: N z;
844: {
845: int w,l,nl,i,j;
846: unsigned int msw;
847: unsigned int *p,*pz;
848:
849: if ( b == 0 ) {
850: copyn(n,PL(n),BD(z)); PL(z) = PL(n); return;
851: }
852: if ( b > 0 ) { /* >> */
853: w = b / BSH; l = PL(n)-w;
854: if ( l <= 0 ) {
855: PL(z) = 0; return;
856: }
857: b %= BSH; p = BD(n)+w;
858: if ( !b ) {
859: PL(z) = l;
860: bcopy(p,BD(z),l*sizeof(unsigned int));
861: return;
862: }
863: msw = p[l-1];
864: for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
865: if ( b >= i ) {
866: l--;
867: if ( !l ) {
868: PL(z) = 0; return;
869: }
870: PL(z) = l; pz = BD(z);
871: for ( j = 0; j < l; j++, p++ )
872: *pz++ = (*(p+1)<<(BSH-b))|(*p>>b);
873: } else {
874: PL(z) = l; pz = BD(z);
875: for ( j = 1; j < l; j++, p++ )
876: *pz++ = (*(p+1)<<(BSH-b))|(*p>>b);
877: *pz = *p>>b;
878: }
879: } else { /* << */
880: b = -b;
881: w = b / BSH; b %= BSH; l = PL(n); p = BD(n);
882: if ( !b ) {
883: nl = l+w; PL(z) = nl;
884: bzero((char *)BD(z),w*sizeof(unsigned int));
885: bcopy(p,BD(z)+w,l*sizeof(unsigned int));
886: return;
887: }
888: msw = p[l-1];
889: for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
890: if ( b + i > BSH ) {
891: nl = l+w+1;
892: PL(z) = nl; pz = BD(z)+w;
893: bzero((char *)BD(z),w*sizeof(unsigned int));
894: *pz++ = *p++<<b;
895: for ( j = 1; j < l; j++, p++ )
896: *pz++ = (*p<<b)|(*(p-1)>>(BSH-b));
897: *pz = *(p-1)>>(BSH-b);
898: } else {
899: nl = l+w;
900: PL(z) = nl; pz = BD(z)+w;
901: bzero((char *)BD(z),w*sizeof(unsigned int));
902: *pz++ = *p++<<b;
903: for ( j = 1; j < l; j++, p++ )
904: *pz++ = (*p<<b)|(*(p-1)>>(BSH-b));
905: }
906: }
907: }
908: #endif
909:
910: void shiftn(n,w,r)
911: N n;
912: int w;
913: N *r;
914: {
915: int l,nl;
916: N z;
917:
918: if ( w == 0 )
919: COPY(n,*r);
920: else if ( w > 0 ) { /* >> */
921: l = PL(n)-w;
922: if ( l <= 0 )
923: *r = 0;
924: else {
925: *r = z = NALLOC(l); INITRC(z); PL(z) = l;
926: bcopy(BD(n)+w,BD(z),l*sizeof(unsigned int));
927: }
928: } else { /* << */
929: w = -w;
930: l = PL(n); nl = l+w;
931: *r = z = NALLOC(nl); INITRC(z); PL(z) = nl;
932: bzero((char *)BD(z),w*sizeof(unsigned int));
933: bcopy(BD(n),BD(z)+w,l*sizeof(unsigned int));
934: }
935: }
936:
937: void randomn(bits,r)
938: int bits;
939: N *r;
940: {
941: int l,i;
942: unsigned int *tb;
943: N t;
944:
945: l = (bits+31)>>5; /* word length */
946: *r = t = NALLOC(l);
947: tb = BD(t);
948: for ( i = 0; i < l; i++ )
949: tb[i] = mt_genrand();
950: if ( bits&31 )
951: tb[l-1] &= (1<<(bits&31))-1;
952: for ( i = l-1; i >= 0 && !tb[i]; i-- );
953: if ( i < 0 )
954: *r = 0;
955: else
956: PL(t) = i+1;
957: }
958:
959: void freen(n)
960: N n;
961: {
962: if ( n && (n != ONEN) )
963: free(n);
964: }
965:
966: int n_bits(n)
967: N n;
968: {
969: unsigned int l,i,t;
970:
971: if ( !n )
972: return 0;
973: l = PL(n); t = BD(n)[l-1];
974: for ( i = 0; t; t>>=1, i++);
975: return i + (l-1)*BSH;
976: }