![[BACK]](/icons/back.gif){kind=icon}
Return to dfsqrt.c CVS log ![[TXT]](/icons/text.gif){kind=icon}
![[DIR]](/icons/dir.gif){kind=icon}
Up to [local] / sys / arch / hppa / spmath|
Return to dfsqrt.c CVS log | Up to [local] / sys / arch / hppa / spmath |
1.1 nbrk 1: /* $OpenBSD: dfsqrt.c,v 1.7 2003/04/10 17:27:58 mickey Exp $ */
2: /*
3: (c) Copyright 1986 HEWLETT-PACKARD COMPANY
4: To anyone who acknowledges that this file is provided "AS IS"
5: without any express or implied warranty:
6: permission to use, copy, modify, and distribute this file
7: for any purpose is hereby granted without fee, provided that
8: the above copyright notice and this notice appears in all
9: copies, and that the name of Hewlett-Packard Company not be
10: used in advertising or publicity pertaining to distribution
11: of the software without specific, written prior permission.
12: Hewlett-Packard Company makes no representations about the
13: suitability of this software for any purpose.
14: */
15: /* @(#)dfsqrt.c: Revision: 1.9.88.1 Date: 93/12/07 15:05:46 */
16:
17: #include "float.h"
18: #include "dbl_float.h"
19:
20: /*
21: * Double Floating-point Square Root
22: */
23:
24: /*ARGSUSED*/
25: int
26: dbl_fsqrt(srcptr, null, dstptr, status)
27: dbl_floating_point *srcptr, *null, *dstptr;
28: unsigned int *status;
29: {
30: register unsigned int srcp1, srcp2, resultp1, resultp2;
31: register unsigned int newbitp1, newbitp2, sump1, sump2;
32: register int src_exponent;
33: register int guardbit = FALSE, even_exponent;
34:
35: Dbl_copyfromptr(srcptr,srcp1,srcp2);
36: /*
37: * check source operand for NaN or infinity
38: */
39: if ((src_exponent = Dbl_exponent(srcp1)) == DBL_INFINITY_EXPONENT) {
40: /*
41: * is signaling NaN?
42: */
43: if (Dbl_isone_signaling(srcp1)) {
44: /* trap if INVALIDTRAP enabled */
45: if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
46: /* make NaN quiet */
47: Set_invalidflag();
48: Dbl_set_quiet(srcp1);
49: }
50: /*
51: * Return quiet NaN or positive infinity.
52: * Fall thru to negative test if negative infinity.
53: */
54: if (Dbl_iszero_sign(srcp1) ||
55: Dbl_isnotzero_mantissa(srcp1,srcp2)) {
56: Dbl_copytoptr(srcp1,srcp2,dstptr);
57: return(NOEXCEPTION);
58: }
59: }
60:
61: /*
62: * check for zero source operand
63: */
64: if (Dbl_iszero_exponentmantissa(srcp1,srcp2)) {
65: Dbl_copytoptr(srcp1,srcp2,dstptr);
66: return(NOEXCEPTION);
67: }
68:
69: /*
70: * check for negative source operand
71: */
72: if (Dbl_isone_sign(srcp1)) {
73: /* trap if INVALIDTRAP enabled */
74: if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
75: /* make NaN quiet */
76: Set_invalidflag();
77: Dbl_makequietnan(srcp1,srcp2);
78: Dbl_copytoptr(srcp1,srcp2,dstptr);
79: return(NOEXCEPTION);
80: }
81:
82: /*
83: * Generate result
84: */
85: if (src_exponent > 0) {
86: even_exponent = Dbl_hidden(srcp1);
87: Dbl_clear_signexponent_set_hidden(srcp1);
88: }
89: else {
90: /* normalize operand */
91: Dbl_clear_signexponent(srcp1);
92: src_exponent++;
93: Dbl_normalize(srcp1,srcp2,src_exponent);
94: even_exponent = src_exponent & 1;
95: }
96: if (even_exponent) {
97: /* exponent is even */
98: /* Add comment here. Explain why odd exponent needs correction */
99: Dbl_leftshiftby1(srcp1,srcp2);
100: }
101: /*
102: * Add comment here. Explain following algorithm.
103: *
104: * Trust me, it works.
105: *
106: */
107: Dbl_setzero(resultp1,resultp2);
108: Dbl_allp1(newbitp1) = 1 << (DBL_P - 32);
109: Dbl_setzero_mantissap2(newbitp2);
110: while (Dbl_isnotzero(newbitp1,newbitp2) && Dbl_isnotzero(srcp1,srcp2)) {
111: Dbl_addition(resultp1,resultp2,newbitp1,newbitp2,sump1,sump2);
112: if(Dbl_isnotgreaterthan(sump1,sump2,srcp1,srcp2)) {
113: Dbl_leftshiftby1(newbitp1,newbitp2);
114: /* update result */
115: Dbl_addition(resultp1,resultp2,newbitp1,newbitp2,
116: resultp1,resultp2);
117: Dbl_subtract(srcp1,srcp2,sump1,sump2,srcp1,srcp2);
118: Dbl_rightshiftby2(newbitp1,newbitp2);
119: }
120: else {
121: Dbl_rightshiftby1(newbitp1,newbitp2);
122: }
123: Dbl_leftshiftby1(srcp1,srcp2);
124: }
125: /* correct exponent for pre-shift */
126: if (even_exponent) {
127: Dbl_rightshiftby1(resultp1,resultp2);
128: }
129:
130: /* check for inexact */
131: if (Dbl_isnotzero(srcp1,srcp2)) {
132: if (!even_exponent & Dbl_islessthan(resultp1,resultp2,srcp1,srcp2)) {
133: Dbl_increment(resultp1,resultp2);
134: }
135: guardbit = Dbl_lowmantissap2(resultp2);
136: Dbl_rightshiftby1(resultp1,resultp2);
137:
138: /* now round result */
139: switch (Rounding_mode()) {
140: case ROUNDPLUS:
141: Dbl_increment(resultp1,resultp2);
142: break;
143: case ROUNDNEAREST:
144: /* stickybit is always true, so guardbit
145: * is enough to determine rounding */
146: if (guardbit) {
147: Dbl_increment(resultp1,resultp2);
148: }
149: break;
150: }
151: /* increment result exponent by 1 if mantissa overflowed */
152: if (Dbl_isone_hiddenoverflow(resultp1)) src_exponent+=2;
153:
154: if (Is_inexacttrap_enabled()) {
155: Dbl_set_exponent(resultp1,
156: ((src_exponent-DBL_BIAS)>>1)+DBL_BIAS);
157: Dbl_copytoptr(resultp1,resultp2,dstptr);
158: return(INEXACTEXCEPTION);
159: }
160: else Set_inexactflag();
161: }
162: else {
163: Dbl_rightshiftby1(resultp1,resultp2);
164: }
165: Dbl_set_exponent(resultp1,((src_exponent-DBL_BIAS)>>1)+DBL_BIAS);
166: Dbl_copytoptr(resultp1,resultp2,dstptr);
167: return(NOEXCEPTION);
168: }