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