1*25c28e83SPiotr Jasiukajtis /* 2*25c28e83SPiotr Jasiukajtis * CDDL HEADER START 3*25c28e83SPiotr Jasiukajtis * 4*25c28e83SPiotr Jasiukajtis * The contents of this file are subject to the terms of the 5*25c28e83SPiotr Jasiukajtis * Common Development and Distribution License (the "License"). 6*25c28e83SPiotr Jasiukajtis * You may not use this file except in compliance with the License. 7*25c28e83SPiotr Jasiukajtis * 8*25c28e83SPiotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9*25c28e83SPiotr Jasiukajtis * or http://www.opensolaris.org/os/licensing. 10*25c28e83SPiotr Jasiukajtis * See the License for the specific language governing permissions 11*25c28e83SPiotr Jasiukajtis * and limitations under the License. 12*25c28e83SPiotr Jasiukajtis * 13*25c28e83SPiotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each 14*25c28e83SPiotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15*25c28e83SPiotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the 16*25c28e83SPiotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying 17*25c28e83SPiotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner] 18*25c28e83SPiotr Jasiukajtis * 19*25c28e83SPiotr Jasiukajtis * CDDL HEADER END 20*25c28e83SPiotr Jasiukajtis */ 21*25c28e83SPiotr Jasiukajtis 22*25c28e83SPiotr Jasiukajtis /* 23*25c28e83SPiotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 24*25c28e83SPiotr Jasiukajtis */ 25*25c28e83SPiotr Jasiukajtis /* 26*25c28e83SPiotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 27*25c28e83SPiotr Jasiukajtis * Use is subject to license terms. 28*25c28e83SPiotr Jasiukajtis */ 29*25c28e83SPiotr Jasiukajtis 30*25c28e83SPiotr Jasiukajtis #if defined(ELFOBJ) 31*25c28e83SPiotr Jasiukajtis #pragma weak powl = __powl 32*25c28e83SPiotr Jasiukajtis #endif 33*25c28e83SPiotr Jasiukajtis 34*25c28e83SPiotr Jasiukajtis #include "libm.h" 35*25c28e83SPiotr Jasiukajtis #include "xpg6.h" /* __xpg6 */ 36*25c28e83SPiotr Jasiukajtis #define _C99SUSv3_pow _C99SUSv3_pow_treats_Inf_as_an_even_int 37*25c28e83SPiotr Jasiukajtis 38*25c28e83SPiotr Jasiukajtis #if defined(__sparc) 39*25c28e83SPiotr Jasiukajtis #define i0 0 40*25c28e83SPiotr Jasiukajtis #define i1 1 41*25c28e83SPiotr Jasiukajtis #define i2 2 42*25c28e83SPiotr Jasiukajtis #define i3 3 43*25c28e83SPiotr Jasiukajtis 44*25c28e83SPiotr Jasiukajtis static const long double zero = 0.0L, one = 1.0L, two = 2.0L; 45*25c28e83SPiotr Jasiukajtis 46*25c28e83SPiotr Jasiukajtis extern const long double _TBL_logl_hi[], _TBL_logl_lo[]; 47*25c28e83SPiotr Jasiukajtis 48*25c28e83SPiotr Jasiukajtis static const long double 49*25c28e83SPiotr Jasiukajtis two113 = 10384593717069655257060992658440192.0L, 50*25c28e83SPiotr Jasiukajtis ln2hi = 6.931471805599453094172319547495844850203e-0001L, 51*25c28e83SPiotr Jasiukajtis ln2lo = 1.667085920830552208890449330400379754169e-0025L, 52*25c28e83SPiotr Jasiukajtis A2 = 6.666666666666666666666666666666091393804e-0001L, 53*25c28e83SPiotr Jasiukajtis A3 = 4.000000000000000000000000407167070220671e-0001L, 54*25c28e83SPiotr Jasiukajtis A4 = 2.857142857142857142730077490612903681164e-0001L, 55*25c28e83SPiotr Jasiukajtis A5 = 2.222222222222242577702836920812882605099e-0001L, 56*25c28e83SPiotr Jasiukajtis A6 = 1.818181816435493395985912667105885828356e-0001L, 57*25c28e83SPiotr Jasiukajtis A7 = 1.538537835211839751112067512805496931725e-0001L, 58*25c28e83SPiotr Jasiukajtis B1 = 6.666666666666666666666666666666666667787e-0001L, 59*25c28e83SPiotr Jasiukajtis B2 = 3.999999999999999999999999999999848524411e-0001L, 60*25c28e83SPiotr Jasiukajtis B3 = 2.857142857142857142857142865084581075070e-0001L, 61*25c28e83SPiotr Jasiukajtis B4 = 2.222222222222222222222010781800643808497e-0001L, 62*25c28e83SPiotr Jasiukajtis B5 = 1.818181818181818185051442171337036403674e-0001L, 63*25c28e83SPiotr Jasiukajtis B6 = 1.538461538461508363540720286292008207673e-0001L, 64*25c28e83SPiotr Jasiukajtis B7 = 1.333333333506731842033180638329317108428e-0001L, 65*25c28e83SPiotr Jasiukajtis B8 = 1.176469984587418890634302788283946761670e-0001L, 66*25c28e83SPiotr Jasiukajtis B9 = 1.053794891561452331722969901564862497132e-0001L; 67*25c28e83SPiotr Jasiukajtis 68*25c28e83SPiotr Jasiukajtis static long double 69*25c28e83SPiotr Jasiukajtis logl_x(long double x, long double *w) { 70*25c28e83SPiotr Jasiukajtis long double f, f1, v, s, z, qn, h, t; 71*25c28e83SPiotr Jasiukajtis int *px = (int *) &x; 72*25c28e83SPiotr Jasiukajtis int *pz = (int *) &z; 73*25c28e83SPiotr Jasiukajtis int i, j, ix, n; 74*25c28e83SPiotr Jasiukajtis 75*25c28e83SPiotr Jasiukajtis n = 0; 76*25c28e83SPiotr Jasiukajtis ix = px[i0]; 77*25c28e83SPiotr Jasiukajtis if (ix > 0x3ffef03f && ix < 0x3fff0820) { /* 65/63 > x > 63/65 */ 78*25c28e83SPiotr Jasiukajtis f = x - one; 79*25c28e83SPiotr Jasiukajtis z = f * f; 80*25c28e83SPiotr Jasiukajtis if (((ix - 0x3fff0000) | px[i1] | px[i2] | px[i3]) == 0) { 81*25c28e83SPiotr Jasiukajtis *w = zero; 82*25c28e83SPiotr Jasiukajtis return (zero); /* log(1)= +0 */ 83*25c28e83SPiotr Jasiukajtis } 84*25c28e83SPiotr Jasiukajtis qn = one / (two + f); 85*25c28e83SPiotr Jasiukajtis s = f * qn; /* |s|<2**-6 */ 86*25c28e83SPiotr Jasiukajtis v = s * s; 87*25c28e83SPiotr Jasiukajtis h = (long double) (2.0 * (double) s); 88*25c28e83SPiotr Jasiukajtis f1 = (long double) ((double) f); 89*25c28e83SPiotr Jasiukajtis t = ((two * (f - h) - h * f1) - h * (f - f1)) * qn + 90*25c28e83SPiotr Jasiukajtis s * (v * (B1 + v * (B2 + v * (B3 + v * (B4 + 91*25c28e83SPiotr Jasiukajtis v * (B5 + v * (B6 + v * (B7 + v * (B8 + v * B9))))))))); 92*25c28e83SPiotr Jasiukajtis s = (long double) ((double) (h + t)); 93*25c28e83SPiotr Jasiukajtis *w = t - (s - h); 94*25c28e83SPiotr Jasiukajtis return (s); 95*25c28e83SPiotr Jasiukajtis } 96*25c28e83SPiotr Jasiukajtis if (ix < 0x00010000) { /* subnormal x */ 97*25c28e83SPiotr Jasiukajtis x *= two113; 98*25c28e83SPiotr Jasiukajtis n = -113; 99*25c28e83SPiotr Jasiukajtis ix = px[i0]; 100*25c28e83SPiotr Jasiukajtis } 101*25c28e83SPiotr Jasiukajtis /* LARGE_N */ 102*25c28e83SPiotr Jasiukajtis n += ((ix + 0x200) >> 16) - 0x3fff; 103*25c28e83SPiotr Jasiukajtis ix = (ix & 0x0000ffff) | 0x3fff0000; /* scale x to [1,2] */ 104*25c28e83SPiotr Jasiukajtis px[i0] = ix; 105*25c28e83SPiotr Jasiukajtis i = ix + 0x200; 106*25c28e83SPiotr Jasiukajtis pz[i0] = i & 0xfffffc00; 107*25c28e83SPiotr Jasiukajtis pz[i1] = pz[i2] = pz[i3] = 0; 108*25c28e83SPiotr Jasiukajtis qn = one / (x + z); 109*25c28e83SPiotr Jasiukajtis f = x - z; 110*25c28e83SPiotr Jasiukajtis s = f * qn; 111*25c28e83SPiotr Jasiukajtis f1 = (long double) ((double) f); 112*25c28e83SPiotr Jasiukajtis h = (long double) (2.0 * (double) s); 113*25c28e83SPiotr Jasiukajtis t = qn * ((two * (f - z * h) - h * f1) - h * (f - f1)); 114*25c28e83SPiotr Jasiukajtis j = (i >> 10) & 0x3f; 115*25c28e83SPiotr Jasiukajtis v = s * s; 116*25c28e83SPiotr Jasiukajtis qn = (long double) n; 117*25c28e83SPiotr Jasiukajtis t += qn * ln2lo + _TBL_logl_lo[j]; 118*25c28e83SPiotr Jasiukajtis t += s * (v * (A2 + v * (A3 + v * (A4 + v * (A5 + v * (A6 + 119*25c28e83SPiotr Jasiukajtis v * A7)))))); 120*25c28e83SPiotr Jasiukajtis v = qn * ln2hi + _TBL_logl_hi[j]; 121*25c28e83SPiotr Jasiukajtis s = h + v; 122*25c28e83SPiotr Jasiukajtis t += (h - (s - v)); 123*25c28e83SPiotr Jasiukajtis z = (long double) ((double) (s + t)); 124*25c28e83SPiotr Jasiukajtis *w = t - (z - s); 125*25c28e83SPiotr Jasiukajtis return (z); 126*25c28e83SPiotr Jasiukajtis } 127*25c28e83SPiotr Jasiukajtis 128*25c28e83SPiotr Jasiukajtis extern const long double _TBL_expl_hi[], _TBL_expl_lo[]; 129*25c28e83SPiotr Jasiukajtis static const long double 130*25c28e83SPiotr Jasiukajtis invln2_32 = 4.616624130844682903551758979206054839765e+1L, 131*25c28e83SPiotr Jasiukajtis ln2_32hi = 2.166084939249829091928849858592451515688e-2L, 132*25c28e83SPiotr Jasiukajtis ln2_32lo = 5.209643502595475652782654157501186731779e-27L, 133*25c28e83SPiotr Jasiukajtis ln2_64 = 1.083042469624914545964425189778400898568e-2L; 134*25c28e83SPiotr Jasiukajtis 135*25c28e83SPiotr Jasiukajtis long double 136*25c28e83SPiotr Jasiukajtis powl(long double x, long double y) { 137*25c28e83SPiotr Jasiukajtis long double z, ax; 138*25c28e83SPiotr Jasiukajtis long double y1, y2, w1, w2; 139*25c28e83SPiotr Jasiukajtis int sbx, sby, j, k, yisint, m; 140*25c28e83SPiotr Jasiukajtis int hx, lx, hy, ly, ahx, ahy; 141*25c28e83SPiotr Jasiukajtis int *pz = (int *) &z; 142*25c28e83SPiotr Jasiukajtis int *px = (int *) &x; 143*25c28e83SPiotr Jasiukajtis int *py = (int *) &y; 144*25c28e83SPiotr Jasiukajtis 145*25c28e83SPiotr Jasiukajtis hx = px[i0]; 146*25c28e83SPiotr Jasiukajtis lx = px[i1] | px[i2] | px[i3]; 147*25c28e83SPiotr Jasiukajtis hy = py[i0]; 148*25c28e83SPiotr Jasiukajtis ly = py[i1] | py[i2] | py[i3]; 149*25c28e83SPiotr Jasiukajtis ahx = hx & ~0x80000000; 150*25c28e83SPiotr Jasiukajtis ahy = hy & ~0x80000000; 151*25c28e83SPiotr Jasiukajtis 152*25c28e83SPiotr Jasiukajtis if ((ahy | ly) == 0) 153*25c28e83SPiotr Jasiukajtis return (one); /* x**+-0 = 1 */ 154*25c28e83SPiotr Jasiukajtis else if (hx == 0x3fff0000 && lx == 0 && 155*25c28e83SPiotr Jasiukajtis (__xpg6 & _C99SUSv3_pow) != 0) 156*25c28e83SPiotr Jasiukajtis return (one); /* C99: 1**anything = 1 */ 157*25c28e83SPiotr Jasiukajtis else if (ahx > 0x7fff0000 || (ahx == 0x7fff0000 && lx != 0) || 158*25c28e83SPiotr Jasiukajtis ahy > 0x7fff0000 || (ahy == 0x7fff0000 && ly != 0)) 159*25c28e83SPiotr Jasiukajtis return (x + y); /* +-NaN return x+y */ 160*25c28e83SPiotr Jasiukajtis /* includes Sun: 1**NaN = NaN */ 161*25c28e83SPiotr Jasiukajtis sbx = (unsigned) hx >> 31; 162*25c28e83SPiotr Jasiukajtis sby = (unsigned) hy >> 31; 163*25c28e83SPiotr Jasiukajtis ax = fabsl(x); 164*25c28e83SPiotr Jasiukajtis /* 165*25c28e83SPiotr Jasiukajtis * determine if y is an odd int when x < 0 166*25c28e83SPiotr Jasiukajtis * yisint = 0 ... y is not an integer 167*25c28e83SPiotr Jasiukajtis * yisint = 1 ... y is an odd int 168*25c28e83SPiotr Jasiukajtis * yisint = 2 ... y is an even int 169*25c28e83SPiotr Jasiukajtis */ 170*25c28e83SPiotr Jasiukajtis yisint = 0; 171*25c28e83SPiotr Jasiukajtis if (sbx) { 172*25c28e83SPiotr Jasiukajtis if (ahy >= 0x40700000) /* if |y|>=2**113 */ 173*25c28e83SPiotr Jasiukajtis yisint = 2; /* even integer y */ 174*25c28e83SPiotr Jasiukajtis else if (ahy >= 0x3fff0000) { 175*25c28e83SPiotr Jasiukajtis k = (ahy >> 16) - 0x3fff; /* exponent */ 176*25c28e83SPiotr Jasiukajtis if (k > 80) { 177*25c28e83SPiotr Jasiukajtis j = ((unsigned) py[i3]) >> (112 - k); 178*25c28e83SPiotr Jasiukajtis if ((j << (112 - k)) == py[i3]) 179*25c28e83SPiotr Jasiukajtis yisint = 2 - (j & 1); 180*25c28e83SPiotr Jasiukajtis } else if (k > 48) { 181*25c28e83SPiotr Jasiukajtis j = ((unsigned) py[i2]) >> (80 - k); 182*25c28e83SPiotr Jasiukajtis if ((j << (80 - k)) == py[i2]) 183*25c28e83SPiotr Jasiukajtis yisint = 2 - (j & 1); 184*25c28e83SPiotr Jasiukajtis } else if (k > 16) { 185*25c28e83SPiotr Jasiukajtis j = ((unsigned) py[i1]) >> (48 - k); 186*25c28e83SPiotr Jasiukajtis if ((j << (48 - k)) == py[i1]) 187*25c28e83SPiotr Jasiukajtis yisint = 2 - (j & 1); 188*25c28e83SPiotr Jasiukajtis } else if (ly == 0) { 189*25c28e83SPiotr Jasiukajtis j = ahy >> (16 - k); 190*25c28e83SPiotr Jasiukajtis if ((j << (16 - k)) == ahy) 191*25c28e83SPiotr Jasiukajtis yisint = 2 - (j & 1); 192*25c28e83SPiotr Jasiukajtis } 193*25c28e83SPiotr Jasiukajtis } 194*25c28e83SPiotr Jasiukajtis } 195*25c28e83SPiotr Jasiukajtis 196*25c28e83SPiotr Jasiukajtis /* special value of y */ 197*25c28e83SPiotr Jasiukajtis if (ly == 0) { 198*25c28e83SPiotr Jasiukajtis if (ahy == 0x7fff0000) { /* y is +-inf */ 199*25c28e83SPiotr Jasiukajtis if (((ahx - 0x3fff0000) | lx) == 0) { 200*25c28e83SPiotr Jasiukajtis if ((__xpg6 & _C99SUSv3_pow) != 0) 201*25c28e83SPiotr Jasiukajtis return (one); 202*25c28e83SPiotr Jasiukajtis /* C99: (-1)**+-inf = 1 */ 203*25c28e83SPiotr Jasiukajtis else 204*25c28e83SPiotr Jasiukajtis return (y - y); 205*25c28e83SPiotr Jasiukajtis /* Sun: (+-1)**+-inf = NaN */ 206*25c28e83SPiotr Jasiukajtis } else if (ahx >= 0x3fff0000) 207*25c28e83SPiotr Jasiukajtis /* (|x|>1)**+,-inf = inf,0 */ 208*25c28e83SPiotr Jasiukajtis return (sby == 0 ? y : zero); 209*25c28e83SPiotr Jasiukajtis else /* (|x|<1)**-,+inf = inf,0 */ 210*25c28e83SPiotr Jasiukajtis return (sby != 0 ? -y : zero); 211*25c28e83SPiotr Jasiukajtis } else if (ahy == 0x3fff0000) { /* y is +-1 */ 212*25c28e83SPiotr Jasiukajtis if (sby != 0) 213*25c28e83SPiotr Jasiukajtis return (one / x); 214*25c28e83SPiotr Jasiukajtis else 215*25c28e83SPiotr Jasiukajtis return (x); 216*25c28e83SPiotr Jasiukajtis } else if (hy == 0x40000000) /* y is 2 */ 217*25c28e83SPiotr Jasiukajtis return (x * x); 218*25c28e83SPiotr Jasiukajtis else if (hy == 0x3ffe0000) { /* y is 0.5 */ 219*25c28e83SPiotr Jasiukajtis if (!((ahx | lx) == 0 || ((ahx - 0x7fff0000) | lx) == 220*25c28e83SPiotr Jasiukajtis 0)) 221*25c28e83SPiotr Jasiukajtis return (sqrtl(x)); 222*25c28e83SPiotr Jasiukajtis } 223*25c28e83SPiotr Jasiukajtis } 224*25c28e83SPiotr Jasiukajtis 225*25c28e83SPiotr Jasiukajtis /* special value of x */ 226*25c28e83SPiotr Jasiukajtis if (lx == 0) { 227*25c28e83SPiotr Jasiukajtis if (ahx == 0x7fff0000 || ahx == 0 || ahx == 0x3fff0000) { 228*25c28e83SPiotr Jasiukajtis /* x is +-0,+-inf,+-1 */ 229*25c28e83SPiotr Jasiukajtis z = ax; 230*25c28e83SPiotr Jasiukajtis if (sby == 1) 231*25c28e83SPiotr Jasiukajtis z = one / z; /* z = 1/|x| if y is negative */ 232*25c28e83SPiotr Jasiukajtis if (sbx == 1) { 233*25c28e83SPiotr Jasiukajtis if (ahx == 0x3fff0000 && yisint == 0) 234*25c28e83SPiotr Jasiukajtis z = zero / zero; 235*25c28e83SPiotr Jasiukajtis /* (-1)**non-int is NaN */ 236*25c28e83SPiotr Jasiukajtis else if (yisint == 1) 237*25c28e83SPiotr Jasiukajtis z = -z; /* (x<0)**odd = -(|x|**odd) */ 238*25c28e83SPiotr Jasiukajtis } 239*25c28e83SPiotr Jasiukajtis return (z); 240*25c28e83SPiotr Jasiukajtis } 241*25c28e83SPiotr Jasiukajtis } 242*25c28e83SPiotr Jasiukajtis 243*25c28e83SPiotr Jasiukajtis /* (x<0)**(non-int) is NaN */ 244*25c28e83SPiotr Jasiukajtis if (sbx == 1 && yisint == 0) 245*25c28e83SPiotr Jasiukajtis return (zero / zero); /* should be volatile */ 246*25c28e83SPiotr Jasiukajtis 247*25c28e83SPiotr Jasiukajtis /* Now ax is finite, y is finite */ 248*25c28e83SPiotr Jasiukajtis /* first compute log(ax) = w1+w2, with 53 bits w1 */ 249*25c28e83SPiotr Jasiukajtis w1 = logl_x(ax, &w2); 250*25c28e83SPiotr Jasiukajtis 251*25c28e83SPiotr Jasiukajtis /* split up y into y1+y2 and compute (y1+y2)*(w1+w2) */ 252*25c28e83SPiotr Jasiukajtis if (ly == 0 || ahy >= 0x43fe0000) { 253*25c28e83SPiotr Jasiukajtis y1 = y * w1; 254*25c28e83SPiotr Jasiukajtis y2 = y * w2; 255*25c28e83SPiotr Jasiukajtis } else { 256*25c28e83SPiotr Jasiukajtis y1 = (long double) ((double) y); 257*25c28e83SPiotr Jasiukajtis y2 = (y - y1) * w1 + y * w2; 258*25c28e83SPiotr Jasiukajtis y1 *= w1; 259*25c28e83SPiotr Jasiukajtis } 260*25c28e83SPiotr Jasiukajtis z = y1 + y2; 261*25c28e83SPiotr Jasiukajtis j = pz[i0]; 262*25c28e83SPiotr Jasiukajtis if ((unsigned) j >= 0xffff0000) { /* NaN or -inf */ 263*25c28e83SPiotr Jasiukajtis if (sbx == 1 && yisint == 1) 264*25c28e83SPiotr Jasiukajtis return (one / z); 265*25c28e83SPiotr Jasiukajtis else 266*25c28e83SPiotr Jasiukajtis return (-one / z); 267*25c28e83SPiotr Jasiukajtis } else if ((j & ~0x80000000) < 0x3fc30000) { /* |x|<2^-60 */ 268*25c28e83SPiotr Jasiukajtis if (sbx == 1 && yisint == 1) 269*25c28e83SPiotr Jasiukajtis return (-one - z); 270*25c28e83SPiotr Jasiukajtis else 271*25c28e83SPiotr Jasiukajtis return (one + z); 272*25c28e83SPiotr Jasiukajtis } else if (j > 0) { 273*25c28e83SPiotr Jasiukajtis if (j > 0x400d0000) { 274*25c28e83SPiotr Jasiukajtis if (sbx == 1 && yisint == 1) 275*25c28e83SPiotr Jasiukajtis return (scalbnl(-one, 20000)); 276*25c28e83SPiotr Jasiukajtis else 277*25c28e83SPiotr Jasiukajtis return (scalbnl(one, 20000)); 278*25c28e83SPiotr Jasiukajtis } 279*25c28e83SPiotr Jasiukajtis k = (int) (invln2_32 * (z + ln2_64)); 280*25c28e83SPiotr Jasiukajtis } else { 281*25c28e83SPiotr Jasiukajtis if ((unsigned) j > 0xc00d0000) { 282*25c28e83SPiotr Jasiukajtis if (sbx == 1 && yisint == 1) 283*25c28e83SPiotr Jasiukajtis return (scalbnl(-one, -20000)); 284*25c28e83SPiotr Jasiukajtis else 285*25c28e83SPiotr Jasiukajtis return (scalbnl(one, -20000)); 286*25c28e83SPiotr Jasiukajtis } 287*25c28e83SPiotr Jasiukajtis k = (int) (invln2_32 * (z - ln2_64)); 288*25c28e83SPiotr Jasiukajtis } 289*25c28e83SPiotr Jasiukajtis j = k & 0x1f; 290*25c28e83SPiotr Jasiukajtis m = k >> 5; 291*25c28e83SPiotr Jasiukajtis { 292*25c28e83SPiotr Jasiukajtis /* rational approximation coeffs for [-(ln2)/64,(ln2)/64] */ 293*25c28e83SPiotr Jasiukajtis long double 294*25c28e83SPiotr Jasiukajtis t1 = 1.666666666666666666666666666660876387437e-1L, 295*25c28e83SPiotr Jasiukajtis t2 = -2.777777777777777777777707812093173478756e-3L, 296*25c28e83SPiotr Jasiukajtis t3 = 6.613756613756613482074280932874221202424e-5L, 297*25c28e83SPiotr Jasiukajtis t4 = -1.653439153392139954169609822742235851120e-6L, 298*25c28e83SPiotr Jasiukajtis t5 = 4.175314851769539751387852116610973796053e-8L; 299*25c28e83SPiotr Jasiukajtis long double t = (long double) k; 300*25c28e83SPiotr Jasiukajtis 301*25c28e83SPiotr Jasiukajtis w1 = (y2 - (t * ln2_32hi - y1)) - t * ln2_32lo; 302*25c28e83SPiotr Jasiukajtis t = w1 * w1; 303*25c28e83SPiotr Jasiukajtis w2 = (w1 - t * (t1 + t * (t2 + t * (t3 + t * (t4 + t * t5))))) - 304*25c28e83SPiotr Jasiukajtis two; 305*25c28e83SPiotr Jasiukajtis z = _TBL_expl_hi[j] - ((_TBL_expl_hi[j] * (w1 + w1)) / w2 - 306*25c28e83SPiotr Jasiukajtis _TBL_expl_lo[j]); 307*25c28e83SPiotr Jasiukajtis } 308*25c28e83SPiotr Jasiukajtis j = m + (pz[i0] >> 16); 309*25c28e83SPiotr Jasiukajtis if (j && (unsigned) j < 0x7fff) 310*25c28e83SPiotr Jasiukajtis pz[i0] += m << 16; 311*25c28e83SPiotr Jasiukajtis else 312*25c28e83SPiotr Jasiukajtis z = scalbnl(z, m); 313*25c28e83SPiotr Jasiukajtis 314*25c28e83SPiotr Jasiukajtis if (sbx == 1 && yisint == 1) 315*25c28e83SPiotr Jasiukajtis z = -z; /* (-ve)**(odd int) */ 316*25c28e83SPiotr Jasiukajtis return (z); 317*25c28e83SPiotr Jasiukajtis } 318*25c28e83SPiotr Jasiukajtis #else 319*25c28e83SPiotr Jasiukajtis #error Unsupported Architecture 320*25c28e83SPiotr Jasiukajtis #endif /* defined(__sparc) */ 321