mpf_root?
hv at crypt.org
hv at crypt.org
Fri Dec 28 00:26:43 CET 2007
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
More information about the gmp-discuss
mailing list