error in floor

Genius 0.4.3
Copyright (c) 1997,1998,1999 Free Software Foundation, Inc.
This is free software with ABSOLUTELY NO WARRANTY.
For details type `warranty'.

genius> floor(73! * 1.0) - 73!
= 0
genius> floor(74! * 1.0) - 74!
= -39306667005013562067473399808
genius> floor(74!) - 74!
= 0
genius> quit

This is the unmodified version.

Basically 74! doesn't fit into a double, so mpf_get_d(74!*1.0) doesn't
return the correct value. FWIW it returns the correct decimal digits, it's
off by a pwoer of 10. (I find that strange to say the least).

FIX: have floor return a mpz_t instead of a double.

The gmp C++ stuff I've written I am currently translating back to C and
I'll pass that along. It has some functions to do this reasonably
correctly.

This example code is untested in C form, but if you want to use it as a
start:

ManI returns the mantissa as an integer.
ExpI returns the power of 2 to multiply the mantissa by to get the
        original number.
Floor does the obvious.
Init is used to take the various fields of an mpz_t and initialize it.
Shift is a wrapper macro for shifting a _signed_ amount.
Cntlz takes an int and tells how many leading (most significant) bits are
        zero. It is a single instruction on many architectures, but a
        somewhat lengthy c routine. gmp contains a definition of it, but
        you have some #define and #include messiness. I'll send a cleaned
        up version with the rest of the C code after its been checked abit
        more.

no warranty, copyright under GPL blah blah:

/* get count_leading_zeros from gmp source. it is a macro that
   calls assembly versions. see longlong.h for details.
 */

int cntlz (mp_limb_t l) { int c; count_leading_zeros(c,l); return c; }

void mpf_manI (mpz_t rop,mpf_srcptr z)
{
        mpz_init_si_ui_limbp (rop,z[0]._mp_size, z[0]._mp_prec + 1,
z[0]._mp_d);
}

mp_exp_t mpf_exp2 (mpf_srcptr z)
{
        return z[0]._mp_size ?
                z[0]._mp_exp * BITS_PER_MP_LIMB - cntlz (z[0]._mp_d[ABS
(z[0]._mp_size) - 1]) - 1 : 0;
}

void mpf_floor (mpz_t rop,mpf_srcptr z)
{
        mpf_manI(rop,z);
        mpz_shift_left(rop,rop,mpf_expI(z));
}

void mpz_init_si_ui_limbp(mpz_t rop,int size,unsigned int alloc,mp_limb_t
* limbs)
{
        rop[0]._mp_size=size;
        rop[0]._mp_alloc=alloc;
        rop[0]._mp_d=malloc(alloc*BYTES_PER_MP_LIMB);
        memcpy(rop[0]._mp_d,limbs,alloc*BYTES_PER_MP_LIMB);
}

void mpz_shift_left (mpz_t rop, mpz_srcptr op1, long shift)
{
        if (shift >= 0)
                mpz_mul_2exp (rop, op1, shift);
        else
                mpz_tdiv_q_2exp (rop, op1, -shift);
}
Received on Sat Jun 05 1999 - 19:50:38 CDT