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