1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
|
/***
*sqrt4000.c - square root for the R4000
*
* Copyright (c) 1991-1991, Microsoft Corporation. All rights reserved.
*
*Purpose:
*
*Revision History:
* 3-31-92 GDP written
*******************************************************************************/
#ifdef R4000
#include <math.h>
#include <trans.h>
static double scale = 4503599627370496.0; // 2^52
static double lgscale = 26; // log2(2^52) / 2
#define SWMASK 0x78
#define INEXACT_MASK 0x4
/***
*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,cw,sw;
double result,savedx;
int scaled = 0;
// mask exceptions, keep user's rounding mode
savedcw = _ctrlfp(ICW, IMCW & ~IMCW_RC);
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);
}
savedx = x;
if (IS_D_DENORM(x)) {
x *= scale;
scaled = 1;
}
sw = _clrfp();
result = _fsqrt(x);
cw = _get_fsr();
_set_fsr(cw & ~SWMASK | sw & SWMASK); // restore all but inexact
if (scaled) {
result = _add_exp(result, -lgscale);
}
if ((cw & INEXACT_MASK) == INEXACT_MASK) {
return _except1(FP_P,OP_SQRT,savedx,result,savedcw);
}
RETURN(savedcw,result);
}
#endif
|
