Annotation of sys/arch/hppa/spmath/sfsqrt.c, Revision 1.1
1.1 ! nbrk 1: /* $OpenBSD: sfsqrt.c,v 1.7 2003/04/10 17:27:59 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: /* @(#)sfsqrt.c: Revision: 1.9.88.1 Date: 93/12/07 15:07:13 */
! 16:
! 17: #include "float.h"
! 18: #include "sgl_float.h"
! 19:
! 20: /*
! 21: * Single Floating-point Square Root
! 22: */
! 23:
! 24: /*ARGSUSED*/
! 25: int
! 26: sgl_fsqrt(srcptr, null, dstptr, status)
! 27: sgl_floating_point *srcptr, *null, *dstptr;
! 28: unsigned int *status;
! 29: {
! 30: register unsigned int src, result;
! 31: register int src_exponent, newbit, sum;
! 32: register int guardbit = FALSE, even_exponent;
! 33:
! 34: src = *srcptr;
! 35: /*
! 36: * check source operand for NaN or infinity
! 37: */
! 38: if ((src_exponent = Sgl_exponent(src)) == SGL_INFINITY_EXPONENT) {
! 39: /*
! 40: * is signaling NaN?
! 41: */
! 42: if (Sgl_isone_signaling(src)) {
! 43: /* trap if INVALIDTRAP enabled */
! 44: if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
! 45: /* make NaN quiet */
! 46: Set_invalidflag();
! 47: Sgl_set_quiet(src);
! 48: }
! 49: /*
! 50: * Return quiet NaN or positive infinity.
! 51: * Fall thru to negative test if negative infinity.
! 52: */
! 53: if (Sgl_iszero_sign(src) || Sgl_isnotzero_mantissa(src)) {
! 54: *dstptr = src;
! 55: return(NOEXCEPTION);
! 56: }
! 57: }
! 58:
! 59: /*
! 60: * check for zero source operand
! 61: */
! 62: if (Sgl_iszero_exponentmantissa(src)) {
! 63: *dstptr = src;
! 64: return(NOEXCEPTION);
! 65: }
! 66:
! 67: /*
! 68: * check for negative source operand
! 69: */
! 70: if (Sgl_isone_sign(src)) {
! 71: /* trap if INVALIDTRAP enabled */
! 72: if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
! 73: /* make NaN quiet */
! 74: Set_invalidflag();
! 75: Sgl_makequietnan(src);
! 76: *dstptr = src;
! 77: return(NOEXCEPTION);
! 78: }
! 79:
! 80: /*
! 81: * Generate result
! 82: */
! 83: if (src_exponent > 0) {
! 84: even_exponent = Sgl_hidden(src);
! 85: Sgl_clear_signexponent_set_hidden(src);
! 86: }
! 87: else {
! 88: /* normalize operand */
! 89: Sgl_clear_signexponent(src);
! 90: src_exponent++;
! 91: Sgl_normalize(src,src_exponent);
! 92: even_exponent = src_exponent & 1;
! 93: }
! 94: if (even_exponent) {
! 95: /* exponent is even */
! 96: /* Add comment here. Explain why odd exponent needs correction */
! 97: Sgl_leftshiftby1(src);
! 98: }
! 99: /*
! 100: * Add comment here. Explain following algorithm.
! 101: *
! 102: * Trust me, it works.
! 103: *
! 104: */
! 105: Sgl_setzero(result);
! 106: newbit = 1 << SGL_P;
! 107: while (newbit && Sgl_isnotzero(src)) {
! 108: Sgl_addition(result,newbit,sum);
! 109: if(sum <= Sgl_all(src)) {
! 110: /* update result */
! 111: Sgl_addition(result,(newbit<<1),result);
! 112: Sgl_subtract(src,sum,src);
! 113: }
! 114: Sgl_rightshiftby1(newbit);
! 115: Sgl_leftshiftby1(src);
! 116: }
! 117: /* correct exponent for pre-shift */
! 118: if (even_exponent) {
! 119: Sgl_rightshiftby1(result);
! 120: }
! 121:
! 122: /* check for inexact */
! 123: if (Sgl_isnotzero(src)) {
! 124: if (!even_exponent & Sgl_islessthan(result,src))
! 125: Sgl_increment(result);
! 126: guardbit = Sgl_lowmantissa(result);
! 127: Sgl_rightshiftby1(result);
! 128:
! 129: /* now round result */
! 130: switch (Rounding_mode()) {
! 131: case ROUNDPLUS:
! 132: Sgl_increment(result);
! 133: break;
! 134: case ROUNDNEAREST:
! 135: /* stickybit is always true, so guardbit
! 136: * is enough to determine rounding */
! 137: if (guardbit) {
! 138: Sgl_increment(result);
! 139: }
! 140: break;
! 141: }
! 142: /* increment result exponent by 1 if mantissa overflowed */
! 143: if (Sgl_isone_hiddenoverflow(result)) src_exponent+=2;
! 144:
! 145: if (Is_inexacttrap_enabled()) {
! 146: Sgl_set_exponent(result,
! 147: ((src_exponent-SGL_BIAS)>>1)+SGL_BIAS);
! 148: *dstptr = result;
! 149: return(INEXACTEXCEPTION);
! 150: }
! 151: else Set_inexactflag();
! 152: }
! 153: else {
! 154: Sgl_rightshiftby1(result);
! 155: }
! 156: Sgl_set_exponent(result,((src_exponent-SGL_BIAS)>>1)+SGL_BIAS);
! 157: *dstptr = result;
! 158: return(NOEXCEPTION);
! 159: }
CVSweb
![[BACK]](/icons/back.gif){kind=icon}
Return to sfsqrt.c CVS log ![[TXT]](/icons/text.gif){kind=icon}
![[DIR]](/icons/dir.gif){kind=icon}
Up to [local] / sys / arch / hppa / spmath