Print this page
5261 libm should stop using synonyms.h
5298 fabs is 0-sized, confuses dis(1) and others
Reviewed by: Josef 'Jeff' Sipek <jeffpc@josefsipek.net>
Approved by: Gordon Ross <gwr@nexenta.com>


  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 #include <sys/isa_defs.h>
  31 #include "libm_synonyms.h"
  32 #include "libm_inlines.h"
  33 
  34 #ifdef _LITTLE_ENDIAN
  35 #define HI(x)   *(1+(int*)x)
  36 #define LO(x)   *(unsigned*)x
  37 #else
  38 #define HI(x)   *(int*)x
  39 #define LO(x)   *(1+(unsigned*)x)
  40 #endif
  41 
  42 #ifdef __RESTRICT
  43 #define restrict _Restrict
  44 #else
  45 #define restrict
  46 #endif
  47 
  48 /* double rhypot(double x, double y)
  49  *
  50  * Method :
  51  *      1. Special cases:
  52  *              x or  y = Inf                                   => 0
  53  *              x or  y = NaN                                   => QNaN
  54  *              x and y = 0                                     => Inf + divide-by-zero
  55  *      2. Computes rhypot(x,y):
  56  *              rhypot(x,y) = m * sqrt(1/(xnm * xnm + ynm * ynm))
  57  *      Where:
  58  *              m = 1/max(|x|,|y|)
  59  *              xnm = x * m
  60  *              ynm = y * m
  61  *
  62  *      Compute 1/(xnm * xnm + ynm * ynm) by simulating
  63  *      muti-precision arithmetic.
  64  *
  65  * Accuracy:
  66  *      Maximum error observed: less than 0.869 ulp after 1.000.000.000
  67  *      results.
  68  */
  69 
  70 #define sqrt __sqrt
  71 
  72 extern double sqrt(double);
  73 
  74 extern double fabs(double);
  75 
  76 static const int __vlibm_TBL_rhypot[] = {
  77 /* i = [0,127]
  78  * TBL[i] = 0x3ff00000 + *(int*)&(1.0 / *(double*)&(0x3ff0000000000000ULL + (i << 45))); */
  79  0x7fe00000,  0x7fdfc07f, 0x7fdf81f8,  0x7fdf4465,
  80  0x7fdf07c1,  0x7fdecc07, 0x7fde9131,  0x7fde573a,
  81  0x7fde1e1e,  0x7fdde5d6, 0x7fddae60,  0x7fdd77b6,
  82  0x7fdd41d4,  0x7fdd0cb5, 0x7fdcd856,  0x7fdca4b3,
  83  0x7fdc71c7,  0x7fdc3f8f, 0x7fdc0e07,  0x7fdbdd2b,
  84  0x7fdbacf9,  0x7fdb7d6c, 0x7fdb4e81,  0x7fdb2036,
  85  0x7fdaf286,  0x7fdac570, 0x7fda98ef,  0x7fda6d01,
  86  0x7fda41a4,  0x7fda16d3, 0x7fd9ec8e,  0x7fd9c2d1,
  87  0x7fd99999,  0x7fd970e4, 0x7fd948b0,  0x7fd920fb,
  88  0x7fd8f9c1,  0x7fd8d301, 0x7fd8acb9,  0x7fd886e5,
  89  0x7fd86186,  0x7fd83c97, 0x7fd81818,  0x7fd7f405,
  90  0x7fd7d05f,  0x7fd7ad22, 0x7fd78a4c,  0x7fd767dc,
  91  0x7fd745d1,  0x7fd72428, 0x7fd702e0,  0x7fd6e1f7,
  92  0x7fd6c16c,  0x7fd6a13c, 0x7fd68168,  0x7fd661ec,
  93  0x7fd642c8,  0x7fd623fa, 0x7fd60581,  0x7fd5e75b,


 411                         itbl1 = ((int*)((char*)__vlibm_TBL_rhypot + iarr1))[0];
 412                         itbl1 -= iexp1;
 413                         HI(&dd1) = itbl1;
 414                         LO(&dd1) = 0;
 415 
 416                         dd1 = dd1 * (DTWO - dd1 * dres1);
 417                         dd1 = dd1 * (DTWO - dd1 * dres1);
 418                         dres1 = dd1 * (DTWO - dd1 * dres1);
 419 
 420                         HI(&res1) = HI(&dres1) & 0xffffff00;
 421                         LO(&res1) = 0;
 422                         res1 += (DONE - res1_hi * res1 - res1_lo * res1) * dres1;
 423                         res1 = sqrt (res1);
 424 
 425                         res1 = scl1 * res1;
 426 
 427                         *pz1 = res1;
 428                 }
 429         }
 430 }
 431 


  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 #include <sys/isa_defs.h>

  31 #include "libm_inlines.h"
  32 
  33 #ifdef _LITTLE_ENDIAN
  34 #define HI(x)   *(1+(int*)x)
  35 #define LO(x)   *(unsigned*)x
  36 #else
  37 #define HI(x)   *(int*)x
  38 #define LO(x)   *(1+(unsigned*)x)
  39 #endif
  40 
  41 #ifdef __RESTRICT
  42 #define restrict _Restrict
  43 #else
  44 #define restrict
  45 #endif
  46 
  47 /* double rhypot(double x, double y)
  48  *
  49  * Method :
  50  *      1. Special cases:
  51  *              x or  y = Inf                                   => 0
  52  *              x or  y = NaN                                   => QNaN
  53  *              x and y = 0                                     => Inf + divide-by-zero
  54  *      2. Computes rhypot(x,y):
  55  *              rhypot(x,y) = m * sqrt(1/(xnm * xnm + ynm * ynm))
  56  *      Where:
  57  *              m = 1/max(|x|,|y|)
  58  *              xnm = x * m
  59  *              ynm = y * m
  60  *
  61  *      Compute 1/(xnm * xnm + ynm * ynm) by simulating
  62  *      muti-precision arithmetic.
  63  *
  64  * Accuracy:
  65  *      Maximum error observed: less than 0.869 ulp after 1.000.000.000
  66  *      results.
  67  */
  68 


  69 extern double sqrt(double);

  70 extern double fabs(double);
  71 
  72 static const int __vlibm_TBL_rhypot[] = {
  73 /* i = [0,127]
  74  * TBL[i] = 0x3ff00000 + *(int*)&(1.0 / *(double*)&(0x3ff0000000000000ULL + (i << 45))); */
  75  0x7fe00000,  0x7fdfc07f, 0x7fdf81f8,  0x7fdf4465,
  76  0x7fdf07c1,  0x7fdecc07, 0x7fde9131,  0x7fde573a,
  77  0x7fde1e1e,  0x7fdde5d6, 0x7fddae60,  0x7fdd77b6,
  78  0x7fdd41d4,  0x7fdd0cb5, 0x7fdcd856,  0x7fdca4b3,
  79  0x7fdc71c7,  0x7fdc3f8f, 0x7fdc0e07,  0x7fdbdd2b,
  80  0x7fdbacf9,  0x7fdb7d6c, 0x7fdb4e81,  0x7fdb2036,
  81  0x7fdaf286,  0x7fdac570, 0x7fda98ef,  0x7fda6d01,
  82  0x7fda41a4,  0x7fda16d3, 0x7fd9ec8e,  0x7fd9c2d1,
  83  0x7fd99999,  0x7fd970e4, 0x7fd948b0,  0x7fd920fb,
  84  0x7fd8f9c1,  0x7fd8d301, 0x7fd8acb9,  0x7fd886e5,
  85  0x7fd86186,  0x7fd83c97, 0x7fd81818,  0x7fd7f405,
  86  0x7fd7d05f,  0x7fd7ad22, 0x7fd78a4c,  0x7fd767dc,
  87  0x7fd745d1,  0x7fd72428, 0x7fd702e0,  0x7fd6e1f7,
  88  0x7fd6c16c,  0x7fd6a13c, 0x7fd68168,  0x7fd661ec,
  89  0x7fd642c8,  0x7fd623fa, 0x7fd60581,  0x7fd5e75b,


 407                         itbl1 = ((int*)((char*)__vlibm_TBL_rhypot + iarr1))[0];
 408                         itbl1 -= iexp1;
 409                         HI(&dd1) = itbl1;
 410                         LO(&dd1) = 0;
 411 
 412                         dd1 = dd1 * (DTWO - dd1 * dres1);
 413                         dd1 = dd1 * (DTWO - dd1 * dres1);
 414                         dres1 = dd1 * (DTWO - dd1 * dres1);
 415 
 416                         HI(&res1) = HI(&dres1) & 0xffffff00;
 417                         LO(&res1) = 0;
 418                         res1 += (DONE - res1_hi * res1 - res1_lo * res1) * dres1;
 419                         res1 = sqrt (res1);
 420 
 421                         res1 = scl1 * res1;
 422 
 423                         *pz1 = res1;
 424                 }
 425         }
 426 }