mpn/generic/rootrem.c

1 /* mpn_rootrem(rootp,remp,ap,an,nth) -- Compute the nth root of {ap,an}, and
26 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
29      This implementation is not optimal when remp == NULL, since the complexity
30      is M(n), whereas it should be M(n/k) on average.
68 	  or a non-zero value iff the remainder is non-zero when remp = NULL.
70    (a) up[un-1] is not zero
73    (d) the operands do not overlap.
105 	 one limb. (In case sp[0]=1, we can deduce the root, but not decide
106 	 whether it is exact or not.) */
117 /* if approx is non-zero, does not compute the final remainder */
123   mp_size_t qn, rn, sn, wn, nl, bn;
125   unsigned long int unb; /* number of significant bits of {up,un} */
126   unsigned long int xnb; /* number of significant bits of the result */
130   int ni, i;
137   /* qp and wp need enough space to store S'^k where S' is an approximate
139      S'=4. But then since we know the number of bits of S in advance, S'
142      fits in un limbs, the number of extra limbs needed is bounded by
162   /* unb is the number of bits of the input U */
165   /* xnb is the number of bits of the root R */
173 				   if we demand u to be normalized  */
182   kk = k * (xnb - 1);		/* number of truncated bits in the input */
183   rn = un - kk / GMP_NUMB_BITS; /* number of limbs of the non-truncated part */
186 				   the non-truncated part is less than 2^k, it
195   b = xnb - 1; /* number of remaining bits to determine in the kth root */
196   ni = 0;
200       sizes[ni] = b;
201       /* if c is the new value of b, this means that we'll go from a root
205 	 if s' >= k*beta, then at most one correction is necessary.
209       if (b >= sizes[ni])
210 	b = sizes[ni] - 1;	/* add just one bit at a time */
211       ni++;
213   sizes[ni] = 0;
214   ASSERT_ALWAYS (ni < GMP_NUMB_BITS + 1);
215   /* We have sizes[0] = b > sizes[1] > ... > sizes[ni] = 0 with
217      Newton iteration will first compute sizes[ni-1] extra bits,
218      then sizes[ni-2], ..., then sizes[0] = b. */
222   for (i = ni; i != 0; i--)
226 		  exactly 1 + sizes[ni] bits.
229 	 kk = number of truncated bits of the input
231       b = sizes[i - 1] - sizes[i]; /* number of bits to compute in that
234       /* Reinsert a low zero limb if we normalized away the entire remainder */
257       /* nl is the number of limbs in U which contain bits [kk,kk+b-1] */
258       nl = 1 + (kk + b - 1) / GMP_NUMB_BITS - (kk / GMP_NUMB_BITS);
259       /* nl  = 1 + floor((kk + b - 1) / GMP_NUMB_BITS)
264 	 thus since nl is an integer:
265 	 nl <= 2 + floor(b/GMP_NUMB_BITS) <= 2 + bn. */
266       /* we have to save rp[bn] up to rp[nl-1], i.e. 1 or 2 limbs */
267       if (nl - 1 > bn)
269       MPN_RSHIFT (cy, rp, up + kk / GMP_NUMB_BITS, nl, kk % GMP_NUMB_BITS);
274       if (nl - 1 > bn)
288       /* now divide {rp, rn} by {wp, wn} to get the low part of the root */
298 	     The quotient needs rn-wn+1 limbs, thus quotient+remainder
299 	     need altogether rn+1 limbs. */
306 	 Note: {rp,rn} is not needed any more since we'll compute it from
353 	 we now have to take another (k-1)*b bits from the input. */
367 	      /* Last iteration: we don't need W anymore. */
375 	      /* W <- S^(k-1) for the next iteration,