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 /* INDENT OFF */
  31 /*
  32  * long double __k_cexpl(long double x, int *n);
  33  * Returns the exponential of x in the form of 2**n * y, y=__k_cexpl(x,&n).
  34  *
  35  *      1. Argument Reduction: given the input x, find r and integer k
  36  *         and j such that
  37  *                   x = (32k+j)*ln2 + r,  |r| <= (1/64)*ln2 .
  38  *
  39  *      2. expl(x) = 2^k * (2^(j/32) + 2^(j/32)*expm1(r))
  40  *         Note:
  41  *         a. expm1(r) = (2r)/(2-R), R = r - r^2*(t1 + t2*r^2)
  42  *         b. 2^(j/32) is represented as
  43  *                      exp2_32_hi[j]+exp2_32_lo[j]
  44  *         where
  45  *              exp2_32_hi[j] = 2^(j/32) rounded
  46  *              exp2_32_lo[j] = 2^(j/32) - exp2_32_hi[j].
  47  *
  48  * Special cases:
  49  *      expl(INF) is INF, expl(NaN) is NaN;
  50  *      expl(-INF)=  0;
  51  *      for finite argument, only expl(0)=1 is exact.
  52  *
  53  * Accuracy:
  54  *      according to an error analysis, the error is always less than
  55  *      an ulp (unit in the last place).
  56  *
  57  * Misc. info.
  58  *      When |x| is really big, say |x| > 1000000, the accuracy
  59  *      is not important because the ultimate result will over or under
  60  *      flow. So we will simply replace n = 1000000 and r = 0.0. For
  61  *      moderate size x, according to an error analysis, the error is
  62  *      always less than 1 ulp (unit in the last place).
  63  *
  64  * Constants:
  65  * Only decimal values are given. We assume that the compiler will convert
  66  * from decimal to binary accurately enough to produce the correct
  67  * hexadecimal values.
  68  */
  69 /* INDENT ON */
  70 
  71 #include "libm.h"               /* __k_cexpl */
  72 #include "complex_wrapper.h"    /* HI_XWORD */
  73 
  74 /* INDENT OFF */
  75 /* ln2/32 = 0.0216608493924982909192885037955680177523593791987579766912713 */
  76 #if defined(__x86)
  77 static const long double
  78                         /* 43 significant bits, 21 trailing zeros */
  79 ln2_32hi = 2.166084939249657281834515742957592010498046875e-2L,
  80 ln2_32lo = 1.7181009433463659920976473789104487579766912713e-15L;
  81 static const long double exp2_32_hi[] = {       /* exp2_32[j] = 2^(j/32) */
  82         1.0000000000000000000000000e+00L,
  83         1.0218971486541166782081522e+00L,
  84         1.0442737824274138402382006e+00L,
  85         1.0671404006768236181297224e+00L,
  86         1.0905077326652576591003302e+00L,
  87         1.1143867425958925362894369e+00L,
  88         1.1387886347566916536971221e+00L,
  89         1.1637248587775775137938619e+00L,
  90         1.1892071150027210666875674e+00L,
  91         1.2152473599804688780476325e+00L,
  92         1.2418578120734840485256747e+00L,
  93         1.2690509571917332224885722e+00L,
  94         1.2968395546510096659215822e+00L,
  95         1.3252366431597412945939118e+00L,
  96         1.3542555469368927282668852e+00L,
  97         1.3839098819638319548151403e+00L,
  98         1.4142135623730950487637881e+00L,
  99         1.4451808069770466200253470e+00L,
 100         1.4768261459394993113155431e+00L,
 101         1.5091644275934227397133885e+00L,
 102         1.5422108254079408235859630e+00L,
 103         1.5759808451078864864006862e+00L,
 104         1.6104903319492543080837174e+00L,
 105         1.6457554781539648445110730e+00L,
 106         1.6817928305074290860378350e+00L,
 107         1.7186192981224779156032914e+00L,
 108         1.7562521603732994831094730e+00L,
 109         1.7947090750031071864148413e+00L,
 110         1.8340080864093424633989166e+00L,
 111         1.8741676341102999013002103e+00L,
 112         1.9152065613971472938202589e+00L,
 113         1.9571441241754002689657438e+00L,
 114 };
 115 static const long double exp2_32_lo[] = {
 116         0.0000000000000000000000000e+00L,
 117         2.6327965667180882569382524e-20L,
 118         8.3765863521895191129661899e-20L,
 119         3.9798705777454504249209575e-20L,
 120         1.0668046596651558640993042e-19L,
 121         1.9376009847285360448117114e-20L,
 122         6.7081819456112953751277576e-21L,
 123         1.9711680502629186462729727e-20L,
 124         2.9932584438449523689104569e-20L,
 125         6.8887754153039109411061914e-20L,
 126         6.8002718741225378942847820e-20L,
 127         6.5846917376975403439742349e-20L,
 128         1.2171958727511372194876001e-20L,
 129         3.5625253228704087115438260e-20L,
 130         3.1129551559077560956309179e-20L,
 131         5.7519192396164779846216492e-20L,
 132         3.7900651177865141593101239e-20L,
 133         1.1659262405698741798080115e-20L,
 134         7.1364385105284695967172478e-20L,
 135         5.2631003710812203588788949e-20L,
 136         2.6328853788732632868460580e-20L,
 137         5.4583950085438242788190141e-20L,
 138         9.5803254376938269960718656e-20L,
 139         7.6837733983874245823512279e-21L,
 140         2.4415965910835093824202087e-20L,
 141         2.6052966871016580981769728e-20L,
 142         2.6876456344632553875309579e-21L,
 143         1.2861930155613700201703279e-20L,
 144         8.8166633394037485606572294e-20L,
 145         2.9788615389580190940837037e-20L,
 146         5.2352341619805098677422139e-20L,
 147         5.2578463064010463732242363e-20L,
 148 };
 149 #else   /* sparc */
 150 static const long double
 151                         /* 0x3FF962E4 2FEFA39E F35793C7 00000000 */
 152 ln2_32hi = 2.166084939249829091928849858592451515688e-2L,
 153 ln2_32lo = 5.209643502595475652782654157501186731779e-27L;
 154 static const long double exp2_32_hi[] = {       /* exp2_32[j] = 2^(j/32) */
 155         1.000000000000000000000000000000000000000e+0000L,
 156         1.021897148654116678234480134783299439782e+0000L,
 157         1.044273782427413840321966478739929008785e+0000L,
 158         1.067140400676823618169521120992809162607e+0000L,
 159         1.090507732665257659207010655760707978993e+0000L,
 160         1.114386742595892536308812956919603067800e+0000L,
 161         1.138788634756691653703830283841511254720e+0000L,
 162         1.163724858777577513813573599092185312343e+0000L,
 163         1.189207115002721066717499970560475915293e+0000L,
 164         1.215247359980468878116520251338798457624e+0000L,
 165         1.241857812073484048593677468726595605511e+0000L,
 166         1.269050957191733222554419081032338004715e+0000L,
 167         1.296839554651009665933754117792451159835e+0000L,
 168         1.325236643159741294629537095498721674113e+0000L,
 169         1.354255546936892728298014740140702804343e+0000L,
 170         1.383909881963831954872659527265192818002e+0000L,
 171         1.414213562373095048801688724209698078570e+0000L,
 172         1.445180806977046620037006241471670905678e+0000L,
 173         1.476826145939499311386907480374049923924e+0000L,
 174         1.509164427593422739766019551033193531420e+0000L,
 175         1.542210825407940823612291862090734841307e+0000L,
 176         1.575980845107886486455270160181905008906e+0000L,
 177         1.610490331949254308179520667357400583459e+0000L,
 178         1.645755478153964844518756724725822445667e+0000L,
 179         1.681792830507429086062250952466429790080e+0000L,
 180         1.718619298122477915629344376456312504516e+0000L,
 181         1.756252160373299483112160619375313221294e+0000L,
 182         1.794709075003107186427703242127781814354e+0000L,
 183         1.834008086409342463487083189588288856077e+0000L,
 184         1.874167634110299901329998949954446534439e+0000L,
 185         1.915206561397147293872611270295830887850e+0000L,
 186         1.957144124175400269018322251626871491190e+0000L,
 187 };
 188 
 189 static const long double exp2_32_lo[] = {
 190         +0.000000000000000000000000000000000000000e+0000L,
 191         +1.805067874203309547455733330545737864651e-0035L,
 192         -9.374520292280427421957567419730832143843e-0035L,
 193         -1.596968447292758770712909630231499971233e-0035L,
 194         +9.112493410125022978511686101672486662119e-0035L,
 195         -6.504228206978548287230374775259388710985e-0035L,
 196         -8.148468844525851137325691767488155323605e-0035L,
 197         -5.066214576721800313372330745142903350963e-0035L,
 198         -1.359830974688816973749875638245919118924e-0035L,
 199         +9.497427635563196470307710566433246597109e-0035L,
 200         -3.283170523176998601615065965333915261932e-0036L,
 201         -5.017235709387190410290186530458428950862e-0035L,
 202         -2.391474797689109171622834301602640139258e-0035L,
 203         -8.350571357633908815298890737944083853080e-0036L,
 204         +7.036756889073265042421737190671876440729e-0035L,
 205         -5.182484853064646457536893018566956189817e-0035L,
 206         +9.422242548621832065692116736394064879758e-0035L,
 207         -3.967500825398862309167306130216418281103e-0035L,
 208         +7.143528991563300614523273615092767243521e-0035L,
 209         +1.159871252867985124246517834100444327747e-0035L,
 210         +4.696933478358115495309739213201874466685e-0035L,
 211         -3.386513175995004710799241984999819165197e-0035L,
 212         -8.587318774298247068868655935103874453522e-0035L,
 213         -9.605951548749350503185499362246069088835e-0035L,
 214         +9.609733932128012784507558697141785813655e-0035L,
 215         +6.378397921440028439244761449780848545957e-0035L,
 216         +7.792430785695864249456461125169277701177e-0035L,
 217         +7.361337767588456524131930836633932195088e-0035L,
 218         -6.472995147913347230035214575612170525266e-0035L,
 219         +8.587474417953698694278798062295229624207e-0035L,
 220         +2.371815422825174835691651228302690977951e-0035L,
 221         -3.026891682096118773004597373421900314256e-0037L,
 222 };
 223 #endif
 224 
 225 static const long double
 226         one = 1.0L,
 227         two = 2.0L,
 228         ln2_64 = 1.083042469624914545964425189778400898568e-2L,
 229         invln2_32 = 4.616624130844682903551758979206054839765e+1L;
 230 
 231 /* rational approximation coeffs for [-(ln2)/64,(ln2)/64] */
 232 static const long double
 233         t1 =  1.666666666666666666666666666660876387437e-1L,
 234         t2 = -2.777777777777777777777707812093173478756e-3L,
 235         t3 =  6.613756613756613482074280932874221202424e-5L,
 236         t4 = -1.653439153392139954169609822742235851120e-6L,
 237         t5 =  4.175314851769539751387852116610973796053e-8L;
 238 /* INDENT ON */
 239 
 240 long double
 241 __k_cexpl(long double x, int *n) {
 242         int hx, ix, j, k;
 243         long double t, r;
 244 
 245         *n = 0;
 246         hx = HI_XWORD(x);
 247         ix = hx & 0x7fffffff;
 248         if (hx >= 0x7fff0000)
 249                 return (x + x); /* NaN of +inf */
 250         if (((unsigned) hx) >= 0xffff0000)
 251                 return (-one / x);      /* NaN or -inf */
 252         if (ix < 0x3fc30000)
 253                 return (one + x);       /* |x|<2^-60 */
 254         if (hx > 0) {
 255                 if (hx > 0x401086a0) {       /* x > 200000 */
 256                         *n = 200000;
 257                         return (one);
 258                 }
 259                 k = (int) (invln2_32 * (x + ln2_64));
 260         } else {
 261                 if (ix > 0x401086a0) {       /* x < -200000 */
 262                         *n = -200000;
 263                         return (one);
 264                 }
 265                 k = (int) (invln2_32 * (x - ln2_64));
 266         }
 267         j = k & 0x1f;
 268         *n = k >> 5;
 269         t = (long double) k;
 270         x = (x - t * ln2_32hi) - t * ln2_32lo;
 271         t = x * x;
 272         r = (x - t * (t1 + t * (t2 + t * (t3 + t * (t4 + t * t5))))) - two;
 273         x = exp2_32_hi[j] - ((exp2_32_hi[j] * (x + x)) / r - exp2_32_lo[j]);
 274         k >>= 5;
 275         if (k > 240) {
 276                 XFSCALE(x, 240);
 277                 *n -= 240;
 278         } else if (k > 0) {
 279                 XFSCALE(x, k);
 280                 *n = 0;
 281         }
 282         return (x);
 283 }