Annotation of OpenXM_contrib2/asir2000/engine/Ngcd.c, Revision 1.1
1.1 ! murao 1: /*
! 2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
! 3: * All rights reserved.
! 4: *
! 5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
! 6: * non-exclusive and royalty-free license to use, copy, modify and
! 7: * redistribute, solely for non-commercial and non-profit purposes, the
! 8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
! 9: * conditions of this Agreement. For the avoidance of doubt, you acquire
! 10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
! 11: * third party developer retains all rights, including but not limited to
! 12: * copyrights, in and to the SOFTWARE.
! 13: *
! 14: * (1) FLL does not grant you a license in any way for commercial
! 15: * purposes. You may use the SOFTWARE only for non-commercial and
! 16: * non-profit purposes only, such as academic, research and internal
! 17: * business use.
! 18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
! 19: * international copyright treaties. If you make copies of the SOFTWARE,
! 20: * with or without modification, as permitted hereunder, you shall affix
! 21: * to all such copies of the SOFTWARE the above copyright notice.
! 22: * (3) An explicit reference to this SOFTWARE and its copyright owner
! 23: * shall be made on your publication or presentation in any form of the
! 24: * results obtained by use of the SOFTWARE.
! 25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
! 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
! 27: * for such modification or the source code of the modified part of the
! 28: * SOFTWARE.
! 29: *
! 30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
! 31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
! 32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
! 33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
! 34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
! 35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
! 36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
! 37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
! 38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
! 39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
! 40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
! 41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
! 42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
! 43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
! 44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
! 45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
! 46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
! 47: *
! 48: */
! 49: /*
! 50: #include "ca.h"
! 51: #include "base.h"
! 52: */
! 53:
! 54: #undef CALL
! 55:
! 56: /**** Machine specific ****/
! 57: #define BIT_WIDTH_OF_INT 32
! 58:
! 59:
! 60: #if BSH == BIT_WIDTH_OF_INT
! 61: #define BMask ((unsigned)(-1))
! 62: #define BMASK BMask
! 63: #endif
! 64:
! 65:
! 66: int igcd_algorithm = 0;
! 67: /* == 0 : Euclid,
! 68: * == 1 : binary,
! 69: * == 2 : bmod,
! 70: * >= 3 : (Weber's accelerated)/(Jebelean's generalized binary) algorithm,
! 71: */
! 72: int igcd_thre_inidiv = 50;
! 73: /*
! 74: * In the non-Euclidean algorithms, if the ratio of the lengths (number
! 75: * of words) of two integers is >= igcd_thre_inidiv, we first perform
! 76: * remainder calculation.
! 77: * If == 0, this remainder calculation is not performed.
! 78: */
! 79: int igcdacc_thre = 10;
! 80: /*
! 81: * In the accelerated algorithm, if the bit-lengths of two integers is
! 82: * > igcdacc_thre, "bmod" reduction is done.
! 83: */
! 84:
! 85: #include "inline.h"
! 86:
! 87: #define TRAILINGZEROS(t,cntr) for(cntr=0;(t&1)==0;t>>=1)cntr++;
! 88:
! 89: #define W_NALLOC(d) ((N)ALLOCA(TRUESIZE(oN,(d)-1,int)))
! 90:
! 91: #define ShouldCompRemInit(n1,n2) (igcd_thre_inidiv != 0 && PL(n1) >= igcd_thre_inidiv*PL(n2))
! 92:
! 93: #define IniDiv(n1,n2) \
! 94: if ( ShouldCompRemInit(n1,n2) ) {\
! 95: N q, r; int w, b; \
! 96: divn(n1,n2,&q,&r); \
! 97: if ( !r ) return(n2); \
! 98: b = trailingzerosn( r, &w ); \
! 99: q = n1; n1 = n2; n2 = q; \
! 100: rshiftn( r, w, b, n2 ); \
! 101: }
! 102:
! 103: /*
! 104: * Binary GCD algorithm by J.Stein
! 105: * [J. Comp. Phys. Vol. 1 (1967), pp. 397-405)]:
! 106: * The right-shift binary algorithm is used.
! 107: */
! 108:
! 109:
! 110: /*
! 111: * subsidiary routines for gcdbinn below.
! 112: */
! 113: static int /* number of bits */ trailingzeros_nbd( /* BD of N */ nbd, pnw )
! 114: int *nbd, *pnw /* number of zero words */;
! 115: {
! 116: #if BSH == BIT_WIDTH_OF_INT
! 117: unsigned
! 118: #endif
! 119: int nw, nb, w;
! 120:
! 121: for ( nw = 0; (w = *nbd) == 0; nbd++ ) nw++;
! 122: TRAILINGZEROS(w,nb);
! 123: *pnw = nw;
! 124: return nb;
! 125: }
! 126:
! 127: #define trailingzerosn(n,pnw) trailingzeros_nbd(BD(n),pnw)
! 128:
! 129: static int /* PL of N */ rshift_nbd( /* BD of N */ nbd, /* PL of N */ nl,
! 130: /* # words */ shw, /* # bits */ shb, /* BD of N */ p )
! 131: #if BSH == BIT_WIDTH_OF_INT
! 132: unsigned
! 133: #endif
! 134: int *nbd, nl, shw, shb, *p;
! 135: {
! 136: unsigned int i, v, w, lshb; /* <---- */
! 137:
! 138: nbd += shw, i = (nl -= shw);
! 139: if ( shb == 0 ) {
! 140: for ( ; nl > 0; nl-- ) *p++ = *nbd++;
! 141: return i;
! 142: } else if ( nl < 2 ) {
! 143: *p = (*nbd) >> shb;
! 144: return 1;
! 145: }
! 146: for ( lshb = BSH - shb, v = *nbd++; --nl > 0; v = w ) {
! 147: w = *nbd++;
! 148: #if BSH == BIT_WIDTH_OF_INT
! 149: *p++ = (v >> shb) | (w << lshb); /********/
! 150: #else
! 151: *p++ = (v >> shb) | ((w << lshb)&BMASK);
! 152: #endif
! 153: }
! 154: if ( (v >>= shb) == 0 ) return( i-1 );
! 155: *p = v;
! 156: return i;
! 157: }
! 158:
! 159: #define rshiftn(ns,shw,shb,nd) (PL(nd)=rshift_nbd(BD(ns),PL(ns),shw,shb,BD(nd)))
! 160: /* nd <= ns << (shb + shw*BSH), returns PL of the result */
! 161:
! 162: #ifdef FULLSET
! 163: static N N_of_i_lshifted_by_wb( i, gw, gb )
! 164: int i, gw, gb;
! 165: /*
! 166: * returns pointer to a new struct (N)(((int)i) >> (gb + gw*BSH))
! 167: */
! 168: {
! 169: unsigned int j, l, *p; /* <---- */
! 170: N n;
! 171:
! 172: j = i >> (BSH - gb);
! 173: #if BSH == BIT_WIDTH_OF_INT
! 174: i = (i << gb); /********/
! 175: #else
! 176: i = (i << gb)&BMASK;
! 177: #endif
! 178: l = j != 0 ? gw + 2 : gw + 1;
! 179: n = NALLOC(l);
! 180: PL(n) = l;
! 181: for ( p = BD(n); gw-- > 0; ) *p++ = 0;
! 182: *p++ = i;
! 183: if ( j != 0 ) *p = j;
! 184: return n;
! 185: }
! 186: #endif /* FULLSET */
! 187:
! 188: /*
! 189: * routines to make a new struct
! 190: * (N)(((BD of N)(b[0],...,b[lb-1])) << (gb + gw*BSH))
! 191: */
! 192: static N N_of_nbd_lshifted_by_wb( /* BD of N */ b, /* PL of N */ lb, gw, gb )
! 193: int *b, lb, gw, gb;
! 194: /*
! 195: * returns pointer to a new struct
! 196: * (N)(((BD of N)(b[0],...,b[lb-1])) << (gb + gw*BSH))
! 197: */
! 198: {
! 199: unsigned int rsh, s, t, *p, l; /* <---- */
! 200: N n;
! 201:
! 202: l = lb + gw;
! 203: if ( gb == 0 ) {
! 204: n = NALLOC(l);
! 205: PL(n) = l;
! 206: for ( p = BD(n); gw-- > 0; ) *p++ = 0;
! 207: while ( lb-- > 0 ) *p++ = *b++;
! 208: return n;
! 209: }
! 210: rsh = BSH - gb; s = b[lb-1];
! 211: if ( (t = s >> rsh) != 0 ) {
! 212: n = NALLOC(l+1);
! 213: PL(n) = l+1;
! 214: (p = BD(n))[l] = t;
! 215: } else {
! 216: n = NALLOC(l);
! 217: PL(n) = l;
! 218: p = BD(n);
! 219: }
! 220: while ( gw-- > 0 ) *p++ = 0;
! 221: #if BSH == BIT_WIDTH_OF_INT
! 222: *p++ = (t = *b++) << gb; /********/
! 223: for ( ; --lb > 0; t = s )
! 224: *p++ = (t >> rsh) | ((s = *b++) << gb); /********/
! 225: #else
! 226: *p++ = ((t = *b++) << gb)&BMASK;
! 227: for ( ; --lb > 0; t = s )
! 228: *p++ = (t >> rsh) | (((s = *b++) << gb)&BMASK);
! 229: #endif
! 230: return n;
! 231: }
! 232:
! 233: #define N_of_n_lshifted_by_wb(a,gw,gb) N_of_nbd_lshifted_by_wb(BD(a),PL(a),gw,gb)
! 234:
! 235: #define SWAP(a,b,Type) { Type temp=a;a=b;b=temp;}
! 236: #define SIGNED_VAL(a,s) ((s)>0?(a):-(a))
! 237:
! 238:
! 239: #ifdef CALL
! 240: static int bw_int32( n )
! 241: unsigned int n; /* <---- */
! 242: {
! 243: int w;
! 244:
! 245: w = 0;
! 246: #if BSH > 32
! 247: if ( n > 0xffffffff ) w += 32, n >>= 32;
! 248: #endif
! 249: if ( n >= 0x10000 ) w += 16, n >>= 16;
! 250: if ( n >= 0x100 ) w += 8, n >>= 8;
! 251: if ( n >= 0x10 ) w += 4, n >>= 4;
! 252: if ( n >= 0x4 ) w += 2, n >>= 2;
! 253: if ( n >= 0x2 ) w += 1, n >>= 1;
! 254: if ( n != 0 ) ++w;
! 255: return w;
! 256: }
! 257: #define BitWidth(n,bw) bw = bw_int32( n )
! 258: #else
! 259:
! 260: #if BSH > 32
! 261: #define BitWidth(n,bw) {\
! 262: unsigned int k = (n); \
! 263: bw = 0; \
! 264: if ( k > 0xffffffff ) bw += 32, k >>= 32; \
! 265: if ( k >= 0x10000 ) bw += 16, k >>= 16; \
! 266: if ( k >= 0x100 ) bw += 8, k >>= 8; \
! 267: if ( k >= 0x10 ) bw += 4, k >>= 4; \
! 268: if ( k >= 0x4 ) bw += 2, k >>= 2; \
! 269: if ( k >= 0x2 ) bw += 1, k >>= 1; \
! 270: if ( k != 0 ) bw++; \
! 271: }
! 272: #else
! 273: #define BitWidth(n,bw) {\
! 274: unsigned int k = (n); \
! 275: bw = 0; \
! 276: if ( k >= 0x10000 ) bw += 16, k >>= 16; \
! 277: if ( k >= 0x100 ) bw += 8, k >>= 8; \
! 278: if ( k >= 0x10 ) bw += 4, k >>= 4; \
! 279: if ( k >= 0x4 ) bw += 2, k >>= 2; \
! 280: if ( k >= 0x2 ) bw += 1, k >>= 1; \
! 281: if ( k != 0 ) bw++; \
! 282: }
! 283: #endif
! 284: #endif
! 285:
! 286: #include "igcdhack.c"
! 287:
! 288: /*
! 289: * Implementation of the binary GCD algorithm for two oN structs
! 290: * (big-integers) in risa.
! 291: *
! 292: * The major operations in the following algorithms are the binary-shifts
! 293: * and the updates of (u, v) by (min(u,v), |u-v|), and are to be open-coded
! 294: * without using routines for oN structures just as in addn() or subn().
! 295: */
! 296:
! 297: static int igcd_binary_2w( u, lu, v, lv, pans )
! 298: #if BSH == BIT_WIDTH_OF_INT
! 299: unsigned
! 300: #endif
! 301: int *u, lu, *v, lv, *pans; /* <---- */
! 302: /* both u[0:lu-1] and v[0:lv-1] are assumed to be odd */
! 303: {
! 304: #if BSH == BIT_WIDTH_OF_INT
! 305: unsigned
! 306: #endif
! 307: int i, h1, l1, h2, l2; /* <---- */
! 308:
! 309: l1 = u[0], l2 = v[0];
! 310: h1 = lu <= 1 ? 0 : u[1];
! 311: h2 = lv <= 1 ? 0 : v[1];
! 312: /**/
! 313: loop: if ( h1 == 0 ) {
! 314: no_hi1: if ( h2 == 0 ) goto one_word;
! 315: no_hi1n:if ( l1 == 1 ) return 0;
! 316: #if BSH == BIT_WIDTH_OF_INT
! 317: if ( l2 == l1 ) {
! 318: for ( l2 = h2; (l2&1) == 0; l2 >>= 1 ) ;
! 319: goto one_word;
! 320: } else if ( l2 < l1 ) {
! 321: l2 -= l1, h2--;
! 322: } else l2 -= l1;
! 323: i = 0; do { l2 >>= 1, i++; } while ( (l2&1) == 0 );
! 324: l2 |= (h2 << (BSH - i));
! 325: #else
! 326: if ( (l2 -= l1) == 0 ) {
! 327: for ( l2 = h2; (l2&1) == 0; l2 >>= 1 ) ;
! 328: goto one_word;
! 329: } else if ( l2 < 0 ) h2--, l2 += BASE;
! 330: i = 0; do { l2 >>= 1, i++; } while ( (l2&1) == 0 );
! 331: l2 |= ((h2 << (BSH - i)) & BMASK);
! 332: #endif
! 333: h2 >>= i;
! 334: goto no_hi1;
! 335: } else if ( h2 == 0 ) {
! 336: no_hi2: if ( l2 == 1 ) return 0;
! 337: #if BSH == BIT_WIDTH_OF_INT
! 338: if ( l1 == l2 ) {
! 339: for ( l1 = h1; (l1&1) == 0; l1 >>= 1 ) ;
! 340: goto one_word;
! 341: } else if ( l1 < l2 ) {
! 342: l1 -= l2, h1--;
! 343: } else l1 -= l2;
! 344: i = 0; do { l1 >>= 1, i++; } while ( (l1&1) == 0 );
! 345: l1 |= (h1 << (BSH - i));
! 346: #else
! 347: if ( (l1 -= l2) == 0 ) {
! 348: for ( l1 = h1; (l1&1) == 0; l1 >>= 1 ) ;
! 349: goto one_word;
! 350: } else if ( l1 < 0 ) h1--, l1 += BASE;
! 351: i = 0; do { l1 >>= 1, i++; } while ( (l1&1) == 0 );
! 352: l1 |= ((h1 << (BSH - i)) & BMASK);
! 353: #endif
! 354: if ( (h1 >>= i) == 0 ) goto one_word;
! 355: goto no_hi2;
! 356: } else if ( l1 == l2 ) {
! 357: if ( h1 == h2 ) {
! 358: pans[0] = l1, pans[1] = h1;
! 359: return 2;
! 360: } else if ( h1 > h2 ) {
! 361: for ( l1 = h1 - h2; (l1&1) == 0; l1 >>= 1 ) ;
! 362: goto no_hi1n;
! 363: } else {
! 364: for ( l2 = h2 - h1; (l2&1) == 0; l2 >>= 1 ) ;
! 365: goto no_hi2;
! 366: }
! 367: } else if ( h1 == h2 ) {
! 368: if ( l1 > l2 ) {
! 369: for ( l1 -= l2; (l1&1) == 0; l1 >>= 1 ) ;
! 370: goto no_hi1n;
! 371: } else {
! 372: for ( l2 -= l1; (l2&1) == 0; l2 >>= 1 ) ;
! 373: goto no_hi2;
! 374: }
! 375: } else if ( h1 > h2 ) {
! 376: h1 -= h2;
! 377: #if BSH == BIT_WIDTH_OF_INT
! 378: if ( l1 < l2 ) l1 -= l2, h1--;
! 379: else l1 -= l2;
! 380: i = 0; do { l1 >>= 1, i++; } while ( (l1&1) == 0 );
! 381: l1 |= (h1 << (BSH - i));
! 382: #else
! 383: if ( (l1 -= l2) < 0 ) h1--, l1 += BASE;
! 384: i = 0; do { l1 >>= 1, i++; } while ( (l1&1) == 0 );
! 385: l1 |= ((h1 << (BSH - i)) & BMASK);
! 386: #endif
! 387: h1 >>= i;
! 388: } else {
! 389: h2 -= h1;
! 390: #if BSH == BIT_WIDTH_OF_INT
! 391: if ( l2 < l1 ) l2 -= l1, h2--;
! 392: else l2 -= l1;
! 393: i = 0; do { l2 >>= 1, i++; } while ( (l2&1) == 0 );
! 394: l2 |= (h2 << (BSH - i));
! 395: #else
! 396: if ( (l2 -= l1) < 0 ) h2--, l2 += BASE;
! 397: i = 0; do { l2 >>= 1, i++; } while ( (l2&1) == 0 );
! 398: l2 |= ((h2 << (BSH - i)) & BMASK);
! 399: #endif
! 400: h2 >>= i;
! 401: }
! 402: goto loop;
! 403: one_word:
! 404: if ( l1 == 1 || l2 == 1 ) return 0;
! 405: else if ( l1 == l2 ) {
! 406: pans[0] = l1;
! 407: return 1;
! 408: }
! 409: one_word_neq:
! 410: if ( l1 > l2 ) {
! 411: l1 -= l2;
! 412: do { l1 >>= 1; } while ( (l1&1) == 0 );
! 413: goto one_word;
! 414: } else {
! 415: l2 -= l1;
! 416: do { l2 >>= 1; } while ( (l2&1) == 0 );
! 417: goto one_word;
! 418: }
! 419: }
! 420:
! 421: static N igcd_binary( n1, n2, nt )
! 422: N n1, n2, nt;
! 423: /* both n1 and n2 are assumed to be odd */
! 424: {
! 425: int l1, *b1, l2, *b2, *bt = BD(nt);
! 426: int l;
! 427:
! 428: if ( (l = cmpn( n1, n2 )) == 0 ) return n1;
! 429: else if ( l < 0 ) { SWAP( n1, n2, N ); }
! 430: IniDiv( n1, n2 );
! 431: if ( UNIN(n2) ) return 0;
! 432: l1 = PL(n1), b1 = BD(n1), l2 = PL(n2), b2 = BD(n2);
! 433: loop:
! 434: #if 0000
! 435: {
! 436: int i, b, w;
! 437:
! 438: printf( "===============\n" );
! 439: for ( i = 0; i < l1; i++ ) printf( "0x%08x ", b1[i] );
! 440: printf( "\n" );
! 441: for ( i = 0; i < l2; i++ ) printf( "0x%08x ", b2[i] );
! 442: printf( "\n" );
! 443: }
! 444: #endif
! 445: if ( l1 <= 2 && l2 <= 2 ) {
! 446: l = igcd_binary_2w( b1, l1, b2, l2, bt );
! 447: if ( l == 0 ) return 0;
! 448: PL(nt) = l;
! 449: return nt;
! 450: }
! 451: /**/
! 452: l = abs_U_V_maxrshift( b1, l1, b2, l2, bt );
! 453: /**/
! 454: if ( l == 0 ) {
! 455: PL(n1) = l1;
! 456: return n1;
! 457: } else if ( l > 0 ) {
! 458: l1 = l;
! 459: SWAP( b1, bt, int * ); SWAP( n1, nt, N );
! 460: } else {
! 461: l2 = -l;
! 462: SWAP( b2, bt, int * ); SWAP( n2, nt, N );
! 463: }
! 464: goto loop;
! 465: }
! 466:
! 467: #define RetTrueGCD(p,gw,gb,nr,l0) \
! 468: if (p==0) { l0: if (gw==0&&gb==0) { *(nr)=ONEN; return; } else p=ONEN; } \
! 469: *(nr) = N_of_n_lshifted_by_wb(p,gw,gb); \
! 470: return;
! 471:
! 472: void gcdbinn( n1, n2, nr )
! 473: N n1, n2, *nr;
! 474: {
! 475: int s1, s2, gw, gb, t1, t2;
! 476: int w1, w2;
! 477: N tn1, tn2, tnt, p;
! 478:
! 479: if ( !n1 ) {
! 480: *nr = n2;
! 481: return;
! 482: } else if ( !n2 ) {
! 483: *nr = n1;
! 484: return;
! 485: }
! 486: s1 = trailingzerosn( n1, &w1 );
! 487: s2 = trailingzerosn( n2, &w2 );
! 488: if ( w1 == w2 ) gw = w1, gb = s1 <= s2 ? s1 : s2;
! 489: else if ( w1 < w2 ) gw = w1, gb = s1;
! 490: else gw = w2, gb = s2;
! 491: /*
! 492: * true GCD must be multiplied by 2^{gw*BSH+gb}.
! 493: */
! 494: t1 = PL(n1) - w1;
! 495: t2 = PL(n2) - w2;
! 496: if ( t1 < t2 ) t1 = t2;
! 497: tn1 = W_NALLOC(t1); tn2 = W_NALLOC(t1); tnt = W_NALLOC(t1);
! 498: rshiftn( n1, w1, s1, tn1 );
! 499: rshiftn( n2, w2, s2, tn2 );
! 500: p = igcd_binary( tn1, tn2, tnt );
! 501: RetTrueGCD( p, gw, gb, nr, L0 )
! 502: }
! 503:
! 504:
! 505: /*
! 506: * The bmod gcd algorithm stated briefly in K.Weber's paper
! 507: * [ACM TOMS, Vol.21, No. 1 (1995), pp. 111-122].
! 508: * It replaces the subtraction (n1 - n2) in the binary algorithm
! 509: * by (n1 - S*n2) with such an S that (n1 - S*n2) \equiv 0 \bmod 2^BSH,
! 510: * which should improve the efficiency when n1 \gg n2.
! 511: */
! 512:
! 513: /* subsidiary routines */
! 514: #if BSH == BIT_WIDTH_OF_INT
! 515: #ifdef CALL
! 516: static int u_div_v_mod_2toBSH( u, v )
! 517: unsigned int u, v;
! 518: /*
! 519: * u/v mod 2^BSH.
! 520: */
! 521: {
! 522: unsigned int i, lsh_i, m;
! 523:
! 524: lsh_i = (sizeof(int) << 3) - 1;
! 525: m = i = 0;
! 526: do {
! 527: if ( u == 0 ) break;
! 528: if ( (u << lsh_i) != 0 ) {
! 529: m += (1 << i);
! 530: u -= (v << i);
! 531: }
! 532: lsh_i--;
! 533: } while ( ++i != BSH );
! 534: return m;
! 535: }
! 536:
! 537: #define Comp_U_div_V_mod_BASE(U,V,R) R = u_div_v_mod_2toBSH(U,V)
! 538: #else /* CALL */
! 539: #define Comp_U_div_V_mod_BASE(U,V,R) {\
! 540: unsigned int u = (U), v = (V), i, lsh; \
! 541: /* U and V are assumed to be odd */ \
! 542: i = R = 1, lsh = (sizeof(int) << 3) - 2; u = (u - v); \
! 543: do { if ( u == 0 ) break; \
! 544: if ( (u << lsh) != 0 ) R += (1 << i), u = (u - (v << i)); \
! 545: i++, lsh--; \
! 546: } while ( i < BSH ); \
! 547: }
! 548: #endif /* CALL */
! 549: #else
! 550: #ifdef CALL
! 551: static int u_div_v_mod_2tos( u, v, s )
! 552: int u, v, s;
! 553: /*
! 554: * u/v mod 2^s.
! 555: */
! 556: {
! 557: int i, lsh_i, mask, m;
! 558:
! 559: mask = (1 << s) - 1;
! 560: lsh_i = (sizeof(int) << 3) - 1;
! 561: m = i = 0;
! 562: u &= mask, v &= mask;
! 563: do {
! 564: if ( u == 0 ) break;
! 565: if ( (u << lsh_i) != 0 ) {
! 566: m += (1 << i);
! 567: u -= (v << i);
! 568: u &= mask;
! 569: }
! 570: lsh_i--;
! 571: } while ( ++i != s );
! 572: return m;
! 573: }
! 574:
! 575: #define Comp_U_div_V_mod_BASE(U,V,R) R = u_div_v_mod_2tos(U,V,BSH)
! 576: #else
! 577: #define Comp_U_div_V_mod_BASE(U,V,R) {\
! 578: int u = (U), v = (V), i, lsh; \
! 579: /* U and V are assumed to be odd */ \
! 580: i = R = 1, lsh = (sizeof(int) << 3) - 2; u = (u - v) & BMASK; \
! 581: do { if ( u == 0 ) break; \
! 582: if ( (u << lsh) != 0 ) R += (1 << i), u = (u - (v << i)) & BMASK; \
! 583: i++, lsh--; \
! 584: } while ( i < BSH ); \
! 585: }
! 586: #endif
! 587: #endif
! 588:
! 589:
! 590: static int bmod_n( nu, nv, na )
! 591: N nu, nv, na;
! 592: /*
! 593: * Computes (u[] \bmod v[]) >> (as much as possible) in r[].
! 594: */
! 595: {
! 596: #if BSH == BIT_WIDTH_OF_INT
! 597: unsigned int *u = BD(nu), *v = BD(nv), *r = BD(na);
! 598: unsigned int *p, a, t, l, v0, vh, bv, v0r;
! 599: int lu = PL(nu), lv = PL(nv), z;
! 600: #else
! 601: int *u = BD(nu), lu = PL(nu), *v = BD(nv), lv = PL(nv),
! 602: *r = BD(na);
! 603: int *p, a, t, l, z, v0, vh, bv, v0r;
! 604: #endif
! 605:
! 606: v0 = v[0];
! 607: if ( lv == 1 ) {
! 608: if ( lu == 1 ) a = u[0] % v0;
! 609: else {
! 610: p = &u[--lu];
! 611: #if BSH == BIT_WIDTH_OF_INT
! 612: a = (*p) % v0, t = (unsigned)(-((int)v0)) % v0;
! 613: #else
! 614: a = (*p) % v0, t = BASE % v0;
! 615: #endif
! 616: for ( ; --lu >= 0; a = l ) {
! 617: --p;
! 618: DMAR(a,t,*p,v0,l);
! 619: /* l <= (a*t + p[0])%v0 */
! 620: }
! 621: }
! 622: if ( a == 0 ) return 0;
! 623: while ( (a&1) == 0 ) a >>= 1;
! 624: *r = a;
! 625: return( PL(na) = 1 );
! 626: }
! 627: Comp_U_div_V_mod_BASE( 1, v0, v0r );
! 628: vh = v[lv -1];
! 629: BitWidth( vh, bv );
! 630: bv--;
! 631: t = 1 << bv;
! 632: l = lv + 1;
! 633: for ( z = -1; lu > l || lu == l && u[lu-1] >= t; z = -z ) {
! 634: #if BSH == BIT_WIDTH_OF_INT
! 635: a = (v0r*u[0]);
! 636: #else
! 637: a = (v0r*u[0])&BMASK;
! 638: #endif
! 639: /**/
! 640: #if 0000
! 641: {
! 642: int i;
! 643: for ( i = 0; i < lu; i++ ) printf( "0x%08x ", u[i] );
! 644: printf( "\n- a=0x%08x, %u*\n", a, a );
! 645: for ( i = 0; i < lv; i++ ) printf( "0x%08x ", v[i] );
! 646: printf( "\n=>\n" );
! 647: }
! 648: #endif
! 649: lu = abs_U_aV_maxrshift( u, lu, a, v, lv, r );
! 650: /**/
! 651: #if 0000
! 652: printf( "***lu=%d\n", lu );
! 653: if ( lu != 0 ) {
! 654: int i;
! 655: for ( i = 0; i < lu; i++ ) printf( "0x%08x ", r[i] );
! 656: printf( "\n" );
! 657: }
! 658: #endif
! 659: if ( lu == 0 ) return 0;
! 660: p = r;
! 661: r = u;
! 662: u = p;
! 663: }
! 664: if ( lu < lv ) goto ret;
! 665: t = u[lu-1];
! 666: if ( lu > lv ) l = BSH;
! 667: else if ( t < vh ) goto ret;
! 668: else l = 0;
! 669: BitWidth( t, a );
! 670: l += (a - bv);
! 671: #if BSH == BIT_WIDTH_OF_INT
! 672: a = (v0r*u[0])&((unsigned)(-1) >> (BSH - l));
! 673: #else
! 674: a = (v0r*u[0])&(BMASK >> (BSH - l));
! 675: #endif
! 676: #if 0000
! 677: {
! 678: int i;
! 679: for ( i = 0; i < lu; i++ ) printf( "0x%08x ", u[i] );
! 680: printf( "\n - a=0x%08x, %u*\n", a, a );
! 681: for ( i = 0; i < lv; i++ ) printf( "0x%08x ", v[i] );
! 682: printf( "\n =>\n" );
! 683: }
! 684: #endif
! 685: /**/
! 686: lu = abs_U_aV_maxrshift( u, lu, a, v, lv, r );
! 687: /**/
! 688: #if 0000
! 689: printf( "::: lu=%d\n", lu );
! 690: if ( lu != 0 ) {
! 691: int i;
! 692: for ( i = 0; i < lu; i++ ) printf( "0x%08x ", r[i] );
! 693: printf( "\n" );
! 694: }
! 695: #endif
! 696: if ( lu == 0 ) return 0;
! 697: z = -z;
! 698: ret: if ( z > 0 ) return( PL(na) = lu );
! 699: PL(nu) = lu;
! 700: return( -lu );
! 701: }
! 702:
! 703:
! 704: static N igcd_bmod( n1, n2, nt )
! 705: N n1, n2, nt;
! 706: /* both n1 and n2 are assumed to be odd */
! 707: {
! 708: int l1, l2;
! 709: int l;
! 710:
! 711: if ( (l = cmpn( n1, n2 )) == 0 ) return n1;
! 712: else if ( l < 0 ) { SWAP( n1, n2, N ); }
! 713: IniDiv( n1, n2 );
! 714: if ( UNIN(n2) ) return 0;
! 715: loop: if ( (l1 = PL(n1)) <= 2 && (l2 = PL(n2)) <= 2 ) {
! 716: l = igcd_binary_2w( BD(n1), l1, BD(n2), l2, BD(nt) );
! 717: if ( l == 0 ) return 0;
! 718: PL(nt) = l;
! 719: return nt;
! 720: }
! 721: /**/
! 722: l = bmod_n( n1, n2, nt );
! 723: /**/
! 724: if ( l == 0 ) return n2;
! 725: else if ( l > 0 ) {
! 726: N tmp = n1;
! 727:
! 728: n1 = n2;
! 729: n2 = nt;
! 730: nt = tmp;
! 731: } else SWAP( n1, n2, N );
! 732: goto loop;
! 733: }
! 734:
! 735: void gcdbmodn( n1, n2, nr )
! 736: N n1, n2, *nr;
! 737: {
! 738: int s1, s2, gw, gb, t1, t2;
! 739: int w1, w2;
! 740: N tn1, tn2, tnt, p;
! 741:
! 742: if ( !n1 ) {
! 743: *nr = n2;
! 744: return;
! 745: } else if ( !n2 ) {
! 746: *nr = n1;
! 747: return;
! 748: }
! 749: s1 = trailingzerosn( n1, &w1 );
! 750: s2 = trailingzerosn( n2, &w2 );
! 751: if ( w1 == w2 ) gw = w1, gb = s1 <= s2 ? s1 : s2;
! 752: else if ( w1 < w2 ) gw = w1, gb = s1;
! 753: else gw = w2, gb = s2;
! 754: /*
! 755: * true GCD must be multiplied by 2^{gw*BSH+gs}.
! 756: */
! 757: t1 = PL(n1) - w1;
! 758: t2 = PL(n2) - w2;
! 759: if ( t1 < t2 ) t1 = t2;
! 760: tn1 = W_NALLOC(t1); tn2 = W_NALLOC(t1); tnt = W_NALLOC(t1);
! 761: rshiftn( n1, w1, s1, tn1 );
! 762: rshiftn( n2, w2, s2, tn2 );
! 763: p = igcd_bmod( tn1, tn2, tnt );
! 764: RetTrueGCD( p, gw, gb, nr, L0 )
! 765: }
! 766:
! 767: /*
! 768: * The accelerated integer GCD algorithm by K.Weber
! 769: * [ACM TOMS, Vol.21, No. 1 (1995), pp. 111-122]:
! 770: */
! 771:
! 772: static int ReducedRatMod( x, y, pn, pd )
! 773: N x, y;
! 774: #if BSH == BIT_WIDTH_OF_INT
! 775: unsigned
! 776: #endif
! 777: int *pn, *pd;
! 778: /*
! 779: * Let m = 2^{2*BSH} = 2*BASE. We assume x, y > 0 and \gcd(x,m)
! 780: * = \gcd(y,m) = 1. This routine computes n and d (resp. returned
! 781: * in *pn and *pd) such that 0 < n, |d| < \sqrt{m} and
! 782: * n*y \equiv x*d \bmod m.
! 783: */
! 784: {
! 785: #if BSH == BIT_WIDTH_OF_INT
! 786: unsigned int n1h, n1l, d1h, d1l, n2h, n2l, d2h, d2l;
! 787: unsigned int th, tl, l1, l2, i, ir;
! 788: int s1, s2;
! 789: #else
! 790: int n1h, n1l, d1h, d1l, n2h, n2l, d2h, d2l;
! 791: int th, tl, l1, l2, i, ir;
! 792: int s1, s2;
! 793: #endif
! 794:
! 795: {
! 796: #if BSH == BIT_WIDTH_OF_INT
! 797: unsigned
! 798: #endif
! 799: int xh, xl, yh, yl, lsh_i;
! 800:
! 801: xl = BD(x)[0];
! 802: xh = PL(x) > 1 ? BD(x)[1] : 0;
! 803: yl = BD(y)[0];
! 804: yh = PL(y) > 1 ? BD(y)[1] : 0;
! 805: #if 0000
! 806: printf( "*** RedRatMod: (0x%08x:0x%08x=%u*2^%d+%u)\n /(0x%08x:0x%08x=%u*2^%d+%u) mod 2^%d\n",
! 807: xh,xl, xh,BSH,xl, yh,yl, yh,BSH,yl, BSH );
! 808: #endif
! 809: Comp_U_div_V_mod_BASE( xl, yl, n2l );
! 810: DM(n2l,yl,th,tl) /* n2l*yl = tl+th*BASE, where tl==xl. */;
! 811: #if BSH == BIT_WIDTH_OF_INT
! 812: xh -= th;
! 813: #else
! 814: if ( xh > th ) xh -= th;
! 815: else xh += (BASE - th);
! 816: #endif
! 817: DM(n2l,yh,th,tl) /* n2l*yh = tl+th*BASE. */;
! 818: #if BSH == BIT_WIDTH_OF_INT
! 819: xh -= tl;
! 820: #else
! 821: if ( xh > tl ) xh -= tl;
! 822: else xh += (BASE - tl);
! 823: #endif
! 824: /* n2h = i = 0, lsh_i = 31;*/
! 825: n2h = i = 0, lsh_i = BIT_WIDTH_OF_INT -1;
! 826: do {
! 827: if ( xh == 0 ) break;
! 828: if ( (xh << lsh_i) != 0 ) {
! 829: n2h += (1 << i);
! 830: #if BSH == BIT_WIDTH_OF_INT
! 831: tl = yl << i;
! 832: xh -= tl;
! 833: #else
! 834: tl = (yl << i)&BMASK;
! 835: if ( xh > tl ) xh -= tl;
! 836: else xh += (BASE - tl);
! 837: #endif
! 838: }
! 839: lsh_i--;
! 840: } while ( ++i != BSH );
! 841: }
! 842: /*
! 843: * n2l + n2h*BASE = x/y mod 2^{2*BSH}.
! 844: */
! 845: #if 0000
! 846: printf( "=====> 0x%08x(%u) + 2^%d*0x%08x(%u)\n", n2l, n2l, BSH, n2h, n2h );
! 847: #endif
! 848: d2h = 0, d2l = 1, s2 = 1;
! 849: #if BSH == BIT_WIDTH_OF_INT
! 850: if ( n2h == 0 ) goto done;
! 851: BitWidth( n2h, l2 );
! 852: if ( l2 == BSH ) {
! 853: d1h = 0, d1l = 1, s1 = -1;
! 854: th = n2h, tl = n2l;
! 855: } else {
! 856: i = BSH - l2;
! 857: d1h = 0, d1l = 1 << i, s1 = -1;
! 858: th = (n2h << i) | (n2l >> l2);
! 859: tl = n2l << i;
! 860: }
! 861: if ( tl == 0 ) n1h = - th, n1l = 0;
! 862: else n1h = BMASK - th, n1l = - tl;
! 863: /**/
! 864: if ( n1h < n2h || (n1h == n2h && n1l < n2l) ) goto swap12;
! 865: BitWidth( n1h, l1 );
! 866: goto sub12;
! 867: #else
! 868: n1h = BASE, n1l = 0, l1 = BSH,
! 869: d1h = d1l = 0, s1 = 0;
! 870: #endif
! 871: /**/
! 872: while ( n2h != 0 ) {
! 873: BitWidth( n2h, l2 );
! 874: sub12: ir = BSH - (i = l1 - l2);
! 875: do {
! 876: if ( i == 0 ) th = n2h, tl = n2l;
! 877: else
! 878: th = (n2h << i) | (n2l >> ir),
! 879: #if BSH == BIT_WIDTH_OF_INT
! 880: tl = n2l << i;
! 881: #else
! 882: tl = (n2l << i) & BMASK;
! 883: #endif
! 884: if ( th > n1h || (th == n1h && tl > n1l) ) goto next_i;
! 885: #if BSH == BIT_WIDTH_OF_INT
! 886: if ( tl > n1l ) n1h--;
! 887: n1l -= tl;
! 888: #else
! 889: if ( tl <= n1l ) n1l -= tl;
! 890: else n1l += (BASE - tl), n1h--;
! 891: #endif
! 892: n1h -= th;
! 893: /* (s1:d1h,d1l) -= ((s2:d2h,d2l) << i); */
! 894: if ( s2 != 0 ) {
! 895: if ( i == 0 ) th = d2h, tl = d2l;
! 896: else
! 897: #if BSH == BIT_WIDTH_OF_INT
! 898: th = (d2h << i) | (d2l >> ir),
! 899: tl = (d2l << i);
! 900: #else
! 901: th = (d2h << i)&BMASK | (d2l >> ir),
! 902: tl = (d2l << i)&BMASK;
! 903: #endif
! 904: if ( s1 == 0 )
! 905: s1 = -s2, d1h = th, d1l = tl;
! 906: else if ( s1 != s2 ) {
! 907: #if BSH == BIT_WIDTH_OF_INT
! 908: d1l += tl,
! 909: d1h += th;
! 910: if ( d1l < tl ) d1h++;
! 911: #else
! 912: tl += d1l;
! 913: d1l = tl&BMASK;
! 914: d1h = (d1h + th + (tl >> BSH))&BMASK;
! 915: #endif
! 916: if ( d1h == 0 && d1l == 0 ) s1 = 0;
! 917: } else if ( d1h > th ) {
! 918: #if BSH == BIT_WIDTH_OF_INT
! 919: if ( d1l < tl ) d1h--;
! 920: d1l -= tl;
! 921: #else
! 922: if ( d1l >= tl ) d1l -= tl;
! 923: else d1l += (BASE - tl), d1h--;
! 924: #endif
! 925: d1h -= th;
! 926: } else if ( d1h == th ) {
! 927: d1h = 0;
! 928: if ( d1l == tl ) s1 = d2h = 0;
! 929: else if ( d1l > tl ) d1l -= tl;
! 930: else d1l = tl - d1l, s1 = -s1;
! 931: } else {
! 932: #if BSH == BIT_WIDTH_OF_INT
! 933: if ( tl < d1l ) th--;
! 934: d1l = tl - d1l;
! 935: #else
! 936: if ( tl >= d1l ) d1l = tl - d1l;
! 937: else d1l = tl + (BASE - d1l), th--;
! 938: #endif
! 939: d1h = th - d1h;
! 940: s1 = -s1;
! 941: }
! 942: }
! 943: next_i: i--, ir++;
! 944: } while ( n1h > n2h || (n1h == n2h && n1l >= n2l) );
! 945: swap12: /* swap 1 and 2 */
! 946: th = n1h, tl = n1l;
! 947: n1h = n2h, n1l = n2l;
! 948: n2h = th, n2l = tl;
! 949: l1 = l2;
! 950: th = d1h, tl = d1l, i = s1;
! 951: d1h = d2h, d1l = d2l, s1 = s2;
! 952: d2h = th, d2l = tl, s2 = i;
! 953: }
! 954: /**/
! 955: done: *pn = n2l, *pd = d2l;
! 956: return s2;
! 957: }
! 958:
! 959: static int igcd_spurious_factor;
! 960:
! 961: #define SaveN(s,d) {\
! 962: int i, l; \
! 963: for ( l = PL(d) = PL(s), i = 0; i < l; i++ ) BD(d)[i] = BD(s)[i]; \
! 964: }
! 965:
! 966: static N igcd_acc( n1, n2, nt )
! 967: N n1, n2, nt;
! 968: /* both n1 and n2 are assumed to be odd */
! 969: {
! 970: int l1, l2, *b1, *b2, bw1, bw2;
! 971: int l;
! 972: int n, d;
! 973: N p, s1, s2;
! 974:
! 975: if ( (l = cmpn( n1, n2 )) == 0 ) return n1;
! 976: else if ( l < 0 ) { SWAP( n1, n2, N ); }
! 977: if ( ShouldCompRemInit(n1,n2) ) {
! 978: int w, b;
! 979:
! 980: divn( n1, n2, &s1, &s2 );
! 981: if ( !s2 ) return n2;
! 982: b = trailingzerosn( s2, &w );
! 983: p = n1; n1 = n2; n2 = p;
! 984: rshiftn( s2, w, b, n2 );
! 985: if ( UNIN(n2) ) return 0;
! 986: l1 = PL(n1);
! 987: if ( !s1 || PL(s1) < l1 ) s1 = NALLOC(l1);
! 988: } else if ( UNIN(n2) ) return 0;
! 989: else {
! 990: s1 = NALLOC(PL(n1));
! 991: s2 = NALLOC(PL(n2));
! 992: }
! 993: SaveN( n1, s1 );
! 994: SaveN( n2, s2 );
! 995: igcd_spurious_factor = 0;
! 996: loop: l1 = PL(n1), l2 = PL(n2);
! 997: if ( l1 <= 2 && l2 <= 2 ) {
! 998: l = igcd_binary_2w( BD(n1), l1, BD(n2), l2, BD(nt) );
! 999: if ( l == 0 ) return 0;
! 1000: PL(nt) = l;
! 1001: SWAP( n2, nt, N );
! 1002: goto ret;
! 1003: }
! 1004: /**/
! 1005: b1 = BD(n1), b2 = BD(n2);
! 1006: BitWidth( b1[l1 -1], bw1 );
! 1007: BitWidth( b2[l2 -1], bw2 );
! 1008: if ( (l1*BSH + bw1) - (l2*BSH + bw2) <= igcdacc_thre ) {
! 1009: l = ReducedRatMod( n1, n2, &n, &d );
! 1010: l = l < 0 ? aUplusbV_maxrshift( n, b2, l2, d, b1, l1, BD(nt) ) :
! 1011: abs_axU_bxV_maxrshift( n, b2, l2, d, b1, l1, BD(nt) );
! 1012: igcd_spurious_factor++;
! 1013: if ( l == 0 ) goto ret;
! 1014: PL(nt) = l;
! 1015: } else {
! 1016: l = bmod_n( n1, n2, nt );
! 1017: if ( l == 0 ) goto ret;
! 1018: else if ( l < 0 ) {
! 1019: SWAP( n1, n2, N );
! 1020: goto loop;
! 1021: }
! 1022: }
! 1023: p = n1;
! 1024: n1 = n2;
! 1025: n2 = nt;
! 1026: nt = p;
! 1027: goto loop;
! 1028: /**/
! 1029: ret: if ( igcd_spurious_factor != 0 && !UNIN(n2) ) {
! 1030: if ( (p = igcd_bmod( n2, s1, n1 )) == 0 ) return 0;
! 1031: if ( (p = igcd_bmod( p, s2, nt )) == 0 ) return 0;
! 1032: return p;
! 1033: } else return n2;
! 1034: }
! 1035:
! 1036: void gcdaccn( n1, n2, nr )
! 1037: N n1, n2, *nr;
! 1038: {
! 1039: int s1, s2, gw, gb, t1, t2;
! 1040: int w1, w2;
! 1041: N tn1, tn2, tnt, p;
! 1042:
! 1043: if ( !n1 ) {
! 1044: *nr = n2;
! 1045: return;
! 1046: } else if ( !n2 ) {
! 1047: *nr = n1;
! 1048: return;
! 1049: }
! 1050: s1 = trailingzerosn( n1, &w1 );
! 1051: s2 = trailingzerosn( n2, &w2 );
! 1052: if ( w1 == w2 ) gw = w1, gb = s1 <= s2 ? s1 : s2;
! 1053: else if ( w1 < w2 ) gw = w1, gb = s1;
! 1054: else gw = w2, gb = s2;
! 1055: /*
! 1056: * true GCD must be multiplied by 2^{gw*BSH+gs}.
! 1057: */
! 1058: t1 = PL(n1) - w1;
! 1059: t2 = PL(n2) - w2;
! 1060: if ( t1 < t2 ) t1 = t2;
! 1061: tn1 = W_NALLOC(t1); tn2 = W_NALLOC(t1); tnt = W_NALLOC(t1);
! 1062: rshiftn( n1, w1, s1, tn1 );
! 1063: rshiftn( n2, w2, s2, tn2 );
! 1064: /**/
! 1065: p = igcd_acc( tn1, tn2, tnt );
! 1066: /**/
! 1067: if ( p == 0 ) goto L0;
! 1068: RetTrueGCD( p, gw, gb, nr, L0 )
! 1069: }
! 1070:
! 1071:
! 1072: /********************************/
! 1073:
! 1074: void gcdn_HMEXT( n1, n2, nr )
! 1075: N n1, n2, *nr;
! 1076: {
! 1077: int b1, b2, w1, w2, gw, gb;
! 1078: int l1, l2;
! 1079: N tn1, tn2, tnt, a;
! 1080:
! 1081: if ( !n1 ) {
! 1082: *nr = n2; return;
! 1083: } else if ( !n2 ) {
! 1084: *nr = n1; return;
! 1085: }
! 1086: b1 = trailingzerosn( n1, &w1 );
! 1087: b2 = trailingzerosn( n2, &w2 );
! 1088: if ( w1 == w2 ) gw = w1, gb = b1 <= b2 ? b1 : b2;
! 1089: else if ( w1 < w2 ) gw = w1, gb = b1;
! 1090: else gw = w2, gb = b2;
! 1091: /*
! 1092: * true GCD must be multiplied by 2^{gw*BSH+gb}.
! 1093: */
! 1094: l1 = PL(n1) - w1;
! 1095: l2 = PL(n2) - w2;
! 1096: if ( l1 < l2 ) l1 = l2;
! 1097: tn1 = W_NALLOC( l1 ); tn2 = W_NALLOC( l1 ); tnt = W_NALLOC( l1 );
! 1098: rshiftn( n1, w1, b1, tn1 );
! 1099: rshiftn( n2, w2, b2, tn2 );
! 1100: /**/
! 1101: if ( igcd_algorithm == 1 ) {
! 1102: a = igcd_binary( tn1, tn2, tnt );
! 1103: } else if ( igcd_algorithm == 2 ) {
! 1104: a = igcd_bmod( tn1, tn2, tnt );
! 1105: } else {
! 1106: a = igcd_acc( tn1, tn2, tnt );
! 1107: if ( igcd_spurious_factor != 0 ) {
! 1108: }
! 1109: }
! 1110: RetTrueGCD( a, gw, gb, nr, L0 )
! 1111: }
! 1112:
! 1113:
! 1114: /**************************/
! 1115: #if 111
! 1116: N maxrshn( n, p )
! 1117: N n;
! 1118: int *p;
! 1119: {
! 1120: int nw, nb, c, l;
! 1121: N new;
! 1122:
! 1123: nb = trailingzerosn( n, &nw );
! 1124: l = PL(n);
! 1125: c = BD(n)[l -1];
! 1126: l -= nw;
! 1127: if ( (c >> nb) == 0 ) l--;
! 1128: new = NALLOC(l);
! 1129: rshiftn( n, nw, nb, new );
! 1130: *p = nb + nw*BSH;
! 1131: return new;
! 1132: }
! 1133: #endif
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to Ngcd.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine