11da177e4SLinus Torvalds /* 21da177e4SLinus Torvalds * Linux/PA-RISC Project (http://www.parisc-linux.org/) 31da177e4SLinus Torvalds * 41da177e4SLinus Torvalds * Floating-point emulation code 51da177e4SLinus Torvalds * Copyright (C) 2001 Hewlett-Packard (Paul Bame) <bame@debian.org> 61da177e4SLinus Torvalds * 71da177e4SLinus Torvalds * This program is free software; you can redistribute it and/or modify 81da177e4SLinus Torvalds * it under the terms of the GNU General Public License as published by 91da177e4SLinus Torvalds * the Free Software Foundation; either version 2, or (at your option) 101da177e4SLinus Torvalds * any later version. 111da177e4SLinus Torvalds * 121da177e4SLinus Torvalds * This program is distributed in the hope that it will be useful, 131da177e4SLinus Torvalds * but WITHOUT ANY WARRANTY; without even the implied warranty of 141da177e4SLinus Torvalds * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 151da177e4SLinus Torvalds * GNU General Public License for more details. 161da177e4SLinus Torvalds * 171da177e4SLinus Torvalds * You should have received a copy of the GNU General Public License 181da177e4SLinus Torvalds * along with this program; if not, write to the Free Software 191da177e4SLinus Torvalds * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 201da177e4SLinus Torvalds */ 211da177e4SLinus Torvalds /* 221da177e4SLinus Torvalds * BEGIN_DESC 231da177e4SLinus Torvalds * 241da177e4SLinus Torvalds * File: 251da177e4SLinus Torvalds * @(#) pa/spmath/dfsqrt.c $Revision: 1.1 $ 261da177e4SLinus Torvalds * 271da177e4SLinus Torvalds * Purpose: 281da177e4SLinus Torvalds * Double Floating-point Square Root 291da177e4SLinus Torvalds * 301da177e4SLinus Torvalds * External Interfaces: 311da177e4SLinus Torvalds * dbl_fsqrt(srcptr,nullptr,dstptr,status) 321da177e4SLinus Torvalds * 331da177e4SLinus Torvalds * Internal Interfaces: 341da177e4SLinus Torvalds * 351da177e4SLinus Torvalds * Theory: 361da177e4SLinus Torvalds * <<please update with a overview of the operation of this file>> 371da177e4SLinus Torvalds * 381da177e4SLinus Torvalds * END_DESC 391da177e4SLinus Torvalds */ 401da177e4SLinus Torvalds 411da177e4SLinus Torvalds 421da177e4SLinus Torvalds #include "float.h" 431da177e4SLinus Torvalds #include "dbl_float.h" 441da177e4SLinus Torvalds 451da177e4SLinus Torvalds /* 461da177e4SLinus Torvalds * Double Floating-point Square Root 471da177e4SLinus Torvalds */ 481da177e4SLinus Torvalds 491da177e4SLinus Torvalds /*ARGSUSED*/ 501da177e4SLinus Torvalds unsigned int 511da177e4SLinus Torvalds dbl_fsqrt( 521da177e4SLinus Torvalds dbl_floating_point *srcptr, 531da177e4SLinus Torvalds unsigned int *nullptr, 541da177e4SLinus Torvalds dbl_floating_point *dstptr, 551da177e4SLinus Torvalds unsigned int *status) 561da177e4SLinus Torvalds { 571da177e4SLinus Torvalds register unsigned int srcp1, srcp2, resultp1, resultp2; 581da177e4SLinus Torvalds register unsigned int newbitp1, newbitp2, sump1, sump2; 591da177e4SLinus Torvalds register int src_exponent; 601da177e4SLinus Torvalds register boolean guardbit = FALSE, even_exponent; 611da177e4SLinus Torvalds 621da177e4SLinus Torvalds Dbl_copyfromptr(srcptr,srcp1,srcp2); 631da177e4SLinus Torvalds /* 641da177e4SLinus Torvalds * check source operand for NaN or infinity 651da177e4SLinus Torvalds */ 661da177e4SLinus Torvalds if ((src_exponent = Dbl_exponent(srcp1)) == DBL_INFINITY_EXPONENT) { 671da177e4SLinus Torvalds /* 681da177e4SLinus Torvalds * is signaling NaN? 691da177e4SLinus Torvalds */ 701da177e4SLinus Torvalds if (Dbl_isone_signaling(srcp1)) { 711da177e4SLinus Torvalds /* trap if INVALIDTRAP enabled */ 721da177e4SLinus Torvalds if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 731da177e4SLinus Torvalds /* make NaN quiet */ 741da177e4SLinus Torvalds Set_invalidflag(); 751da177e4SLinus Torvalds Dbl_set_quiet(srcp1); 761da177e4SLinus Torvalds } 771da177e4SLinus Torvalds /* 781da177e4SLinus Torvalds * Return quiet NaN or positive infinity. 797022672eSSimon Arlott * Fall through to negative test if negative infinity. 801da177e4SLinus Torvalds */ 811da177e4SLinus Torvalds if (Dbl_iszero_sign(srcp1) || 821da177e4SLinus Torvalds Dbl_isnotzero_mantissa(srcp1,srcp2)) { 831da177e4SLinus Torvalds Dbl_copytoptr(srcp1,srcp2,dstptr); 841da177e4SLinus Torvalds return(NOEXCEPTION); 851da177e4SLinus Torvalds } 861da177e4SLinus Torvalds } 871da177e4SLinus Torvalds 881da177e4SLinus Torvalds /* 891da177e4SLinus Torvalds * check for zero source operand 901da177e4SLinus Torvalds */ 911da177e4SLinus Torvalds if (Dbl_iszero_exponentmantissa(srcp1,srcp2)) { 921da177e4SLinus Torvalds Dbl_copytoptr(srcp1,srcp2,dstptr); 931da177e4SLinus Torvalds return(NOEXCEPTION); 941da177e4SLinus Torvalds } 951da177e4SLinus Torvalds 961da177e4SLinus Torvalds /* 971da177e4SLinus Torvalds * check for negative source operand 981da177e4SLinus Torvalds */ 991da177e4SLinus Torvalds if (Dbl_isone_sign(srcp1)) { 1001da177e4SLinus Torvalds /* trap if INVALIDTRAP enabled */ 1011da177e4SLinus Torvalds if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 1021da177e4SLinus Torvalds /* make NaN quiet */ 1031da177e4SLinus Torvalds Set_invalidflag(); 1041da177e4SLinus Torvalds Dbl_makequietnan(srcp1,srcp2); 1051da177e4SLinus Torvalds Dbl_copytoptr(srcp1,srcp2,dstptr); 1061da177e4SLinus Torvalds return(NOEXCEPTION); 1071da177e4SLinus Torvalds } 1081da177e4SLinus Torvalds 1091da177e4SLinus Torvalds /* 1101da177e4SLinus Torvalds * Generate result 1111da177e4SLinus Torvalds */ 1121da177e4SLinus Torvalds if (src_exponent > 0) { 1131da177e4SLinus Torvalds even_exponent = Dbl_hidden(srcp1); 1141da177e4SLinus Torvalds Dbl_clear_signexponent_set_hidden(srcp1); 1151da177e4SLinus Torvalds } 1161da177e4SLinus Torvalds else { 1171da177e4SLinus Torvalds /* normalize operand */ 1181da177e4SLinus Torvalds Dbl_clear_signexponent(srcp1); 1191da177e4SLinus Torvalds src_exponent++; 1201da177e4SLinus Torvalds Dbl_normalize(srcp1,srcp2,src_exponent); 1211da177e4SLinus Torvalds even_exponent = src_exponent & 1; 1221da177e4SLinus Torvalds } 1231da177e4SLinus Torvalds if (even_exponent) { 1241da177e4SLinus Torvalds /* exponent is even */ 1251da177e4SLinus Torvalds /* Add comment here. Explain why odd exponent needs correction */ 1261da177e4SLinus Torvalds Dbl_leftshiftby1(srcp1,srcp2); 1271da177e4SLinus Torvalds } 1281da177e4SLinus Torvalds /* 1291da177e4SLinus Torvalds * Add comment here. Explain following algorithm. 1301da177e4SLinus Torvalds * 1311da177e4SLinus Torvalds * Trust me, it works. 1321da177e4SLinus Torvalds * 1331da177e4SLinus Torvalds */ 1341da177e4SLinus Torvalds Dbl_setzero(resultp1,resultp2); 1351da177e4SLinus Torvalds Dbl_allp1(newbitp1) = 1 << (DBL_P - 32); 1361da177e4SLinus Torvalds Dbl_setzero_mantissap2(newbitp2); 1371da177e4SLinus Torvalds while (Dbl_isnotzero(newbitp1,newbitp2) && Dbl_isnotzero(srcp1,srcp2)) { 1381da177e4SLinus Torvalds Dbl_addition(resultp1,resultp2,newbitp1,newbitp2,sump1,sump2); 1391da177e4SLinus Torvalds if(Dbl_isnotgreaterthan(sump1,sump2,srcp1,srcp2)) { 1401da177e4SLinus Torvalds Dbl_leftshiftby1(newbitp1,newbitp2); 1411da177e4SLinus Torvalds /* update result */ 1421da177e4SLinus Torvalds Dbl_addition(resultp1,resultp2,newbitp1,newbitp2, 1431da177e4SLinus Torvalds resultp1,resultp2); 1441da177e4SLinus Torvalds Dbl_subtract(srcp1,srcp2,sump1,sump2,srcp1,srcp2); 1451da177e4SLinus Torvalds Dbl_rightshiftby2(newbitp1,newbitp2); 1461da177e4SLinus Torvalds } 1471da177e4SLinus Torvalds else { 1481da177e4SLinus Torvalds Dbl_rightshiftby1(newbitp1,newbitp2); 1491da177e4SLinus Torvalds } 1501da177e4SLinus Torvalds Dbl_leftshiftby1(srcp1,srcp2); 1511da177e4SLinus Torvalds } 1521da177e4SLinus Torvalds /* correct exponent for pre-shift */ 1531da177e4SLinus Torvalds if (even_exponent) { 1541da177e4SLinus Torvalds Dbl_rightshiftby1(resultp1,resultp2); 1551da177e4SLinus Torvalds } 1561da177e4SLinus Torvalds 1571da177e4SLinus Torvalds /* check for inexact */ 1581da177e4SLinus Torvalds if (Dbl_isnotzero(srcp1,srcp2)) { 1591da177e4SLinus Torvalds if (!even_exponent && Dbl_islessthan(resultp1,resultp2,srcp1,srcp2)) { 1601da177e4SLinus Torvalds Dbl_increment(resultp1,resultp2); 1611da177e4SLinus Torvalds } 1621da177e4SLinus Torvalds guardbit = Dbl_lowmantissap2(resultp2); 1631da177e4SLinus Torvalds Dbl_rightshiftby1(resultp1,resultp2); 1641da177e4SLinus Torvalds 1651da177e4SLinus Torvalds /* now round result */ 1661da177e4SLinus Torvalds switch (Rounding_mode()) { 1671da177e4SLinus Torvalds case ROUNDPLUS: 1681da177e4SLinus Torvalds Dbl_increment(resultp1,resultp2); 1691da177e4SLinus Torvalds break; 1701da177e4SLinus Torvalds case ROUNDNEAREST: 1711da177e4SLinus Torvalds /* stickybit is always true, so guardbit 1721da177e4SLinus Torvalds * is enough to determine rounding */ 1731da177e4SLinus Torvalds if (guardbit) { 1741da177e4SLinus Torvalds Dbl_increment(resultp1,resultp2); 1751da177e4SLinus Torvalds } 1761da177e4SLinus Torvalds break; 1771da177e4SLinus Torvalds } 1781da177e4SLinus Torvalds /* increment result exponent by 1 if mantissa overflowed */ 1791da177e4SLinus Torvalds if (Dbl_isone_hiddenoverflow(resultp1)) src_exponent+=2; 1801da177e4SLinus Torvalds 1811da177e4SLinus Torvalds if (Is_inexacttrap_enabled()) { 1821da177e4SLinus Torvalds Dbl_set_exponent(resultp1, 1831da177e4SLinus Torvalds ((src_exponent-DBL_BIAS)>>1)+DBL_BIAS); 1841da177e4SLinus Torvalds Dbl_copytoptr(resultp1,resultp2,dstptr); 1851da177e4SLinus Torvalds return(INEXACTEXCEPTION); 1861da177e4SLinus Torvalds } 1871da177e4SLinus Torvalds else Set_inexactflag(); 1881da177e4SLinus Torvalds } 1891da177e4SLinus Torvalds else { 1901da177e4SLinus Torvalds Dbl_rightshiftby1(resultp1,resultp2); 1911da177e4SLinus Torvalds } 1921da177e4SLinus Torvalds Dbl_set_exponent(resultp1,((src_exponent-DBL_BIAS)>>1)+DBL_BIAS); 1931da177e4SLinus Torvalds Dbl_copytoptr(resultp1,resultp2,dstptr); 1941da177e4SLinus Torvalds return(NOEXCEPTION); 1951da177e4SLinus Torvalds } 196