Changeset ca113cf in mainline for uspace/lib/math/generic/sqrt.c


Ignore:
Timestamp:
2021-06-18T16:06:10Z (3 years ago)
Author:
Maurizio Lombardi <mlombard@…>
Children:
f23dbf4
Parents:
197d59d
git-author:
Maurizio Lombardi <mlombard@…> (2018-11-25 16:51:03)
git-committer:
Maurizio Lombardi <mlombard@…> (2021-06-18 16:06:10)
Message:

math: sync sqrt() to FreeBSD 11.2

File:
1 edited

Legend:

Unmodified
Added
Removed
  • uspace/lib/math/generic/sqrt.c

    r197d59d rca113cf  
    9292#include "internal.h"
    9393
    94 static const double tiny = 1.0e-300;
     94static  const double    one     = 1.0, tiny=1.0e-300;
    9595
    9696double sqrt(double x)
     
    101101        uint32_t r,t1,s1,ix1,q1;
    102102
    103         EXTRACT_WORDS(ix0, ix1, x);
    104 
    105         /* take care of Inf and NaN */
    106         if ((ix0&0x7ff00000) == 0x7ff00000) {
    107                 return x*x + x;  /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */
    108         }
    109         /* take care of zero */
    110         if (ix0 <= 0) {
    111                 if (((ix0&~sign)|ix1) == 0)
    112                         return x;  /* sqrt(+-0) = +-0 */
    113                 if (ix0 < 0)
    114                         return (x-x)/(x-x);  /* sqrt(-ve) = sNaN */
    115         }
    116         /* normalize x */
    117         m = ix0>>20;
    118         if (m == 0) {  /* subnormal x */
    119                 while (ix0 == 0) {
    120                         m -= 21;
    121                         ix0 |= (ix1>>11);
    122                         ix1 <<= 21;
    123                 }
    124                 for (i=0; (ix0&0x00100000) == 0; i++)
    125                         ix0<<=1;
    126                 m -= i - 1;
    127                 ix0 |= ix1>>(32-i);
    128                 ix1 <<= i;
    129         }
    130         m -= 1023;    /* unbias exponent */
     103        EXTRACT_WORDS(ix0,ix1,x);
     104
     105    /* take care of Inf and NaN */
     106        if((ix0&0x7ff00000)==0x7ff00000) {                     
     107            return x*x+x;               /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
     108                                           sqrt(-inf)=sNaN */
     109        }
     110    /* take care of zero */
     111        if(ix0<=0) {
     112            if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
     113            else if(ix0<0)
     114                return (x-x)/(x-x);             /* sqrt(-ve) = sNaN */
     115        }
     116    /* normalize x */
     117        m = (ix0>>20);
     118        if(m==0) {                              /* subnormal x */
     119            while(ix0==0) {
     120                m -= 21;
     121                ix0 |= (ix1>>11); ix1 <<= 21;
     122            }
     123            for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1;
     124            m -= i-1;
     125            ix0 |= (ix1>>(32-i));
     126            ix1 <<= i;
     127        }
     128        m -= 1023;      /* unbias exponent */
    131129        ix0 = (ix0&0x000fffff)|0x00100000;
    132         if (m & 1) {  /* odd m, double x to make it even */
    133                 ix0 += ix0 + ((ix1&sign)>>31);
    134                 ix1 += ix1;
    135         }
    136         m >>= 1;      /* m = [m/2] */
    137 
    138         /* generate sqrt(x) bit by bit */
     130        if(m&1){        /* odd m, double x to make it even */
     131            ix0 += ix0 + ((ix1&sign)>>31);
     132            ix1 += ix1;
     133        }
     134        m >>= 1;        /* m = [m/2] */
     135
     136    /* generate sqrt(x) bit by bit */
    139137        ix0 += ix0 + ((ix1&sign)>>31);
    140138        ix1 += ix1;
    141         q = q1 = s0 = s1 = 0;  /* [q,q1] = sqrt(x) */
    142         r = 0x00200000;        /* r = moving bit from right to left */
    143 
    144         while (r != 0) {
    145                 t = s0 + r;
    146                 if (t <= ix0) {
    147                         s0   = t + r;
    148                         ix0 -= t;
    149                         q   += r;
    150                 }
    151                 ix0 += ix0 + ((ix1&sign)>>31);
    152                 ix1 += ix1;
    153                 r >>= 1;
     139        q = q1 = s0 = s1 = 0;   /* [q,q1] = sqrt(x) */
     140        r = 0x00200000;         /* r = moving bit from right to left */
     141
     142        while(r!=0) {
     143            t = s0+r;
     144            if(t<=ix0) {
     145                s0   = t+r;
     146                ix0 -= t;
     147                q   += r;
     148            }
     149            ix0 += ix0 + ((ix1&sign)>>31);
     150            ix1 += ix1;
     151            r>>=1;
    154152        }
    155153
    156154        r = sign;
    157         while (r != 0) {
    158                 t1 = s1 + r;
    159                 t  = s0;
    160                 if (t < ix0 || (t == ix0 && t1 <= ix1)) {
    161                         s1 = t1 + r;
    162                         if ((t1&sign) == (uint32_t)sign && (s1&sign) == 0)
    163                                 s0++;
    164                         ix0 -= t;
    165                         if (ix1 < t1)
    166                                 ix0--;
    167                         ix1 -= t1;
    168                         q1 += r;
    169                 }
    170                 ix0 += ix0 + ((ix1&sign)>>31);
    171                 ix1 += ix1;
    172                 r >>= 1;
    173         }
    174 
    175         /* use floating add to find out rounding direction */
    176         if ((ix0|ix1) != 0) {
    177                 z = 1.0 - tiny; /* raise inexact flag */
    178                 if (z >= 1.0) {
    179                         z = 1.0 + tiny;
    180                         if (q1 == (uint32_t)0xffffffff) {
    181                                 q1 = 0;
    182                                 q++;
    183                         } else if (z > 1.0) {
    184                                 if (q1 == (uint32_t)0xfffffffe)
    185                                         q++;
    186                                 q1 += 2;
    187                         } else
    188                                 q1 += q1 & 1;
    189                 }
    190         }
    191         ix0 = (q>>1) + 0x3fe00000;
    192         ix1 = q1>>1;
    193         if (q&1)
    194                 ix1 |= sign;
    195         ix0 += m << 20;
    196         INSERT_WORDS(z, ix0, ix1);
     155        while(r!=0) {
     156            t1 = s1+r;
     157            t  = s0;
     158            if((t<ix0)||((t==ix0)&&(t1<=ix1))) {
     159                s1  = t1+r;
     160                if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1;
     161                ix0 -= t;
     162                if (ix1 < t1) ix0 -= 1;
     163                ix1 -= t1;
     164                q1  += r;
     165            }
     166            ix0 += ix0 + ((ix1&sign)>>31);
     167            ix1 += ix1;
     168            r>>=1;
     169        }
     170
     171    /* use floating add to find out rounding direction */
     172        if((ix0|ix1)!=0) {
     173            z = one-tiny; /* trigger inexact flag */
     174            if (z>=one) {
     175                z = one+tiny;
     176                if (q1==(uint32_t)0xffffffff) { q1=0; q += 1;}
     177                else if (z>one) {
     178                    if (q1==(uint32_t)0xfffffffe) q+=1;
     179                    q1+=2;
     180                } else
     181                    q1 += (q1&1);
     182            }
     183        }
     184        ix0 = (q>>1)+0x3fe00000;
     185        ix1 =  q1>>1;
     186        if ((q&1)==1) ix1 |= sign;
     187        ix0 += (m <<20);
     188        INSERT_WORDS(z,ix0,ix1);
    197189        return z;
    198190}
Note: See TracChangeset for help on using the changeset viewer.