include/impl/sqrt.hpp
/**
 * @brief This ieee754 conforming sqrt rounds correctly.
 *
 * @param x
 * @return constexpr double
 */
@brief This ieee754 conforming sqrt rounds correctly.
@param x
@return constexpr double
@brief This ieee754 conforming sqrt rounds correctly . @param x @return constexpr double
constexpr double __ieee754_sqrt(double x) {
    if (isnan(x) || isinf(x)) { return x * x + x; }
    if (x == 0) {
        return x; 
    } else if (x < 0) {
        return (x - x) / (x - x);
    }
    constexpr int sign = (int)0x80000000;
    constexpr double ONE = 1.;
    constexpr double TINY = 1.e-300;
    int ix0 = f64_high(x);
    unsigned ix1 = f64_low(x);
    int m = ix0 >> 20;
    if (m == 0) { 
        while (ix0 == 0) {
            m -= 21;
            ix0 |= (ix1 >> 11);
            ix1 <<= 21;
        }
        int i = 0;
        for (i = 0; (ix0 & 0x00100000) == 0; i++) ix0 <<= 1;
        m -= i - 1;
        ix0 |= (ix1 >> (32 - i));
        ix1 <<= i;
    }
    m -= 1023;
    ix0 = (ix0 & 0x000fffff) | 0x00100000;
    if (m & 1) { 
        ix0 += ix0 + ((ix1 & sign) >> 31);
        ix1 += ix1;
    }
    m >>= 1;
    ix0 += ix0 + ((ix1 & sign) >> 31);
    ix1 += ix1;
    int q{}, s0{}, t{};
    unsigned q1{}, s1{}, t1{};
    q = q1 = s0 = s1 = 0;    
    unsigned r = 0x00200000; 
    while (r != 0) {
        t = s0 + r;
        if (t <= ix0) {
            s0 = t + r;
            ix0 -= t;
            q += r;
        }
        ix0 += ix0 + ((ix1 & sign) >> 31);
        ix1 += ix1;
        r >>= 1;
    }
    r = sign;
    while (r != 0) {
        t1 = s1 + r;
        t = s0;
        if ((t < ix0) || ((t == ix0) && (t1 <= ix1))) {
            s1 = t1 + r;
            if (((t1 & sign) == sign) && (s1 & sign) == 0) s0 += 1;
            ix0 -= t;
            if (ix1 < t1) ix0 -= 1;
            ix1 -= t1;
            q1 += r;
        }
        ix0 += ix0 + ((ix1 & sign) >> 31);
        ix1 += ix1;
        r >>= 1;
    }
    if ((ix0 | ix1) != 0) {
        double z = ONE - TINY; 
        if (z >= ONE) {
            z = ONE + TINY;
            if (q1 == (unsigned)0xffffffff) {
                q1 = 0;
                q += 1;
            } else if (z > ONE) {
                if (q1 == (unsigned)0xfffffffe) q += 1;
                q1 += 2;
            } else
                q1 += (q1 & 1);
        }
    }
    ix0 = (q >> 1) + 0x3fe00000;
    ix1 = q1 >> 1;
    if ((q & 1) == 1) ix1 |= sign;
    ix0 += (m << 20);
    return f64_from_bits(ix0, ix1);
}
constexpr double __ieee754_sqrt ( double x ) { if ( isnan ( x ) || isinf ( x ) ) { return x * x + x ; } if ( x == 0 ) { return x ; } else if ( x < 0 ) { return ( x - x ) / ( x - x ) ; } constexpr int sign = ( int ) 0x80000000 ; constexpr double ONE = 1. ; constexpr double TINY = 1.e-300 ; int ix0 = f64_high ( x ) ; unsigned ix1 = f64_low ( x ) ; int m = ix0 >> 20 ; if ( m == 0 ) { while ( ix0 == 0 ) { m -= 21 ; ix0 |= ( ix1 >> 11 ) ; ix1 <<= 21 ; } int i = 0 ; for ( i = 0 ; ( ix0 & 0x00100000 ) == 0 ; i ++ ) ix0 <<= 1 ; m -= i - 1 ; ix0 |= ( ix1 >> ( 32 - i ) ) ; ix1 <<= i ; } m -= 1023 ; ix0 = ( ix0 & 0x000fffff ) | 0x00100000 ; if ( m & 1 ) { ix0 += ix0 + ( ( ix1 & sign ) >> 31 ) ; ix1 += ix1 ; } m >>= 1 ; ix0 += ix0 + ( ( ix1 & sign ) >> 31 ) ; ix1 += ix1 ; int q { } , s0 { } , t { } ; unsigned q1 { } , s1 { } , t1 { } ; q = q1 = s0 = s1 = 0 ; unsigned r = 0x00200000 ; while ( r != 0 ) { t = s0 + r ; if ( t <= ix0 ) { s0 = t + r ; ix0 -= t ; q += r ; } ix0 += ix0 + ( ( ix1 & sign ) >> 31 ) ; ix1 += ix1 ; r >>= 1 ; } r = sign ; while ( r != 0 ) { t1 = s1 + r ; t = s0 ; if ( ( t < ix0 ) || ( ( t == ix0 ) && ( t1 <= ix1 ) ) ) { s1 = t1 + r ; if ( ( ( t1 & sign ) == sign ) && ( s1 & sign ) == 0 ) s0 += 1 ; ix0 -= t ; if ( ix1 < t1 ) ix0 -= 1 ; ix1 -= t1 ; q1 += r ; } ix0 += ix0 + ( ( ix1 & sign ) >> 31 ) ; ix1 += ix1 ; r >>= 1 ; } if ( ( ix0 | ix1 ) != 0 ) { double z = ONE - TINY ; if ( z >= ONE ) { z = ONE + TINY ; if ( q1 == ( unsigned ) 0xffffffff ) { q1 = 0 ; q += 1 ; } else if ( z > ONE ) { if ( q1 == ( unsigned ) 0xfffffffe ) q += 1 ; q1 += 2 ; } else q1 += ( q1 & 1 ) ; } } ix0 = ( q >> 1 ) + 0x3fe00000 ; ix1 = q1 >> 1 ; if ( ( q & 1 ) == 1 ) ix1 |= sign ; ix0 += ( m << 20 ) ; return f64_from_bits ( ix0 , ix1 ) ; }
@brief This ieee754 conforming sqrt rounds correctly.
@brief This ieee754 conforming sqrt rounds correctly .
/* $sqrt(\pm 0) = \pm 0$ */ /* subnormal x */ /* odd m, double x to make it even */ /* generate sqrt(x) bit by bit */ /* [q,q1] = sqrt(x) */ /* r = moving bit from right to left */ /* use floating add to find out rounding direction */ /* trigger inexact flag */