![[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
! 26: * e-mail at risa-admin@flab.fujitsu.co.jp of the detailed specification
! 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: *
! 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/N.c,v 1.1.1.1 1999/12/03 07:39:08 noro Exp $
! 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:
686: extern int igcd_algorithm;
687:
688: void gcdEuclidn();
689:
690: void gcdn(n1,n2,nr)
691: N n1,n2,*nr;
692: {
693: N m1,m2,g;
694:
695: if ( !igcd_algorithm )
696: gcdEuclidn(n1,n2,nr);
697: else {
698: if ( !n1 )
699: *nr = n2;
700: else if ( !n2 )
701: *nr = n1;
702: else {
703: n32ton27(n1,&m1);
704: n32ton27(n2,&m2);
705: gcdBinary_27n(m1,m2,&g);
706: n27ton32(g,nr);
707: }
708: }
709: }
710:
711: void gcdEuclidn(n1,n2,nr)
712: N n1,n2,*nr;
713: {
714: N m1,m2,q,r;
715: unsigned int i1,i2,ir;
716:
717: if ( !n1 )
718: COPY(n2,*nr);
719: else if ( !n2 )
720: COPY(n1,*nr);
721: else {
722: if ( PL(n1) > PL(n2) ) {
723: COPY(n1,m1); COPY(n2,m2);
724: } else {
725: COPY(n1,m2); COPY(n2,m1);
726: }
727: while ( PL(m1) > 1 ) {
728: divn(m1,m2,&q,&r); FREEN(m1); FREEN(q);
729: if ( !r ) {
730: *nr = m2;
731: return;
732: } else {
733: m1 = m2; m2 = r;
734: }
735: }
736: for ( i1 = BD(m1)[0], i2 = BD(m2)[0]; ir = i1 % i2; i1 = i2, i2 = ir );
737: if ( i2 == 1 )
738: COPY(ONEN,*nr);
739: else {
740: *nr = r = NALLOC(1); INITRC(r); PL(r) = 1; BD(r)[0] = i2;
741: }
742: }
743: }
744:
745: int cmpn(n1,n2)
746: N n1,n2;
747: {
748: int i;
749: unsigned int *m1,*m2;
750:
751: if ( !n1 )
752: if ( !n2 )
753: return 0;
754: else
755: return -1;
756: else if ( !n2 )
757: return 1;
758: else if ( PL(n1) > PL(n2) )
759: return 1;
760: else if ( PL(n1) < PL(n2) )
761: return -1;
762: else {
763: for ( i = PL(n1)-1, m1 = BD(n1)+i, m2 = BD(n2)+i;
764: i >= 0; i--, m1--, m2-- )
765: if ( *m1 > *m2 )
766: return 1;
767: else if ( *m1 < *m2 )
768: return -1;
769: return 0;
770: }
771: }
772:
773: void bshiftn(n,b,r)
774: N n;
775: int b;
776: N *r;
777: {
778: int w,l,nl,i,j;
779: N z;
780: unsigned int msw;
781: unsigned int *p,*pz;
782:
783: if ( b == 0 ) {
784: COPY(n,*r); return;
785: }
786: if ( b > 0 ) { /* >> */
787: w = b / BSH; l = PL(n)-w;
788: if ( l <= 0 ) {
789: *r = 0; return;
790: }
791: b %= BSH; p = BD(n)+w;
792: if ( !b ) {
793: *r = z = NALLOC(l); INITRC(z); PL(z) = l;
794: bcopy(p,BD(z),l*sizeof(unsigned int));
795: return;
796: }
797: msw = p[l-1];
798: for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
799: if ( b >= i ) {
800: l--;
801: if ( !l ) {
802: *r = 0; return;
803: }
804: *r = z = NALLOC(l); INITRC(z); PL(z) = l; pz = BD(z);
805: for ( j = 0; j < l; j++, p++ )
806: *pz++ = (*(p+1)<<(BSH-b))|(*p>>b);
807: } else {
808: *r = z = NALLOC(l); INITRC(z); PL(z) = l; pz = BD(z);
809: for ( j = 1; j < l; j++, p++ )
810: *pz++ = (*(p+1)<<(BSH-b))|(*p>>b);
811: *pz = *p>>b;
812: }
813: } else { /* << */
814: b = -b;
815: w = b / BSH; b %= BSH; l = PL(n); p = BD(n);
816: if ( !b ) {
817: nl = l+w; *r = z = NALLOC(nl); INITRC(z); PL(z) = nl;
818: bzero((char *)BD(z),w*sizeof(unsigned int));
819: bcopy(p,BD(z)+w,l*sizeof(unsigned int));
820: return;
821: }
822: msw = p[l-1];
823: for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
824: if ( b + i > BSH ) {
825: nl = l+w+1;
826: *r = z = NALLOC(nl); INITRC(z); PL(z) = nl; pz = BD(z)+w;
827: bzero((char *)BD(z),w*sizeof(unsigned int));
828: *pz++ = *p++<<b;
829: for ( j = 1; j < l; j++, p++ )
830: *pz++ = (*p<<b)|(*(p-1)>>(BSH-b));
831: *pz = *(p-1)>>(BSH-b);
832: } else {
833: nl = l+w;
834: *r = z = NALLOC(nl); INITRC(z); PL(z) = nl; pz = BD(z)+w;
835: bzero((char *)BD(z),w*sizeof(unsigned int));
836: *pz++ = *p++<<b;
837: for ( j = 1; j < l; j++, p++ )
838: *pz++ = (*p<<b)|(*(p-1)>>(BSH-b));
839: }
840: }
841: }
842:
843: #if 0
844: void _bshiftn(n,b,z)
845: N n;
846: int b;
847: N z;
848: {
849: int w,l,nl,i,j;
850: unsigned int msw;
851: unsigned int *p,*pz;
852:
853: if ( b == 0 ) {
854: copyn(n,PL(n),BD(z)); PL(z) = PL(n); return;
855: }
856: if ( b > 0 ) { /* >> */
857: w = b / BSH; l = PL(n)-w;
858: if ( l <= 0 ) {
859: PL(z) = 0; return;
860: }
861: b %= BSH; p = BD(n)+w;
862: if ( !b ) {
863: PL(z) = l;
864: bcopy(p,BD(z),l*sizeof(unsigned int));
865: return;
866: }
867: msw = p[l-1];
868: for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
869: if ( b >= i ) {
870: l--;
871: if ( !l ) {
872: PL(z) = 0; return;
873: }
874: PL(z) = l; pz = BD(z);
875: for ( j = 0; j < l; j++, p++ )
876: *pz++ = (*(p+1)<<(BSH-b))|(*p>>b);
877: } else {
878: PL(z) = l; pz = BD(z);
879: for ( j = 1; j < l; j++, p++ )
880: *pz++ = (*(p+1)<<(BSH-b))|(*p>>b);
881: *pz = *p>>b;
882: }
883: } else { /* << */
884: b = -b;
885: w = b / BSH; b %= BSH; l = PL(n); p = BD(n);
886: if ( !b ) {
887: nl = l+w; PL(z) = nl;
888: bzero((char *)BD(z),w*sizeof(unsigned int));
889: bcopy(p,BD(z)+w,l*sizeof(unsigned int));
890: return;
891: }
892: msw = p[l-1];
893: for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
894: if ( b + i > BSH ) {
895: nl = l+w+1;
896: PL(z) = nl; pz = BD(z)+w;
897: bzero((char *)BD(z),w*sizeof(unsigned int));
898: *pz++ = *p++<<b;
899: for ( j = 1; j < l; j++, p++ )
900: *pz++ = (*p<<b)|(*(p-1)>>(BSH-b));
901: *pz = *(p-1)>>(BSH-b);
902: } else {
903: nl = l+w;
904: PL(z) = nl; pz = BD(z)+w;
905: bzero((char *)BD(z),w*sizeof(unsigned int));
906: *pz++ = *p++<<b;
907: for ( j = 1; j < l; j++, p++ )
908: *pz++ = (*p<<b)|(*(p-1)>>(BSH-b));
909: }
910: }
911: }
912: #endif
913:
914: void shiftn(n,w,r)
915: N n;
916: int w;
917: N *r;
918: {
919: int l,nl;
920: N z;
921:
922: if ( w == 0 )
923: COPY(n,*r);
924: else if ( w > 0 ) { /* >> */
925: l = PL(n)-w;
926: if ( l <= 0 )
927: *r = 0;
928: else {
929: *r = z = NALLOC(l); INITRC(z); PL(z) = l;
930: bcopy(BD(n)+w,BD(z),l*sizeof(unsigned int));
931: }
932: } else { /* << */
933: w = -w;
934: l = PL(n); nl = l+w;
935: *r = z = NALLOC(nl); INITRC(z); PL(z) = nl;
936: bzero((char *)BD(z),w*sizeof(unsigned int));
937: bcopy(BD(n),BD(z)+w,l*sizeof(unsigned int));
938: }
939: }
940:
941: void randomn(bits,r)
942: int bits;
943: N *r;
944: {
945: int l,i;
946: unsigned int *tb;
947: N t;
948:
949: l = (bits+31)>>5; /* word length */
950: *r = t = NALLOC(l);
951: tb = BD(t);
952: for ( i = 0; i < l; i++ )
953: tb[i] = mt_genrand();
954: if ( bits&31 )
955: tb[l-1] &= (1<<(bits&31))-1;
956: for ( i = l-1; i >= 0 && !tb[i]; i-- );
957: if ( i < 0 )
958: *r = 0;
959: else
960: PL(t) = i+1;
961: }
962:
963: void freen(n)
964: N n;
965: {
966: if ( n && (n != ONEN) )
967: free(n);
968: }
969:
970: int n_bits(n)
971: N n;
972: {
973: unsigned int l,i,t;
974:
975: if ( !n )
976: return 0;
977: l = PL(n); t = BD(n)[l-1];
978: for ( i = 0; t; t>>=1, i++);
979: return i + (l-1)*BSH;
980: }