mpf_root?

[hv at crypt.org]

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