mpf_root?

[Paul Leyland]

I am clearly missing something, which may not be too surprising as I now
have a nasty attack of flu with all the unpleasant symptoms that
implies.

I would use double precision floating point to get an initial estimate
thus:

double dc, dp = p, dq = q, dn = n;
dc = exp ((log (dq) - log (dp)) / dn) - 1.0;

followed by a single iteration of Newton-Raphson in integer arithmetic
to polish the result.

(Hope the C above is ok.  I took the liberty of inverting p/q to remove
the need for the reciprocal.  I further assumed that p and q are both >=
1. If this is not the case, you run into all sorts of nasty special
cases.


Paul


On Thu, 2007-12-27 at 23:26, hv at crypt.org wrote:
> Paul Zimmermann <Paul.Zimmermann at loria.fr> wrote:
> :> Hi, I'm trying to convert some code from PARI/GP to C, using gmp-4.2.2.
> :> The code requires only (big) integer arithmetic, except for one calculation:
> :>   c = floor(1 / (sqrtn(p / q, n) - 1))
> :> .. in which q < p < 2^128 and 3 <= n <= 8.
> :> 
> :> Now I have a fair amount of experience in integer calculations, but much
> :> less in fp maths: in the above form I would need the equivalent of an
> :> mpf_root() function, which of course isn't in the library.
> :> 
> :> I think I can see how to converge on the result using the inequality:
> :>   p . c^n <= q . (c+1)^n
> :> using integer arithmetic only, effectively doing a binary chop (though
> :> I need to think whether I can find a way to calculate the upper bound
> :> without another set of trial calculations). Can anyone suggest a better
> :> way?
> :
> :a possible alternative is to use mpfr (mpfr.org), which extends mpf:
> :
> : -- Function: int mpfr_cbrt (mpfr_t ROP, mpfr_t OP, mp_rnd_t RND)
> : -- Function: int mpfr_root (mpfr_t ROP, mpfr_t OP, unsigned long int
> :          K, mp_rnd_t RND)
> :     Set ROP to the cubic root (resp. the Kth root) of OP rounded in
> :     the direction RND.  An odd (resp. even) root of a negative number
> :     (including -Inf) returns a negative number (resp. NaN).  The Kth
> :     root of -0 is defined to be -0, whatever the parity of K.
> 
> Thanks, I'll note that for next time. For the current project I've come
> to the conclusion that I simply don't know how to ensure sufficient
> accuracy, so have decided to go for the integer approach (which is not
> disastrously slow).
> 
> (This was brought home to me when a difference between results reported
> by the PARI/GP code and the C/GMP version turned out to be because
> floor(1/(sqrtn(21^3/20^3, 3)-1)) returned 19 in PARI, which caused it
> to miss one of the solutions it was supposed to be counting. Replacing
> (p/q) with ((p+1)/q) would be sufficient to fix that for this case, but
> I would need to do some maths to determine whether it would guarantee me
> correct results for the full range of the numbers I'm dealing with.)
> 
> Hugo
> _______________________________________________
> gmp-discuss mailing list
> gmp-discuss at swox.com
> https://gmplib.org/mailman/listinfo/gmp-discuss
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 189 bytes
Desc: This is a digitally signed message part
Url : http://gmplib.org/list-archives/gmp-discuss/attachments/20071228/e298adca/attachment.bin