|
Packit |
5c3484 |
/* mpn_get_d -- limbs to double conversion.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY. THEY'RE ALMOST
|
|
Packit |
5c3484 |
CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
|
|
Packit |
5c3484 |
FUTURE GNU MP RELEASES.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
Copyright 2003, 2004, 2007, 2009, 2010, 2012 Free Software Foundation, Inc.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
This file is part of the GNU MP Library.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
The GNU MP Library is free software; you can redistribute it and/or modify
|
|
Packit |
5c3484 |
it under the terms of either:
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
* the GNU Lesser General Public License as published by the Free
|
|
Packit |
5c3484 |
Software Foundation; either version 3 of the License, or (at your
|
|
Packit |
5c3484 |
option) any later version.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
or
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
* the GNU General Public License as published by the Free Software
|
|
Packit |
5c3484 |
Foundation; either version 2 of the License, or (at your option) any
|
|
Packit |
5c3484 |
later version.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
or both in parallel, as here.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
The GNU MP Library is distributed in the hope that it will be useful, but
|
|
Packit |
5c3484 |
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
|
|
Packit |
5c3484 |
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
Packit |
5c3484 |
for more details.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
You should have received copies of the GNU General Public License and the
|
|
Packit |
5c3484 |
GNU Lesser General Public License along with the GNU MP Library. If not,
|
|
Packit |
5c3484 |
see https://www.gnu.org/licenses/. */
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#include "gmp.h"
|
|
Packit |
5c3484 |
#include "gmp-impl.h"
|
|
Packit |
5c3484 |
#include "longlong.h"
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#ifndef _GMP_IEEE_FLOATS
|
|
Packit |
5c3484 |
#define _GMP_IEEE_FLOATS 0
|
|
Packit |
5c3484 |
#endif
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* To force use of the generic C code for testing, put
|
|
Packit |
5c3484 |
"#define _GMP_IEEE_FLOATS 0" at this point. */
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* In alpha gcc prior to 3.4, signed DI comparisons involving constants are
|
|
Packit |
5c3484 |
rearranged from "x < n" to "x+(-n) < 0", which is of course hopelessly
|
|
Packit |
5c3484 |
wrong if that addition overflows.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
The workaround here avoids this bug by ensuring n is not a literal constant.
|
|
Packit |
5c3484 |
Note that this is alpha specific. The offending transformation is/was in
|
|
Packit |
5c3484 |
alpha.c alpha_emit_conditional_branch() under "We want to use cmpcc/bcc".
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
Bizarrely, this happens also with Cray cc on alphaev5-cray-unicosmk2.0.6.X,
|
|
Packit |
5c3484 |
and has the same solution. Don't know why or how. */
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#if HAVE_HOST_CPU_FAMILY_alpha \
|
|
Packit |
5c3484 |
&& ((defined (__GNUC__) && ! __GMP_GNUC_PREREQ(3,4)) \
|
|
Packit |
5c3484 |
|| defined (_CRAY))
|
|
Packit |
5c3484 |
static volatile const long CONST_1024 = 1024;
|
|
Packit |
5c3484 |
static volatile const long CONST_NEG_1023 = -1023;
|
|
Packit |
5c3484 |
static volatile const long CONST_NEG_1022_SUB_53 = -1022 - 53;
|
|
Packit |
5c3484 |
#else
|
|
Packit |
5c3484 |
#define CONST_1024 (1024)
|
|
Packit |
5c3484 |
#define CONST_NEG_1023 (-1023)
|
|
Packit |
5c3484 |
#define CONST_NEG_1022_SUB_53 (-1022 - 53)
|
|
Packit |
5c3484 |
#endif
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Return the value {ptr,size}*2^exp, and negative if sign<0. Must have
|
|
Packit |
5c3484 |
size>=1, and a non-zero high limb ptr[size-1].
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
When we know the fp format, the result is truncated towards zero. This is
|
|
Packit |
5c3484 |
consistent with other gmp conversions, like mpz_set_f or mpz_set_q, and is
|
|
Packit |
5c3484 |
easy to implement and test.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
When we do not know the format, such truncation seems much harder. One
|
|
Packit |
5c3484 |
would need to defeat any rounding mode, including round-up.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
It's felt that GMP is not primarily concerned with hardware floats, and
|
|
Packit |
5c3484 |
really isn't enhanced by getting involved with hardware rounding modes
|
|
Packit |
5c3484 |
(which could even be some weird unknown style), so something unambiguous and
|
|
Packit |
5c3484 |
straightforward is best.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
The IEEE code below is the usual case, it knows either a 32-bit or 64-bit
|
|
Packit |
5c3484 |
limb and is done with shifts and masks. The 64-bit case in particular
|
|
Packit |
5c3484 |
should come out nice and compact.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
The generic code used to work one bit at a time, which was not only slow,
|
|
Packit |
5c3484 |
but implicitly relied upon denorms for intermediates, since the lowest bits'
|
|
Packit |
5c3484 |
weight of a perfectly valid fp number underflows in non-denorm. Therefore,
|
|
Packit |
5c3484 |
the generic code now works limb-per-limb, initially creating a number x such
|
|
Packit |
5c3484 |
that 1 <= x <= BASE. (BASE is reached only as result of rounding.) Then
|
|
Packit |
5c3484 |
x's exponent is scaled with explicit code (not ldexp to avoid libm
|
|
Packit |
5c3484 |
dependency). It is a tap-dance to avoid underflow or overflow, beware!
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
Traps:
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
Hardware traps for overflow to infinity, underflow to zero, or unsupported
|
|
Packit |
5c3484 |
denorms may or may not be taken. The IEEE code works bitwise and so
|
|
Packit |
5c3484 |
probably won't trigger them, the generic code works by float operations and
|
|
Packit |
5c3484 |
so probably will. This difference might be thought less than ideal, but
|
|
Packit |
5c3484 |
again its felt straightforward code is better than trying to get intimate
|
|
Packit |
5c3484 |
with hardware exceptions (of perhaps unknown nature).
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
Not done:
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
mpz_get_d in the past handled size==1 with a cast limb->double. This might
|
|
Packit |
5c3484 |
still be worthwhile there (for up to the mantissa many bits), but for
|
|
Packit |
5c3484 |
mpn_get_d here, the cost of applying "exp" to the resulting exponent would
|
|
Packit |
5c3484 |
probably use up any benefit a cast may have over bit twiddling. Also, if
|
|
Packit |
5c3484 |
the exponent is pushed into denorm range then bit twiddling is the only
|
|
Packit |
5c3484 |
option, to ensure the desired truncation is obtained.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
Other:
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
For reference, note that HPPA 8000, 8200, 8500 and 8600 trap FCNV,UDW,DBL
|
|
Packit |
5c3484 |
to the kernel for values >= 2^63. This makes it slow, and worse the kernel
|
|
Packit |
5c3484 |
Linux (what versions?) apparently uses untested code in its trap handling
|
|
Packit |
5c3484 |
routines, and gets the sign wrong. We don't use such a limb-to-double
|
|
Packit |
5c3484 |
cast, neither in the IEEE or generic code. */
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#undef FORMAT_RECOGNIZED
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
double
|
|
Packit |
5c3484 |
mpn_get_d (mp_srcptr up, mp_size_t size, mp_size_t sign, long exp)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
int lshift, nbits;
|
|
Packit |
5c3484 |
mp_limb_t x, mhi, mlo;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
ASSERT (size >= 0);
|
|
Packit |
5c3484 |
ASSERT_MPN (up, size);
|
|
Packit |
5c3484 |
ASSERT (size == 0 || up[size-1] != 0);
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
if (size == 0)
|
|
Packit |
5c3484 |
return 0.0;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Adjust exp to a radix point just above {up,size}, guarding against
|
|
Packit |
5c3484 |
overflow. After this exp can of course be reduced to anywhere within
|
|
Packit |
5c3484 |
the {up,size} region without underflow. */
|
|
Packit |
5c3484 |
if (UNLIKELY ((unsigned long) (GMP_NUMB_BITS * size)
|
|
Packit |
5c3484 |
> ((unsigned long) LONG_MAX - exp)))
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
#if _GMP_IEEE_FLOATS
|
|
Packit |
5c3484 |
goto ieee_infinity;
|
|
Packit |
5c3484 |
#endif
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* generic */
|
|
Packit |
5c3484 |
exp = LONG_MAX;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
else
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
exp += GMP_NUMB_BITS * size;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#if _GMP_IEEE_FLOATS
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
union ieee_double_extract u;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
up += size;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#if GMP_LIMB_BITS == 64
|
|
Packit |
5c3484 |
mlo = up[-1];
|
|
Packit |
5c3484 |
count_leading_zeros (lshift, mlo);
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
exp -= (lshift - GMP_NAIL_BITS) + 1;
|
|
Packit |
5c3484 |
mlo <<= lshift;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
nbits = GMP_LIMB_BITS - lshift;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
if (nbits < 53 && size > 1)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
x = up[-2];
|
|
Packit |
5c3484 |
x <<= GMP_NAIL_BITS;
|
|
Packit |
5c3484 |
x >>= nbits;
|
|
Packit |
5c3484 |
mlo |= x;
|
|
Packit |
5c3484 |
nbits += GMP_NUMB_BITS;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
if (LIMBS_PER_DOUBLE >= 3 && nbits < 53 && size > 2)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
x = up[-3];
|
|
Packit |
5c3484 |
x <<= GMP_NAIL_BITS;
|
|
Packit |
5c3484 |
x >>= nbits;
|
|
Packit |
5c3484 |
mlo |= x;
|
|
Packit |
5c3484 |
nbits += GMP_NUMB_BITS;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
mhi = mlo >> (32 + 11);
|
|
Packit |
5c3484 |
mlo = mlo >> 11; /* later implicitly truncated to 32 bits */
|
|
Packit |
5c3484 |
#endif
|
|
Packit |
5c3484 |
#if GMP_LIMB_BITS == 32
|
|
Packit |
5c3484 |
x = *--up;
|
|
Packit |
5c3484 |
count_leading_zeros (lshift, x);
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
exp -= (lshift - GMP_NAIL_BITS) + 1;
|
|
Packit |
5c3484 |
x <<= lshift;
|
|
Packit |
5c3484 |
mhi = x >> 11;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
if (lshift < 11) /* FIXME: never true if NUMB < 20 bits */
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
/* All 20 bits in mhi */
|
|
Packit |
5c3484 |
mlo = x << 21;
|
|
Packit |
5c3484 |
/* >= 1 bit in mlo */
|
|
Packit |
5c3484 |
nbits = GMP_LIMB_BITS - lshift - 21;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
else
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
if (size > 1)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
nbits = GMP_LIMB_BITS - lshift;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
x = *--up, size--;
|
|
Packit |
5c3484 |
x <<= GMP_NAIL_BITS;
|
|
Packit |
5c3484 |
mhi |= x >> nbits >> 11;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
mlo = x << GMP_LIMB_BITS - nbits - 11;
|
|
Packit |
5c3484 |
nbits = nbits + 11 - GMP_NAIL_BITS;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
else
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
mlo = 0;
|
|
Packit |
5c3484 |
goto done;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Now all needed bits in mhi have been accumulated. Add bits to mlo. */
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
if (LIMBS_PER_DOUBLE >= 2 && nbits < 32 && size > 1)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
x = up[-1];
|
|
Packit |
5c3484 |
x <<= GMP_NAIL_BITS;
|
|
Packit |
5c3484 |
x >>= nbits;
|
|
Packit |
5c3484 |
mlo |= x;
|
|
Packit |
5c3484 |
nbits += GMP_NUMB_BITS;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
if (LIMBS_PER_DOUBLE >= 3 && nbits < 32 && size > 2)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
x = up[-2];
|
|
Packit |
5c3484 |
x <<= GMP_NAIL_BITS;
|
|
Packit |
5c3484 |
x >>= nbits;
|
|
Packit |
5c3484 |
mlo |= x;
|
|
Packit |
5c3484 |
nbits += GMP_NUMB_BITS;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
if (LIMBS_PER_DOUBLE >= 4 && nbits < 32 && size > 3)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
x = up[-3];
|
|
Packit |
5c3484 |
x <<= GMP_NAIL_BITS;
|
|
Packit |
5c3484 |
x >>= nbits;
|
|
Packit |
5c3484 |
mlo |= x;
|
|
Packit |
5c3484 |
nbits += GMP_NUMB_BITS;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
done:;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#endif
|
|
Packit |
5c3484 |
if (UNLIKELY (exp >= CONST_1024))
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
/* overflow, return infinity */
|
|
Packit |
5c3484 |
ieee_infinity:
|
|
Packit |
5c3484 |
mhi = 0;
|
|
Packit |
5c3484 |
mlo = 0;
|
|
Packit |
5c3484 |
exp = 1024;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
else if (UNLIKELY (exp <= CONST_NEG_1023))
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
int rshift;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
if (LIKELY (exp <= CONST_NEG_1022_SUB_53))
|
|
Packit |
5c3484 |
return 0.0; /* denorm underflows to zero */
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
rshift = -1022 - exp;
|
|
Packit |
5c3484 |
ASSERT (rshift > 0 && rshift < 53);
|
|
Packit |
5c3484 |
#if GMP_LIMB_BITS > 53
|
|
Packit |
5c3484 |
mlo >>= rshift;
|
|
Packit |
5c3484 |
mhi = mlo >> 32;
|
|
Packit |
5c3484 |
#else
|
|
Packit |
5c3484 |
if (rshift >= 32)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
mlo = mhi;
|
|
Packit |
5c3484 |
mhi = 0;
|
|
Packit |
5c3484 |
rshift -= 32;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
lshift = GMP_LIMB_BITS - rshift;
|
|
Packit |
5c3484 |
mlo = (mlo >> rshift) | (rshift == 0 ? 0 : mhi << lshift);
|
|
Packit |
5c3484 |
mhi >>= rshift;
|
|
Packit |
5c3484 |
#endif
|
|
Packit |
5c3484 |
exp = -1023;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
u.s.manh = mhi;
|
|
Packit |
5c3484 |
u.s.manl = mlo;
|
|
Packit |
5c3484 |
u.s.exp = exp + 1023;
|
|
Packit |
5c3484 |
u.s.sig = (sign < 0);
|
|
Packit |
5c3484 |
return u.d;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
#define FORMAT_RECOGNIZED 1
|
|
Packit |
5c3484 |
#endif
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#if HAVE_DOUBLE_VAX_D
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
union double_extract u;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
up += size;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
mhi = up[-1];
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
count_leading_zeros (lshift, mhi);
|
|
Packit |
5c3484 |
exp -= lshift;
|
|
Packit |
5c3484 |
mhi <<= lshift;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
mlo = 0;
|
|
Packit |
5c3484 |
if (size > 1)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
mlo = up[-2];
|
|
Packit |
5c3484 |
if (lshift != 0)
|
|
Packit |
5c3484 |
mhi += mlo >> (GMP_LIMB_BITS - lshift);
|
|
Packit |
5c3484 |
mlo <<= lshift;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
if (size > 2 && lshift > 8)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
x = up[-3];
|
|
Packit |
5c3484 |
mlo += x >> (GMP_LIMB_BITS - lshift);
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
if (UNLIKELY (exp >= 128))
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
/* overflow, return maximum number */
|
|
Packit |
5c3484 |
mhi = 0xffffffff;
|
|
Packit |
5c3484 |
mlo = 0xffffffff;
|
|
Packit |
5c3484 |
exp = 127;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
else if (UNLIKELY (exp < -128))
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
return 0.0; /* underflows to zero */
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
u.s.man3 = mhi >> 24; /* drop msb, since implicit */
|
|
Packit |
5c3484 |
u.s.man2 = mhi >> 8;
|
|
Packit |
5c3484 |
u.s.man1 = (mhi << 8) + (mlo >> 24);
|
|
Packit |
5c3484 |
u.s.man0 = mlo >> 8;
|
|
Packit |
5c3484 |
u.s.exp = exp + 128;
|
|
Packit |
5c3484 |
u.s.sig = sign < 0;
|
|
Packit |
5c3484 |
return u.d;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
#define FORMAT_RECOGNIZED 1
|
|
Packit |
5c3484 |
#endif
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#if ! FORMAT_RECOGNIZED
|
|
Packit |
5c3484 |
{ /* Non-IEEE or strange limb size, do something generic. */
|
|
Packit |
5c3484 |
mp_size_t i;
|
|
Packit |
5c3484 |
double d, weight;
|
|
Packit |
5c3484 |
unsigned long uexp;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* First generate an fp number disregarding exp, instead keeping things
|
|
Packit |
5c3484 |
within the numb base factor from 1, which should prevent overflow and
|
|
Packit |
5c3484 |
underflow even for the most exponent limited fp formats. The
|
|
Packit |
5c3484 |
termination criteria should be refined, since we now include too many
|
|
Packit |
5c3484 |
limbs. */
|
|
Packit |
5c3484 |
weight = 1/MP_BASE_AS_DOUBLE;
|
|
Packit |
5c3484 |
d = up[size - 1];
|
|
Packit |
5c3484 |
for (i = size - 2; i >= 0; i--)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
d += up[i] * weight;
|
|
Packit |
5c3484 |
weight /= MP_BASE_AS_DOUBLE;
|
|
Packit |
5c3484 |
if (weight == 0)
|
|
Packit |
5c3484 |
break;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Now apply exp. */
|
|
Packit |
5c3484 |
exp -= GMP_NUMB_BITS;
|
|
Packit |
5c3484 |
if (exp > 0)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
weight = 2.0;
|
|
Packit |
5c3484 |
uexp = exp;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
else
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
weight = 0.5;
|
|
Packit |
5c3484 |
uexp = 1 - (unsigned long) (exp + 1);
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
#if 1
|
|
Packit |
5c3484 |
/* Square-and-multiply exponentiation. */
|
|
Packit |
5c3484 |
if (uexp & 1)
|
|
Packit |
5c3484 |
d *= weight;
|
|
Packit |
5c3484 |
while (uexp >>= 1)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
weight *= weight;
|
|
Packit |
5c3484 |
if (uexp & 1)
|
|
Packit |
5c3484 |
d *= weight;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
#else
|
|
Packit |
5c3484 |
/* Plain exponentiation. */
|
|
Packit |
5c3484 |
while (uexp > 0)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
d *= weight;
|
|
Packit |
5c3484 |
uexp--;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
#endif
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
return sign >= 0 ? d : -d;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
#endif
|
|
Packit |
5c3484 |
}
|