Annotation of OpenXM_contrib2/asir2000/engine/N.c, Revision 1.12
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.12 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/N.c,v 1.11 2015/08/29 04:15:04 fujimoto Exp $
1.2 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include "base.h"
52:
1.11 fujimoto 53: #if (defined(_M_IX86) || defined(i386)) && !defined(__MINGW32__)
1.5 noro 54: void addn(N n1,N n2,N *nr)
1.1 noro 55: {
1.12 ! noro 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);
1.1 noro 73:
1.10 ohara 74: #if defined(_M_IX86)
1.12 ! 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: }
1.1 noro 99: #else
1.12 ! noro 100: asm volatile("\
! 101: pushl %%ebx;\
! 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;\
! 119: movl %%eax,%0;\
! 120: popl %%ebx"\
! 121: :"=m"(c)\
! 122: :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
! 123: :"eax","ecx","edx","esi","edi");
1.1 noro 124: #endif
1.12 ! noro 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: }
1.1 noro 135: }
136:
1.5 noro 137: int subn(N n1,N n2,N *nr)
1.1 noro 138: {
1.12 ! noro 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);
1.1 noro 179:
1.10 ohara 180: #if defined(_M_IX86)
1.12 ! 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: }
1.1 noro 205: #else
1.12 ! noro 206: asm volatile("\
! 207: pushl %%ebx;\
! 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;\
! 225: movl %%eax,%0;\
! 226: popl %%ebx"\
! 227: :"=m"(br)\
! 228: :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
! 229: :"eax","ecx","edx","esi","edi");
1.1 noro 230: #endif
1.12 ! noro 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: }
1.1 noro 243: }
244:
1.5 noro 245: void _addn(N n1,N n2,N nr)
1.1 noro 246: {
1.12 ! noro 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);
1.1 noro 263:
1.10 ohara 264: #if defined(_M_IX86)
1.12 ! 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: }
1.1 noro 289: #else
1.12 ! noro 290: asm volatile("\
! 291: pushl %%ebx;\
! 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;\
! 309: movl %%eax,%0;\
! 310: popl %%ebx"\
! 311: :"=m"(c)\
! 312: :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
! 313: :"eax","ecx","edx","esi","edi");
1.1 noro 314: #endif
1.12 ! noro 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: }
1.1 noro 325: }
326:
1.5 noro 327: int _subn(N n1,N n2,N nr)
1.1 noro 328: {
1.12 ! noro 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);
1.1 noro 368:
1.10 ohara 369: #if defined(_M_IX86)
1.12 ! 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: }
1.1 noro 394: #else
1.12 ! noro 395: asm volatile("\
! 396: pushl %%ebx;\
! 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;\
! 414: movl %%eax,%0;\
! 415: popl %%ebx"\
! 416: :"=m"(br)\
! 417: :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
! 418: :"eax","ecx","edx","esi","edi");
1.1 noro 419: #endif
1.12 ! noro 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: }
1.1 noro 432: }
433: #else
434:
1.5 noro 435: void addn(N n1,N n2,N *nr)
1.1 noro 436: {
1.12 ! noro 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: }
1.1 noro 474: }
475:
1.5 noro 476: int subn(N n1,N n2,N *nr)
1.1 noro 477: {
1.12 ! noro 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: }
1.1 noro 532: }
533:
1.5 noro 534: void _addn(N n1,N n2,N nr)
1.1 noro 535: {
1.12 ! noro 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: }
1.1 noro 571: }
572:
1.5 noro 573: int _subn(N n1,N n2,N nr)
1.1 noro 574: {
1.12 ! noro 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: }
1.1 noro 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: {
1.12 ! noro 635: int i;
! 636: unsigned int c,tmp;
1.1 noro 637:
1.12 ! noro 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;
1.1 noro 656: }
657:
1.5 noro 658: void pwrn(N n,int e,N *nr)
1.1 noro 659: {
1.12 ! noro 660: N nw,nw1;
1.1 noro 661:
1.12 ! noro 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);
1.1 noro 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: {
1.12 ! noro 683: N g,t;
1.7 noro 684:
1.12 ! noro 685: gcdn(n1,n2,&g);
! 686: divsn(n1,g,&t);
! 687: muln(t,n2,nr);
1.7 noro 688: }
1.1 noro 689:
1.5 noro 690: void gcdn(N n1,N n2,N *nr)
1.1 noro 691: {
1.12 ! noro 692: if ( !igcd_algorithm )
! 693: gcdEuclidn(n1,n2,nr);
! 694: else {
! 695: gcdn_HMEXT(n1,n2,nr);
! 696: }
1.1 noro 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: {
1.12 ! noro 703: N m1,m2,q,r;
! 704: unsigned int i1,i2,ir;
1.1 noro 705:
1.12 ! noro 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: }
! 725: for ( i1 = BD(m1)[0], i2 = BD(m2)[0]; ir = i1 % i2; i2 = ir )
! 726: i1 = i2;
! 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: }
1.1 noro 733: }
734:
1.5 noro 735: int cmpn(N n1,N n2)
1.1 noro 736: {
1.12 ! noro 737: int i;
! 738: unsigned int *m1,*m2;
1.1 noro 739:
1.12 ! noro 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: }
1.1 noro 760: }
761:
1.5 noro 762: void bshiftn(N n,int b,N *r)
1.1 noro 763: {
1.12 ! noro 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: }
1.1 noro 827: }
828:
829: #if 0
1.5 noro 830: void _bshiftn(N n,int b,N z)
1.1 noro 831: {
1.12 ! noro 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: }
1.1 noro 894: }
895: #endif
896:
1.5 noro 897: void shiftn(N n,int w,N *r)
1.1 noro 898: {
1.12 ! noro 899: int l,nl;
! 900: N z;
1.1 noro 901:
1.12 ! noro 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: }
1.1 noro 919: }
920:
1.5 noro 921: void randomn(int bits,N *r)
1.1 noro 922: {
1.12 ! noro 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;
1.1 noro 939: }
940:
1.5 noro 941: void freen(N n)
1.1 noro 942: {
1.12 ! noro 943: if ( n && (n != ONEN) )
! 944: free(n);
1.1 noro 945: }
946:
1.6 noro 947: /* accepts Z */
1.5 noro 948: int n_bits(N n)
1.1 noro 949: {
1.12 ! noro 950: unsigned int i,t;
! 951: int l;
1.1 noro 952:
1.12 ! noro 953: if ( !n )
! 954: return 0;
! 955: l = PL(n);
! 956: if ( l < 0 ) l = -l;
! 957: t = BD(n)[l-1];
! 958: for ( i = 0; t; t>>=1, i++);
! 959: return i + (l-1)*BSH;
1.1 noro 960: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[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