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