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));
}
}
\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
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
\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,
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
if (!init)
{
- real_arithmetic (&halfthree, PLUS_EXPR, &dconst1, &dconsthalf);
+ 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, &dconsthalf);
- 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, &dconsthalf, &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
