diff options
Diffstat (limited to 'private/fp32/tran/alpha/sqrt.c')
| -rw-r--r-- | private/fp32/tran/alpha/sqrt.c | 192 |
1 files changed, 192 insertions, 0 deletions
diff --git a/private/fp32/tran/alpha/sqrt.c b/private/fp32/tran/alpha/sqrt.c new file mode 100644 index 000000000..8cff2f7bf --- /dev/null +++ b/private/fp32/tran/alpha/sqrt.c @@ -0,0 +1,192 @@ +/*** +*sqrt.c - square root +* +* Copyright (c) 1991-1991, Microsoft Corporation. All rights reserved. +* +*Purpose: +* +*Revision History: +* 8-15-91 GDP written +* 1-29-91 GDP Kahan's algorithm for final rounding +* 3-11-92 GDP new interval and initial approximation +* 2-22-93 TVB Adapted for Alpha AXP +* +*******************************************************************************/ +#ifdef _ALPHA_ + +// +// This Alpha-specific version of sqrt is identical to portable version +// except for instances where chopped floating point operations are performed +// implicitly by setting the rounding mode to chopped. For Alpha, these +// instances are replaced with calls to assembler routines that explicitly +// perform the chopped operations. This is necessary because there is no +// compiler support yet for generating floating point instructions using the +// Alpha dynamic rounding mode. +// + +#include <math.h> +#include <trans.h> + +extern double _addtc(double x, double y); +extern double _divtc(double x, double y); +extern double _multc(double x, double y); +extern double _subtc(double x, double y); + +// +// Coefficients for initial approximation (Hart & al) +// + +static double p00 = .2592768763e+0; +static double p01 = .1052021187e+1; +static double p02 = -.3163221431e+0; + + +/*** +*double sqrt(double x) - square root +* +*Purpose: +* Compute the square root of a number. +* This function should be provided by the underlying +* hardware (IEEE spec). +*Entry: +* +*Exit: +* +*Exceptions: +* I P +*******************************************************************************/ +double sqrt(double x) +{ + unsigned int savedcw, sw; + double result,t; + unsigned int stat,rc; + + savedcw = _ctrlfp(ICW, IMCW); + + if (IS_D_SPECIAL(x)){ + switch (_sptype(x)) { + case T_PINF: + RETURN(savedcw, x); + case T_QNAN: + return _handle_qnan1(OP_SQRT, x, savedcw); + case T_SNAN: + return _except1(FP_I,OP_SQRT,x,QNAN_SQRT,savedcw); + } + /* -INF will be handled in the x<0 case */ + } + if (x < 0.0) { + return _except1(FP_I, OP_SQRT, x, QNAN_SQRT,savedcw); + } + + if (x == 0.0) { + RETURN (savedcw, x); + } + + + result = _fsqrt(x); + + + // + // Kahan's algorithm + // + + t = _divtc(x, result); + + // + // Multiply back to see if division was exact. + // Compare using subtraction to avoid invalid exceptions. + // + + if (_subtc(x, _multc(t, result)) == 0.0) { + // exact + if (t == result) { + RETURN(savedcw, result); + } + else { + // t = t-1 + if (*D_LO(t) == 0) { + (*D_HI(t)) --; + } + (*D_LO(t)) --; + } + + } + + rc = savedcw & IMCW_RC; + if (rc == IRC_UP || rc == IRC_NEAR) { + // t = t+1 + (*D_LO(t)) ++; + if (*D_LO(t) == 0) { + (*D_HI(t)) ++; + } + if (rc == IRC_UP) { + // y = y+1 + (*D_LO(t)) ++; + if (*D_LO(t) == 0) { + (*D_HI(t)) ++; + } + } + } + + result = _multc(0.5, _addtc(t, result)); + + _set_statfp(ISW_INEXACT); // update status word + RETURN_INEXACT1(OP_SQRT, x, result, savedcw); +} + + + +/*** +* _fsqrt - non IEEE conforming square root +* +*Purpose: +* compute a square root of a normal number without performing +* IEEE rounding. The argument is a finite number (no NaN or INF) +* +*Entry: +* +*Exit: +* +*Exceptions: +* +*******************************************************************************/ + +double _fsqrt(double x) +{ + double f,y,result; + int n; + + f = _decomp(x,&n); + + if (n & 0x1) { + // n is odd + n++; + f *= 0.5; + } + + // + // approximation for sqrt in the interval [.25, 1] + // (Computer Approximationsn, Hart & al.) + // gives more than 7 bits of accuracy + // + + y = p00 + f * (p01 + f * p02); + + y += f / y; + y *= 0.5; + + y += f / y; + y *= 0.5; + + y += f / y; + y *= 0.5; + + n >>= 1; + result = _add_exp(y,n); + + return result; +} + + + +#endif |