cheng/mpir
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
1c2f0506f5
mpir/mpn/generic/rootrem.c
jasonmoxham eb134da0ea C89 changes
2009-09-03 18:18:35 +00:00

993 lines
30 KiB
C
Raw Blame History

/* mpn_rootrem
Copyright 2009 Jason Moxham
This file is part of the MPIR Library.
The MPIR Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation; either version 2.1 of the License, or (at
your option) any later version.
The MPIR Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the MPIR Library; see the file COPYING.LIB. If not, write
to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
Boston, MA 02110-1301, USA.
*/
// HERE IS A COPY OF THE OLD GMP ROOTREM WHICH WE RENAMED MPN-ROOTREM_BASECASE
// WE USE THIS FOR SMALL SIZES
// AND OF THE COURSE THE OLD GMP BOILERPLATE
/* mpn_rootrem(rootp,remp,ap,an,nth) -- Compute the nth root of {ap,an}, and
store the truncated integer part at rootp and the remainder at remp.
THE FUNCTIONS IN THIS FILE ARE INTERNAL FUNCTIONS WITH MUTABLE
INTERFACES. IT IS ONLY SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.
IN FACT, IT IS ALMOST GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A
FUTURE GNU MP RELEASE.
Copyright 2002, 2005 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.
The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library; see the file COPYING.LIB. If not, write
to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
Boston, MA 02110-1301, USA. */
/*
We use Newton's method to compute the root of a:
n
f(x) := x - a
n - 1
f'(x) := x n
n-1 n-1 n-1
x - a/x a/x - x a/x + (n-1)x
new x = x - f(x)/f'(x) = x - ---------- = x + --------- = --------------
n n n
*/
#include "mpir.h"
#include "gmp-impl.h"
#include "longlong.h"
static mp_size_t
mpn_rootrem_basecase (mp_ptr rootp, mp_ptr remp,mp_srcptr up, mp_size_t un, mp_limb_t nth)
{
mp_ptr pp, qp, xp;
mp_size_t pn, xn, qn;
unsigned long int unb, xnb, bit;
unsigned int cnt;
mp_size_t i;
unsigned long int n_valid_bits, adj;
TMP_DECL;
TMP_MARK;
/* The extra factor 1.585 = log(3)/log(2) here is for the worst case
overestimate of the root, i.e., when the code rounds a root that is
2+epsilon to 3, and then powers this to a potentially huge power. We
could generalize the code for detecting root=1 a few lines below to deal
with xnb <= k, for some small k. For example, when xnb <= 2, meaning
the root should be 1, 2, or 3, we could replace this factor by the much
smaller log(5)/log(4). */
#define PP_ALLOC (2 + (mp_size_t) (un*1.585))
pp = TMP_ALLOC_LIMBS (PP_ALLOC);
count_leading_zeros (cnt, up[un - 1]);
unb = un * GMP_NUMB_BITS - cnt + GMP_NAIL_BITS;
xnb = (unb - 1) / nth + 1;
if (xnb == 1)
{
if (remp == 0)
remp = pp;
mpn_sub_1 (remp, up, un, (mp_limb_t) 1);
MPN_NORMALIZE (remp, un);
rootp[0] = 1;
TMP_FREE;
return un;
}
xn = (xnb + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
qp = TMP_ALLOC_LIMBS (PP_ALLOC);
xp = TMP_ALLOC_LIMBS (xn + 1);
/* Set initial root to only ones. This is an overestimate of the actual root
by less than a factor of 2. */
for (i = 0; i < xn; i++)
xp[i] = GMP_NUMB_MAX;
xp[xnb / GMP_NUMB_BITS] = ((mp_limb_t) 1 << (xnb % GMP_NUMB_BITS)) - 1;
/* Improve the initial approximation, one bit at a time. Keep the
approximations >= root(U,nth). */
bit = xnb - 2;
n_valid_bits = 0;
for (i = 0; (nth >> i) != 0; i++)
{
mp_limb_t xl = xp[bit / GMP_NUMB_BITS];
xp[bit / GMP_NUMB_BITS] = xl ^ (mp_limb_t) 1 << bit % GMP_NUMB_BITS;
pn = mpn_pow_1 (pp, xp, xn, nth, qp);
ASSERT_ALWAYS (pn < PP_ALLOC);
/* If the new root approximation is too small, restore old value. */
if (! (un < pn || (un == pn && mpn_cmp (up, pp, pn) < 0)))
xp[bit / GMP_NUMB_BITS] = xl; /* restore old value */
n_valid_bits += 1;
if (bit == 0)
goto done;
bit--;
}
adj = n_valid_bits - 1;
/* Newton loop. Converges downwards towards root(U,nth). Currently we use
full precision from iteration 1. Clearly, we should use just n_valid_bits
of precision in each step, and thus save most of the computations. */
while (n_valid_bits <= xnb)
{
mp_limb_t cy;
pn = mpn_pow_1 (pp, xp, xn, nth - 1, qp);
ASSERT_ALWAYS (pn < PP_ALLOC);
qp[xn - 1] = 0; /* pad quotient to make it always xn limbs */
mpn_tdiv_qr (qp, pp, (mp_size_t) 0, up, un, pp, pn); /* junk remainder */
cy = mpn_addmul_1 (qp, xp, xn, nth - 1);
if (un - pn == xn)
{
cy += qp[xn];
if (cy == nth)
{
for (i = xn - 1; i >= 0; i--)
qp[i] = GMP_NUMB_MAX;
cy = nth - 1;
}
}
qp[xn] = cy;
qn = xn + (cy != 0);
mpn_divrem_euclidean_qr_1 (xp, qp, qn, nth);
n_valid_bits = n_valid_bits * 2 - adj;
}
/* The computed result might be one unit too large. Adjust as necessary. */
done:
pn = mpn_pow_1 (pp, xp, xn, nth, qp);
ASSERT_ALWAYS (pn < PP_ALLOC);
if (un < pn || (un == pn && mpn_cmp (up, pp, pn) < 0))
{
mpn_decr_u (xp, 1);
pn = mpn_pow_1 (pp, xp, xn, nth, qp);
ASSERT_ALWAYS (pn < PP_ALLOC);
ASSERT_ALWAYS (! (un < pn || (un == pn && mpn_cmp (up, pp, pn) < 0)));
}
if (remp == 0)
remp = pp;
mpn_sub (remp, up, un, pp, pn);
MPN_NORMALIZE (remp, un);
MPN_COPY (rootp, xp, xn);
TMP_FREE;
return un;
}
// HERE IS THE NEW CODE
/*
TODO
For large k we can calulate x^k faster as a float ie exp(k*ln(x)) or x^(1/k)=exp(ln(x)/k)
rather than doing it bitwise , round up all the truncation to the next limb , this should save
quite a lot of shifts , don't know how much this will save (if any) in practice
The powering is now a base2 left to right binary expansion , we could the usual sliding base 2^k
expansion , although the most common roots are small so this is not likely to give us much in the common case
As most roots are for small k , we can do the powering via an optimized addition chain , ie some sort of
table lookup
Merge this reciprocal with our reciprocal used in our barratt (and/or newton division)
Currently we calc x^(1/k) as (x^(-1/k))^(-1/1)
or (x^(-1/1))^(-1/k)
could also try x(x^(-1/k)^(k-1)) (*)
or (x^(-1/a))^(-1/b) where k=ab
this last one is SLOWER as high k is fast as so make out computation as small as poss as fast as poss
So (*) is the only alternative , which I guess is only faster for small k ???
Rewrite in term of mpf (or similar) like it was when I started , but I lost it , will make the code
below much clearer and smaller.
multrunc can use high half mul
if k<496 (32 bit cpus) then nroot_vsmall can be further reduced for a nroot_vvsmall
change signed long etc to mp_size_t ? mainly for MSVC
At the moment we have just one threshold , need a separate one for each k , and some sort of rule for large k
*/
/* Algortihms from "Detecting Perfect Powers in Essentially Linear Time" ,
Daniel J Bernstein http://cr.yp.to/papers.html */
// define this to 1 to test the nroot_small code
#define TESTSMALL 0
// if k<=floor((2^(GMP_LIMB_BITS-1)-33)/66) && k<=2^(GMP_LIMB_BITS-4) then can call vsmall
// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ is always the smallest
#define NROOT_VSMALL_MIN (((((mp_limb_t)1)<<(GMP_LIMB_BITS-1))-33)/66)
// shiftrights requires an extra gmp_numb_bits
#define shiftright(x,xn,c) \
do{(xn)=(xn)-(c)/GMP_NUMB_BITS; \
if((c)%GMP_NUMB_BITS!=0) \
{mpn_rshift((x),(x)+(c)/GMP_NUMB_BITS,(xn),(c)%GMP_NUMB_BITS); \
if((x)[(xn)-1]==0)(xn)--;} \
else{if((c)/GMP_NUMB_BITS!=0)MPN_COPY_INCR((x),(x)+(c)/GMP_NUMB_BITS,(xn));} \
}while(0)
// shiftrights requires an extra gmp_numb_bits
#define shiftrights(x,xn,y,yn,c) \
do{(xn)=(yn)-(c)/GMP_NUMB_BITS; \
if((c)%GMP_NUMB_BITS!=0) \
{mpn_rshift((x),(y)+(c)/GMP_NUMB_BITS,(xn),(c)%GMP_NUMB_BITS); \
if((x)[(xn)-1]==0)(xn)--;} \
else{MPN_COPY_INCR((x),(y)+(c)/GMP_NUMB_BITS,(xn));} \
}while(0)
#define shiftleft(x,xn,c) \
do{mp_limb_t __t; \
if((c)%GMP_NUMB_BITS!=0) \
{__t=mpn_lshift((x)+(c)/GMP_NUMB_BITS,(x),(xn),(c)%GMP_NUMB_BITS); \
(xn)=(xn)+(c)/GMP_NUMB_BITS; \
if(__t!=0){(x)[(xn)]=__t;(xn)++;}} \
else \
{if((c)/GMP_NUMB_BITS!=0) \
{MPN_COPY_DECR((x)+(c)/GMP_NUMB_BITS,(x),(xn)); \
(xn)=(xn)+(c)/GMP_NUMB_BITS;}} \
if((c)/GMP_NUMB_BITS!=0)MPN_ZERO((x),(c)/GMP_NUMB_BITS); \
}while(0)
#define shiftlefts(x,xn,y,yn,c) \
do{mp_limb_t __t; \
if((c)%GMP_NUMB_BITS!=0) \
{__t=mpn_lshift((x)+(c)/GMP_NUMB_BITS,(y),(yn),(c)%GMP_NUMB_BITS); \
(xn)=(yn)+(c)/GMP_NUMB_BITS; \
if(__t!=0){(x)[(xn)]=__t;(xn)++;}} \
else \
{MPN_COPY_DECR((x)+(c)/GMP_NUMB_BITS,(y),(yn)); \
(xn)=(yn)+(c)/GMP_NUMB_BITS;} \
if((c)/GMP_NUMB_BITS!=0)MPN_ZERO((x),(c)/GMP_NUMB_BITS); \
}while(0)
#define mul_ui(x,xn,k) do{mp_limb_t __t;__t=mpn_mul_1((x),(x),(xn),(k));if(__t!=0){(x)[(xn)]=__t;(xn)++;}}while(0)
// tdiv_q_ui requires an extra gmp_numb_bits
#define tdiv_q_ui(x,xn,k) do{mpn_divrem_1((x),0,(x),(xn),(k));if((x)[(xn)-1]==0)(xn)--;}while(0)
// bigmultrunc requires an extra gmp_numb_bits
#define bigmultrunc(xv,xn,xp,yv,yn,yp,B) \
do{signed long __f;mp_limb_t __t; \
(xp)+=(yp); \
if((xn)>=(yn)){__t=mpn_mul(t1,(xv),(xn),(yv),(yn));} \
else{__t=mpn_mul(t1,(yv),(yn),(xv),(xn));} \
t1n=(xn)+(yn);if(__t==0)t1n--; \
__f=sizetwo(t1,t1n);__f=__f-(B); \
if(__f>0) \
{shiftrights((xv),(xn),t1,t1n,__f);(xp)+=__f;} \
else{MPN_COPY_INCR((xv),t1,t1n);(xn)=t1n;} \
}while(0)
// bigsqrtrunc requires an extra gmp_numb_bits
#define bigsqrtrunc(xv,xn,xp,B) \
do{signed long __f; \
(xp)+=(xp); \
mpn_sqr_n(t1,(xv),(xn)); \
t1n=(xn)*2;if(t1[t1n-1]==0)t1n--; \
__f=sizetwo(t1,t1n);__f=__f-(B); \
if(__f>0) \
{shiftrights((xv),(xn),t1,t1n,__f);(xp)+=__f;} \
else{MPN_COPY_INCR((xv),t1,t1n);(xn)=t1n;} \
}while(0)
// must have y>z value wise
#define subtract(x,xn,y,yn,z,zn) \
do{mpn_sub((x),(y),(yn),(z),(zn));/* no carry */ \
(xn)=(yn);while((x)[(xn)-1]==0)(xn)--; \
}while(0)
// returns ceil(lg(x)) where x!=0
signed long
clg (unsigned long x)
{
mp_limb_t t;
ASSERT (x != 0);
#if BITS_PER_ULONG<=GMP_LIMB_BITS
if (x == 1)
return 0;
count_leading_zeros (t, (mp_limb_t) (x - 1));
return GMP_LIMB_BITS - t;
#endif
#if BITS_PER_ULONG>GMP_LIMB_BITS
#error FIXME
#endif
}
// returns sizeinbase(x,2) x!=0
static inline signed long
sizetwo (mp_srcptr x, mp_size_t xn)
{
signed long r;
count_leading_zeros (r, x[xn - 1]);
return xn * GMP_NUMB_BITS + GMP_NAIL_BITS - r;
}
// returns sizeinbase(x-1,2) and returns 0 if x=1
static inline signed long
sizetwom1 (mp_srcptr x, mp_size_t xn)
{
signed long r, i;
mp_limb_t v;
ASSERT (xn > 1 || (xn == 1 && x[0] != 0));
if (xn == 1 && x[0] == 1)
return 0;
r = sizetwo (x, xn);
i = xn - 1;
v = x[i];
if ((v & (v - 1)) != 0)
return r;
for (i--; i >= 0; i--)
if (x[i] != 0)
return r;
return r - 1;
}
/* Algorithm B
Calculates Z such that Z*(1-2^(-b)) < Y^(-1/k) < Z*(1+2^(-b))
ie a b bit approximation the reciprocal of the kth root of Y
where Z,Y>0 are real , 1<=b<=3+ceil(lg(k)) is an int , k>=1 is an int
Z={z,zn}*2^zp where zp is the return value , and zn is modified
Y={y,yn}*2^yp where {y,yn}>=2 and leading limb of {y,yn} is not zero
{z,zn} requires space for GMP_LIMB_BITS+4 bits
{z,zn} and {y,yn} must be completly distinct
*/
static signed long
nroot_small (mp_ptr z, mp_size_t * zn, mp_srcptr y, mp_size_t yn,
signed long yp, mp_limb_t b, mp_limb_t k)
{
signed long zp, f, j, g, B, t2p, t3p, t4p;
int ret;
mp_limb_t mask;
mp_size_t t1n, t2n, t3n, t4n;
mp_limb_t t1[BITS_TO_LIMBS (2 * (GMP_LIMB_BITS + 8) + GMP_NUMB_BITS)];
mp_limb_t t2[BITS_TO_LIMBS (GMP_LIMB_BITS + 8 + GMP_NUMB_BITS)];
mp_limb_t t3[BITS_TO_LIMBS (GMP_LIMB_BITS + 8 + 2 + GMP_NUMB_BITS)];
mp_limb_t t4[BITS_TO_LIMBS (GMP_LIMB_BITS + 8 + GMP_NUMB_BITS)];
ASSERT (k != 0);
ASSERT (b >= 1);
ASSERT (b <= (mp_limb_t) (clg (k) + 3)); // bit counts are maximums , ie can have less
ASSERT (yn > 1 || (yn == 1 && y[0] >= 2));
ASSERT (y[yn - 1] != 0);
g = sizetwom1 (y, yn);
g = g + yp;
g = -g;
if (g >= 0)
{
g = g / k;
}
else
{
g = -((k - 1 - g) / k);
}
B = 66 * (2 * k + 1);
B = clg (B);
ASSERT (B <= GMP_LIMB_BITS + 8);
ASSERT (b + 1 <= GMP_LIMB_BITS + 4);
f = sizetwo (y, yn);
if (f > B)
{
shiftrights (t4, t4n, y, yn, f - B);
t4p = yp + f - B;
}
else
{
MPN_COPY_INCR (t4, y, yn);
t4n = yn;
t4p = yp;
} // t4 has B bits+numb space
*zn = 1;
z[0] = 3;
zp = g - 1; // z has 2 bits
for (j = 1; (unsigned long) j < b; j++)
{
f = sizetwo (z, *zn);
if (f > B)
{
shiftrights (t2, t2n, z, *zn, f - B);
t2p = zp + f - B;
}
else
{
MPN_COPY_INCR (t2, z, *zn);
t2n = *zn;
t2p = zp;
} // t2 has B bits+numb space
if (k != 1)
{
MPN_COPY_INCR (t3, t2, t2n);
t3n = t2n;
t3p = t2p; // t3 has B bits
mask = (((mp_limb_t) 1) << (GMP_LIMB_BITS - 1));
while ((mask & k) == 0)
mask >>= 1;
mask >>= 1;
for (; mask != 0; mask >>= 1)
{
bigsqrtrunc (t2, t2n, t2p, B); // t2 has B bits+numb space , t1 has 2*B bits+numb space
if ((k & mask) != 0)
{
bigmultrunc (t2, t2n, t2p, t3, t3n, t3p, B);
}
}
} // t2 has B bits+numb space , t1 has 2*B bits+numb space
bigmultrunc (t2, t2n, t2p, t4, t4n, t4p, B); // t2 has B bits+numb space , t1 has 2*B bits+numb space
ret = 0;
f = sizetwo (t2, t2n);
if (f - 1 <= 8 - (t2p + 10))
ret = 1;
if (f - 1 >= 10 - (t2p + 10))
ret = 0;
if (f - 1 == 9 - (t2p + 10))
{ // so 512 <= t2.2^(t2p+10) < 1024
if (t2p + 10 >= 0)
{
shiftlefts (t3, t3n, t2, t2n, t2p + 10);
} // t3 has 10 bits
else
{
shiftrights (t3, t3n, t2, t2n, -t2p - 10);
} // t3 has 10 bits+numb space
if (t3n == 1 && t3[0] <= 993)
ret = 1;
}
if (ret != 0)
{
shiftleft (z, *zn, zp - (g - j - 1)); // z has j+2 bits
{
mp_limb_t __t;
__t = mpn_add_1 (z, z, *zn, 1);
if (__t != 0)
{
z[*zn] = 1;
(*zn)++;
}
}
zp = g - j - 1;
continue;
}
f = sizetwom1 (t2, t2n);
if (f + t2p >= 1)
{
shiftleft (z, *zn, zp - (g - j - 1));
mpn_sub_1 (z, z, *zn, 1);
zp = g - j - 1;
}
} // z has j+2 bits
return zp;
} // z has b+1 bits
/* Algorithm B for k<=NROOT_VSMALL_MIN=(((((mp_limb_t)1)<<(GMP_LIMB_BITS-1))-33)/66)
Calculates Z such that Z*(1-2^(-b)) < Y^(-1/k) < Z*(1+2^(-b))
ie a b bit approximation the reciprocal of the kth root of Y
where Z,Y>0 are real , 1<=b<=3+ceil(lg(k)) is an int , k>=1 is an int
Z=z[0]*2^zp where zp is the return value
Y={y,yn}*2^yp where {y,yn}>=2 and leading limb of {y,yn} is not zero
{z,1} and {y,yn} must be completly distinct
Note : the restriction on k allows calculations to be less than limb sized
assumes GMP_LIMB_BITS>=10
*/
static signed long
nroot_vsmall (mp_ptr z, mp_srcptr y, mp_size_t yn, signed long yp,
mp_limb_t b, mp_limb_t k)
{
signed long f1, zp, f, j, g, B, t1p, t2p, t3p;
int ret;
mp_limb_t t1, t2, t3, qh, ql, mask;
ASSERT (k != 0);
ASSERT (b >= 1);
ASSERT (b <= (mp_limb_t) (clg (k) + 3));
ASSERT (yn > 1 || (yn == 1 && y[0] >= 2));
ASSERT (y[yn - 1] != 0);
ASSERT (GMP_LIMB_BITS >= 10);
ASSERT (k <= NROOT_VSMALL_MIN);
g = sizetwom1 (y, yn);
B = 66 * (2 * k + 1);
B = clg (B);
ASSERT (B <= GMP_LIMB_BITS);
ASSERT (b <= GMP_LIMB_BITS - 1);
#if GMP_NAIL_BITS==0
t3p = yp;
t3 = y[yn - 1];
count_leading_zeros (f1, t3);
f = yn * GMP_NUMB_BITS + GMP_NAIL_BITS - f1; //related to g(internally)
f1 = GMP_LIMB_BITS - f1;
if (f1 >= B)
{
t3 >>= f1 - B;
t3p += f - B;
}
else
{
if (yn != 1)
{
t3 = (t3 << (B - f1)) | ((y[yn - 2]) >> (GMP_LIMB_BITS - B + f1));
t3p += f - B;
}
}
#endif
#if GMP_NAIL_BITS!=0
#if GMP_NUMB_BITS*2 < GMP_LIMB_BITS
#error not supported
#endif
f = sizetwo (y, yn);
if (f > B)
{
mp_limb_t t3t[2];
mp_size_t t3n;
t3p = yp + f - B;
shiftrights (t3t, t3n, y, yn, f - B);
t3 = t3t[0];
}
else
{
t3p = yp;
if (f <= GMP_NUMB_BITS)
{
t3 = y[0];
}
else
{
t3 = (y[0] | (y[1] << (GMP_NUMB_BITS)));
}
}
#endif
g = g + yp;
g = -g;
if (g >= 0)
{
g = g / k;
}
else
{
g = -((k - 1 - g) / k);
}
z[0] = 3;
zp = g - 1;
for (j = 1; (unsigned long) j < b; j++)
{
count_leading_zeros (f, z[0]);
f = GMP_LIMB_BITS - f;
if (f > B)
{
t1 = (z[0] >> (f - B));
t1p = zp + f - B;
}
else
{
t1 = z[0];
t1p = zp;
}
if (k != 1)
{
t2 = t1;
t2p = t1p;
mask = (((mp_limb_t) 1) << (GMP_LIMB_BITS - 1));
while ((mask & k) == 0)
mask >>= 1;
mask >>= 1;
for (; mask != 0; mask >>= 1)
{
umul_ppmm (qh, ql, t1, t1);
t1p += t1p;
if (qh == 0)
{ //count_leading_zeros(f,ql);f=GMP_LIMB_BITS-f;f=f-B;if(f>0){t1=(ql>>f);t1p+=f;}else{t1=ql;}
t1 = ql; // be lazy
}
else
{
count_leading_zeros (f, qh);
f = 2 * GMP_LIMB_BITS - f;
f = f - B;
t1p += f; //only need these cases when B>=16
if (f < GMP_LIMB_BITS)
{
t1 = (ql >> f);
t1 |= (qh << (GMP_LIMB_BITS - f));
}
else
{
t1 = (qh >> (f - GMP_LIMB_BITS));
}
}
if ((k & mask) != 0)
{
umul_ppmm (qh, ql, t1, t2);
t1p += t2p;
if (qh == 0)
{ //count_leading_zeros(f,ql);f=GMP_LIMB_BITS-f;f=f-B;if(f>0){t1=(ql>>f);t1p+=f;}else{t1=ql;}
t1 = ql; // be lazy
}
else
{
count_leading_zeros (f, qh);
f = 2 * GMP_LIMB_BITS - f;
f = f - B;
t1p += f;
if (f < GMP_LIMB_BITS)
{
t1 = (ql >> f);
t1 |= (qh << (GMP_LIMB_BITS - f));
}
else
{
t1 = (qh >> (f - GMP_LIMB_BITS));
}
}
}
}
}
umul_ppmm (qh, ql, t1, t3);
t1p += t3p;
if (qh == 0)
{
count_leading_zeros (f, ql);
f = GMP_LIMB_BITS - f;
f = f - B;
if (f > 0)
{
t1 = (ql >> f);
t1p += f;
}
else
{
t1 = ql;
}
// dont be lazy here as it could screw up the compairison below
}
else
{
count_leading_zeros (f, qh);
f = 2 * GMP_LIMB_BITS - f;
f = f - B;
t1p += f;
if (f < GMP_LIMB_BITS)
{
t1 = (ql >> f);
t1 |= (qh << (GMP_LIMB_BITS - f));
}
else
{
t1 = (qh >> (f - GMP_LIMB_BITS));
}
}
ret = 0;
ASSERT (t1 != 0);
count_leading_zeros (f, t1);
f = GMP_LIMB_BITS - f;
if (f - 1 <= 8 - (t1p + 10))
ret = 1;
if (f - 1 >= 10 - (t1p + 10))
ret = 0;
if (f - 1 == 9 - (t1p + 10))
{ // so 512 <= t1.2^(t1p+10) < 1024
if (t1p + 10 >= 0)
{
t2 = (t1 << (t1p + 10));
}
else
{
t2 = (t1 >> (-t1p - 10));
}
if (t2 <= 993)
ret = 1;
}
if (ret != 0)
{
z[0] = (z[0] << (zp - (g - j - 1)));
z[0]++;
zp = g - j - 1;
continue;
}
if (t1 == 1)
{
f = 0;
}
else
{
count_leading_zeros (f, t1 - 1);
f = GMP_LIMB_BITS - f;
}
if (f + t1p >= 1)
{
z[0] = (z[0] << (zp - (g - j - 1)));
z[0]--;
zp = g - j - 1;
}
}
return zp;
}
/* Algorithm N
Calculates Z such that Z*(1-2^(-b)) < Y^(-1/k) < Z*(1+2^(-b))
ie a b bit approximation the reciprocal of the kth root of Y
where Z,Y>0 are real , b>=1 is an int , k>=1 is an int
Z={z,zn}*2^zp where zp is the return value , and zn is modified
Y={y,yn}*2^yp where {y,yn}>=2 and leading limb of {y,yn} is not zero
z satisfies 1 <= z < 2^(b+7)
zp satisfies -lg(Y)/k-b-7-lg(3/2) < zp < -lg(Y)/k+1
{z,zn} and {y,yn} and temps t1,t2,t3 must be completly distinct
z requires b+6+GMP_NUMB_BITS+max(1,clgk)
t1 requires max( 2*b+12+GMP_NUMB_BITS , b+6+clg(k+1) )
t2 requires b+6+GMP_NUMB_BITS
t3 requires b+6+GMP_NUMB_BITS
*/
static signed long
nroot (mp_ptr z, mp_size_t * zn, mp_srcptr y, mp_size_t yn, signed long yp,
mp_limb_t b, mp_limb_t k, signed long clgk, mp_ptr t1, mp_ptr t2,
mp_ptr t3)
{ mp_size_t t1n, t2n, t3n; mp_limb_t mask, kpow2, k1pow2;signed long t1p, zp, t2p, t3p, f, bd, bs[GMP_LIMB_BITS * 2], c; // FIXME how many
ASSERT (k != 0);ASSERT (yn > 1 || (yn == 1 && y[0] >= 2));ASSERT (y[yn - 1] != 0);
bs[0] = b; // bit counts are maximums , ie can have less
for (c = 0;; c++)
{ if (bs[c] <= 3 + clgk)
break;
bs[c + 1] = 1 + (bs[c] + clgk) / 2; } // so bs[c]<=3+clgk
#if GMP_LIMB_BITS>=10 && TESTSMALL==0
if (k <= NROOT_VSMALL_MIN)
{ zp = nroot_vsmall (z, y, yn, yp, bs[c], k); *zn = 1; }
else
{ zp = nroot_small (z, (mp_limb_t *) zn, y, yn, yp, bs[c], k); }
#endif
#if GMP_LIMB_BITS<10 || TESTSMALL==1
zp = nroot_small (z, zn, y, yn, yp, bs[c], k);
#endif // bs[1]=1+floor((b+clgk)/2) max bd=b+6 // z has bs[c]+1 bits
kpow2 = 0;k1pow2 = 0; // shortcut for div,mul to a shift instead
if (POW2_P(k)){count_leading_zeros (kpow2, k);kpow2 = GMP_LIMB_BITS - kpow2;} // k=2^(kpow2-1)
if (POW2_P(k+1)){count_leading_zeros (k1pow2, k + 1);k1pow2 = GMP_LIMB_BITS - k1pow2;} // k+1=2^(k1pow2-1)
for (; c != 0; c--)
{ bd = 2 * bs[c] + 4 - clgk;
f = sizetwo (z, *zn); // is this trunc ever going to do something real?
if (f > bd){ shiftright (z, *zn, f - bd); zp = zp + f - bd;} // z has bd bits + numb space
MPN_COPY_INCR (t3, z, *zn); t3n = *zn;t3p = zp;
mask = (((mp_limb_t) 1) << (GMP_LIMB_BITS - 1)); // t3 has bd bits
while ((mask & (k + 1)) == 0)mask >>= 1;
for (mask >>= 1; mask != 0; mask >>= 1)
{ bigsqrtrunc (t3, t3n, t3p, bd); // t3 has bd bits + numb space t1 has 2*bd bits + numb space
if (((k + 1) & mask) != 0){bigmultrunc (t3, t3n, t3p, z, *zn, zp, bd);}}// t3 has bd bits + numb space t1 has 2*bd bits + numb space
if (k1pow2){shiftleft (z, *zn, k1pow2 - 1);}else{mul_ui (z, *zn, k + 1);} // z has bd+clg(k+1) bits
f = sizetwo (y, yn);
if (f > bd){ shiftrights (t2, t2n, y, yn, f - bd);t2p = yp + f - bd;} // t2 has bd bits + numb space
else{ MPN_COPY_INCR (t2, y, yn);t2n = yn;t2p = yp;} // this case may not happen if this is only called by mpn_root
bigmultrunc (t3, t3n, t3p, t2, t2n, t2p, bd); // t3 has bd bits + numb space t1 has 2*bd bits + numb space
if (zp <= t3p) // which branch depends on yp ????? and only want the top bd+clgk bits exactly
{ shiftlefts (t1, t1n, t3, t3n, t3p - zp); // t1 has bd+clg(k+1) bits
subtract (t1, t1n, z, *zn, t1, t1n);
t1p = zp;} // t1 has bd+clg(k+1) bits
else
{ ASSERT(zp - t3p + sizetwo (z, *zn) <= 2 * b + 12 + GMP_NUMB_BITS);// not allocated enough mem
shiftlefts (t1, t1n, z, *zn, zp - t3p); // t1 has 2*b+12+numb
subtract (t1, t1n, t1, t1n, t3, t3n);
t1p = t3p;} // t1 has 2*b+12+numb
f = sizetwo (t1, t1n);
if (f >= bd + clgk){shiftrights (z, *zn, t1, t1n, f - bd - clgk);}
else{shiftlefts (z, *zn, t1, t1n, bd + clgk - f);} // z has bd+clgk bits + numb space
zp = t1p + f - bd - clgk;
if (kpow2){shiftright (z, *zn, kpow2 - 1);}else{tdiv_q_ui (z, *zn, k);}
} // z has bd+1 bits + numb space (maybe prove just bd bits ?)
return zp;
} // z has b+7 bits
/* same as Algorithm N but for k=1
Calculates Z such that Z*(1-2^(-b)) < Y^(-1/k) < Z*(1+2^(-b))
ie a b bit approximation the reciprocal of the kth root of Y
where Z,Y>0 are real , b>=1 is an int , k>=1 is an int
Z={z,zn}*2^zp where zp is the return value , and zn is modified
Y={y,yn}*2^yp where {y,yn}>=2 and leading limb of {y,yn} is not zero
and z satisfies 2^b <= z <= 2^(b+1)
and zp satisfies zp=-sizetwo(y,yn)-b-yp
{z,zn} and {y,yn} and temps t1,t2 must be completly distinct
z requires 2+floor(((sizetwo(y,yn)+b+1)/GMP_NUMB_BITS)-yn limbs
t1 requires 1+floor((sizetwo(y,yn)+b+1)/GMP_NUMB_BITS) limbs
t2 requires yn limbs
*/
static signed long
finv_fast (mp_ptr z, int *zn, mp_srcptr y, mp_size_t yn, signed long yp,
unsigned long b, mp_ptr t1, mp_ptr t2)
{
signed long c;
signed long zp;
mp_size_t t1n;
c = sizetwo (y, yn) + b;
MPN_COPY_INCR (t1, y, yn);
t1n = yn; // t1 has yn limbs
MPN_ZERO(t1 + t1n, (c + 1) / GMP_NUMB_BITS + 1 - t1n); // t1 has 1+floor((c+1)/numb) limbs
t1[(c + 1) / GMP_NUMB_BITS] = (((mp_limb_t) 1) << ((c + 1) % GMP_NUMB_BITS)); // t1 has 1+floor((c+1)/numb) limbs
t1n = (c + 1) / GMP_NUMB_BITS + 1;
ASSERT (y[yn - 1] != 0);
mpn_tdiv_qr (z, t2, 0, t1, t1n, y, yn); //bdivmod could be faster // z has 2+floor((c+1)/numb)-yn t2 has yn limbs
*zn = t1n - yn + 1;
while (*zn != 0 && z[*zn - 1] == 0)
(*zn)--;
shiftright (z, *zn, 1);
zp = -c - yp;
return zp;
}
/* calculates X and R such that X^k<=Y and (X+1)^k>Y
where X={x,xn} Y={y,yn} R={r,rn} , only calculates R if r!=NULL
R satisfies R < (X+1)^k-X^k
X satisfies X^k <= Y
X needs ceil(yn/k) limb space
R needs yn limb space if r!=0
return sizeof remainder if r!=0
*/
mp_size_t mpn_rootrem(mp_ptr xp, mp_ptr r, mp_srcptr y,mp_size_t yn, mp_limb_t k)
{
unsigned long b, clgk;
signed long d, tp, zp;
mpz_t t4, t3;
mp_ptr x,t1,t2;
mp_size_t t2n,xn,rn;
mp_limb_t val;mp_size_t pos, bit;
if(BELOW_THRESHOLD(yn,ROOTREM_THRESHOLD))return mpn_rootrem_basecase(xp,r,y,yn,k);
d = 8; // any d>=1 will do , for testing to its limits use d=1 TUNEME
b = sizetwo (y, yn);
b = (b + k - 1) / k + 2 + d;
clgk = clg (k);
x=__GMP_ALLOCATE_FUNC_LIMBS(BITS_TO_LIMBS(b+7+GMP_NUMB_BITS));
t1=__GMP_ALLOCATE_FUNC_LIMBS(BITS_TO_LIMBS (2 * b + 12 + GMP_NUMB_BITS));
t2=__GMP_ALLOCATE_FUNC_LIMBS(BITS_TO_LIMBS (b + 6 + clgk + 1 + GMP_NUMB_BITS));
mpz_init2 (t3, b + 6 + GMP_NUMB_BITS * 2);
mpz_init2 (t4, b + 6 + GMP_NUMB_BITS);
zp = nroot (t2, &t2n, y, yn, 0, b, k, clgk, t1, PTR (t3), PTR (t4));
/* 1 <= t2 < 2^(b+7) -lg(Y)/k-b-7-lg(3/2) < zp < -lg(Y)/k+1 where Y={y,yn} */
tp = finv_fast (PTR (t3), &SIZ (t3), t2, t2n, zp, b, t1, PTR (t4)); // t3 is our approx root
/* 2^b <= t3 <= 2^(b+1) tp=-sizetwo(t2,t2n)-b-zp */
ASSERT (tp <= -d - 1);
pos = (-tp - d - 1 + 1) / GMP_NUMB_BITS;
bit = (-tp - d - 1 + 1) % GMP_NUMB_BITS;
val = (((mp_limb_t) 1) << bit);
mpn_sub_1 (PTR (t3) + pos, PTR (t3) + pos, SIZ (t3) - pos, val);
if (PTR (t3)[SIZ (t3) - 1] == 0)
SIZ (t3)--;
shiftrights (PTR (t4), SIZ (t4), PTR (t3), SIZ (t3), -tp);
if (mpn_add_1 (PTR (t3) + pos, PTR (t3) + pos, SIZ (t3) - pos, val))
{
PTR (t3)[SIZ (t3)] = 1;
SIZ (t3)++;
}
pos = (-tp - d - 1) / GMP_NUMB_BITS;
bit = (-tp - d - 1) % GMP_NUMB_BITS;
val = (((mp_limb_t) 1) << bit);
if (mpn_add_1 (PTR (t3) + pos, PTR (t3) + pos, SIZ (t3) - pos, val))
{
PTR (t3)[SIZ (t3)] = 1;
SIZ (t3)++;
}
shiftright (PTR (t3), SIZ (t3), -tp);
if (mpz_cmp (t4, t3) == 0)
{
xn = SIZ (t3);
MPN_COPY_INCR (x, PTR (t3), xn);
if (r != 0)
{
mpz_pow_ui (t4, t3, k);
mpn_sub (r, y, yn, PTR (t4), SIZ (t4)); /* no carry */
rn = yn;
while (rn != 0 && r[rn - 1] == 0)rn--;
}
mpz_clear (t4);
mpz_clear (t3);
MPN_COPY(xp,x,(yn+k-1)/k);
__GMP_FREE_FUNC_LIMBS(x,BITS_TO_LIMBS(b+7+GMP_NUMB_BITS));
__GMP_FREE_FUNC_LIMBS(t1,BITS_TO_LIMBS(2*b+12+GMP_NUMB_BITS));
__GMP_FREE_FUNC_LIMBS(t2,BITS_TO_LIMBS(b+6+clgk+1+GMP_NUMB_BITS));
return rn;
}
mpz_pow_ui (t4, t3, k);
if (SIZ (t4) > yn || (SIZ (t4) == yn && mpn_cmp (PTR (t4), y, yn) > 0))
{
mpz_sub_ui (t3, t3, 1);
xn = SIZ (t3);
MPN_COPY_INCR (x, PTR (t3), xn);
if (r != 0)
{
mpz_pow_ui (t4, t3, k);
mpn_sub (r, y, yn, PTR (t4), SIZ (t4)); /* no carry */
rn = yn;
while (rn != 0 && r[rn - 1] == 0) rn--;
}
mpz_clear (t4);
mpz_clear (t3);
MPN_COPY(xp,x,(yn+k-1)/k);
__GMP_FREE_FUNC_LIMBS(x,BITS_TO_LIMBS(b+7+GMP_NUMB_BITS));
__GMP_FREE_FUNC_LIMBS(t1,BITS_TO_LIMBS(2*b+12+GMP_NUMB_BITS));
__GMP_FREE_FUNC_LIMBS(t2,BITS_TO_LIMBS(b+6+clgk+1+GMP_NUMB_BITS));
return rn;
}
xn = SIZ (t3);
MPN_COPY_INCR (x, PTR (t3), xn);
if (r != 0)
{
mpn_sub (r, y, yn, PTR (t4), SIZ (t4)); /* no carry */
rn = yn;
while (rn != 0 && r[rn - 1] == 0) rn--;
}
mpz_clear (t4);
mpz_clear (t3);
MPN_COPY(xp,x,(yn+k-1)/k);
__GMP_FREE_FUNC_LIMBS(x,BITS_TO_LIMBS(b+7+GMP_NUMB_BITS));
__GMP_FREE_FUNC_LIMBS(t1,BITS_TO_LIMBS(2*b+12+GMP_NUMB_BITS));
__GMP_FREE_FUNC_LIMBS(t2,BITS_TO_LIMBS(b+6+clgk+1+GMP_NUMB_BITS));
return rn;}
Reference in New Issue View Git Blame Copy Permalink