version 1.3, 2001/10/09 01:36:11 |
version 1.4, 2018/03/29 01:32:51 |
|
|
/* |
/* |
* $OpenXM: OpenXM_contrib2/asir2000/engine/Ngcd.c,v 1.2 2000/12/21 02:56:31 murao Exp $ |
* $OpenXM: OpenXM_contrib2/asir2000/engine/Ngcd.c,v 1.3 2001/10/09 01:36:11 noro Exp $ |
*/ |
*/ |
/* |
/* |
#include "ca.h" |
#include "ca.h" |
|
|
|
|
|
|
int igcd_algorithm = 0; |
int igcd_algorithm = 0; |
/* == 0 : Euclid, |
/* == 0 : Euclid, |
* == 1 : binary, |
* == 1 : binary, |
* == 2 : bmod, |
* == 2 : bmod, |
* >= 3 : (Weber's accelerated)/(Jebelean's generalized binary) algorithm, |
* >= 3 : (Weber's accelerated)/(Jebelean's generalized binary) algorithm, |
*/ |
*/ |
int igcd_thre_inidiv = 50; |
int igcd_thre_inidiv = 50; |
/* |
/* |
Line 49 int igcdacc_thre = 10; |
|
Line 49 int igcdacc_thre = 10; |
|
#define ShouldCompRemInit(n1,n2) (igcd_thre_inidiv != 0 && PL(n1) >= igcd_thre_inidiv*PL(n2)) |
#define ShouldCompRemInit(n1,n2) (igcd_thre_inidiv != 0 && PL(n1) >= igcd_thre_inidiv*PL(n2)) |
|
|
#define IniDiv(n1,n2) \ |
#define IniDiv(n1,n2) \ |
if ( ShouldCompRemInit(n1,n2) ) {\ |
if ( ShouldCompRemInit(n1,n2) ) {\ |
N q, r; int w, b; \ |
N q, r; int w, b; \ |
divn(n1,n2,&q,&r); \ |
divn(n1,n2,&q,&r); \ |
if ( !r ) return(n2); \ |
if ( !r ) return(n2); \ |
b = trailingzerosn( r, &w ); \ |
b = trailingzerosn( r, &w ); \ |
q = n1; n1 = n2; n2 = q; \ |
q = n1; n1 = n2; n2 = q; \ |
rshiftn( r, w, b, n2 ); \ |
rshiftn( r, w, b, n2 ); \ |
} |
} |
|
|
/* |
/* |
* Binary GCD algorithm by J.Stein |
* Binary GCD algorithm by J.Stein |
* [J. Comp. Phys. Vol. 1 (1967), pp. 397-405)]: |
* [J. Comp. Phys. Vol. 1 (1967), pp. 397-405)]: |
* The right-shift binary algorithm is used. |
* The right-shift binary algorithm is used. |
*/ |
*/ |
|
|
|
|
/* |
/* |
* subsidiary routines for gcdbinn below. |
* subsidiary routines for gcdbinn below. |
*/ |
*/ |
static int /* number of bits */ trailingzeros_nbd( /* BD of N */ nbd, pnw ) |
static int /* number of bits */ trailingzeros_nbd( /* BD of N */ nbd, pnw ) |
int *nbd, *pnw /* number of zero words */; |
int *nbd, *pnw /* number of zero words */; |
{ |
{ |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
unsigned |
unsigned |
#endif |
#endif |
int nw, nb, w; |
int nw, nb, w; |
|
|
for ( nw = 0; (w = *nbd) == 0; nbd++ ) nw++; |
for ( nw = 0; (w = *nbd) == 0; nbd++ ) nw++; |
TRAILINGZEROS(w,nb); |
TRAILINGZEROS(w,nb); |
*pnw = nw; |
*pnw = nw; |
return nb; |
return nb; |
} |
} |
|
|
#define trailingzerosn(n,pnw) trailingzeros_nbd(BD(n),pnw) |
#define trailingzerosn(n,pnw) trailingzeros_nbd(BD(n),pnw) |
|
|
static int /* PL of N */ rshift_nbd( /* BD of N */ nbd, /* PL of N */ nl, |
static int /* PL of N */ rshift_nbd( /* BD of N */ nbd, /* PL of N */ nl, |
/* # words */ shw, /* # bits */ shb, /* BD of N */ p ) |
/* # words */ shw, /* # bits */ shb, /* BD of N */ p ) |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
unsigned |
unsigned |
#endif |
#endif |
int *nbd, nl, shw, shb, *p; |
int *nbd, nl, shw, shb, *p; |
{ |
{ |
unsigned int i, v, w, lshb; /* <---- */ |
unsigned int i, v, w, lshb; /* <---- */ |
|
|
nbd += shw, i = (nl -= shw); |
nbd += shw, i = (nl -= shw); |
if ( shb == 0 ) { |
if ( shb == 0 ) { |
for ( ; nl > 0; nl-- ) *p++ = *nbd++; |
for ( ; nl > 0; nl-- ) *p++ = *nbd++; |
return i; |
return i; |
} else if ( nl < 2 ) { |
} else if ( nl < 2 ) { |
*p = (*nbd) >> shb; |
*p = (*nbd) >> shb; |
return 1; |
return 1; |
} |
} |
for ( lshb = BSH - shb, v = *nbd++; --nl > 0; v = w ) { |
for ( lshb = BSH - shb, v = *nbd++; --nl > 0; v = w ) { |
w = *nbd++; |
w = *nbd++; |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
*p++ = (v >> shb) | (w << lshb); /********/ |
*p++ = (v >> shb) | (w << lshb); /********/ |
#else |
#else |
*p++ = (v >> shb) | ((w << lshb)&BMASK); |
*p++ = (v >> shb) | ((w << lshb)&BMASK); |
#endif |
#endif |
} |
} |
if ( (v >>= shb) == 0 ) return( i-1 ); |
if ( (v >>= shb) == 0 ) return( i-1 ); |
*p = v; |
*p = v; |
return i; |
return i; |
} |
} |
|
|
#define rshiftn(ns,shw,shb,nd) (PL(nd)=rshift_nbd(BD(ns),PL(ns),shw,shb,BD(nd))) |
#define rshiftn(ns,shw,shb,nd) (PL(nd)=rshift_nbd(BD(ns),PL(ns),shw,shb,BD(nd))) |
Line 120 int *nbd, nl, shw, shb, *p; |
|
Line 120 int *nbd, nl, shw, shb, *p; |
|
#ifdef FULLSET |
#ifdef FULLSET |
static N N_of_i_lshifted_by_wb( i, gw, gb ) |
static N N_of_i_lshifted_by_wb( i, gw, gb ) |
int i, gw, gb; |
int i, gw, gb; |
/* |
/* |
* returns pointer to a new struct (N)(((int)i) >> (gb + gw*BSH)) |
* returns pointer to a new struct (N)(((int)i) >> (gb + gw*BSH)) |
*/ |
*/ |
{ |
{ |
unsigned int j, l, *p; /* <---- */ |
unsigned int j, l, *p; /* <---- */ |
N n; |
N n; |
|
|
j = i >> (BSH - gb); |
j = i >> (BSH - gb); |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
i = (i << gb); /********/ |
i = (i << gb); /********/ |
#else |
#else |
i = (i << gb)&BMASK; |
i = (i << gb)&BMASK; |
#endif |
#endif |
l = j != 0 ? gw + 2 : gw + 1; |
l = j != 0 ? gw + 2 : gw + 1; |
n = NALLOC(l); |
n = NALLOC(l); |
PL(n) = l; |
PL(n) = l; |
for ( p = BD(n); gw-- > 0; ) *p++ = 0; |
for ( p = BD(n); gw-- > 0; ) *p++ = 0; |
*p++ = i; |
*p++ = i; |
if ( j != 0 ) *p = j; |
if ( j != 0 ) *p = j; |
return n; |
return n; |
} |
} |
#endif /* FULLSET */ |
#endif /* FULLSET */ |
|
|
/* |
/* |
* routines to make a new struct |
* routines to make a new struct |
* (N)(((BD of N)(b[0],...,b[lb-1])) << (gb + gw*BSH)) |
* (N)(((BD of N)(b[0],...,b[lb-1])) << (gb + gw*BSH)) |
*/ |
*/ |
static N N_of_nbd_lshifted_by_wb( /* BD of N */ b, /* PL of N */ lb, gw, gb ) |
static N N_of_nbd_lshifted_by_wb( /* BD of N */ b, /* PL of N */ lb, gw, gb ) |
int *b, lb, gw, gb; |
int *b, lb, gw, gb; |
/* |
/* |
* returns pointer to a new struct |
* returns pointer to a new struct |
* (N)(((BD of N)(b[0],...,b[lb-1])) << (gb + gw*BSH)) |
* (N)(((BD of N)(b[0],...,b[lb-1])) << (gb + gw*BSH)) |
*/ |
*/ |
{ |
{ |
unsigned int rsh, s, t, *p, l; /* <---- */ |
unsigned int rsh, s, t, *p, l; /* <---- */ |
N n; |
N n; |
|
|
l = lb + gw; |
l = lb + gw; |
if ( gb == 0 ) { |
if ( gb == 0 ) { |
n = NALLOC(l); |
n = NALLOC(l); |
PL(n) = l; |
PL(n) = l; |
for ( p = BD(n); gw-- > 0; ) *p++ = 0; |
for ( p = BD(n); gw-- > 0; ) *p++ = 0; |
while ( lb-- > 0 ) *p++ = *b++; |
while ( lb-- > 0 ) *p++ = *b++; |
return n; |
return n; |
} |
} |
rsh = BSH - gb; s = b[lb-1]; |
rsh = BSH - gb; s = b[lb-1]; |
if ( (t = s >> rsh) != 0 ) { |
if ( (t = s >> rsh) != 0 ) { |
n = NALLOC(l+1); |
n = NALLOC(l+1); |
PL(n) = l+1; |
PL(n) = l+1; |
(p = BD(n))[l] = t; |
(p = BD(n))[l] = t; |
} else { |
} else { |
n = NALLOC(l); |
n = NALLOC(l); |
PL(n) = l; |
PL(n) = l; |
p = BD(n); |
p = BD(n); |
} |
} |
while ( gw-- > 0 ) *p++ = 0; |
while ( gw-- > 0 ) *p++ = 0; |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
*p++ = (t = *b++) << gb; /********/ |
*p++ = (t = *b++) << gb; /********/ |
for ( ; --lb > 0; t = s ) |
for ( ; --lb > 0; t = s ) |
*p++ = (t >> rsh) | ((s = *b++) << gb); /********/ |
*p++ = (t >> rsh) | ((s = *b++) << gb); /********/ |
#else |
#else |
*p++ = ((t = *b++) << gb)&BMASK; |
*p++ = ((t = *b++) << gb)&BMASK; |
for ( ; --lb > 0; t = s ) |
for ( ; --lb > 0; t = s ) |
*p++ = (t >> rsh) | (((s = *b++) << gb)&BMASK); |
*p++ = (t >> rsh) | (((s = *b++) << gb)&BMASK); |
#endif |
#endif |
return n; |
return n; |
} |
} |
|
|
#define N_of_n_lshifted_by_wb(a,gw,gb) N_of_nbd_lshifted_by_wb(BD(a),PL(a),gw,gb) |
#define N_of_n_lshifted_by_wb(a,gw,gb) N_of_nbd_lshifted_by_wb(BD(a),PL(a),gw,gb) |
Line 196 int *b, lb, gw, gb; |
|
Line 196 int *b, lb, gw, gb; |
|
|
|
#ifdef CALL |
#ifdef CALL |
static int bw_int32( n ) |
static int bw_int32( n ) |
unsigned int n; /* <---- */ |
unsigned int n; /* <---- */ |
{ |
{ |
int w; |
int w; |
|
|
w = 0; |
w = 0; |
#if BSH > 32 |
#if BSH > 32 |
if ( n > 0xffffffff ) w += 32, n >>= 32; |
if ( n > 0xffffffff ) w += 32, n >>= 32; |
#endif |
#endif |
if ( n >= 0x10000 ) w += 16, n >>= 16; |
if ( n >= 0x10000 ) w += 16, n >>= 16; |
if ( n >= 0x100 ) w += 8, n >>= 8; |
if ( n >= 0x100 ) w += 8, n >>= 8; |
if ( n >= 0x10 ) w += 4, n >>= 4; |
if ( n >= 0x10 ) w += 4, n >>= 4; |
if ( n >= 0x4 ) w += 2, n >>= 2; |
if ( n >= 0x4 ) w += 2, n >>= 2; |
if ( n >= 0x2 ) w += 1, n >>= 1; |
if ( n >= 0x2 ) w += 1, n >>= 1; |
if ( n != 0 ) ++w; |
if ( n != 0 ) ++w; |
return w; |
return w; |
} |
} |
#define BitWidth(n,bw) bw = bw_int32( n ) |
#define BitWidth(n,bw) bw = bw_int32( n ) |
#else |
#else |
|
|
#if BSH > 32 |
#if BSH > 32 |
#define BitWidth(n,bw) {\ |
#define BitWidth(n,bw) {\ |
unsigned int k = (n); \ |
unsigned int k = (n); \ |
bw = 0; \ |
bw = 0; \ |
if ( k > 0xffffffff ) bw += 32, k >>= 32; \ |
if ( k > 0xffffffff ) bw += 32, k >>= 32; \ |
if ( k >= 0x10000 ) bw += 16, k >>= 16; \ |
if ( k >= 0x10000 ) bw += 16, k >>= 16; \ |
if ( k >= 0x100 ) bw += 8, k >>= 8; \ |
if ( k >= 0x100 ) bw += 8, k >>= 8; \ |
if ( k >= 0x10 ) bw += 4, k >>= 4; \ |
if ( k >= 0x10 ) bw += 4, k >>= 4; \ |
if ( k >= 0x4 ) bw += 2, k >>= 2; \ |
if ( k >= 0x4 ) bw += 2, k >>= 2; \ |
if ( k >= 0x2 ) bw += 1, k >>= 1; \ |
if ( k >= 0x2 ) bw += 1, k >>= 1; \ |
if ( k != 0 ) bw++; \ |
if ( k != 0 ) bw++; \ |
} |
} |
#else |
#else |
#define BitWidth(n,bw) {\ |
#define BitWidth(n,bw) {\ |
unsigned int k = (n); \ |
unsigned int k = (n); \ |
bw = 0; \ |
bw = 0; \ |
if ( k >= 0x10000 ) bw += 16, k >>= 16; \ |
if ( k >= 0x10000 ) bw += 16, k >>= 16; \ |
if ( k >= 0x100 ) bw += 8, k >>= 8; \ |
if ( k >= 0x100 ) bw += 8, k >>= 8; \ |
if ( k >= 0x10 ) bw += 4, k >>= 4; \ |
if ( k >= 0x10 ) bw += 4, k >>= 4; \ |
if ( k >= 0x4 ) bw += 2, k >>= 2; \ |
if ( k >= 0x4 ) bw += 2, k >>= 2; \ |
if ( k >= 0x2 ) bw += 1, k >>= 1; \ |
if ( k >= 0x2 ) bw += 1, k >>= 1; \ |
if ( k != 0 ) bw++; \ |
if ( k != 0 ) bw++; \ |
} |
} |
#endif |
#endif |
#endif |
#endif |
Line 244 unsigned int n; /* <---- */ |
|
Line 244 unsigned int n; /* <---- */ |
|
#include "igcdhack.c" |
#include "igcdhack.c" |
|
|
/* |
/* |
* Implementation of the binary GCD algorithm for two oN structs |
* Implementation of the binary GCD algorithm for two oN structs |
* (big-integers) in risa. |
* (big-integers) in risa. |
* |
* |
* The major operations in the following algorithms are the binary-shifts |
* The major operations in the following algorithms are the binary-shifts |
* and the updates of (u, v) by (min(u,v), |u-v|), and are to be open-coded |
* and the updates of (u, v) by (min(u,v), |u-v|), and are to be open-coded |
* without using routines for oN structures just as in addn() or subn(). |
* without using routines for oN structures just as in addn() or subn(). |
*/ |
*/ |
|
|
static int igcd_binary_2w( u, lu, v, lv, pans ) |
static int igcd_binary_2w( u, lu, v, lv, pans ) |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
unsigned |
unsigned |
#endif |
#endif |
int *u, lu, *v, lv, *pans; /* <---- */ |
int *u, lu, *v, lv, *pans; /* <---- */ |
/* both u[0:lu-1] and v[0:lv-1] are assumed to be odd */ |
/* both u[0:lu-1] and v[0:lv-1] are assumed to be odd */ |
{ |
{ |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
unsigned |
unsigned |
#endif |
#endif |
int i, h1, l1, h2, l2; /* <---- */ |
int i, h1, l1, h2, l2; /* <---- */ |
|
|
l1 = u[0], l2 = v[0]; |
l1 = u[0], l2 = v[0]; |
h1 = lu <= 1 ? 0 : u[1]; |
h1 = lu <= 1 ? 0 : u[1]; |
h2 = lv <= 1 ? 0 : v[1]; |
h2 = lv <= 1 ? 0 : v[1]; |
/**/ |
/**/ |
loop: if ( h1 == 0 ) { |
loop: if ( h1 == 0 ) { |
no_hi1: if ( h2 == 0 ) goto one_word; |
no_hi1: if ( h2 == 0 ) goto one_word; |
no_hi1n:if ( l1 == 1 ) return 0; |
no_hi1n:if ( l1 == 1 ) return 0; |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
if ( l2 == l1 ) { |
if ( l2 == l1 ) { |
for ( l2 = h2; (l2&1) == 0; l2 >>= 1 ) ; |
for ( l2 = h2; (l2&1) == 0; l2 >>= 1 ) ; |
goto one_word; |
goto one_word; |
} else if ( l2 < l1 ) { |
} else if ( l2 < l1 ) { |
l2 -= l1, h2--; |
l2 -= l1, h2--; |
} else l2 -= l1; |
} else l2 -= l1; |
i = 0; do { l2 >>= 1, i++; } while ( (l2&1) == 0 ); |
i = 0; do { l2 >>= 1, i++; } while ( (l2&1) == 0 ); |
l2 |= (h2 << (BSH - i)); |
l2 |= (h2 << (BSH - i)); |
#else |
#else |
if ( (l2 -= l1) == 0 ) { |
if ( (l2 -= l1) == 0 ) { |
for ( l2 = h2; (l2&1) == 0; l2 >>= 1 ) ; |
for ( l2 = h2; (l2&1) == 0; l2 >>= 1 ) ; |
goto one_word; |
goto one_word; |
} else if ( l2 < 0 ) h2--, l2 += BASE; |
} else if ( l2 < 0 ) h2--, l2 += BASE; |
i = 0; do { l2 >>= 1, i++; } while ( (l2&1) == 0 ); |
i = 0; do { l2 >>= 1, i++; } while ( (l2&1) == 0 ); |
l2 |= ((h2 << (BSH - i)) & BMASK); |
l2 |= ((h2 << (BSH - i)) & BMASK); |
#endif |
#endif |
h2 >>= i; |
h2 >>= i; |
goto no_hi1; |
goto no_hi1; |
} else if ( h2 == 0 ) { |
} else if ( h2 == 0 ) { |
no_hi2: if ( l2 == 1 ) return 0; |
no_hi2: if ( l2 == 1 ) return 0; |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
if ( l1 == l2 ) { |
if ( l1 == l2 ) { |
for ( l1 = h1; (l1&1) == 0; l1 >>= 1 ) ; |
for ( l1 = h1; (l1&1) == 0; l1 >>= 1 ) ; |
goto one_word; |
goto one_word; |
} else if ( l1 < l2 ) { |
} else if ( l1 < l2 ) { |
l1 -= l2, h1--; |
l1 -= l2, h1--; |
} else l1 -= l2; |
} else l1 -= l2; |
i = 0; do { l1 >>= 1, i++; } while ( (l1&1) == 0 ); |
i = 0; do { l1 >>= 1, i++; } while ( (l1&1) == 0 ); |
l1 |= (h1 << (BSH - i)); |
l1 |= (h1 << (BSH - i)); |
#else |
#else |
if ( (l1 -= l2) == 0 ) { |
if ( (l1 -= l2) == 0 ) { |
for ( l1 = h1; (l1&1) == 0; l1 >>= 1 ) ; |
for ( l1 = h1; (l1&1) == 0; l1 >>= 1 ) ; |
goto one_word; |
goto one_word; |
} else if ( l1 < 0 ) h1--, l1 += BASE; |
} else if ( l1 < 0 ) h1--, l1 += BASE; |
i = 0; do { l1 >>= 1, i++; } while ( (l1&1) == 0 ); |
i = 0; do { l1 >>= 1, i++; } while ( (l1&1) == 0 ); |
l1 |= ((h1 << (BSH - i)) & BMASK); |
l1 |= ((h1 << (BSH - i)) & BMASK); |
#endif |
#endif |
if ( (h1 >>= i) == 0 ) goto one_word; |
if ( (h1 >>= i) == 0 ) goto one_word; |
goto no_hi2; |
goto no_hi2; |
} else if ( l1 == l2 ) { |
} else if ( l1 == l2 ) { |
if ( h1 == h2 ) { |
if ( h1 == h2 ) { |
pans[0] = l1, pans[1] = h1; |
pans[0] = l1, pans[1] = h1; |
return 2; |
return 2; |
} else if ( h1 > h2 ) { |
} else if ( h1 > h2 ) { |
for ( l1 = h1 - h2; (l1&1) == 0; l1 >>= 1 ) ; |
for ( l1 = h1 - h2; (l1&1) == 0; l1 >>= 1 ) ; |
goto no_hi1n; |
goto no_hi1n; |
} else { |
} else { |
for ( l2 = h2 - h1; (l2&1) == 0; l2 >>= 1 ) ; |
for ( l2 = h2 - h1; (l2&1) == 0; l2 >>= 1 ) ; |
goto no_hi2; |
goto no_hi2; |
} |
} |
} else if ( h1 == h2 ) { |
} else if ( h1 == h2 ) { |
if ( l1 > l2 ) { |
if ( l1 > l2 ) { |
for ( l1 -= l2; (l1&1) == 0; l1 >>= 1 ) ; |
for ( l1 -= l2; (l1&1) == 0; l1 >>= 1 ) ; |
goto no_hi1n; |
goto no_hi1n; |
} else { |
} else { |
for ( l2 -= l1; (l2&1) == 0; l2 >>= 1 ) ; |
for ( l2 -= l1; (l2&1) == 0; l2 >>= 1 ) ; |
goto no_hi2; |
goto no_hi2; |
} |
} |
} else if ( h1 > h2 ) { |
} else if ( h1 > h2 ) { |
h1 -= h2; |
h1 -= h2; |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
if ( l1 < l2 ) l1 -= l2, h1--; |
if ( l1 < l2 ) l1 -= l2, h1--; |
else l1 -= l2; |
else l1 -= l2; |
i = 0; do { l1 >>= 1, i++; } while ( (l1&1) == 0 ); |
i = 0; do { l1 >>= 1, i++; } while ( (l1&1) == 0 ); |
l1 |= (h1 << (BSH - i)); |
l1 |= (h1 << (BSH - i)); |
#else |
#else |
if ( (l1 -= l2) < 0 ) h1--, l1 += BASE; |
if ( (l1 -= l2) < 0 ) h1--, l1 += BASE; |
i = 0; do { l1 >>= 1, i++; } while ( (l1&1) == 0 ); |
i = 0; do { l1 >>= 1, i++; } while ( (l1&1) == 0 ); |
l1 |= ((h1 << (BSH - i)) & BMASK); |
l1 |= ((h1 << (BSH - i)) & BMASK); |
#endif |
#endif |
h1 >>= i; |
h1 >>= i; |
} else { |
} else { |
h2 -= h1; |
h2 -= h1; |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
if ( l2 < l1 ) l2 -= l1, h2--; |
if ( l2 < l1 ) l2 -= l1, h2--; |
else l2 -= l1; |
else l2 -= l1; |
i = 0; do { l2 >>= 1, i++; } while ( (l2&1) == 0 ); |
i = 0; do { l2 >>= 1, i++; } while ( (l2&1) == 0 ); |
l2 |= (h2 << (BSH - i)); |
l2 |= (h2 << (BSH - i)); |
#else |
#else |
if ( (l2 -= l1) < 0 ) h2--, l2 += BASE; |
if ( (l2 -= l1) < 0 ) h2--, l2 += BASE; |
i = 0; do { l2 >>= 1, i++; } while ( (l2&1) == 0 ); |
i = 0; do { l2 >>= 1, i++; } while ( (l2&1) == 0 ); |
l2 |= ((h2 << (BSH - i)) & BMASK); |
l2 |= ((h2 << (BSH - i)) & BMASK); |
#endif |
#endif |
h2 >>= i; |
h2 >>= i; |
} |
} |
goto loop; |
goto loop; |
one_word: |
one_word: |
if ( l1 == 1 || l2 == 1 ) return 0; |
if ( l1 == 1 || l2 == 1 ) return 0; |
else if ( l1 == l2 ) { |
else if ( l1 == l2 ) { |
pans[0] = l1; |
pans[0] = l1; |
return 1; |
return 1; |
} |
} |
one_word_neq: |
one_word_neq: |
if ( l1 > l2 ) { |
if ( l1 > l2 ) { |
l1 -= l2; |
l1 -= l2; |
do { l1 >>= 1; } while ( (l1&1) == 0 ); |
do { l1 >>= 1; } while ( (l1&1) == 0 ); |
goto one_word; |
goto one_word; |
} else { |
} else { |
l2 -= l1; |
l2 -= l1; |
do { l2 >>= 1; } while ( (l2&1) == 0 ); |
do { l2 >>= 1; } while ( (l2&1) == 0 ); |
goto one_word; |
goto one_word; |
} |
} |
} |
} |
|
|
static N igcd_binary( n1, n2, nt ) |
static N igcd_binary( n1, n2, nt ) |
N n1, n2, nt; |
N n1, n2, nt; |
/* both n1 and n2 are assumed to be odd */ |
/* both n1 and n2 are assumed to be odd */ |
{ |
{ |
int l1, *b1, l2, *b2, *bt = BD(nt); |
int l1, *b1, l2, *b2, *bt = BD(nt); |
int l; |
int l; |
|
|
if ( (l = cmpn( n1, n2 )) == 0 ) return n1; |
if ( (l = cmpn( n1, n2 )) == 0 ) return n1; |
else if ( l < 0 ) { SWAP( n1, n2, N ); } |
else if ( l < 0 ) { SWAP( n1, n2, N ); } |
IniDiv( n1, n2 ); |
IniDiv( n1, n2 ); |
if ( UNIN(n2) ) return 0; |
if ( UNIN(n2) ) return 0; |
l1 = PL(n1), b1 = BD(n1), l2 = PL(n2), b2 = BD(n2); |
l1 = PL(n1), b1 = BD(n1), l2 = PL(n2), b2 = BD(n2); |
loop: |
loop: |
#if 0000 |
#if 0000 |
{ |
{ |
|
|
printf( "\n" ); |
printf( "\n" ); |
} |
} |
#endif |
#endif |
if ( l1 <= 2 && l2 <= 2 ) { |
if ( l1 <= 2 && l2 <= 2 ) { |
l = igcd_binary_2w( b1, l1, b2, l2, bt ); |
l = igcd_binary_2w( b1, l1, b2, l2, bt ); |
if ( l == 0 ) return 0; |
if ( l == 0 ) return 0; |
PL(nt) = l; |
PL(nt) = l; |
return nt; |
return nt; |
} |
} |
/**/ |
/**/ |
l = abs_U_V_maxrshift( b1, l1, b2, l2, bt ); |
l = abs_U_V_maxrshift( b1, l1, b2, l2, bt ); |
/**/ |
/**/ |
if ( l == 0 ) { |
if ( l == 0 ) { |
PL(n1) = l1; |
PL(n1) = l1; |
return n1; |
return n1; |
} else if ( l > 0 ) { |
} else if ( l > 0 ) { |
l1 = l; |
l1 = l; |
SWAP( b1, bt, int * ); SWAP( n1, nt, N ); |
SWAP( b1, bt, int * ); SWAP( n1, nt, N ); |
} else { |
} else { |
l2 = -l; |
l2 = -l; |
SWAP( b2, bt, int * ); SWAP( n2, nt, N ); |
SWAP( b2, bt, int * ); SWAP( n2, nt, N ); |
} |
} |
goto loop; |
goto loop; |
} |
} |
|
|
#define RetTrueGCD(p,gw,gb,nr,l0) \ |
#define RetTrueGCD(p,gw,gb,nr,l0) \ |
|
|
void gcdbinn( n1, n2, nr ) |
void gcdbinn( n1, n2, nr ) |
N n1, n2, *nr; |
N n1, n2, *nr; |
{ |
{ |
int s1, s2, gw, gb, t1, t2; |
int s1, s2, gw, gb, t1, t2; |
int w1, w2; |
int w1, w2; |
N tn1, tn2, tnt, p; |
N tn1, tn2, tnt, p; |
|
|
if ( !n1 ) { |
if ( !n1 ) { |
*nr = n2; |
*nr = n2; |
return; |
return; |
} else if ( !n2 ) { |
} else if ( !n2 ) { |
*nr = n1; |
*nr = n1; |
return; |
return; |
} |
} |
s1 = trailingzerosn( n1, &w1 ); |
s1 = trailingzerosn( n1, &w1 ); |
s2 = trailingzerosn( n2, &w2 ); |
s2 = trailingzerosn( n2, &w2 ); |
if ( w1 == w2 ) gw = w1, gb = s1 <= s2 ? s1 : s2; |
if ( w1 == w2 ) gw = w1, gb = s1 <= s2 ? s1 : s2; |
else if ( w1 < w2 ) gw = w1, gb = s1; |
else if ( w1 < w2 ) gw = w1, gb = s1; |
else gw = w2, gb = s2; |
else gw = w2, gb = s2; |
/* |
/* |
* true GCD must be multiplied by 2^{gw*BSH+gb}. |
* true GCD must be multiplied by 2^{gw*BSH+gb}. |
*/ |
*/ |
t1 = PL(n1) - w1; |
t1 = PL(n1) - w1; |
t2 = PL(n2) - w2; |
t2 = PL(n2) - w2; |
if ( t1 < t2 ) t1 = t2; |
if ( t1 < t2 ) t1 = t2; |
tn1 = W_NALLOC(t1); tn2 = W_NALLOC(t1); tnt = W_NALLOC(t1); |
tn1 = W_NALLOC(t1); tn2 = W_NALLOC(t1); tnt = W_NALLOC(t1); |
rshiftn( n1, w1, s1, tn1 ); |
rshiftn( n1, w1, s1, tn1 ); |
rshiftn( n2, w2, s2, tn2 ); |
rshiftn( n2, w2, s2, tn2 ); |
p = igcd_binary( tn1, tn2, tnt ); |
p = igcd_binary( tn1, tn2, tnt ); |
RetTrueGCD( p, gw, gb, nr, L0 ) |
RetTrueGCD( p, gw, gb, nr, L0 ) |
} |
} |
|
|
|
|
/* |
/* |
* The bmod gcd algorithm stated briefly in K.Weber's paper |
* The bmod gcd algorithm stated briefly in K.Weber's paper |
* [ACM TOMS, Vol.21, No. 1 (1995), pp. 111-122]. |
* [ACM TOMS, Vol.21, No. 1 (1995), pp. 111-122]. |
* It replaces the subtraction (n1 - n2) in the binary algorithm |
* It replaces the subtraction (n1 - n2) in the binary algorithm |
* by (n1 - S*n2) with such an S that (n1 - S*n2) \equiv 0 \bmod 2^BSH, |
* by (n1 - S*n2) with such an S that (n1 - S*n2) \equiv 0 \bmod 2^BSH, |
* which should improve the efficiency when n1 \gg n2. |
* which should improve the efficiency when n1 \gg n2. |
*/ |
*/ |
|
|
/* subsidiary routines */ |
/* subsidiary routines */ |
|
|
#ifdef CALL |
#ifdef CALL |
static int u_div_v_mod_2toBSH( u, v ) |
static int u_div_v_mod_2toBSH( u, v ) |
unsigned int u, v; |
unsigned int u, v; |
/* |
/* |
* u/v mod 2^BSH. |
* u/v mod 2^BSH. |
*/ |
*/ |
{ |
{ |
unsigned int i, lsh_i, m; |
unsigned int i, lsh_i, m; |
|
|
lsh_i = (sizeof(int) << 3) - 1; |
lsh_i = (sizeof(int) << 3) - 1; |
m = i = 0; |
m = i = 0; |
do { |
do { |
if ( u == 0 ) break; |
if ( u == 0 ) break; |
if ( (u << lsh_i) != 0 ) { |
if ( (u << lsh_i) != 0 ) { |
m += (1 << i); |
m += (1 << i); |
u -= (v << i); |
u -= (v << i); |
} |
} |
lsh_i--; |
lsh_i--; |
} while ( ++i != BSH ); |
} while ( ++i != BSH ); |
return m; |
return m; |
} |
} |
|
|
#define Comp_U_div_V_mod_BASE(U,V,R) R = u_div_v_mod_2toBSH(U,V) |
#define Comp_U_div_V_mod_BASE(U,V,R) R = u_div_v_mod_2toBSH(U,V) |
#else /* CALL */ |
#else /* CALL */ |
#define Comp_U_div_V_mod_BASE(U,V,R) {\ |
#define Comp_U_div_V_mod_BASE(U,V,R) {\ |
unsigned int u = (U), v = (V), i, lsh; \ |
unsigned int u = (U), v = (V), i, lsh; \ |
/* U and V are assumed to be odd */ \ |
/* U and V are assumed to be odd */ \ |
i = R = 1, lsh = (sizeof(int) << 3) - 2; u = (u - v); \ |
i = R = 1, lsh = (sizeof(int) << 3) - 2; u = (u - v); \ |
do { if ( u == 0 ) break; \ |
do { if ( u == 0 ) break; \ |
if ( (u << lsh) != 0 ) R += (1 << i), u = (u - (v << i)); \ |
if ( (u << lsh) != 0 ) R += (1 << i), u = (u - (v << i)); \ |
i++, lsh--; \ |
i++, lsh--; \ |
} while ( i < BSH ); \ |
} while ( i < BSH ); \ |
} |
} |
#endif /* CALL */ |
#endif /* CALL */ |
#else |
#else |
#ifdef CALL |
#ifdef CALL |
static int u_div_v_mod_2tos( u, v, s ) |
static int u_div_v_mod_2tos( u, v, s ) |
int u, v, s; |
int u, v, s; |
/* |
/* |
* u/v mod 2^s. |
* u/v mod 2^s. |
*/ |
*/ |
{ |
{ |
int i, lsh_i, mask, m; |
int i, lsh_i, mask, m; |
|
|
mask = (1 << s) - 1; |
mask = (1 << s) - 1; |
lsh_i = (sizeof(int) << 3) - 1; |
lsh_i = (sizeof(int) << 3) - 1; |
m = i = 0; |
m = i = 0; |
u &= mask, v &= mask; |
u &= mask, v &= mask; |
do { |
do { |
if ( u == 0 ) break; |
if ( u == 0 ) break; |
if ( (u << lsh_i) != 0 ) { |
if ( (u << lsh_i) != 0 ) { |
m += (1 << i); |
m += (1 << i); |
u -= (v << i); |
u -= (v << i); |
u &= mask; |
u &= mask; |
} |
} |
lsh_i--; |
lsh_i--; |
} while ( ++i != s ); |
} while ( ++i != s ); |
return m; |
return m; |
} |
} |
|
|
#define Comp_U_div_V_mod_BASE(U,V,R) R = u_div_v_mod_2tos(U,V,BSH) |
#define Comp_U_div_V_mod_BASE(U,V,R) R = u_div_v_mod_2tos(U,V,BSH) |
#else |
#else |
#define Comp_U_div_V_mod_BASE(U,V,R) {\ |
#define Comp_U_div_V_mod_BASE(U,V,R) {\ |
int u = (U), v = (V), i, lsh; \ |
int u = (U), v = (V), i, lsh; \ |
/* U and V are assumed to be odd */ \ |
/* U and V are assumed to be odd */ \ |
i = R = 1, lsh = (sizeof(int) << 3) - 2; u = (u - v) & BMASK; \ |
i = R = 1, lsh = (sizeof(int) << 3) - 2; u = (u - v) & BMASK; \ |
do { if ( u == 0 ) break; \ |
do { if ( u == 0 ) break; \ |
if ( (u << lsh) != 0 ) R += (1 << i), u = (u - (v << i)) & BMASK; \ |
if ( (u << lsh) != 0 ) R += (1 << i), u = (u - (v << i)) & BMASK; \ |
i++, lsh--; \ |
i++, lsh--; \ |
} while ( i < BSH ); \ |
} while ( i < BSH ); \ |
} |
} |
#endif |
#endif |
#endif |
#endif |
|
|
|
|
static int bmod_n( nu, nv, na ) |
static int bmod_n( nu, nv, na ) |
N nu, nv, na; |
N nu, nv, na; |
/* |
/* |
* Computes (u[] \bmod v[]) >> (as much as possible) in r[]. |
* Computes (u[] \bmod v[]) >> (as much as possible) in r[]. |
*/ |
*/ |
{ |
{ |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
unsigned int *u = BD(nu), *v = BD(nv), *r = BD(na); |
unsigned int *u = BD(nu), *v = BD(nv), *r = BD(na); |
unsigned int *p, a, t, l, v0, vh, bv, v0r; |
unsigned int *p, a, t, l, v0, vh, bv, v0r; |
int lu = PL(nu), lv = PL(nv), z; |
int lu = PL(nu), lv = PL(nv), z; |
#else |
#else |
int *u = BD(nu), lu = PL(nu), *v = BD(nv), lv = PL(nv), |
int *u = BD(nu), lu = PL(nu), *v = BD(nv), lv = PL(nv), |
*r = BD(na); |
*r = BD(na); |
int *p, a, t, l, z, v0, vh, bv, v0r; |
int *p, a, t, l, z, v0, vh, bv, v0r; |
#endif |
#endif |
|
|
v0 = v[0]; |
v0 = v[0]; |
if ( lv == 1 ) { |
if ( lv == 1 ) { |
if ( lu == 1 ) a = u[0] % v0; |
if ( lu == 1 ) a = u[0] % v0; |
else { |
else { |
p = &u[--lu]; |
p = &u[--lu]; |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
a = (*p) % v0, t = (unsigned)(-((int)v0)) % v0; |
a = (*p) % v0, t = (unsigned)(-((int)v0)) % v0; |
#else |
#else |
a = (*p) % v0, t = BASE % v0; |
a = (*p) % v0, t = BASE % v0; |
#endif |
#endif |
for ( ; --lu >= 0; a = l ) { |
for ( ; --lu >= 0; a = l ) { |
--p; |
--p; |
DMAR(a,t,*p,v0,l); |
DMAR(a,t,*p,v0,l); |
/* l <= (a*t + p[0])%v0 */ |
/* l <= (a*t + p[0])%v0 */ |
} |
} |
} |
} |
if ( a == 0 ) return 0; |
if ( a == 0 ) return 0; |
while ( (a&1) == 0 ) a >>= 1; |
while ( (a&1) == 0 ) a >>= 1; |
*r = a; |
*r = a; |
return( PL(na) = 1 ); |
return( PL(na) = 1 ); |
} |
} |
Comp_U_div_V_mod_BASE( 1, v0, v0r ); |
Comp_U_div_V_mod_BASE( 1, v0, v0r ); |
vh = v[lv -1]; |
vh = v[lv -1]; |
BitWidth( vh, bv ); |
BitWidth( vh, bv ); |
bv--; |
bv--; |
t = 1 << bv; |
t = 1 << bv; |
l = lv + 1; |
l = lv + 1; |
for ( z = -1; lu > l || lu == l && u[lu-1] >= t; z = -z ) { |
for ( z = -1; lu > l || lu == l && u[lu-1] >= t; z = -z ) { |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
a = (v0r*u[0]); |
a = (v0r*u[0]); |
#else |
#else |
a = (v0r*u[0])&BMASK; |
a = (v0r*u[0])&BMASK; |
#endif |
#endif |
/**/ |
/**/ |
#if 0000 |
#if 0000 |
{ |
{ |
int i; |
int i; |
|
|
printf( "\n=>\n" ); |
printf( "\n=>\n" ); |
} |
} |
#endif |
#endif |
lu = abs_U_aV_maxrshift( u, lu, a, v, lv, r ); |
lu = abs_U_aV_maxrshift( u, lu, a, v, lv, r ); |
/**/ |
/**/ |
#if 0000 |
#if 0000 |
printf( "***lu=%d\n", lu ); |
printf( "***lu=%d\n", lu ); |
if ( lu != 0 ) { |
if ( lu != 0 ) { |
Line 614 if ( lu != 0 ) { |
|
Line 614 if ( lu != 0 ) { |
|
printf( "\n" ); |
printf( "\n" ); |
} |
} |
#endif |
#endif |
if ( lu == 0 ) return 0; |
if ( lu == 0 ) return 0; |
p = r; |
p = r; |
r = u; |
r = u; |
u = p; |
u = p; |
} |
} |
if ( lu < lv ) goto ret; |
if ( lu < lv ) goto ret; |
t = u[lu-1]; |
t = u[lu-1]; |
if ( lu > lv ) l = BSH; |
if ( lu > lv ) l = BSH; |
else if ( t < vh ) goto ret; |
else if ( t < vh ) goto ret; |
else l = 0; |
else l = 0; |
BitWidth( t, a ); |
BitWidth( t, a ); |
l += (a - bv); |
l += (a - bv); |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
a = (v0r*u[0])&((unsigned)(-1) >> (BSH - l)); |
a = (v0r*u[0])&((unsigned)(-1) >> (BSH - l)); |
#else |
#else |
a = (v0r*u[0])&(BMASK >> (BSH - l)); |
a = (v0r*u[0])&(BMASK >> (BSH - l)); |
#endif |
#endif |
#if 0000 |
#if 0000 |
{ |
{ |
Line 640 if ( lu != 0 ) { |
|
Line 640 if ( lu != 0 ) { |
|
printf( "\n =>\n" ); |
printf( "\n =>\n" ); |
} |
} |
#endif |
#endif |
/**/ |
/**/ |
lu = abs_U_aV_maxrshift( u, lu, a, v, lv, r ); |
lu = abs_U_aV_maxrshift( u, lu, a, v, lv, r ); |
/**/ |
/**/ |
#if 0000 |
#if 0000 |
printf( "::: lu=%d\n", lu ); |
printf( "::: lu=%d\n", lu ); |
if ( lu != 0 ) { |
if ( lu != 0 ) { |
Line 651 if ( lu != 0 ) { |
|
Line 651 if ( lu != 0 ) { |
|
printf( "\n" ); |
printf( "\n" ); |
} |
} |
#endif |
#endif |
if ( lu == 0 ) return 0; |
if ( lu == 0 ) return 0; |
z = -z; |
z = -z; |
ret: if ( z > 0 ) return( PL(na) = lu ); |
ret: if ( z > 0 ) return( PL(na) = lu ); |
PL(nu) = lu; |
PL(nu) = lu; |
return( -lu ); |
return( -lu ); |
} |
} |
|
|
|
|
static N igcd_bmod( n1, n2, nt ) |
static N igcd_bmod( n1, n2, nt ) |
N n1, n2, nt; |
N n1, n2, nt; |
/* both n1 and n2 are assumed to be odd */ |
/* both n1 and n2 are assumed to be odd */ |
{ |
{ |
int l1, l2; |
int l1, l2; |
int l; |
int l; |
|
|
if ( (l = cmpn( n1, n2 )) == 0 ) return n1; |
if ( (l = cmpn( n1, n2 )) == 0 ) return n1; |
else if ( l < 0 ) { SWAP( n1, n2, N ); } |
else if ( l < 0 ) { SWAP( n1, n2, N ); } |
IniDiv( n1, n2 ); |
IniDiv( n1, n2 ); |
if ( UNIN(n2) ) return 0; |
if ( UNIN(n2) ) return 0; |
loop: if ( (l1 = PL(n1)) <= 2 && (l2 = PL(n2)) <= 2 ) { |
loop: if ( (l1 = PL(n1)) <= 2 && (l2 = PL(n2)) <= 2 ) { |
l = igcd_binary_2w( BD(n1), l1, BD(n2), l2, BD(nt) ); |
l = igcd_binary_2w( BD(n1), l1, BD(n2), l2, BD(nt) ); |
if ( l == 0 ) return 0; |
if ( l == 0 ) return 0; |
PL(nt) = l; |
PL(nt) = l; |
return nt; |
return nt; |
} |
} |
/**/ |
/**/ |
l = bmod_n( n1, n2, nt ); |
l = bmod_n( n1, n2, nt ); |
/**/ |
/**/ |
if ( l == 0 ) return n2; |
if ( l == 0 ) return n2; |
else if ( l > 0 ) { |
else if ( l > 0 ) { |
N tmp = n1; |
N tmp = n1; |
|
|
n1 = n2; |
n1 = n2; |
n2 = nt; |
n2 = nt; |
nt = tmp; |
nt = tmp; |
} else SWAP( n1, n2, N ); |
} else SWAP( n1, n2, N ); |
goto loop; |
goto loop; |
} |
} |
|
|
void gcdbmodn( n1, n2, nr ) |
void gcdbmodn( n1, n2, nr ) |
N n1, n2, *nr; |
N n1, n2, *nr; |
{ |
{ |
int s1, s2, gw, gb, t1, t2; |
int s1, s2, gw, gb, t1, t2; |
int w1, w2; |
int w1, w2; |
N tn1, tn2, tnt, p; |
N tn1, tn2, tnt, p; |
|
|
if ( !n1 ) { |
if ( !n1 ) { |
*nr = n2; |
*nr = n2; |
return; |
return; |
} else if ( !n2 ) { |
} else if ( !n2 ) { |
*nr = n1; |
*nr = n1; |
return; |
return; |
} |
} |
s1 = trailingzerosn( n1, &w1 ); |
s1 = trailingzerosn( n1, &w1 ); |
s2 = trailingzerosn( n2, &w2 ); |
s2 = trailingzerosn( n2, &w2 ); |
if ( w1 == w2 ) gw = w1, gb = s1 <= s2 ? s1 : s2; |
if ( w1 == w2 ) gw = w1, gb = s1 <= s2 ? s1 : s2; |
else if ( w1 < w2 ) gw = w1, gb = s1; |
else if ( w1 < w2 ) gw = w1, gb = s1; |
else gw = w2, gb = s2; |
else gw = w2, gb = s2; |
/* |
/* |
* true GCD must be multiplied by 2^{gw*BSH+gs}. |
* true GCD must be multiplied by 2^{gw*BSH+gs}. |
*/ |
*/ |
t1 = PL(n1) - w1; |
t1 = PL(n1) - w1; |
t2 = PL(n2) - w2; |
t2 = PL(n2) - w2; |
if ( t1 < t2 ) t1 = t2; |
if ( t1 < t2 ) t1 = t2; |
tn1 = W_NALLOC(t1); tn2 = W_NALLOC(t1); tnt = W_NALLOC(t1); |
tn1 = W_NALLOC(t1); tn2 = W_NALLOC(t1); tnt = W_NALLOC(t1); |
rshiftn( n1, w1, s1, tn1 ); |
rshiftn( n1, w1, s1, tn1 ); |
rshiftn( n2, w2, s2, tn2 ); |
rshiftn( n2, w2, s2, tn2 ); |
p = igcd_bmod( tn1, tn2, tnt ); |
p = igcd_bmod( tn1, tn2, tnt ); |
RetTrueGCD( p, gw, gb, nr, L0 ) |
RetTrueGCD( p, gw, gb, nr, L0 ) |
} |
} |
|
|
/* |
/* |
* The accelerated integer GCD algorithm by K.Weber |
* The accelerated integer GCD algorithm by K.Weber |
* [ACM TOMS, Vol.21, No. 1 (1995), pp. 111-122]: |
* [ACM TOMS, Vol.21, No. 1 (1995), pp. 111-122]: |
*/ |
*/ |
|
|
static int ReducedRatMod( x, y, pn, pd ) |
static int ReducedRatMod( x, y, pn, pd ) |
|
|
unsigned |
unsigned |
#endif |
#endif |
int *pn, *pd; |
int *pn, *pd; |
/* |
/* |
* Let m = 2^{2*BSH} = 2*BASE. We assume x, y > 0 and \gcd(x,m) |
* Let m = 2^{2*BSH} = 2*BASE. We assume x, y > 0 and \gcd(x,m) |
* = \gcd(y,m) = 1. This routine computes n and d (resp. returned |
* = \gcd(y,m) = 1. This routine computes n and d (resp. returned |
* in *pn and *pd) such that 0 < n, |d| < \sqrt{m} and |
* in *pn and *pd) such that 0 < n, |d| < \sqrt{m} and |
* n*y \equiv x*d \bmod m. |
* n*y \equiv x*d \bmod m. |
*/ |
*/ |
{ |
{ |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
unsigned int n1h, n1l, d1h, d1l, n2h, n2l, d2h, d2l; |
unsigned int n1h, n1l, d1h, d1l, n2h, n2l, d2h, d2l; |
unsigned int th, tl, l1, l2, i, ir; |
unsigned int th, tl, l1, l2, i, ir; |
int s1, s2; |
int s1, s2; |
#else |
#else |
int n1h, n1l, d1h, d1l, n2h, n2l, d2h, d2l; |
int n1h, n1l, d1h, d1l, n2h, n2l, d2h, d2l; |
int th, tl, l1, l2, i, ir; |
int th, tl, l1, l2, i, ir; |
int s1, s2; |
int s1, s2; |
#endif |
#endif |
|
|
{ |
{ |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
unsigned |
unsigned |
#endif |
#endif |
int xh, xl, yh, yl, lsh_i; |
int xh, xl, yh, yl, lsh_i; |
|
|
xl = BD(x)[0]; |
xl = BD(x)[0]; |
xh = PL(x) > 1 ? BD(x)[1] : 0; |
xh = PL(x) > 1 ? BD(x)[1] : 0; |
yl = BD(y)[0]; |
yl = BD(y)[0]; |
yh = PL(y) > 1 ? BD(y)[1] : 0; |
yh = PL(y) > 1 ? BD(y)[1] : 0; |
#if 0000 |
#if 0000 |
printf( "*** RedRatMod: (0x%08x:0x%08x=%u*2^%d+%u)\n /(0x%08x:0x%08x=%u*2^%d+%u) mod 2^%d\n", |
printf( "*** RedRatMod: (0x%08x:0x%08x=%u*2^%d+%u)\n /(0x%08x:0x%08x=%u*2^%d+%u) mod 2^%d\n", |
xh,xl, xh,BSH,xl, yh,yl, yh,BSH,yl, BSH ); |
xh,xl, xh,BSH,xl, yh,yl, yh,BSH,yl, BSH ); |
#endif |
#endif |
Comp_U_div_V_mod_BASE( xl, yl, n2l ); |
Comp_U_div_V_mod_BASE( xl, yl, n2l ); |
DM(n2l,yl,th,tl) /* n2l*yl = tl+th*BASE, where tl==xl. */; |
DM(n2l,yl,th,tl) /* n2l*yl = tl+th*BASE, where tl==xl. */; |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
xh -= th; |
xh -= th; |
#else |
#else |
if ( xh > th ) xh -= th; |
if ( xh > th ) xh -= th; |
else xh += (BASE - th); |
else xh += (BASE - th); |
#endif |
#endif |
DM(n2l,yh,th,tl) /* n2l*yh = tl+th*BASE. */; |
DM(n2l,yh,th,tl) /* n2l*yh = tl+th*BASE. */; |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
xh -= tl; |
xh -= tl; |
#else |
#else |
if ( xh > tl ) xh -= tl; |
if ( xh > tl ) xh -= tl; |
else xh += (BASE - tl); |
else xh += (BASE - tl); |
#endif |
#endif |
/* n2h = i = 0, lsh_i = 31;*/ |
/* n2h = i = 0, lsh_i = 31;*/ |
n2h = i = 0, lsh_i = BIT_WIDTH_OF_INT -1; |
n2h = i = 0, lsh_i = BIT_WIDTH_OF_INT -1; |
do { |
do { |
if ( xh == 0 ) break; |
if ( xh == 0 ) break; |
if ( (xh << lsh_i) != 0 ) { |
if ( (xh << lsh_i) != 0 ) { |
n2h += (1 << i); |
n2h += (1 << i); |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
tl = yl << i; |
tl = yl << i; |
xh -= tl; |
xh -= tl; |
#else |
#else |
tl = (yl << i)&BMASK; |
tl = (yl << i)&BMASK; |
if ( xh > tl ) xh -= tl; |
if ( xh > tl ) xh -= tl; |
else xh += (BASE - tl); |
else xh += (BASE - tl); |
#endif |
#endif |
} |
} |
lsh_i--; |
lsh_i--; |
} while ( ++i != BSH ); |
} while ( ++i != BSH ); |
} |
} |
/* |
/* |
* n2l + n2h*BASE = x/y mod 2^{2*BSH}. |
* n2l + n2h*BASE = x/y mod 2^{2*BSH}. |
*/ |
*/ |
#if 0000 |
#if 0000 |
printf( "=====> 0x%08x(%u) + 2^%d*0x%08x(%u)\n", n2l, n2l, BSH, n2h, n2h ); |
printf( "=====> 0x%08x(%u) + 2^%d*0x%08x(%u)\n", n2l, n2l, BSH, n2h, n2h ); |
#endif |
#endif |
d2h = 0, d2l = 1, s2 = 1; |
d2h = 0, d2l = 1, s2 = 1; |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
if ( n2h == 0 ) goto done; |
if ( n2h == 0 ) goto done; |
BitWidth( n2h, l2 ); |
BitWidth( n2h, l2 ); |
if ( l2 == BSH ) { |
if ( l2 == BSH ) { |
d1h = 0, d1l = 1, s1 = -1; |
d1h = 0, d1l = 1, s1 = -1; |
th = n2h, tl = n2l; |
th = n2h, tl = n2l; |
} else { |
} else { |
i = BSH - l2; |
i = BSH - l2; |
d1h = 0, d1l = 1 << i, s1 = -1; |
d1h = 0, d1l = 1 << i, s1 = -1; |
th = (n2h << i) | (n2l >> l2); |
th = (n2h << i) | (n2l >> l2); |
tl = n2l << i; |
tl = n2l << i; |
} |
} |
if ( tl == 0 ) n1h = - th, n1l = 0; |
if ( tl == 0 ) n1h = - th, n1l = 0; |
else n1h = BMASK - th, n1l = - tl; |
else n1h = BMASK - th, n1l = - tl; |
/**/ |
/**/ |
if ( n1h < n2h || (n1h == n2h && n1l < n2l) ) goto swap12; |
if ( n1h < n2h || (n1h == n2h && n1l < n2l) ) goto swap12; |
BitWidth( n1h, l1 ); |
BitWidth( n1h, l1 ); |
goto sub12; |
goto sub12; |
#else |
#else |
n1h = BASE, n1l = 0, l1 = BSH, |
n1h = BASE, n1l = 0, l1 = BSH, |
d1h = d1l = 0, s1 = 0; |
d1h = d1l = 0, s1 = 0; |
#endif |
#endif |
/**/ |
/**/ |
while ( n2h != 0 ) { |
while ( n2h != 0 ) { |
BitWidth( n2h, l2 ); |
BitWidth( n2h, l2 ); |
sub12: ir = BSH - (i = l1 - l2); |
sub12: ir = BSH - (i = l1 - l2); |
do { |
do { |
if ( i == 0 ) th = n2h, tl = n2l; |
if ( i == 0 ) th = n2h, tl = n2l; |
else |
else |
th = (n2h << i) | (n2l >> ir), |
th = (n2h << i) | (n2l >> ir), |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
tl = n2l << i; |
tl = n2l << i; |
#else |
#else |
tl = (n2l << i) & BMASK; |
tl = (n2l << i) & BMASK; |
#endif |
#endif |
if ( th > n1h || (th == n1h && tl > n1l) ) goto next_i; |
if ( th > n1h || (th == n1h && tl > n1l) ) goto next_i; |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
if ( tl > n1l ) n1h--; |
if ( tl > n1l ) n1h--; |
n1l -= tl; |
n1l -= tl; |
#else |
#else |
if ( tl <= n1l ) n1l -= tl; |
if ( tl <= n1l ) n1l -= tl; |
else n1l += (BASE - tl), n1h--; |
else n1l += (BASE - tl), n1h--; |
#endif |
#endif |
n1h -= th; |
n1h -= th; |
/* (s1:d1h,d1l) -= ((s2:d2h,d2l) << i); */ |
/* (s1:d1h,d1l) -= ((s2:d2h,d2l) << i); */ |
if ( s2 != 0 ) { |
if ( s2 != 0 ) { |
if ( i == 0 ) th = d2h, tl = d2l; |
if ( i == 0 ) th = d2h, tl = d2l; |
else |
else |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
th = (d2h << i) | (d2l >> ir), |
th = (d2h << i) | (d2l >> ir), |
tl = (d2l << i); |
tl = (d2l << i); |
#else |
#else |
th = (d2h << i)&BMASK | (d2l >> ir), |
th = (d2h << i)&BMASK | (d2l >> ir), |
tl = (d2l << i)&BMASK; |
tl = (d2l << i)&BMASK; |
#endif |
#endif |
if ( s1 == 0 ) |
if ( s1 == 0 ) |
s1 = -s2, d1h = th, d1l = tl; |
s1 = -s2, d1h = th, d1l = tl; |
else if ( s1 != s2 ) { |
else if ( s1 != s2 ) { |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
d1l += tl, |
d1l += tl, |
d1h += th; |
d1h += th; |
if ( d1l < tl ) d1h++; |
if ( d1l < tl ) d1h++; |
#else |
#else |
tl += d1l; |
tl += d1l; |
d1l = tl&BMASK; |
d1l = tl&BMASK; |
d1h = (d1h + th + (tl >> BSH))&BMASK; |
d1h = (d1h + th + (tl >> BSH))&BMASK; |
#endif |
#endif |
if ( d1h == 0 && d1l == 0 ) s1 = 0; |
if ( d1h == 0 && d1l == 0 ) s1 = 0; |
} else if ( d1h > th ) { |
} else if ( d1h > th ) { |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
if ( d1l < tl ) d1h--; |
if ( d1l < tl ) d1h--; |
d1l -= tl; |
d1l -= tl; |
#else |
#else |
if ( d1l >= tl ) d1l -= tl; |
if ( d1l >= tl ) d1l -= tl; |
else d1l += (BASE - tl), d1h--; |
else d1l += (BASE - tl), d1h--; |
#endif |
#endif |
d1h -= th; |
d1h -= th; |
} else if ( d1h == th ) { |
} else if ( d1h == th ) { |
d1h = 0; |
d1h = 0; |
if ( d1l == tl ) s1 = d2h = 0; |
if ( d1l == tl ) s1 = d2h = 0; |
else if ( d1l > tl ) d1l -= tl; |
else if ( d1l > tl ) d1l -= tl; |
else d1l = tl - d1l, s1 = -s1; |
else d1l = tl - d1l, s1 = -s1; |
} else { |
} else { |
#if BSH == BIT_WIDTH_OF_INT |
#if BSH == BIT_WIDTH_OF_INT |
if ( tl < d1l ) th--; |
if ( tl < d1l ) th--; |
d1l = tl - d1l; |
d1l = tl - d1l; |
#else |
#else |
if ( tl >= d1l ) d1l = tl - d1l; |
if ( tl >= d1l ) d1l = tl - d1l; |
else d1l = tl + (BASE - d1l), th--; |
else d1l = tl + (BASE - d1l), th--; |
#endif |
#endif |
d1h = th - d1h; |
d1h = th - d1h; |
s1 = -s1; |
s1 = -s1; |
} |
} |
} |
} |
next_i: i--, ir++; |
next_i: i--, ir++; |
} while ( n1h > n2h || (n1h == n2h && n1l >= n2l) ); |
} while ( n1h > n2h || (n1h == n2h && n1l >= n2l) ); |
swap12: /* swap 1 and 2 */ |
swap12: /* swap 1 and 2 */ |
th = n1h, tl = n1l; |
th = n1h, tl = n1l; |
n1h = n2h, n1l = n2l; |
n1h = n2h, n1l = n2l; |
n2h = th, n2l = tl; |
n2h = th, n2l = tl; |
l1 = l2; |
l1 = l2; |
th = d1h, tl = d1l, i = s1; |
th = d1h, tl = d1l, i = s1; |
d1h = d2h, d1l = d2l, s1 = s2; |
d1h = d2h, d1l = d2l, s1 = s2; |
d2h = th, d2l = tl, s2 = i; |
d2h = th, d2l = tl, s2 = i; |
} |
} |
/**/ |
/**/ |
done: *pn = n2l, *pd = d2l; |
done: *pn = n2l, *pd = d2l; |
return s2; |
return s2; |
} |
} |
|
|
static int igcd_spurious_factor; |
static int igcd_spurious_factor; |
|
|
#define SaveN(s,d) {\ |
#define SaveN(s,d) {\ |
int i, l; \ |
int i, l; \ |
for ( l = PL(d) = PL(s), i = 0; i < l; i++ ) BD(d)[i] = BD(s)[i]; \ |
for ( l = PL(d) = PL(s), i = 0; i < l; i++ ) BD(d)[i] = BD(s)[i]; \ |
} |
} |
|
|
static N igcd_acc( n1, n2, nt ) |
static N igcd_acc( n1, n2, nt ) |
N n1, n2, nt; |
N n1, n2, nt; |
/* both n1 and n2 are assumed to be odd */ |
/* both n1 and n2 are assumed to be odd */ |
{ |
{ |
int l1, l2, *b1, *b2, bw1, bw2; |
int l1, l2, *b1, *b2, bw1, bw2; |
int l; |
int l; |
int n, d; |
int n, d; |
N p, s1, s2; |
N p, s1, s2; |
|
|
if ( (l = cmpn( n1, n2 )) == 0 ) return n1; |
if ( (l = cmpn( n1, n2 )) == 0 ) return n1; |
else if ( l < 0 ) { SWAP( n1, n2, N ); } |
else if ( l < 0 ) { SWAP( n1, n2, N ); } |
if ( ShouldCompRemInit(n1,n2) ) { |
if ( ShouldCompRemInit(n1,n2) ) { |
int w, b; |
int w, b; |
|
|
divn( n1, n2, &s1, &s2 ); |
divn( n1, n2, &s1, &s2 ); |
if ( !s2 ) return n2; |
if ( !s2 ) return n2; |
b = trailingzerosn( s2, &w ); |
b = trailingzerosn( s2, &w ); |
p = n1; n1 = n2; n2 = p; |
p = n1; n1 = n2; n2 = p; |
rshiftn( s2, w, b, n2 ); |
rshiftn( s2, w, b, n2 ); |
if ( UNIN(n2) ) return 0; |
if ( UNIN(n2) ) return 0; |
l1 = PL(n1); |
l1 = PL(n1); |
if ( !s1 || PL(s1) < l1 ) s1 = NALLOC(l1); |
if ( !s1 || PL(s1) < l1 ) s1 = NALLOC(l1); |
} else if ( UNIN(n2) ) return 0; |
} else if ( UNIN(n2) ) return 0; |
else { |
else { |
s1 = NALLOC(PL(n1)); |
s1 = NALLOC(PL(n1)); |
s2 = NALLOC(PL(n2)); |
s2 = NALLOC(PL(n2)); |
} |
} |
SaveN( n1, s1 ); |
SaveN( n1, s1 ); |
SaveN( n2, s2 ); |
SaveN( n2, s2 ); |
igcd_spurious_factor = 0; |
igcd_spurious_factor = 0; |
loop: l1 = PL(n1), l2 = PL(n2); |
loop: l1 = PL(n1), l2 = PL(n2); |
if ( l1 <= 2 && l2 <= 2 ) { |
if ( l1 <= 2 && l2 <= 2 ) { |
l = igcd_binary_2w( BD(n1), l1, BD(n2), l2, BD(nt) ); |
l = igcd_binary_2w( BD(n1), l1, BD(n2), l2, BD(nt) ); |
if ( l == 0 ) return 0; |
if ( l == 0 ) return 0; |
PL(nt) = l; |
PL(nt) = l; |
SWAP( n2, nt, N ); |
SWAP( n2, nt, N ); |
goto ret; |
goto ret; |
} |
} |
/**/ |
/**/ |
b1 = BD(n1), b2 = BD(n2); |
b1 = BD(n1), b2 = BD(n2); |
BitWidth( b1[l1 -1], bw1 ); |
BitWidth( b1[l1 -1], bw1 ); |
BitWidth( b2[l2 -1], bw2 ); |
BitWidth( b2[l2 -1], bw2 ); |
if ( (l1*BSH + bw1) - (l2*BSH + bw2) <= igcdacc_thre ) { |
if ( (l1*BSH + bw1) - (l2*BSH + bw2) <= igcdacc_thre ) { |
l = ReducedRatMod( n1, n2, &n, &d ); |
l = ReducedRatMod( n1, n2, &n, &d ); |
l = l < 0 ? aUplusbV_maxrshift( n, b2, l2, d, b1, l1, BD(nt) ) : |
l = l < 0 ? aUplusbV_maxrshift( n, b2, l2, d, b1, l1, BD(nt) ) : |
abs_axU_bxV_maxrshift( n, b2, l2, d, b1, l1, BD(nt) ); |
abs_axU_bxV_maxrshift( n, b2, l2, d, b1, l1, BD(nt) ); |
igcd_spurious_factor++; |
igcd_spurious_factor++; |
if ( l == 0 ) goto ret; |
if ( l == 0 ) goto ret; |
PL(nt) = l; |
PL(nt) = l; |
} else { |
} else { |
l = bmod_n( n1, n2, nt ); |
l = bmod_n( n1, n2, nt ); |
if ( l == 0 ) goto ret; |
if ( l == 0 ) goto ret; |
else if ( l < 0 ) { |
else if ( l < 0 ) { |
SWAP( n1, n2, N ); |
SWAP( n1, n2, N ); |
goto loop; |
goto loop; |
} |
} |
} |
} |
p = n1; |
p = n1; |
n1 = n2; |
n1 = n2; |
n2 = nt; |
n2 = nt; |
nt = p; |
nt = p; |
goto loop; |
goto loop; |
/**/ |
/**/ |
ret: if ( igcd_spurious_factor != 0 && !UNIN(n2) ) { |
ret: if ( igcd_spurious_factor != 0 && !UNIN(n2) ) { |
if ( (p = igcd_bmod( n2, s1, n1 )) == 0 ) return 0; |
if ( (p = igcd_bmod( n2, s1, n1 )) == 0 ) return 0; |
if ( (p = igcd_bmod( p, s2, nt )) == 0 ) return 0; |
if ( (p = igcd_bmod( p, s2, nt )) == 0 ) return 0; |
return p; |
return p; |
} else return n2; |
} else return n2; |
} |
} |
|
|
void gcdaccn( n1, n2, nr ) |
void gcdaccn( n1, n2, nr ) |
N n1, n2, *nr; |
N n1, n2, *nr; |
{ |
{ |
int s1, s2, gw, gb, t1, t2; |
int s1, s2, gw, gb, t1, t2; |
int w1, w2; |
int w1, w2; |
N tn1, tn2, tnt, p; |
N tn1, tn2, tnt, p; |
|
|
if ( !n1 ) { |
if ( !n1 ) { |
*nr = n2; |
*nr = n2; |
return; |
return; |
} else if ( !n2 ) { |
} else if ( !n2 ) { |
*nr = n1; |
*nr = n1; |
return; |
return; |
} |
} |
s1 = trailingzerosn( n1, &w1 ); |
s1 = trailingzerosn( n1, &w1 ); |
s2 = trailingzerosn( n2, &w2 ); |
s2 = trailingzerosn( n2, &w2 ); |
if ( w1 == w2 ) gw = w1, gb = s1 <= s2 ? s1 : s2; |
if ( w1 == w2 ) gw = w1, gb = s1 <= s2 ? s1 : s2; |
else if ( w1 < w2 ) gw = w1, gb = s1; |
else if ( w1 < w2 ) gw = w1, gb = s1; |
else gw = w2, gb = s2; |
else gw = w2, gb = s2; |
/* |
/* |
* true GCD must be multiplied by 2^{gw*BSH+gs}. |
* true GCD must be multiplied by 2^{gw*BSH+gs}. |
*/ |
*/ |
t1 = PL(n1) - w1; |
t1 = PL(n1) - w1; |
t2 = PL(n2) - w2; |
t2 = PL(n2) - w2; |
if ( t1 < t2 ) t1 = t2; |
if ( t1 < t2 ) t1 = t2; |
tn1 = W_NALLOC(t1); tn2 = W_NALLOC(t1); tnt = W_NALLOC(t1); |
tn1 = W_NALLOC(t1); tn2 = W_NALLOC(t1); tnt = W_NALLOC(t1); |
rshiftn( n1, w1, s1, tn1 ); |
rshiftn( n1, w1, s1, tn1 ); |
rshiftn( n2, w2, s2, tn2 ); |
rshiftn( n2, w2, s2, tn2 ); |
/**/ |
/**/ |
p = igcd_acc( tn1, tn2, tnt ); |
p = igcd_acc( tn1, tn2, tnt ); |
/**/ |
/**/ |
if ( p == 0 ) goto L0; |
if ( p == 0 ) goto L0; |
RetTrueGCD( p, gw, gb, nr, L0 ) |
RetTrueGCD( p, gw, gb, nr, L0 ) |
} |
} |
|
|
|
|
|
|
void gcdn_HMEXT( n1, n2, nr ) |
void gcdn_HMEXT( n1, n2, nr ) |
N n1, n2, *nr; |
N n1, n2, *nr; |
{ |
{ |
int b1, b2, w1, w2, gw, gb; |
int b1, b2, w1, w2, gw, gb; |
int l1, l2; |
int l1, l2; |
N tn1, tn2, tnt, a; |
N tn1, tn2, tnt, a; |
|
|
if ( !n1 ) { |
if ( !n1 ) { |
*nr = n2; return; |
*nr = n2; return; |
} else if ( !n2 ) { |
} else if ( !n2 ) { |
*nr = n1; return; |
*nr = n1; return; |
} |
} |
b1 = trailingzerosn( n1, &w1 ); |
b1 = trailingzerosn( n1, &w1 ); |
b2 = trailingzerosn( n2, &w2 ); |
b2 = trailingzerosn( n2, &w2 ); |
if ( w1 == w2 ) gw = w1, gb = b1 <= b2 ? b1 : b2; |
if ( w1 == w2 ) gw = w1, gb = b1 <= b2 ? b1 : b2; |
else if ( w1 < w2 ) gw = w1, gb = b1; |
else if ( w1 < w2 ) gw = w1, gb = b1; |
else gw = w2, gb = b2; |
else gw = w2, gb = b2; |
/* |
/* |
* true GCD must be multiplied by 2^{gw*BSH+gb}. |
* true GCD must be multiplied by 2^{gw*BSH+gb}. |
*/ |
*/ |
l1 = PL(n1) - w1; |
l1 = PL(n1) - w1; |
l2 = PL(n2) - w2; |
l2 = PL(n2) - w2; |
if ( l1 < l2 ) l1 = l2; |
if ( l1 < l2 ) l1 = l2; |
tn1 = W_NALLOC( l1 ); tn2 = W_NALLOC( l1 ); tnt = W_NALLOC( l1 ); |
tn1 = W_NALLOC( l1 ); tn2 = W_NALLOC( l1 ); tnt = W_NALLOC( l1 ); |
rshiftn( n1, w1, b1, tn1 ); |
rshiftn( n1, w1, b1, tn1 ); |
rshiftn( n2, w2, b2, tn2 ); |
rshiftn( n2, w2, b2, tn2 ); |
/**/ |
/**/ |
if ( igcd_algorithm == 1 ) { |
if ( igcd_algorithm == 1 ) { |
a = igcd_binary( tn1, tn2, tnt ); |
a = igcd_binary( tn1, tn2, tnt ); |
} else if ( igcd_algorithm == 2 ) { |
} else if ( igcd_algorithm == 2 ) { |
a = igcd_bmod( tn1, tn2, tnt ); |
a = igcd_bmod( tn1, tn2, tnt ); |
} else { |
} else { |
a = igcd_acc( tn1, tn2, tnt ); |
a = igcd_acc( tn1, tn2, tnt ); |
if ( igcd_spurious_factor != 0 ) { |
if ( igcd_spurious_factor != 0 ) { |
} |
} |
} |
} |
RetTrueGCD( a, gw, gb, nr, L0 ) |
RetTrueGCD( a, gw, gb, nr, L0 ) |
} |
} |
|
|
|
|
Line 1075 N maxrshn( n, p ) |
|
Line 1075 N maxrshn( n, p ) |
|
N n; |
N n; |
int *p; |
int *p; |
{ |
{ |
int nw, nb, c, l; |
int nw, nb, c, l; |
N new; |
N new; |
|
|
nb = trailingzerosn( n, &nw ); |
nb = trailingzerosn( n, &nw ); |
l = PL(n); |
l = PL(n); |
c = BD(n)[l -1]; |
c = BD(n)[l -1]; |
l -= nw; |
l -= nw; |
if ( (c >> nb) == 0 ) l--; |
if ( (c >> nb) == 0 ) l--; |
new = NALLOC(l); |
new = NALLOC(l); |
rshiftn( n, nw, nb, new ); |
rshiftn( n, nw, nb, new ); |
*p = nb + nw*BSH; |
*p = nb + nw*BSH; |
return new; |
return new; |
} |
} |
#endif |
#endif |