__x2 = (uintmax_t) __uh * __vl; \
__x3 = (uintmax_t) __uh * __vh; \
\
- __x1 += __ll_highpart (__x0);/* this can't give carry */ \
- __x1 += __x2; /* but this indeed can */ \
- if (__x1 < __x2) /* did we get it? */ \
- __x3 += __ll_B; /* yes, add it in the proper pos. */ \
+ __x1 += __ll_highpart (__x0);/* This can't give carry. */ \
+ __x1 += __x2; /* But this indeed can. */ \
+ if (__x1 < __x2) /* Did we get it? */ \
+ __x3 += __ll_B; /* Yes, add it in the proper pos. */ \
\
(w1) = __x3 + __ll_highpart (__x1); \
(w0) = (__x1 << W_TYPE_SIZE / 2) + __ll_lowpart (__x0); \
#endif
#if !defined udiv_qrnnd || defined UDIV_NEEDS_NORMALIZATION
-/* Define our own, not needing normalization. This function is
- currently not performance critical, so keep it simple. Similar to
- the mod macro below. */
+/* Define our own, not needing normalization. This function is
+ currently not performance critical, so keep it simple. Similar to
+ the mod macro below. */
# undef udiv_qrnnd
# define udiv_qrnnd(q, r, n1, n0, d) \
do { \
static uintmax_t _GL_ATTRIBUTE_CONST
gcd_odd (uintmax_t a, uintmax_t b)
{
- if ( (b & 1) == 0)
+ if ((b & 1) == 0)
{
uintmax_t t = b;
b = a;
{
assert (b0 & 1);
- if ( (a0 | a1) == 0)
+ if ((a0 | a1) == 0)
{
*r1 = b1;
return b0;
/* Prove primality or run probabilistic tests. */
static bool flag_prove_primality = PROVE_PRIMALITY;
-/* Number of Miller-Rabin tests to run when not proving primality. */
+/* Number of Miller-Rabin tests to run when not proving primality. */
#define MR_REPS 25
static void
/* Trial division with odd primes uses the following trick.
- Let p be an odd prime, and B = 2^{W_TYPE_SIZE}. For simplicity,
+ Let p be an odd prime, and B = 2^{W_TYPE_SIZE}. For simplicity,
consider the case t < B (this is the second loop below).
From our tables we get
binv = p^{-1} (mod B)
- lim = floor ( (B-1) / p ).
+ lim = floor ((B-1) / p).
- First assume that t is a multiple of p, t = q * p. Then 0 <= q <= lim
+ First assume that t is a multiple of p, t = q * p. Then 0 <= q <= lim
(and all quotients in this range occur for some t).
Then t = q * p is true also (mod B), and p is invertible we get
q = t * binv (mod B).
- Next, assume that t is *not* divisible by p. Since multiplication
+ Next, assume that t is *not* divisible by p. Since multiplication
by binv (mod B) is a one-to-one mapping,
t * binv (mod B) > lim,
_q0 = (u0) * _di; \
if ((u1) >= (d)) \
{ \
- uintmax_t _p1, _p0 _GL_UNUSED; \
+ uintmax_t _p1, _p0 _GL_UNUSED; \
umul_ppmm (_p1, _p0, _q0, d); \
(q1) = ((u1) - _p1) * _di; \
(q0) = _q0; \
} \
} while (0)
-/* x B (mod n). */
+/* x B (mod n). */
#define redcify(r_prim, r, n) \
do { \
- uintmax_t _redcify_q _GL_UNUSED; \
+ uintmax_t _redcify_q _GL_UNUSED; \
udiv_qrnnd (_redcify_q, r_prim, r, 0, n); \
} while (0)
-/* x B^2 (mod n). Requires x > 0, n1 < B/2 */
+/* x B^2 (mod n). Requires x > 0, n1 < B/2. */
#define redcify2(r1, r0, x, n1, n0) \
do { \
uintmax_t _r1, _r0, _i; \
} while (0)
/* Modular two-word multiplication, r = a * b mod m, with mi = m^(-1) mod B.
- Both a and b must be in redc form, the result will be in redc form too. */
+ Both a and b must be in redc form, the result will be in redc form too. */
static inline uintmax_t
mulredc (uintmax_t a, uintmax_t b, uintmax_t m, uintmax_t mi)
{
/* Modular two-word multiplication, r = a * b mod m, with mi = m^(-1) mod B.
Both a and b must be in redc form, the result will be in redc form too.
- For performance reasons, the most significant bit of m must be clear. */
+ For performance reasons, the most significant bit of m must be clear. */
static uintmax_t
mulredc2 (uintmax_t *r1p,
uintmax_t a1, uintmax_t a0, uintmax_t b1, uintmax_t b0,
{
uintmax_t r1, r0, q, p1, p0 _GL_UNUSED, t1, t0, s1, s0;
mi = -mi;
- assert ( (a1 >> (W_TYPE_SIZE - 1)) == 0);
- assert ( (b1 >> (W_TYPE_SIZE - 1)) == 0);
- assert ( (m1 >> (W_TYPE_SIZE - 1)) == 0);
+ assert ((a1 >> (W_TYPE_SIZE - 1)) == 0);
+ assert ((b1 >> (W_TYPE_SIZE - 1)) == 0);
+ assert ((m1 >> (W_TYPE_SIZE - 1)) == 0);
/* First compute a0 * <b1, b0> B^{-1}
+-----+
{
uintmax_t y = powm (b, q, n, ni, one);
- uintmax_t nm1 = n - one; /* -1, but in redc representation. */
+ uintmax_t nm1 = n - one; /* -1, but in redc representation. */
if (y == one || y == nm1)
return true;
if (n <= 1)
return false;
- /* We have already casted out small primes. */
+ /* We have already casted out small primes. */
if (n < (uintmax_t) FIRST_OMITTED_PRIME * FIRST_OMITTED_PRIME)
return true;
}
else
{
- /* After enough Miller-Rabin runs, be content. */
+ /* After enough Miller-Rabin runs, be content. */
is_prime = (r == MR_REPS - 1);
}
redcify2 (one[1], one[0], 1, n1, n0);
addmod2 (a_prim[1], a_prim[0], one[1], one[0], one[1], one[0], n1, n0);
- /* FIXME: Use scalars or pointers in arguments? Some consistency needed. */
+ /* FIXME: Use scalars or pointers in arguments? Some consistency needed. */
na[0] = n0;
na[1] = n1;
}
for (unsigned int i = 0; i < factors.nfactors && is_prime; i++)
{
- /* FIXME: We always have the factor 2. Do we really need to
- handle it here? We have done the same powering as part
- of millerrabin. */
+ /* FIXME: We always have the factor 2. Do we really need to
+ handle it here? We have done the same powering as part
+ of millerrabin. */
if (factors.p[i] == 2)
rsh2 (e[1], e[0], nm1[1], nm1[0], 1);
else
}
else
{
- /* After enough Miller-Rabin runs, be content. */
+ /* After enough Miller-Rabin runs, be content. */
is_prime = (r == MR_REPS - 1);
}
if (mpz_cmp_ui (n, 1) <= 0)
return false;
- /* We have already casted out small primes. */
+ /* We have already casted out small primes. */
if (mpz_cmp_ui (n, (long) FIRST_OMITTED_PRIME * FIRST_OMITTED_PRIME) < 0)
return true;
}
else
{
- /* After enough Miller-Rabin runs, be content. */
+ /* After enough Miller-Rabin runs, be content. */
is_prime = (r == MR_REPS - 1);
}
{
assert (a < n);
- binv (ni, n); /* FIXME: when could we use old 'ni' value? */
+ binv (ni, n); /* FIXME: when could we use old 'ni' value? */
for (;;)
{
if (g1 == 0)
{
- /* The found factor is one word, and > 1. */
+ /* The found factor is one word, and > 1. */
divexact_21 (n1, n0, n1, n0, g0); /* n = n / g */
if (!prime_p (g0))
return;
}
- binv (ginv, g0); /* Compute n = n / g. Since the result will */
- n0 = ginv * n0; /* fit one word, we can compute the quotient */
- n1 = 0; /* modulo B, ignoring the high divisor word. */
+ /* Compute n = n / g. Since the result will fit one word,
+ we can compute the quotient modulo B, ignoring the high
+ divisor word. */
+ binv (ginv, g0);
+ n0 = ginv * n0;
+ n1 = 0;
if (!prime2_p (g1, g0))
factor_using_pollard_rho2 (g1, g0, a + 1, factors);
#if USE_SQUFOF
/* FIXME: Maybe better to use an iteration converging to 1/sqrt(n)? If
- algorithm is replaced, consider also returning the remainder. */
+ algorithm is replaced, consider also returning the remainder. */
static uintmax_t _GL_ATTRIBUTE_CONST
isqrt (uintmax_t n)
{
count_leading_zeros (c, n);
- /* Make x > sqrt(n). This will be invariant through the loop. */
+ /* Make x > sqrt(n). This will be invariant through the loop. */
x = (uintmax_t) 1 << ((W_TYPE_SIZE + 1 - c) / 2);
for (;;)
unsigned int shift;
uintmax_t x;
- /* Ensures the remainder fits in an uintmax_t. */
+ /* Ensures the remainder fits in an uintmax_t. */
assert (nh < ((uintmax_t) 1 << (W_TYPE_SIZE - 2)));
if (nh == 0)
count_leading_zeros (shift, nh);
shift &= ~1;
- /* Make x > sqrt(n) */
- x = isqrt ( (nh << shift) + (nl >> (W_TYPE_SIZE - shift))) + 1;
+ /* Make x > sqrt (n). */
+ x = isqrt ((nh << shift) + (nl >> (W_TYPE_SIZE - shift))) + 1;
x <<= (W_TYPE_SIZE - shift) / 2;
- /* Do we need more than one iteration? */
+ /* Do we need more than one iteration? */
for (;;)
{
uintmax_t r _GL_UNUSED;
}
}
-/* MAGIC[N] has a bit i set iff i is a quadratic residue mod N. */
+/* MAGIC[N] has a bit i set iff i is a quadratic residue mod N. */
# define MAGIC64 0x0202021202030213ULL
# define MAGIC63 0x0402483012450293ULL
# define MAGIC65 0x218a019866014613ULL
# define MAGIC11 0x23b
-/* Return the square root if the input is a square, otherwise 0. */
+/* Return the square root if the input is a square, otherwise 0. */
static uintmax_t _GL_ATTRIBUTE_CONST
is_square (uintmax_t x)
{
- /* Uses the tests suggested by Cohen. Excludes 99% of the non-squares before
- computing the square root. */
+ /* Uses the tests suggested by Cohen. Excludes 99% of the non-squares before
+ computing the square root. */
if (((MAGIC64 >> (x & 63)) & 1)
&& ((MAGIC63 >> (x % 63)) & 1)
- /* Both 0 and 64 are squares mod (65) */
+ /* Both 0 and 64 are squares mod (65). */
&& ((MAGIC65 >> ((x % 65) & 63)) & 1)
&& ((MAGIC11 >> (x % 11) & 1)))
{
return 0;
}
-/* invtab[i] = floor(0x10000 / (0x100 + i) */
+/* invtab[i] = floor (0x10000 / (0x100 + i) */
static const unsigned short invtab[0x81] =
{
0x200,
_mask = -(uintmax_t) (_r >= (d)); \
(r) = _r - (_mask & (d)); \
(q) = _q - _mask; \
- assert ( (q) * (d) + (r) == u); \
+ assert ((q) * (d) + (r) == u); \
} \
else \
{ \
} \
} while (0)
-/* Notes: Example N = 22117019. After first phase we find Q1 = 6314, Q
+/* Notes: Example N = 22117019. After first phase we find Q1 = 6314, Q
= 3025, P = 1737, representing F_{18} = (-6314, 2* 1737, 3025),
with 3025 = 55^2.
#if USE_SQUFOF
/* Return true on success. Expected to fail only for numbers
- >= 2^{2*W_TYPE_SIZE - 2}, or close to that limit. */
+ >= 2^{2*W_TYPE_SIZE - 2}, or close to that limit. */
static bool
factor_using_squfof (uintmax_t n1, uintmax_t n0, struct factors *factors)
{
assert (mu * n0 % 4 == 3);
/* In the notation of the paper, with mu * n == 3 (mod 4), we
- get \Delta = 4 mu * n, and the paper's \mu is 2 mu. As far as
+ get \Delta = 4 mu * n, and the paper's \mu is 2 mu. As far as
I understand it, the necessary bound is 4 \mu^3 < n, or 32
mu^3 < n.
105, we get a trivial factor, from the square 38809 = 197^2,
without any corresponding Q earlier in the iteration.
- Requiring 64 mu^3 < n seems sufficient. */
+ Requiring 64 mu^3 < n seems sufficient. */
if (n1 == 0)
{
if ((uintmax_t) mu*mu*mu >= n0 / 64)
Q1 = 1;
P = S;
- /* Square root remainder fits in one word, so ignore high part. */
+ /* Square root remainder fits in one word, so ignore high part. */
Q = Dl - P*P;
- /* FIXME: When can this differ from floor(sqrt(2 sqrt(D)))? */
+ /* FIXME: When can this differ from floor (sqrt (2 * sqrt (D)))? */
L = isqrt (2*S);
B = 2*L;
L1 = mu * 2 * L;
/* The form is (+/- Q1, 2P, -/+ Q), of discriminant 4 (P^2 + Q Q1) =
- 4 D. */
+ 4 D. */
for (i = 0; i <= B; i++)
{
{
uintmax_t g = Q;
- if ( (Q & 1) == 0)
+ if ((Q & 1) == 0)
g /= 2;
g /= gcd_odd (g, mu);
}
/* I think the difference can be either sign, but mod
- 2^W_TYPE_SIZE arithmetic should be fine. */
+ 2^W_TYPE_SIZE arithmetic should be fine. */
t = Q1 + q * (P - P1);
Q1 = Q;
Q = t;
P = P1;
- if ( (i & 1) == 0)
+ if ((i & 1) == 0)
{
uintmax_t r = is_square (Q);
if (r)
if (queue[j].Q == r)
{
if (r == 1)
- /* Traversed entire cycle. */
+ /* Traversed entire cycle. */
goto next_multiplier;
- /* Need the absolute value for divisibility test. */
+ /* Need the absolute value for divisibility test. */
if (P >= queue[j].P)
t = P - queue[j].P;
else
if (t % r == 0)
{
/* Delete entries up to and including entry
- j, which matched. */
+ j, which matched. */
memmove (queue, queue + j + 1,
(qpos - j - 1) * sizeof (queue[0]));
qpos -= (j + 1);
}
/* We have found a square form, which should give a
- factor. */
+ factor. */
Q1 = r;
- assert (S >= P); /* What signs are possible? */
+ assert (S >= P); /* What signs are possible? */
P += r * ((S - P) / r);
/* Note: Paper says (N - P*P) / Q1, that seems incorrect
- for the case D = 2N. */
+ for the case D = 2N. */
/* Compute Q = (D - P*P) / Q1, but we need double
- precision. */
+ precision. */
uintmax_t hi, lo;
umul_ppmm (hi, lo, P, P);
sub_ddmmss (hi, lo, Dh, Dl, hi, lo);
Step 4a in the algorithm description says q <--
floor([S+P]/\hat Q), but looking at the equations
in Sec. 3.1, it should be q <-- floor([S+P] / Q).
- (In this code, \hat Q is Q1). */
+ (In this code, \hat Q is Q1). */
div_smallq (q, rem, S+P, Q);
P1 = S - rem; /* P1 = q*Q - P */
P = P1;
}
- if ( (Q & 1) == 0)
+ if ((Q & 1) == 0)
Q /= 2;
Q /= gcd_odd (Q, mu);
projects / thirdparty / coreutils.git / commitdiff
