1 /*
   2  * CDDL HEADER START
   3  *
   4  * The contents of this file are subject to the terms of the
   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21 
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 /*
  26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27  * Use is subject to license terms.
  28  */
  29 
  30 #pragma weak sqrtl = __sqrtl
  31 
  32 #include "libm.h"
  33 #include "longdouble.h"
  34 
  35 extern int __swapTE(int);
  36 extern int __swapEX(int);
  37 extern enum fp_direction_type __swapRD(enum fp_direction_type);
  38 
  39 /*
  40  * in struct longdouble, msw consists of
  41  *      unsigned short  sgn:1;
  42  *      unsigned short  exp:15;
  43  *      unsigned short  frac1:16;
  44  */
  45 
  46 #ifdef __LITTLE_ENDIAN
  47 
  48 /* array indices used to access words within a double */
  49 #define HIWORD  1
  50 #define LOWORD  0
  51 
  52 /* structure used to access words within a quad */
  53 union longdouble {
  54         struct {
  55                 unsigned int    frac4;
  56                 unsigned int    frac3;
  57                 unsigned int    frac2;
  58                 unsigned int    msw;
  59         } l;
  60         long double     d;
  61 };
  62 
  63 /* default NaN returned for sqrt(neg) */
  64 static const union longdouble
  65         qnan = { 0xffffffff, 0xffffffff, 0xffffffff, 0x7fffffff };
  66 
  67 /* signalling NaN used to raise invalid */
  68 static const union {
  69         unsigned u[2];
  70         double d;
  71 } snan = { 0, 0x7ff00001 };
  72 
  73 #else
  74 
  75 /* array indices used to access words within a double */
  76 #define HIWORD  0
  77 #define LOWORD  1
  78 
  79 /* structure used to access words within a quad */
  80 union longdouble {
  81         struct {
  82                 unsigned int    msw;
  83                 unsigned int    frac2;
  84                 unsigned int    frac3;
  85                 unsigned int    frac4;
  86         } l;
  87         long double     d;
  88 };
  89 
  90 /* default NaN returned for sqrt(neg) */
  91 static const union longdouble
  92         qnan = { 0x7fffffff, 0xffffffff, 0xffffffff, 0xffffffff };
  93 
  94 /* signalling NaN used to raise invalid */
  95 static const union {
  96         unsigned u[2];
  97         double d;
  98 } snan = { 0x7ff00001, 0 };
  99 
 100 #endif /* __LITTLE_ENDIAN */
 101 
 102 
 103 static const double
 104         zero = 0.0,
 105         half = 0.5,
 106         one = 1.0,
 107         huge = 1.0e300,
 108         tiny = 1.0e-300,
 109         two36 = 6.87194767360000000000e+10,
 110         two30 = 1.07374182400000000000e+09,
 111         two6 = 6.40000000000000000000e+01,
 112         two4 = 1.60000000000000000000e+01,
 113         twom18 = 3.81469726562500000000e-06,
 114         twom28 = 3.72529029846191406250e-09,
 115         twom42 = 2.27373675443232059479e-13,
 116         twom60 = 8.67361737988403547206e-19,
 117         twom62 = 2.16840434497100886801e-19,
 118         twom66 = 1.35525271560688054251e-20,
 119         twom90 = 8.07793566946316088742e-28,
 120         twom113 = 9.62964972193617926528e-35,
 121         twom124 = 4.70197740328915003187e-38;
 122 
 123 
 124 /*
 125 *       Extract the exponent and normalized significand (represented as
 126 *       an array of five doubles) from a finite, nonzero quad.
 127 */
 128 static int
 129 __q_unpack( const union longdouble *x, double *s )
 130 {
 131         union {
 132                 double                  d;
 133                 unsigned int    l[2];
 134         } u;
 135         double                  b;
 136         unsigned int    lx, w[3];
 137         int                             ex;
 138 
 139         /* get the normalized significand and exponent */
 140         ex = (int) ( ( x->l.msw & 0x7fffffff ) >> 16 );
 141         lx = x->l.msw & 0xffff;
 142         if ( ex )
 143         {
 144                 lx |= 0x10000;
 145                 w[0] = x->l.frac2;
 146                 w[1] = x->l.frac3;
 147                 w[2] = x->l.frac4;
 148         }
 149         else
 150         {
 151                 if ( lx | ( x->l.frac2 & 0xfffe0000 ) )
 152                 {
 153                         w[0] = x->l.frac2;
 154                         w[1] = x->l.frac3;
 155                         w[2] = x->l.frac4;
 156                         ex = 1;
 157                 }
 158                 else if ( x->l.frac2 | ( x->l.frac3 & 0xfffe0000 ) )
 159                 {
 160                         lx = x->l.frac2;
 161                         w[0] = x->l.frac3;
 162                         w[1] = x->l.frac4;
 163                         w[2] = 0;
 164                         ex = -31;
 165                 }
 166                 else if ( x->l.frac3 | ( x->l.frac4 & 0xfffe0000 ) )
 167                 {
 168                         lx = x->l.frac3;
 169                         w[0] = x->l.frac4;
 170                         w[1] = w[2] = 0;
 171                         ex = -63;
 172                 }
 173                 else
 174                 {
 175                         lx = x->l.frac4;
 176                         w[0] = w[1] = w[2] = 0;
 177                         ex = -95;
 178                 }
 179                 while ( ( lx & 0x10000 ) == 0 )
 180                 {
 181                         lx = ( lx << 1 ) | ( w[0] >> 31 );
 182                         w[0] = ( w[0] << 1 ) | ( w[1] >> 31 );
 183                         w[1] = ( w[1] << 1 ) | ( w[2] >> 31 );
 184                         w[2] <<= 1;
 185                         ex--;
 186                 }
 187         }
 188 
 189         /* extract the significand into five doubles */
 190         u.l[HIWORD] = 0x42300000;
 191         u.l[LOWORD] = 0;
 192         b = u.d;
 193         u.l[LOWORD] = lx;
 194         s[0] = u.d - b;
 195 
 196         u.l[HIWORD] = 0x40300000;
 197         u.l[LOWORD] = 0;
 198         b = u.d;
 199         u.l[LOWORD] = w[0] & 0xffffff00;
 200         s[1] = u.d - b;
 201 
 202         u.l[HIWORD] = 0x3e300000;
 203         u.l[LOWORD] = 0;
 204         b = u.d;
 205         u.l[HIWORD] |= w[0] & 0xff;
 206         u.l[LOWORD] = w[1] & 0xffff0000;
 207         s[2] = u.d - b;
 208 
 209         u.l[HIWORD] = 0x3c300000;
 210         u.l[LOWORD] = 0;
 211         b = u.d;
 212         u.l[HIWORD] |= w[1] & 0xffff;
 213         u.l[LOWORD] = w[2] & 0xff000000;
 214         s[3] = u.d - b;
 215 
 216         u.l[HIWORD] = 0x3c300000;
 217         u.l[LOWORD] = 0;
 218         b = u.d;
 219         u.l[LOWORD] = w[2] & 0xffffff;
 220         s[4] = u.d - b;
 221 
 222         return ex - 0x3fff;
 223 }
 224 
 225 
 226 /*
 227 *       Pack an exponent and array of three doubles representing a finite,
 228 *       nonzero number into a quad.  Assume the sign is already there and
 229 *       the rounding mode has been fudged accordingly.
 230 */
 231 static void
 232 __q_pack( const double *z, int exp, enum fp_direction_type rm,
 233         union longdouble *x, int *inexact )
 234 {
 235         union {
 236                 double                  d;
 237                 unsigned int    l[2];
 238         } u;
 239         double                  s[3], t, t2;
 240         unsigned int    msw, frac2, frac3, frac4;
 241 
 242         /* bias exponent and strip off integer bit */
 243         exp += 0x3fff;
 244         s[0] = z[0] - one;
 245         s[1] = z[1];
 246         s[2] = z[2];
 247 
 248         /*
 249          * chop the significand to obtain the fraction;
 250          * use round-to-minus-infinity to ensure chopping
 251          */
 252         (void) __swapRD( fp_negative );
 253 
 254         /* extract the first eighty bits of fraction */
 255         t = s[1] + s[2];
 256         u.d = two36 + ( s[0] + t );
 257         msw = u.l[LOWORD];
 258         s[0] -= ( u.d - two36 );
 259 
 260         u.d = two4 + ( s[0] + t );
 261         frac2 = u.l[LOWORD];
 262         s[0] -= ( u.d - two4 );
 263 
 264         u.d = twom28 + ( s[0] + t );
 265         frac3 = u.l[LOWORD];
 266         s[0] -= ( u.d - twom28 );
 267 
 268         /* condense the remaining fraction; errors here won't matter */
 269         t = s[0] + s[1];
 270         s[1] = ( ( s[0] - t ) + s[1] ) + s[2];
 271         s[0] = t;
 272 
 273         /* get the last word of fraction */
 274         u.d = twom60 + ( s[0] + s[1] );
 275         frac4 = u.l[LOWORD];
 276         s[0] -= ( u.d - twom60 );
 277 
 278         /*
 279          * keep track of what's left for rounding; note that
 280          * t2 will be non-negative due to rounding mode
 281          */
 282         t = s[0] + s[1];
 283         t2 = ( s[0] - t ) + s[1];
 284 
 285         if ( t != zero )
 286         {
 287                 *inexact = 1;
 288 
 289                 /* decide whether to round the fraction up */
 290                 if ( rm == fp_positive || ( rm == fp_nearest && ( t > twom113 ||
 291                         ( t == twom113 && ( t2 != zero || frac4 & 1 ) ) ) ) )
 292                 {
 293                         /* round up and renormalize if necessary */
 294                         if ( ++frac4 == 0 )
 295                                 if ( ++frac3 == 0 )
 296                                         if ( ++frac2 == 0 )
 297                                                 if ( ++msw == 0x10000 )
 298                                                 {
 299                                                         msw = 0;
 300                                                         exp++;
 301                                                 }
 302                 }
 303         }
 304 
 305         /* assemble the result */
 306         x->l.msw |= msw | ( exp << 16 );
 307         x->l.frac2 = frac2;
 308         x->l.frac3 = frac3;
 309         x->l.frac4 = frac4;
 310 }
 311 
 312 
 313 /*
 314 *       Compute the square root of x and place the TP result in s.
 315 */
 316 static void
 317 __q_tp_sqrt( const double *x, double *s )
 318 {
 319         double  c, rr, r[3], tt[3], t[5];
 320 
 321         /* approximate the divisor for the Newton iteration */
 322         c = sqrt( ( x[0] + x[1] ) + x[2] );
 323         rr = half / c;
 324 
 325         /* compute the first five "digits" of the square root */
 326         t[0] = ( c + two30 ) - two30;
 327         tt[0] = t[0] + t[0];
 328         r[0] = ( ( x[0] - t[0] * t[0] ) + x[1] ) + x[2];
 329 
 330         t[1] = ( rr * ( r[0] + x[3] ) + two6 ) - two6;
 331         tt[1] = t[1] + t[1];
 332         r[0] -= tt[0] * t[1];
 333         r[1] = x[3] - t[1] * t[1];
 334         c = ( r[1] + twom18 ) - twom18;
 335         r[0] += c;
 336         r[1] = ( r[1] - c ) + x[4];
 337 
 338         t[2] = ( rr * ( r[0] + r[1] ) + twom18 ) - twom18;
 339         tt[2] = t[2] + t[2];
 340         r[0] -= tt[0] * t[2];
 341         r[1] -= tt[1] * t[2];
 342         c = ( r[1] + twom42 ) - twom42;
 343         r[0] += c;
 344         r[1] = ( r[1] - c ) - t[2] * t[2];
 345 
 346         t[3] = ( rr * ( r[0] + r[1] ) + twom42 ) - twom42;
 347         r[0] = ( ( r[0] - tt[0] * t[3] ) + r[1] ) - tt[1] * t[3];
 348         r[1] = -tt[2] * t[3];
 349         c = ( r[1] + twom90 ) - twom90;
 350         r[0] += c;
 351         r[1] = ( r[1] - c ) - t[3] * t[3];
 352 
 353         t[4] = ( rr * ( r[0] + r[1] ) + twom66 ) - twom66;
 354 
 355         /* here we just need to get the sign of the remainder */
 356         c = ( ( ( ( ( r[0] - tt[0] * t[4] ) - tt[1] * t[4] ) + r[1] )
 357                 - tt[2] * t[4] ) - ( t[3] + t[3] ) * t[4] ) - t[4] * t[4];
 358 
 359         /* reduce to three doubles */
 360         t[0] += t[1];
 361         t[1] = t[2] + t[3];
 362         t[2] = t[4];
 363 
 364         /* if the third term might lie on a rounding boundary, perturb it */
 365         if ( c != zero && t[2] == ( twom62 + t[2] ) - twom62 )
 366         {
 367                 if ( c < zero )
 368                         t[2] -= twom124;
 369                 else
 370                         t[2] += twom124;
 371         }
 372 
 373         /* condense the square root */
 374         c = t[1] + t[2];
 375         t[2] += ( t[1] - c );
 376         t[1] = c;
 377         c = t[0] + t[1];
 378         s[1] = t[1] + ( t[0] - c );
 379         s[0] = c;
 380         if ( s[1] == zero )
 381         {
 382                 c = s[0] + t[2];
 383                 s[1] = t[2] + ( s[0] - c );
 384                 s[0] = c;
 385                 s[2] = zero;
 386         }
 387         else
 388         {
 389                 c = s[1] + t[2];
 390                 s[2] = t[2] + ( s[1] - c );
 391                 s[1] = c;
 392         }
 393 }
 394 
 395 
 396 long double
 397 sqrtl( long double ldx )
 398 {
 399         union   longdouble              x;
 400         volatile double                 t;
 401         double                                  xx[5], zz[3];
 402         enum fp_direction_type  rm;
 403         int                             ex, inexact, exc, traps;
 404 
 405         /* clear cexc */
 406         t = zero;
 407         t -= zero;
 408 
 409         /* check for zero operand */
 410         x.d = ldx;
 411         if ( !( ( x.l.msw & 0x7fffffff ) | x.l.frac2 | x.l.frac3 | x.l.frac4 ) )
 412                 return ldx;
 413 
 414         /* handle nan and inf cases */
 415         if ( ( x.l.msw & 0x7fffffff ) >= 0x7fff0000 )
 416         {
 417                 if ( ( x.l.msw & 0xffff ) | x.l.frac2 | x.l.frac3 | x.l.frac4 )
 418                 {
 419                         if ( !( x.l.msw & 0x8000 ) )
 420                         {
 421                                 /* snan, signal invalid */
 422                                 t += snan.d;
 423                         }
 424                         x.l.msw |= 0x8000;
 425                         return x.d;
 426                 }
 427                 if ( x.l.msw & 0x80000000 )
 428                 {
 429                         /* sqrt(-inf), signal invalid */
 430                         t = -one;
 431                         t = sqrt( t );
 432                         return qnan.d;
 433                 }
 434                 /* sqrt(inf), return inf */
 435                 return x.d;
 436         }
 437 
 438         /* handle negative numbers */
 439         if ( x.l.msw & 0x80000000 )
 440         {
 441                 t = -one;
 442                 t = sqrt( t );
 443                 return qnan.d;
 444         }
 445 
 446         /* now x is finite, positive */
 447 
 448         traps = __swapTE( 0 );
 449         exc = __swapEX( 0 );
 450         rm = __swapRD( fp_nearest );
 451 
 452         ex = __q_unpack( &x, xx );
 453         if ( ex & 1 )
 454         {
 455                 /* make exponent even */
 456                 xx[0] += xx[0];
 457                 xx[1] += xx[1];
 458                 xx[2] += xx[2];
 459                 xx[3] += xx[3];
 460                 xx[4] += xx[4];
 461                 ex--;
 462         }
 463         __q_tp_sqrt( xx, zz );
 464 
 465         /* put everything together */
 466         x.l.msw = 0;
 467         inexact = 0;
 468         __q_pack( zz, ex >> 1, rm, &x, &inexact );
 469 
 470         (void) __swapRD( rm );
 471         (void) __swapEX( exc );
 472         (void) __swapTE( traps );
 473         if ( inexact )
 474         {
 475                 t = huge;
 476                 t += tiny;
 477         }
 478         return x.d;
 479 }