Annotation of sys/arch/hppa/spmath/dfsqrt.c, Revision 1.1
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: }
CVSweb
![[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