X-Git-Url: http://git.shadowcat.co.uk/gitweb/gitweb.cgi?a=blobdiff_plain;f=numeric.c;h=b193e0c2e4a30c621c4a074542c67f06635449d9;hb=439ba5457a8422144686c1df300aa1dde218dbfd;hp=d2144164cd99cfafa979c6b63176785570395ee7;hpb=a36244b753e9b2747b70a6f37a7bcd86eaa61bcb;p=p5sagit%2Fp5-mst-13.2.git diff --git a/numeric.c b/numeric.c index d214416..b193e0c 100644 --- a/numeric.c +++ b/numeric.c @@ -727,10 +727,6 @@ S_mulexp10(NV value, I32 exponent) if (exponent == 0) return value; - else if (exponent < 0) { - negative = 1; - exponent = -exponent; - } /* On OpenVMS VAX we by default use the D_FLOAT double format, * and that format does not have *easy* capabilities [1] for @@ -743,11 +739,6 @@ S_mulexp10(NV value, I32 exponent) * * [1] Trying to establish a condition handler to trap floating point * exceptions is not a good idea. */ -#if defined(VMS) && !defined(__IEEE_FP) && defined(NV_MAX_10_EXP) - if (!negative && - (log10(value) + exponent) >= (NV_MAX_10_EXP)) - return NV_MAX; -#endif /* In UNICOS and in certain Cray models (such as T90) there is no * IEEE fp, and no way at all from C to catch fp overflows gracefully. @@ -756,18 +747,37 @@ S_mulexp10(NV value, I32 exponent) * disable *all* floating point interrupts, a little bit too large * a hammer. Therefore we need to catch potential overflows before * it's too late. */ -#if defined(_UNICOS) && defined(NV_MAX_10_EXP) - if (!negative && - (log10(value) + exponent) >= NV_MAX_10_EXP) - return NV_MAX; + +#if ((defined(VMS) && !defined(__IEEE_FP)) || defined(_UNICOS)) && defined(NV_MAX_10_EXP) + STMT_START { + NV exp_v = log10(value); + if (exponent >= NV_MAX_10_EXP || exponent + exp_v >= NV_MAX_10_EXP) + return NV_MAX; + if (exponent < 0) { + if (-(exponent + exp_v) >= NV_MAX_10_EXP) + return 0.0; + while (-exponent >= NV_MAX_10_EXP) { + /* combination does not overflow, but 10^(-exponent) does */ + value /= 10; + ++exponent; + } + } + } STMT_END; #endif + if (exponent < 0) { + negative = 1; + exponent = -exponent; + } for (bit = 1; exponent; bit <<= 1) { if (exponent & bit) { exponent ^= bit; result *= power; + /* Floating point exceptions are supposed to be turned off, + * but if we're obviously done, don't risk another iteration. + */ + if (exponent == 0) break; } - /* Floating point exceptions are supposed to be turned off. */ power *= power; } return negative ? value / result : value * result;