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 /* 31 * expl(x) 32 * Table driven method 33 * Written by K.C. Ng, November 1988. 34 * Algorithm : 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 * _TBL_expl_hi[j]+_TBL_expl_lo[j] 44 * where 45 * _TBL_expl_hi[j] = 2^(j/32) rounded 46 * _TBL_expl_lo[j] = 2^(j/32) - _TBL_expl_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 * For 113 bit long double 59 * if x > 1.135652340629414394949193107797076342845e+4 60 * then expl(x) overflow; 61 * if x < -1.143346274333629787883724384345262150341e+4 62 * then expl(x) underflow 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 70 #pragma weak __expl = expl 71 72 #include "libm.h" 73 74 extern const long double _TBL_expl_hi[], _TBL_expl_lo[]; 75 76 static const long double 77 one = 1.0L, 78 two = 2.0L, 79 ln2_64 = 1.083042469624914545964425189778400898568e-2L, 80 ovflthreshold = 1.135652340629414394949193107797076342845e+4L, 81 unflthreshold = -1.143346274333629787883724384345262150341e+4L, 82 invln2_32 = 4.616624130844682903551758979206054839765e+1L, 83 ln2_32hi = 2.166084939249829091928849858592451515688e-2L, 84 ln2_32lo = 5.209643502595475652782654157501186731779e-27L; 85 86 /* rational approximation coeffs for [-(ln2)/64,(ln2)/64] */ 87 static const long double 88 t1 = 1.666666666666666666666666666660876387437e-1L, 89 t2 = -2.777777777777777777777707812093173478756e-3L, 90 t3 = 6.613756613756613482074280932874221202424e-5L, 91 t4 = -1.653439153392139954169609822742235851120e-6L, 92 t5 = 4.175314851769539751387852116610973796053e-8L; 93 94 long double 95 expl(long double x) { 96 int *px = (int *) &x, ix, j, k, m; 97 long double t, r; 98 99 ix = px[0]; /* high word of x */ 100 if (ix >= 0x7fff0000) 101 return (x + x); /* NaN of +inf */ 102 if (((unsigned) ix) >= 0xffff0000) 103 return (-one / x); /* NaN or -inf */ 104 if ((ix & 0x7fffffff) < 0x3fc30000) { 105 if ((int) x < 1) 106 return (one + x); /* |x|<2^-60 */ 107 } 108 if (ix > 0) { 109 if (x > ovflthreshold) 110 return (scalbnl(x, 20000)); 111 k = (int) (invln2_32 * (x + ln2_64)); 112 } else { 113 if (x < unflthreshold) 114 return (scalbnl(-x, -40000)); 115 k = (int) (invln2_32 * (x - ln2_64)); 116 } 117 j = k&0x1f; 118 m = k>>5; 119 t = (long double) k; 120 x = (x - t * ln2_32hi) - t * ln2_32lo; 121 t = x * x; 122 r = (x - t * (t1 + t * (t2 + t * (t3 + t * (t4 + t * t5))))) - two; 123 x = _TBL_expl_hi[j] - ((_TBL_expl_hi[j] * (x + x)) / r - 124 _TBL_expl_lo[j]); 125 return (scalbnl(x, m)); 126 } 127