![[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.5 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/N.c,v 1.4 2000/12/21 02:51:45 murao 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.1 noro 672:
1.5 ! noro 673: void gcdn(N n1,N n2,N *nr)
1.1 noro 674: {
675: if ( !igcd_algorithm )
676: gcdEuclidn(n1,n2,nr);
677: else {
1.4 murao 678: gcdn_HMEXT(n1,n2,nr);
1.1 noro 679: }
680: }
681:
1.4 murao 682: #include "Ngcd.c"
683:
1.5 ! noro 684: void gcdEuclidn(N n1,N n2,N *nr)
1.1 noro 685: {
686: N m1,m2,q,r;
687: unsigned int i1,i2,ir;
688:
689: if ( !n1 )
690: COPY(n2,*nr);
691: else if ( !n2 )
692: COPY(n1,*nr);
693: else {
694: if ( PL(n1) > PL(n2) ) {
695: COPY(n1,m1); COPY(n2,m2);
696: } else {
697: COPY(n1,m2); COPY(n2,m1);
698: }
699: while ( PL(m1) > 1 ) {
700: divn(m1,m2,&q,&r); FREEN(m1); FREEN(q);
701: if ( !r ) {
702: *nr = m2;
703: return;
704: } else {
705: m1 = m2; m2 = r;
706: }
707: }
1.4 murao 708: for ( i1 = BD(m1)[0], i2 = BD(m2)[0]; ir = i1 % i2; i2 = ir )
709: i1 = i2;
1.1 noro 710: if ( i2 == 1 )
711: COPY(ONEN,*nr);
712: else {
713: *nr = r = NALLOC(1); INITRC(r); PL(r) = 1; BD(r)[0] = i2;
714: }
715: }
716: }
717:
1.5 ! noro 718: int cmpn(N n1,N n2)
1.1 noro 719: {
720: int i;
721: unsigned int *m1,*m2;
722:
723: if ( !n1 )
724: if ( !n2 )
725: return 0;
726: else
727: return -1;
728: else if ( !n2 )
729: return 1;
730: else if ( PL(n1) > PL(n2) )
731: return 1;
732: else if ( PL(n1) < PL(n2) )
733: return -1;
734: else {
735: for ( i = PL(n1)-1, m1 = BD(n1)+i, m2 = BD(n2)+i;
736: i >= 0; i--, m1--, m2-- )
737: if ( *m1 > *m2 )
738: return 1;
739: else if ( *m1 < *m2 )
740: return -1;
741: return 0;
742: }
743: }
744:
1.5 ! noro 745: void bshiftn(N n,int b,N *r)
1.1 noro 746: {
747: int w,l,nl,i,j;
748: N z;
749: unsigned int msw;
750: unsigned int *p,*pz;
751:
752: if ( b == 0 ) {
753: COPY(n,*r); return;
754: }
755: if ( b > 0 ) { /* >> */
756: w = b / BSH; l = PL(n)-w;
757: if ( l <= 0 ) {
758: *r = 0; return;
759: }
760: b %= BSH; p = BD(n)+w;
761: if ( !b ) {
762: *r = z = NALLOC(l); INITRC(z); PL(z) = l;
763: bcopy(p,BD(z),l*sizeof(unsigned int));
764: return;
765: }
766: msw = p[l-1];
767: for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
768: if ( b >= i ) {
769: l--;
770: if ( !l ) {
771: *r = 0; return;
772: }
773: *r = z = NALLOC(l); INITRC(z); PL(z) = l; pz = BD(z);
774: for ( j = 0; j < l; j++, p++ )
775: *pz++ = (*(p+1)<<(BSH-b))|(*p>>b);
776: } else {
777: *r = z = NALLOC(l); INITRC(z); PL(z) = l; pz = BD(z);
778: for ( j = 1; j < l; j++, p++ )
779: *pz++ = (*(p+1)<<(BSH-b))|(*p>>b);
780: *pz = *p>>b;
781: }
782: } else { /* << */
783: b = -b;
784: w = b / BSH; b %= BSH; l = PL(n); p = BD(n);
785: if ( !b ) {
786: nl = l+w; *r = z = NALLOC(nl); INITRC(z); PL(z) = nl;
787: bzero((char *)BD(z),w*sizeof(unsigned int));
788: bcopy(p,BD(z)+w,l*sizeof(unsigned int));
789: return;
790: }
791: msw = p[l-1];
792: for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
793: if ( b + i > BSH ) {
794: nl = l+w+1;
795: *r = z = NALLOC(nl); INITRC(z); PL(z) = nl; pz = BD(z)+w;
796: bzero((char *)BD(z),w*sizeof(unsigned int));
797: *pz++ = *p++<<b;
798: for ( j = 1; j < l; j++, p++ )
799: *pz++ = (*p<<b)|(*(p-1)>>(BSH-b));
800: *pz = *(p-1)>>(BSH-b);
801: } else {
802: nl = l+w;
803: *r = z = NALLOC(nl); INITRC(z); PL(z) = nl; pz = BD(z)+w;
804: bzero((char *)BD(z),w*sizeof(unsigned int));
805: *pz++ = *p++<<b;
806: for ( j = 1; j < l; j++, p++ )
807: *pz++ = (*p<<b)|(*(p-1)>>(BSH-b));
808: }
809: }
810: }
811:
812: #if 0
1.5 ! noro 813: void _bshiftn(N n,int b,N z)
1.1 noro 814: {
815: int w,l,nl,i,j;
816: unsigned int msw;
817: unsigned int *p,*pz;
818:
819: if ( b == 0 ) {
820: copyn(n,PL(n),BD(z)); PL(z) = PL(n); return;
821: }
822: if ( b > 0 ) { /* >> */
823: w = b / BSH; l = PL(n)-w;
824: if ( l <= 0 ) {
825: PL(z) = 0; return;
826: }
827: b %= BSH; p = BD(n)+w;
828: if ( !b ) {
829: PL(z) = l;
830: bcopy(p,BD(z),l*sizeof(unsigned int));
831: return;
832: }
833: msw = p[l-1];
834: for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
835: if ( b >= i ) {
836: l--;
837: if ( !l ) {
838: PL(z) = 0; return;
839: }
840: PL(z) = l; pz = BD(z);
841: for ( j = 0; j < l; j++, p++ )
842: *pz++ = (*(p+1)<<(BSH-b))|(*p>>b);
843: } else {
844: PL(z) = l; pz = BD(z);
845: for ( j = 1; j < l; j++, p++ )
846: *pz++ = (*(p+1)<<(BSH-b))|(*p>>b);
847: *pz = *p>>b;
848: }
849: } else { /* << */
850: b = -b;
851: w = b / BSH; b %= BSH; l = PL(n); p = BD(n);
852: if ( !b ) {
853: nl = l+w; PL(z) = nl;
854: bzero((char *)BD(z),w*sizeof(unsigned int));
855: bcopy(p,BD(z)+w,l*sizeof(unsigned int));
856: return;
857: }
858: msw = p[l-1];
859: for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
860: if ( b + i > BSH ) {
861: nl = l+w+1;
862: PL(z) = nl; pz = BD(z)+w;
863: bzero((char *)BD(z),w*sizeof(unsigned int));
864: *pz++ = *p++<<b;
865: for ( j = 1; j < l; j++, p++ )
866: *pz++ = (*p<<b)|(*(p-1)>>(BSH-b));
867: *pz = *(p-1)>>(BSH-b);
868: } else {
869: nl = l+w;
870: PL(z) = nl; pz = BD(z)+w;
871: bzero((char *)BD(z),w*sizeof(unsigned int));
872: *pz++ = *p++<<b;
873: for ( j = 1; j < l; j++, p++ )
874: *pz++ = (*p<<b)|(*(p-1)>>(BSH-b));
875: }
876: }
877: }
878: #endif
879:
1.5 ! noro 880: void shiftn(N n,int w,N *r)
1.1 noro 881: {
882: int l,nl;
883: N z;
884:
885: if ( w == 0 )
886: COPY(n,*r);
887: else if ( w > 0 ) { /* >> */
888: l = PL(n)-w;
889: if ( l <= 0 )
890: *r = 0;
891: else {
892: *r = z = NALLOC(l); INITRC(z); PL(z) = l;
893: bcopy(BD(n)+w,BD(z),l*sizeof(unsigned int));
894: }
895: } else { /* << */
896: w = -w;
897: l = PL(n); nl = l+w;
898: *r = z = NALLOC(nl); INITRC(z); PL(z) = nl;
899: bzero((char *)BD(z),w*sizeof(unsigned int));
900: bcopy(BD(n),BD(z)+w,l*sizeof(unsigned int));
901: }
902: }
903:
1.5 ! noro 904: void randomn(int bits,N *r)
1.1 noro 905: {
906: int l,i;
907: unsigned int *tb;
908: N t;
909:
910: l = (bits+31)>>5; /* word length */
911: *r = t = NALLOC(l);
912: tb = BD(t);
913: for ( i = 0; i < l; i++ )
914: tb[i] = mt_genrand();
915: if ( bits&31 )
916: tb[l-1] &= (1<<(bits&31))-1;
917: for ( i = l-1; i >= 0 && !tb[i]; i-- );
918: if ( i < 0 )
919: *r = 0;
920: else
921: PL(t) = i+1;
922: }
923:
1.5 ! noro 924: void freen(N n)
1.1 noro 925: {
926: if ( n && (n != ONEN) )
927: free(n);
928: }
929:
1.5 ! noro 930: int n_bits(N n)
1.1 noro 931: {
932: unsigned int l,i,t;
933:
934: if ( !n )
935: return 0;
936: l = PL(n); t = BD(n)[l-1];
937: for ( i = 0; t; t>>=1, i++);
938: return i + (l-1)*BSH;
939: }