const REAL_VALUE_TYPE *));
static void normalize PARAMS ((REAL_VALUE_TYPE *));
-static void do_add PARAMS ((REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
+static bool do_add PARAMS ((REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
const REAL_VALUE_TYPE *, int));
-static void do_multiply PARAMS ((REAL_VALUE_TYPE *,
+static bool do_multiply PARAMS ((REAL_VALUE_TYPE *,
const REAL_VALUE_TYPE *,
const REAL_VALUE_TYPE *));
-static void do_divide PARAMS ((REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
+static bool do_divide PARAMS ((REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
const REAL_VALUE_TYPE *));
static int do_compare PARAMS ((const REAL_VALUE_TYPE *,
const REAL_VALUE_TYPE *, int));
memset (r, 0, sizeof (*r));
r->class = rvc_nan;
r->sign = sign;
- r->sig[SIGSZ-1] = SIG_MSB >> 1;
+ r->canonical = 1;
}
static inline void
memset (r, 0, sizeof (*r));
r->class = rvc_nan;
r->sign = sign;
- r->sig[SIGSZ-1] = SIG_MSB >> 2;
+ r->signalling = 1;
+ r->canonical = 1;
}
static inline void
}
}
\f
-/* Return R = A + (SUBTRACT_P ? -B : B). */
+/* Calculate R = A + (SUBTRACT_P ? -B : B). Return true if the
+ result may be inexact due to a loss of precision. */
-static void
+static bool
do_add (r, a, b, subtract_p)
REAL_VALUE_TYPE *r;
const REAL_VALUE_TYPE *a, *b;
case CLASS2 (rvc_zero, rvc_zero):
/* -0 + -0 = -0, -0 - +0 = -0; all other cases yield +0. */
get_zero (r, sign & !subtract_p);
- return;
+ return false;
case CLASS2 (rvc_zero, rvc_normal):
case CLASS2 (rvc_zero, rvc_inf):
/* R + Inf = Inf. */
*r = *b;
r->sign = sign ^ subtract_p;
- return;
+ return false;
case CLASS2 (rvc_normal, rvc_zero):
case CLASS2 (rvc_inf, rvc_zero):
case CLASS2 (rvc_inf, rvc_normal):
/* Inf + R = Inf. */
*r = *a;
- return;
+ return false;
case CLASS2 (rvc_inf, rvc_inf):
if (subtract_p)
else
/* Inf + Inf = Inf. */
*r = *a;
- return;
+ return false;
case CLASS2 (rvc_normal, rvc_normal):
break;
{
*r = *a;
r->sign = sign;
- return;
+ return true;
}
inexact |= sticky_rshift_significand (&t, b, dexp);
if (++exp > MAX_EXP)
{
get_inf (r, sign);
- return;
+ return true;
}
}
}
r->sign = 0;
else
r->sig[0] |= inexact;
+
+ return inexact;
}
-/* Return R = A * B. */
+/* Calculate R = A * B. Return true if the result may be inexact. */
-static void
+static bool
do_multiply (r, a, b)
REAL_VALUE_TYPE *r;
const REAL_VALUE_TYPE *a, *b;
REAL_VALUE_TYPE u, t, *rr;
unsigned int i, j, k;
int sign = a->sign ^ b->sign;
+ bool inexact = false;
switch (CLASS2 (a->class, b->class))
{
case CLASS2 (rvc_normal, rvc_zero):
/* +-0 * ANY = 0 with appropriate sign. */
get_zero (r, sign);
- return;
+ return false;
case CLASS2 (rvc_zero, rvc_nan):
case CLASS2 (rvc_normal, rvc_nan):
/* ANY * NaN = NaN. */
*r = *b;
r->sign = sign;
- return;
+ return false;
case CLASS2 (rvc_nan, rvc_zero):
case CLASS2 (rvc_nan, rvc_normal):
/* NaN * ANY = NaN. */
*r = *a;
r->sign = sign;
- return;
+ return false;
case CLASS2 (rvc_zero, rvc_inf):
case CLASS2 (rvc_inf, rvc_zero):
/* 0 * Inf = NaN */
get_canonical_qnan (r, sign);
- return;
+ return false;
case CLASS2 (rvc_inf, rvc_inf):
case CLASS2 (rvc_normal, rvc_inf):
case CLASS2 (rvc_inf, rvc_normal):
/* Inf * Inf = Inf, R * Inf = Inf */
- overflow:
get_inf (r, sign);
- return;
+ return false;
case CLASS2 (rvc_normal, rvc_normal):
break;
+ (b->exp - (1-j)*(HOST_BITS_PER_LONG/2)));
if (exp > MAX_EXP)
- goto overflow;
+ {
+ get_inf (r, sign);
+ return true;
+ }
if (exp < -MAX_EXP)
- /* Would underflow to zero, which we shouldn't bother adding. */
- continue;
+ {
+ /* Would underflow to zero, which we shouldn't bother adding. */
+ inexact = true;
+ continue;
+ }
u.class = rvc_normal;
u.sign = 0;
}
normalize (&u);
- do_add (rr, rr, &u, 0);
+ inexact |= do_add (rr, rr, &u, 0);
}
}
rr->sign = sign;
if (rr != r)
*r = t;
+
+ return inexact;
}
-/* Return R = A / B. */
+/* Calculate R = A / B. Return true if the result may be inexact. */
-static void
+static bool
do_divide (r, a, b)
REAL_VALUE_TYPE *r;
const REAL_VALUE_TYPE *a, *b;
case CLASS2 (rvc_inf, rvc_inf):
/* Inf / Inf = NaN. */
get_canonical_qnan (r, sign);
- return;
+ return false;
case CLASS2 (rvc_zero, rvc_normal):
case CLASS2 (rvc_zero, rvc_inf):
/* 0 / ANY = 0. */
case CLASS2 (rvc_normal, rvc_inf):
/* R / Inf = 0. */
- underflow:
get_zero (r, sign);
- return;
+ return false;
case CLASS2 (rvc_normal, rvc_zero):
/* R / 0 = Inf. */
case CLASS2 (rvc_inf, rvc_zero):
/* Inf / 0 = Inf. */
get_inf (r, sign);
- return;
+ return false;
case CLASS2 (rvc_zero, rvc_nan):
case CLASS2 (rvc_normal, rvc_nan):
/* ANY / NaN = NaN. */
*r = *b;
r->sign = sign;
- return;
+ return false;
case CLASS2 (rvc_nan, rvc_zero):
case CLASS2 (rvc_nan, rvc_normal):
/* NaN / ANY = NaN. */
*r = *a;
r->sign = sign;
- return;
+ return false;
case CLASS2 (rvc_inf, rvc_normal):
/* Inf / R = Inf. */
- overflow:
get_inf (r, sign);
- return;
+ return false;
case CLASS2 (rvc_normal, rvc_normal):
break;
exp = a->exp - b->exp + 1;
if (exp > MAX_EXP)
- goto overflow;
+ {
+ get_inf (r, sign);
+ return true;
+ }
if (exp < -MAX_EXP)
- goto underflow;
+ {
+ get_zero (r, sign);
+ return true;
+ }
rr->exp = exp;
inexact = div_significands (rr, a, b);
if (rr != r)
*r = t;
+
+ return inexact;
}
/* Return a tri-state comparison of A vs B. Return NAN_RESULT if
/* Compare two floating-point objects for bitwise identity. */
-extern bool
+bool
real_identical (a, b)
const REAL_VALUE_TYPE *a, *b;
{
{
case rvc_zero:
case rvc_inf:
- break;
+ return true;
case rvc_normal:
if (a->exp != b->exp)
return false;
- /* FALLTHRU */
+ break;
+
case rvc_nan:
- for (i = 0; i < SIGSZ; ++i)
- if (a->sig[i] != b->sig[i])
- return false;
+ if (a->signalling != b->signalling)
+ return false;
+ /* The significand is ignored for canonical NaNs. */
+ if (a->canonical || b->canonical)
+ return a->canonical == b->canonical;
break;
default:
abort ();
}
+ for (i = 0; i < SIGSZ; ++i)
+ if (a->sig[i] != b->sig[i])
+ return false;
+
return true;
}
case rvc_normal:
if (r->exp <= 0)
goto underflow;
+ /* Only force overflow for unsigned overflow. Signed overflow is
+ undefined, so it doesn't matter what we return, and some callers
+ expect to be able to use this routine for both signed and
+ unsigned conversions. */
if (r->exp > HOST_BITS_PER_WIDE_INT)
goto overflow;
exp = r->exp;
if (exp <= 0)
goto underflow;
+ /* Only force overflow for unsigned overflow. Signed overflow is
+ undefined, so it doesn't matter what we return, and some callers
+ expect to be able to use this routine for both signed and
+ unsigned conversions. */
if (exp > 2*HOST_BITS_PER_WIDE_INT)
goto overflow;
/* Shift the significand into place such that the bits
are in the most significant bits for the format. */
- lshift_significand (r, r, SIGNIFICAND_BITS - fmt->p);
+ lshift_significand (r, r, SIGNIFICAND_BITS - fmt->pnan);
/* Our MSB is always unset for NaNs. */
r->sig[SIGSZ-1] &= ~SIG_MSB;
/* Force quiet or signalling NaN. */
- if (quiet)
- r->sig[SIGSZ-1] |= SIG_MSB >> 1;
- else
- r->sig[SIGSZ-1] &= ~(SIG_MSB >> 1);
-
- /* Force at least one bit of the significand set. */
- for (d = 0; d < SIGSZ; ++d)
- if (r->sig[d])
- break;
- if (d == SIGSZ)
- r->sig[SIGSZ-1] |= SIG_MSB >> 2;
-
- /* Our intermediate format forces QNaNs to have MSB-1 set.
- If the target format has QNaNs with the top bit unset,
- mirror the output routines and invert the top two bits. */
- if (!fmt->qnan_msb_set)
- r->sig[SIGSZ-1] ^= (SIG_MSB >> 1) | (SIG_MSB >> 2);
+ r->signalling = !quiet;
}
return true;
}
+/* Fills R with the largest finite value representable in mode MODE.
+ If SIGN is nonzero, R is set to the most negative finite value. */
+
+void
+real_maxval (r, sign, mode)
+ REAL_VALUE_TYPE *r;
+ int sign;
+ enum machine_mode mode;
+{
+ const struct real_format *fmt;
+ int np2;
+
+ fmt = real_format_for_mode[mode - QFmode];
+ if (fmt == NULL)
+ abort ();
+
+ r->class = rvc_normal;
+ r->sign = sign;
+ r->signalling = 0;
+ r->canonical = 0;
+ r->exp = fmt->emax * fmt->log2_b;
+
+ np2 = SIGNIFICAND_BITS - fmt->p * fmt->log2_b;
+ memset (r->sig, -1, SIGSZ * sizeof (unsigned long));
+ clear_significand_below (r, np2);
+}
+
/* Fills R with 2**N. */
void
case rvc_nan:
clear_significand_below (r, np2);
-
- /* If we've cleared the entire significand, we need one bit
- set for this to continue to be a NaN. */
- for (i = 0; i < SIGSZ; ++i)
- if (r->sig[i])
- break;
- if (i == SIGSZ)
- r->sig[SIGSZ-1] = SIG_MSB >> 2;
return;
case rvc_normal:
{
case rvc_zero:
case rvc_inf:
- break;
+ return h;
case rvc_normal:
h |= r->exp << 3;
- /* FALLTHRU */
+ break;
case rvc_nan:
- if (sizeof(unsigned long) > sizeof(unsigned int))
- for (i = 0; i < SIGSZ; ++i)
- {
- unsigned long s = r->sig[i];
- h ^= s ^ (s >> (HOST_BITS_PER_LONG / 2));
- }
- else
- for (i = 0; i < SIGSZ; ++i)
- h ^= r->sig[i];
+ if (r->signalling)
+ h ^= (unsigned int)-1;
+ if (r->canonical)
+ return h;
break;
default:
abort ();
}
+ if (sizeof(unsigned long) > sizeof(unsigned int))
+ for (i = 0; i < SIGSZ; ++i)
+ {
+ unsigned long s = r->sig[i];
+ h ^= s ^ (s >> (HOST_BITS_PER_LONG / 2));
+ }
+ else
+ for (i = 0; i < SIGSZ; ++i)
+ h ^= r->sig[i];
+
return h;
}
\f
case rvc_nan:
if (fmt->has_nans)
{
+ if (r->canonical)
+ sig = 0;
+ if (r->signalling == fmt->qnan_msb_set)
+ sig &= ~(1 << 22);
+ else
+ sig |= 1 << 22;
+ /* We overload qnan_msb_set here: it's only clear for
+ mips_ieee_single, which wants all mantissa bits but the
+ quiet/signalling one set in canonical NaNs (at least
+ Quiet ones). */
+ if (r->canonical && !fmt->qnan_msb_set)
+ sig |= (1 << 22) - 1;
+ else if (sig == 0)
+ sig = 1 << 21;
+
image |= 255 << 23;
image |= sig;
- if (!fmt->qnan_msb_set)
- image ^= 1 << 23 | 1 << 22;
}
else
image |= 0x7fffffff;
{
r->class = rvc_nan;
r->sign = sign;
- if (!fmt->qnan_msb_set)
- image ^= (SIG_MSB >> 1 | SIG_MSB >> 2);
+ r->signalling = (((image >> (HOST_BITS_PER_LONG - 2)) & 1)
+ ^ fmt->qnan_msb_set);
r->sig[SIGSZ-1] = image;
}
else
2,
1,
24,
+ 24,
-125,
128,
31,
true
};
+const struct real_format mips_single_format =
+ {
+ encode_ieee_single,
+ decode_ieee_single,
+ 2,
+ 1,
+ 24,
+ 24,
+ -125,
+ 128,
+ 31,
+ true,
+ true,
+ true,
+ true,
+ false
+ };
+
\f
/* IEEE double-precision format. */
case rvc_nan:
if (fmt->has_nans)
{
+ if (r->canonical)
+ sig_hi = sig_lo = 0;
+ if (r->signalling == fmt->qnan_msb_set)
+ sig_hi &= ~(1 << 19);
+ else
+ sig_hi |= 1 << 19;
+ /* We overload qnan_msb_set here: it's only clear for
+ mips_ieee_single, which wants all mantissa bits but the
+ quiet/signalling one set in canonical NaNs (at least
+ Quiet ones). */
+ if (r->canonical && !fmt->qnan_msb_set)
+ {
+ sig_hi |= (1 << 19) - 1;
+ sig_lo = 0xffffffff;
+ }
+ else if (sig_hi == 0 && sig_lo == 0)
+ sig_hi = 1 << 18;
+
image_hi |= 2047 << 20;
image_hi |= sig_hi;
- if (!fmt->qnan_msb_set)
- image_hi ^= 1 << 19 | 1 << 18;
image_lo = sig_lo;
}
else
{
r->class = rvc_nan;
r->sign = sign;
+ r->signalling = ((image_hi >> 30) & 1) ^ fmt->qnan_msb_set;
if (HOST_BITS_PER_LONG == 32)
{
r->sig[SIGSZ-1] = image_hi;
}
else
r->sig[SIGSZ-1] = (image_hi << 31 << 1) | image_lo;
-
- if (!fmt->qnan_msb_set)
- r->sig[SIGSZ-1] ^= (SIG_MSB >> 1 | SIG_MSB >> 2);
}
else
{
2,
1,
53,
+ 53,
-1021,
1024,
63,
true
};
+const struct real_format mips_double_format =
+ {
+ encode_ieee_double,
+ decode_ieee_double,
+ 2,
+ 1,
+ 53,
+ 53,
+ -1021,
+ 1024,
+ 63,
+ true,
+ true,
+ true,
+ true,
+ false
+ };
+
\f
/* IEEE extended double precision format. This comes in three
- flavours: Intel's as a 12 byte image, Intel's as a 16 byte image,
+ flavors: Intel's as a 12 byte image, Intel's as a 16 byte image,
and Motorola's. */
static void encode_ieee_extended PARAMS ((const struct real_format *fmt,
sig_hi = sig_lo >> 31 >> 1;
sig_lo &= 0xffffffff;
}
- if (!fmt->qnan_msb_set)
- sig_hi ^= 1 << 30 | 1 << 29;
+ if (r->signalling == fmt->qnan_msb_set)
+ sig_hi &= ~(1 << 30);
+ else
+ sig_hi |= 1 << 30;
+ if ((sig_hi & 0x7fffffff) == 0 && sig_lo == 0)
+ sig_hi = 1 << 29;
/* Intel requires the explicit integer bit to be set, otherwise
it considers the value a "pseudo-nan". Motorola docs say it
{
r->class = rvc_nan;
r->sign = sign;
+ r->signalling = ((sig_hi >> 30) & 1) ^ fmt->qnan_msb_set;
if (HOST_BITS_PER_LONG == 32)
{
r->sig[SIGSZ-1] = sig_hi;
}
else
r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
-
- if (!fmt->qnan_msb_set)
- r->sig[SIGSZ-1] ^= (SIG_MSB >> 1 | SIG_MSB >> 2);
}
else
{
2,
1,
64,
+ 64,
-16382,
16384,
95,
2,
1,
64,
+ 64,
-16381,
16384,
79,
2,
1,
64,
+ 64,
-16381,
16384,
79,
true
};
+/* The following caters to i386 systems that set the rounding precision
+ to 53 bits instead of 64, e.g. FreeBSD. */
+const struct real_format ieee_extended_intel_96_round_53_format =
+ {
+ encode_ieee_extended,
+ decode_ieee_extended,
+ 2,
+ 1,
+ 53,
+ 53,
+ -16381,
+ 16384,
+ 79,
+ true,
+ true,
+ true,
+ true,
+ true
+ };
\f
/* IBM 128-bit extended precision format: a pair of IEEE double precision
numbers whose sum is equal to the extended precision value. The number
static void
encode_ibm_extended (fmt, buf, r)
- const struct real_format *fmt ATTRIBUTE_UNUSED;
+ const struct real_format *fmt;
long *buf;
const REAL_VALUE_TYPE *r;
{
REAL_VALUE_TYPE u, v;
+ const struct real_format *base_fmt;
+
+ base_fmt = fmt->qnan_msb_set ? &ieee_double_format : &mips_double_format;
switch (r->class)
{
case rvc_inf:
case rvc_nan:
/* Both doubles set to Inf / NaN. */
- encode_ieee_double (&ieee_double_format, &buf[0], r);
+ encode_ieee_double (base_fmt, &buf[0], r);
buf[2] = buf[0];
buf[3] = buf[1];
return;
{
/* v = remainder containing additional 53 bits of significand. */
do_add (&v, r, &u, 1);
- round_for_format (&ieee_double_format, &v);
+ round_for_format (base_fmt, &v);
}
- round_for_format (&ieee_double_format, &u);
+ round_for_format (base_fmt, &u);
- encode_ieee_double (&ieee_double_format, &buf[0], &u);
- encode_ieee_double (&ieee_double_format, &buf[2], &v);
+ encode_ieee_double (base_fmt, &buf[0], &u);
+ encode_ieee_double (base_fmt, &buf[2], &v);
break;
default:
const long *buf;
{
REAL_VALUE_TYPE u, v;
+ const struct real_format *base_fmt;
- decode_ieee_double (&ieee_double_format, &u, &buf[0]);
+ base_fmt = fmt->qnan_msb_set ? &ieee_double_format : &mips_double_format;
+ decode_ieee_double (base_fmt, &u, &buf[0]);
if (u.class != rvc_zero && u.class != rvc_inf && u.class != rvc_nan)
{
- decode_ieee_double (&ieee_double_format, &v, &buf[2]);
+ decode_ieee_double (base_fmt, &v, &buf[2]);
do_add (r, &u, &v, 0);
}
else
2,
1,
53 + 53,
+ 53,
-1021 + 53,
1024,
-1,
true
};
+const struct real_format mips_extended_format =
+ {
+ encode_ibm_extended,
+ decode_ibm_extended,
+ 2,
+ 1,
+ 53 + 53,
+ 53,
+ -1021 + 53,
+ 1024,
+ -1,
+ true,
+ true,
+ true,
+ true,
+ false
+ };
+
\f
/* IEEE quad precision format. */
{
image3 |= 32767 << 16;
- if (HOST_BITS_PER_LONG == 32)
+ if (r->canonical)
+ {
+ /* Don't use bits from the significand. The
+ initialization above is right. */
+ }
+ else if (HOST_BITS_PER_LONG == 32)
{
image0 = u.sig[0];
image1 = u.sig[1];
image0 &= 0xffffffff;
image2 &= 0xffffffff;
}
-
- if (!fmt->qnan_msb_set)
- image3 ^= 1 << 15 | 1 << 14;
+ if (r->signalling == fmt->qnan_msb_set)
+ image3 &= ~0x8000;
+ else
+ image3 |= 0x8000;
+ /* We overload qnan_msb_set here: it's only clear for
+ mips_ieee_single, which wants all mantissa bits but the
+ quiet/signalling one set in canonical NaNs (at least
+ Quiet ones). */
+ if (r->canonical && !fmt->qnan_msb_set)
+ {
+ image3 |= 0x7fff;
+ image2 = image1 = image0 = 0xffffffff;
+ }
+ else if (((image3 & 0xffff) | image2 | image1 | image0) == 0)
+ image3 |= 0x4000;
}
else
{
{
r->class = rvc_nan;
r->sign = sign;
+ r->signalling = ((image3 >> 15) & 1) ^ fmt->qnan_msb_set;
if (HOST_BITS_PER_LONG == 32)
{
r->sig[1] = (image3 << 31 << 1) | image2;
}
lshift_significand (r, r, SIGNIFICAND_BITS - 113);
-
- if (!fmt->qnan_msb_set)
- r->sig[SIGSZ-1] ^= (SIG_MSB >> 1 | SIG_MSB >> 2);
}
else
{
2,
1,
113,
+ 113,
-16381,
16384,
127,
true,
true
};
+
+const struct real_format mips_quad_format =
+ {
+ encode_ieee_quad,
+ decode_ieee_quad,
+ 2,
+ 1,
+ 113,
+ 113,
+ -16381,
+ 16384,
+ 127,
+ true,
+ true,
+ true,
+ true,
+ false
+ };
\f
/* Descriptions of VAX floating point formats can be found beginning at
2,
1,
24,
+ 24,
-127,
127,
15,
2,
1,
56,
+ 56,
-127,
127,
15,
2,
1,
53,
+ 53,
-1023,
1023,
15,
16,
4,
6,
+ 6,
-64,
63,
31,
16,
4,
14,
+ 14,
-64,
63,
63,
2,
1,
24,
+ 24,
-126,
128,
-1,
2,
1,
32,
+ 32,
-126,
128,
-1,
2,
1,
SIGNIFICAND_BITS - 2,
+ SIGNIFICAND_BITS - 2,
-MAX_EXP,
MAX_EXP,
-1,
const REAL_VALUE_TYPE *x;
{
static REAL_VALUE_TYPE halfthree;
- static REAL_VALUE_TYPE half;
static bool init = false;
REAL_VALUE_TYPE h, t, i;
int iter, exp;
if (!init)
{
- real_arithmetic (&half, RDIV_EXPR, &dconst1, &dconst2);
- real_arithmetic (&halfthree, PLUS_EXPR, &dconst1, &half);
+ do_add (&halfthree, &dconst1, &dconsthalf, 0);
init = true;
}
for (iter = 0; iter < 16; iter++)
{
/* i(n+1) = i(n) * (1.5 - 0.5*i(n)*i(n)*x). */
- real_arithmetic (&t, MULT_EXPR, x, &i);
- real_arithmetic (&h, MULT_EXPR, &t, &i);
- real_arithmetic (&t, MULT_EXPR, &h, &half);
- real_arithmetic (&h, MINUS_EXPR, &halfthree, &t);
- real_arithmetic (&t, MULT_EXPR, &i, &h);
+ do_multiply (&t, x, &i);
+ do_multiply (&h, &t, &i);
+ do_multiply (&t, &h, &dconsthalf);
+ do_add (&h, &halfthree, &t, 1);
+ do_multiply (&t, &i, &h);
/* Check for early convergence. */
if (iter >= 6 && real_identical (&i, &t))
}
/* Final iteration: r = i*x + 0.5*i*x*(1.0 - i*(i*x)). */
- real_arithmetic (&t, MULT_EXPR, x, &i);
- real_arithmetic (&h, MULT_EXPR, &t, &i);
- real_arithmetic (&i, MINUS_EXPR, &dconst1, &h);
- real_arithmetic (&h, MULT_EXPR, &t, &i);
- real_arithmetic (&i, MULT_EXPR, &half, &h);
- real_arithmetic (&h, PLUS_EXPR, &t, &i);
+ do_multiply (&t, x, &i);
+ do_multiply (&h, &t, &i);
+ do_add (&i, &dconst1, &h, 1);
+ do_multiply (&h, &t, &i);
+ do_multiply (&i, &dconsthalf, &h);
+ do_add (&h, &t, &i, 0);
/* ??? We need a Tuckerman test to get the last bit. */
return true;
}
+/* Calculate X raised to the integer exponent N in mode MODE and store
+ the result in R. Return true if the result may be inexact due to
+ loss of precision. The algorithm is the classic "left-to-right binary
+ method" described in section 4.6.3 of Donald Knuth's "Seminumerical
+ Algorithms", "The Art of Computer Programming", Volume 2. */
+
+bool
+real_powi (r, mode, x, n)
+ REAL_VALUE_TYPE *r;
+ enum machine_mode mode;
+ const REAL_VALUE_TYPE *x;
+ HOST_WIDE_INT n;
+{
+ unsigned HOST_WIDE_INT bit;
+ REAL_VALUE_TYPE t;
+ bool inexact = false;
+ bool init = false;
+ bool neg;
+ int i;
+
+ if (n == 0)
+ {
+ *r = dconst1;
+ return false;
+ }
+ else if (n < 0)
+ {
+ /* Don't worry about overflow, from now on n is unsigned. */
+ neg = true;
+ n = -n;
+ }
+ else
+ neg = false;
+
+ t = *x;
+ bit = (unsigned HOST_WIDE_INT) 1 << (HOST_BITS_PER_WIDE_INT - 1);
+ for (i = 0; i < HOST_BITS_PER_WIDE_INT; i++)
+ {
+ if (init)
+ {
+ inexact |= do_multiply (&t, &t, &t);
+ if (n & bit)
+ inexact |= do_multiply (&t, &t, x);
+ }
+ else if (n & bit)
+ init = true;
+ bit >>= 1;
+ }
+
+ if (neg)
+ inexact |= do_divide (&t, &dconst1, &t);
+
+ real_convert (r, mode, &t);
+ return inexact;
+}
+
OSDN Git Service
