mirror of
https://git.kernel.org/pub/scm/linux/kernel/git/stable/linux.git
synced 2024-09-12 13:55:32 +00:00
![Thomas Gleixner](/assets/img/avatar_default.png)
Based on 1 normalized pattern(s): this program is free software you can redistribute it and or modify it under the terms of the gnu general public license as published by the free software foundation either version 2 or at your option any later version this program is distributed in the hope that it will be useful but without any warranty without even the implied warranty of merchantability or fitness for a particular purpose see the gnu general public license for more details you should have received a copy of the gnu general public license along with this program if not write to the free software foundation inc 59 temple place suite 330 boston ma 02111 1307 usa extracted by the scancode license scanner the SPDX license identifier GPL-2.0-or-later has been chosen to replace the boilerplate/reference in 42 file(s). Signed-off-by: Thomas Gleixner <tglx@linutronix.de> Reviewed-by: Richard Fontana <rfontana@redhat.com> Reviewed-by: Allison Randal <allison@lohutok.net> Reviewed-by: Kate Stewart <kstewart@linuxfoundation.org> Cc: linux-spdx@vger.kernel.org Link: https://lkml.kernel.org/r/20190524100845.259718220@linutronix.de Signed-off-by: Greg Kroah-Hartman <gregkh@linuxfoundation.org>
182 lines
4.7 KiB
C
182 lines
4.7 KiB
C
// SPDX-License-Identifier: GPL-2.0-or-later
|
|
/*
|
|
* Linux/PA-RISC Project (http://www.parisc-linux.org/)
|
|
*
|
|
* Floating-point emulation code
|
|
* Copyright (C) 2001 Hewlett-Packard (Paul Bame) <bame@debian.org>
|
|
*/
|
|
/*
|
|
* BEGIN_DESC
|
|
*
|
|
* File:
|
|
* @(#) pa/spmath/dfsqrt.c $Revision: 1.1 $
|
|
*
|
|
* Purpose:
|
|
* Double Floating-point Square Root
|
|
*
|
|
* External Interfaces:
|
|
* dbl_fsqrt(srcptr,nullptr,dstptr,status)
|
|
*
|
|
* Internal Interfaces:
|
|
*
|
|
* Theory:
|
|
* <<please update with a overview of the operation of this file>>
|
|
*
|
|
* END_DESC
|
|
*/
|
|
|
|
|
|
#include "float.h"
|
|
#include "dbl_float.h"
|
|
|
|
/*
|
|
* Double Floating-point Square Root
|
|
*/
|
|
|
|
/*ARGSUSED*/
|
|
unsigned int
|
|
dbl_fsqrt(
|
|
dbl_floating_point *srcptr,
|
|
unsigned int *nullptr,
|
|
dbl_floating_point *dstptr,
|
|
unsigned int *status)
|
|
{
|
|
register unsigned int srcp1, srcp2, resultp1, resultp2;
|
|
register unsigned int newbitp1, newbitp2, sump1, sump2;
|
|
register int src_exponent;
|
|
register boolean guardbit = FALSE, even_exponent;
|
|
|
|
Dbl_copyfromptr(srcptr,srcp1,srcp2);
|
|
/*
|
|
* check source operand for NaN or infinity
|
|
*/
|
|
if ((src_exponent = Dbl_exponent(srcp1)) == DBL_INFINITY_EXPONENT) {
|
|
/*
|
|
* is signaling NaN?
|
|
*/
|
|
if (Dbl_isone_signaling(srcp1)) {
|
|
/* trap if INVALIDTRAP enabled */
|
|
if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
|
|
/* make NaN quiet */
|
|
Set_invalidflag();
|
|
Dbl_set_quiet(srcp1);
|
|
}
|
|
/*
|
|
* Return quiet NaN or positive infinity.
|
|
* Fall through to negative test if negative infinity.
|
|
*/
|
|
if (Dbl_iszero_sign(srcp1) ||
|
|
Dbl_isnotzero_mantissa(srcp1,srcp2)) {
|
|
Dbl_copytoptr(srcp1,srcp2,dstptr);
|
|
return(NOEXCEPTION);
|
|
}
|
|
}
|
|
|
|
/*
|
|
* check for zero source operand
|
|
*/
|
|
if (Dbl_iszero_exponentmantissa(srcp1,srcp2)) {
|
|
Dbl_copytoptr(srcp1,srcp2,dstptr);
|
|
return(NOEXCEPTION);
|
|
}
|
|
|
|
/*
|
|
* check for negative source operand
|
|
*/
|
|
if (Dbl_isone_sign(srcp1)) {
|
|
/* trap if INVALIDTRAP enabled */
|
|
if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
|
|
/* make NaN quiet */
|
|
Set_invalidflag();
|
|
Dbl_makequietnan(srcp1,srcp2);
|
|
Dbl_copytoptr(srcp1,srcp2,dstptr);
|
|
return(NOEXCEPTION);
|
|
}
|
|
|
|
/*
|
|
* Generate result
|
|
*/
|
|
if (src_exponent > 0) {
|
|
even_exponent = Dbl_hidden(srcp1);
|
|
Dbl_clear_signexponent_set_hidden(srcp1);
|
|
}
|
|
else {
|
|
/* normalize operand */
|
|
Dbl_clear_signexponent(srcp1);
|
|
src_exponent++;
|
|
Dbl_normalize(srcp1,srcp2,src_exponent);
|
|
even_exponent = src_exponent & 1;
|
|
}
|
|
if (even_exponent) {
|
|
/* exponent is even */
|
|
/* Add comment here. Explain why odd exponent needs correction */
|
|
Dbl_leftshiftby1(srcp1,srcp2);
|
|
}
|
|
/*
|
|
* Add comment here. Explain following algorithm.
|
|
*
|
|
* Trust me, it works.
|
|
*
|
|
*/
|
|
Dbl_setzero(resultp1,resultp2);
|
|
Dbl_allp1(newbitp1) = 1 << (DBL_P - 32);
|
|
Dbl_setzero_mantissap2(newbitp2);
|
|
while (Dbl_isnotzero(newbitp1,newbitp2) && Dbl_isnotzero(srcp1,srcp2)) {
|
|
Dbl_addition(resultp1,resultp2,newbitp1,newbitp2,sump1,sump2);
|
|
if(Dbl_isnotgreaterthan(sump1,sump2,srcp1,srcp2)) {
|
|
Dbl_leftshiftby1(newbitp1,newbitp2);
|
|
/* update result */
|
|
Dbl_addition(resultp1,resultp2,newbitp1,newbitp2,
|
|
resultp1,resultp2);
|
|
Dbl_subtract(srcp1,srcp2,sump1,sump2,srcp1,srcp2);
|
|
Dbl_rightshiftby2(newbitp1,newbitp2);
|
|
}
|
|
else {
|
|
Dbl_rightshiftby1(newbitp1,newbitp2);
|
|
}
|
|
Dbl_leftshiftby1(srcp1,srcp2);
|
|
}
|
|
/* correct exponent for pre-shift */
|
|
if (even_exponent) {
|
|
Dbl_rightshiftby1(resultp1,resultp2);
|
|
}
|
|
|
|
/* check for inexact */
|
|
if (Dbl_isnotzero(srcp1,srcp2)) {
|
|
if (!even_exponent && Dbl_islessthan(resultp1,resultp2,srcp1,srcp2)) {
|
|
Dbl_increment(resultp1,resultp2);
|
|
}
|
|
guardbit = Dbl_lowmantissap2(resultp2);
|
|
Dbl_rightshiftby1(resultp1,resultp2);
|
|
|
|
/* now round result */
|
|
switch (Rounding_mode()) {
|
|
case ROUNDPLUS:
|
|
Dbl_increment(resultp1,resultp2);
|
|
break;
|
|
case ROUNDNEAREST:
|
|
/* stickybit is always true, so guardbit
|
|
* is enough to determine rounding */
|
|
if (guardbit) {
|
|
Dbl_increment(resultp1,resultp2);
|
|
}
|
|
break;
|
|
}
|
|
/* increment result exponent by 1 if mantissa overflowed */
|
|
if (Dbl_isone_hiddenoverflow(resultp1)) src_exponent+=2;
|
|
|
|
if (Is_inexacttrap_enabled()) {
|
|
Dbl_set_exponent(resultp1,
|
|
((src_exponent-DBL_BIAS)>>1)+DBL_BIAS);
|
|
Dbl_copytoptr(resultp1,resultp2,dstptr);
|
|
return(INEXACTEXCEPTION);
|
|
}
|
|
else Set_inexactflag();
|
|
}
|
|
else {
|
|
Dbl_rightshiftby1(resultp1,resultp2);
|
|
}
|
|
Dbl_set_exponent(resultp1,((src_exponent-DBL_BIAS)>>1)+DBL_BIAS);
|
|
Dbl_copytoptr(resultp1,resultp2,dstptr);
|
|
return(NOEXCEPTION);
|
|
}
|