Attachment 54635 Details for Bug 82654 – logl.c

logl.c

logl.c (text/x-c++src), 13.06 KB, created by Steven G. Kargl on 2005-09-13 01:30:27 UTC

(hide)

Description:

Filename:

MIME Type:

Creator: Steven G. Kargl

Created: 2005-09-13 01:30:27 UTC

Size: 13.06 KB

patch

obsolete

>/*-
> * Copyright (c) 2005, Steven G. Kargl
> * All rights reserved.
> *
> * Redistribution and use in source and binary forms, with or without
> * modification, are permitted provided that the following conditions
> * are met:
> * 1. Redistributions of source code must retain the above copyright
> *    notice unmodified, this list of conditions, and the following
> *    disclaimer.
> * 2. Redistributions in binary form must reproduce the above copyright
> *    notice, this list of conditions and the following disclaimer in the
> *    documentation and/or other materials provided with the distribution.
> *
> * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
> * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
> * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
> * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
> * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
> * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
> * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
> * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
> * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
> * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
> */
>
>/*
> *  Compute the natural logarithm of x by decomposing x into its base 2
> *  representation.
> *
> *  log(x) = log(F * 2**n) such that F = [1/2,1).
> *         = log(F) + n * log(2)
> *
> *  Let F = f_n + e such that f_n is interval endpoint.
> *
> *  log(F) = log(f_n + e)
> *         = log[f_n (1 + e / f_n)]
> *         = log(f_n) + log(1 + e / f_n)
> *
> *  The log(f_n) are tabulated below.  Let d = e / f_n and use
> *
> *  log(1 + d) = d - (1/2)*d**2 + (1/3)*d**3  - (1/4)*d**4 + ... + (1/9)*d**9
> *             = p(d)
> *
> *  We have d in [0,1/interval width].  This is mapped into an interval
> *  of t in [-1,1] and p(d) is converted into a polynomial p(t).  Next,
> *  p(t) is written as an expansion in Chebyshev polynomials, T_n(t), on
> *  the interval [-1,1].  Finally, T_n(t) are replaced by their polynomials
> *  in t and this is then reduced to an 8th order polynomial.
> *
> *  log(x) = n * log(2) + log(f_n) + p(t)
> *
> */
>
>#include <math.h>
>
>#define LOGL2  0x.b17217f7d1cf79acp0L  /* 6.931471805599453094e-1L */
>
>#define INTERVALS 128
>
>static long double f[INTERVALS][2] = {
>   { 0x.8000000000000000p0L, -0x.b17217f7d1cf79acp0L },  /* 128/256 */
>   { 0x.8100000000000000p0L, -0x.af74155120c9011cp0L },  /* 129/256 */
>   { 0x.8200000000000000p0L, -0x.ad7a02e1b24efd32p0L },  /* 130/256 */
>   { 0x.8300000000000000p0L, -0x.ab83d135dc633301p0L },  /* 131/256 */
>   { 0x.8400000000000000p0L, -0x.a991713433c2b999p0L },  /* 132/256 */
>   { 0x.8500000000000000p0L, -0x.a7a2d41ad270c9d7p0L },  /* 133/256 */
>   { 0x.8600000000000000p0L, -0x.a5b7eb7cb860fb89p0L },  /* 134/256 */
>   { 0x.8700000000000000p0L, -0x.a3d0a93f45169a4bp0L },  /* 135/256 */
>   { 0x.8800000000000000p0L, -0x.a1ecff97c91e267bp0L },  /* 136/256 */
>   { 0x.8900000000000000p0L, -0x.a00ce1092e5498c3p0L },  /* 137/256 */
>   { 0x.8a00000000000000p0L, -0x.9e304061b5fda919p0L },  /* 138/256 */
>   { 0x.8b00000000000000p0L, -0x.9c5710b8cbb73a43p0L },  /* 139/256 */
>   { 0x.8c00000000000000p0L, -0x.9a81456cec642e10p0L },  /* 140/256 */
>   { 0x.8d00000000000000p0L, -0x.98aed221a03458b6p0L },  /* 141/256 */
>   { 0x.8e00000000000000p0L, -0x.96dfaabd86fa1647p0L },  /* 142/256 */
>   { 0x.8f00000000000000p0L, -0x.9513c36876083696p0L },  /* 143/256 */
>   { 0x.9000000000000000p0L, -0x.934b1089a6dc93c2p0L },  /* 144/256 */
>   { 0x.9100000000000000p0L, -0x.918586c5f5e4bf02p0L },  /* 145/256 */
>   { 0x.9200000000000000p0L, -0x.8fc31afe30b2c6dfp0L },  /* 146/256 */
>   { 0x.9300000000000000p0L, -0x.8e03c24d7300395ap0L },  /* 147/256 */
>   { 0x.9400000000000000p0L, -0x.8c47720791e53314p0L },  /* 148/256 */
>   { 0x.9500000000000000p0L, -0x.8a8e1fb794b09134p0L },  /* 149/256 */
>   { 0x.9600000000000000p0L, -0x.88d7c11e3ad53cdcp0L },  /* 150/256 */
>   { 0x.9700000000000000p0L, -0x.87244c308e670a66p0L },  /* 151/256 */
>   { 0x.9800000000000000p0L, -0x.8573b71682a7d21bp0L },  /* 152/256 */
>   { 0x.9900000000000000p0L, -0x.83c5f8299e2b4091p0L },  /* 153/256 */
>   { 0x.9a00000000000000p0L, -0x.821b05f3b01d6774p0L },  /* 154/256 */
>   { 0x.9b00000000000000p0L, -0x.8072d72d903d588cp0L },  /* 155/256 */
>   { 0x.9c00000000000000p0L, -0x.7ecd62bde92210bfp0L },  /* 156/256 */
>   { 0x.9d00000000000000p0L, -0x.7d2a9fb80c64b379p0L },  /* 157/256 */
>   { 0x.9e00000000000000p0L, -0x.7b8a855ad04f93fap0L },  /* 158/256 */
>   { 0x.9f00000000000000p0L, -0x.79ed0b0f76b5cd58p0L },  /* 159/256 */
>   { 0x.a000000000000000p0L, -0x.785228689c9b3653p0L },  /* 160/256 */
>   { 0x.a100000000000000p0L, -0x.76b9d521325856f5p0L },  /* 161/256 */
>   { 0x.a200000000000000p0L, -0x.7524091b7be9add8p0L },  /* 162/256 */
>   { 0x.a300000000000000p0L, -0x.7390bc60191d0d08p0L },  /* 163/256 */
>   { 0x.a400000000000000p0L, -0x.71ffe71d15532491p0L },  /* 164/256 */
>   { 0x.a500000000000000p0L, -0x.707181a4fe8e7640p0L },  /* 165/256 */
>   { 0x.a600000000000000p0L, -0x.6ee5846e038bec2ep0L },  /* 166/256 */
>   { 0x.a700000000000000p0L, -0x.6d5be81118a42549p0L },  /* 167/256 */
>   { 0x.a800000000000000p0L, -0x.6bd4a5492337419dp0L },  /* 168/256 */
>   { 0x.a900000000000000p0L, -0x.6a4fb4f22b678dbdp0L },  /* 169/256 */
>   { 0x.aa00000000000000p0L, -0x.68cd100893e9e323p0L },  /* 170/256 */
>   { 0x.ab00000000000000p0L, -0x.674cafa857b4ec31p0L },  /* 171/256 */
>   { 0x.ac00000000000000p0L, -0x.65ce8d0c4d5ab73bp0L },  /* 172/256 */
>   { 0x.ad00000000000000p0L, -0x.6452a18d6fda2653p0L },  /* 173/256 */
>   { 0x.ae00000000000000p0L, -0x.62d8e6a22cb7d28fp0L },  /* 174/256 */
>   { 0x.af00000000000000p0L, -0x.616155ddb72feab8p0L },  /* 175/256 */
>   { 0x.b000000000000000p0L, -0x.5febe8ef60546fb8p0L },  /* 176/256 */
>   { 0x.b100000000000000p0L, -0x.5e7899a1f3ecf63fp0L },  /* 177/256 */
>   { 0x.b200000000000000p0L, -0x.5d0761db19eec585p0L },  /* 178/256 */
>   { 0x.b300000000000000p0L, -0x.5b983b9abc65c859p0L },  /* 179/256 */
>   { 0x.b400000000000000p0L, -0x.5a2b20fa71a8506ap0L },  /* 180/256 */
>   { 0x.b500000000000000p0L, -0x.58c00c2ceab124eep0L },  /* 181/256 */
>   { 0x.b600000000000000p0L, -0x.5756f77d657cbe9bp0L },  /* 182/256 */
>   { 0x.b700000000000000p0L, -0x.55efdd4f2347eb7bp0L },  /* 183/256 */
>   { 0x.b800000000000000p0L, -0x.548ab81ce28f5f38p0L },  /* 184/256 */
>   { 0x.b900000000000000p0L, -0x.532782785cb0efbbp0L },  /* 185/256 */
>   { 0x.ba00000000000000p0L, -0x.51c63709c7106c19p0L },  /* 186/256 */
>   { 0x.bb00000000000000p0L, -0x.5066d08f57a31c87p0L },  /* 187/256 */
>   { 0x.bc00000000000000p0L, -0x.4f0949dcccc60ed5p0L },  /* 188/256 */
>   { 0x.bd00000000000000p0L, -0x.4dad9ddaf8445bb3p0L },  /* 189/256 */
>   { 0x.be00000000000000p0L, -0x.4c53c7874d738ec3p0L },  /* 190/256 */
>   { 0x.bf00000000000000p0L, -0x.4afbc1f3724d4e7dp0L },  /* 191/256 */
>   { 0x.c000000000000000p0L, -0x.49a58844d36e49e1p0L },  /* 192/256 */
>   { 0x.c100000000000000p0L, -0x.485115b43ae350fcp0L },  /* 193/256 */
>   { 0x.c200000000000000p0L, -0x.46fe658d69ae5377p0L },  /* 194/256 */
>   { 0x.c300000000000000p0L, -0x.45ad732eb3edcd67p0L },  /* 195/256 */
>   { 0x.c400000000000000p0L, -0x.445e3a089f91ef79p0L },  /* 196/256 */
>   { 0x.c500000000000000p0L, -0x.4310b59d858b8c46p0L },  /* 197/256 */
>   { 0x.c600000000000000p0L, -0x.41c4e181356189cep0L },  /* 198/256 */
>   { 0x.c700000000000000p0L, -0x.407ab9589b1a43ddp0L },  /* 199/256 */
>   { 0x.c800000000000000p0L, -0x.3f3238d96766f2fbp0L },  /* 200/256 */
>   { 0x.c900000000000000p0L, -0x.3deb5bc9b9ffcbbep0L },  /* 201/256 */
>   { 0x.ca00000000000000p0L, -0x.3ca61dffce202424p0L },  /* 202/256 */
>   { 0x.cb00000000000000p0L, -0x.3b627b61a912806bp0L },  /* 203/256 */
>   { 0x.cc00000000000000p0L, -0x.3a206fe4cabcf6b0p0L },  /* 204/256 */
>   { 0x.cd00000000000000p0L, -0x.38dff78de01ee139p0L },  /* 205/256 */
>   { 0x.ce00000000000000p0L, -0x.37a10e7077b15a1ep0L },  /* 206/256 */
>   { 0x.cf00000000000000p0L, -0x.3663b0aeb79c794ep0L },  /* 207/256 */
>   { 0x.d000000000000000p0L, -0x.3527da7915b3c6dep0L },  /* 208/256 */
>   { 0x.d100000000000000p0L, -0x.33ed880e112cc827p0L },  /* 209/256 */
>   { 0x.d200000000000000p0L, -0x.32b4b5b9ee02fe45p0L },  /* 210/256 */
>   { 0x.d300000000000000p0L, -0x.317d5fd671fd1855p0L },  /* 211/256 */
>   { 0x.d400000000000000p0L, -0x.304782caa3478377p0L },  /* 212/256 */
>   { 0x.d500000000000000p0L, -0x.2f131b0a8898e67cp0L },  /* 213/256 */
>   { 0x.d600000000000000p0L, -0x.2de02516ead57739p0L },  /* 214/256 */
>   { 0x.d700000000000000p0L, -0x.2cae9d7d182673e3p0L },  /* 215/256 */
>   { 0x.d800000000000000p0L, -0x.2b7e80d6a87b63f7p0L },  /* 216/256 */
>   { 0x.d900000000000000p0L, -0x.2a4fcbc9436b19f4p0L },  /* 217/256 */
>   { 0x.da00000000000000p0L, -0x.29227b06676ac1bdp0L },  /* 218/256 */
>   { 0x.db00000000000000p0L, -0x.27f68b4b32519714p0L },  /* 219/256 */
>   { 0x.dc00000000000000p0L, -0x.26cbf9602b202c5fp0L },  /* 220/256 */
>   { 0x.dd00000000000000p0L, -0x.25a2c2190d0273aep0L },  /* 221/256 */
>   { 0x.de00000000000000p0L, -0x.247ae25493840349p0L },  /* 222/256 */
>   { 0x.df00000000000000p0L, -0x.235456fc47ee53c7p0L },  /* 223/256 */
>   { 0x.e000000000000000p0L, -0x.222f1d044fc8f7bcp0L },  /* 224/256 */
>   { 0x.e100000000000000p0L, -0x.210b316b3c740d11p0L },  /* 225/256 */
>   { 0x.e200000000000000p0L, -0x.1fe89139dbd56595p0L },  /* 226/256 */
>   { 0x.e300000000000000p0L, -0x.1ec739830a111fccp0L },  /* 227/256 */
>   { 0x.e400000000000000p0L, -0x.1da727638446a250p0L },  /* 228/256 */
>   { 0x.e500000000000000p0L, -0x.1c885801bc4b2369p0L },  /* 229/256 */
>   { 0x.e600000000000000p0L, -0x.1b6ac88dad5b1be0p0L },  /* 230/256 */
>   { 0x.e700000000000000p0L, -0x.1a4e7640b1bc37a9p0L },  /* 231/256 */
>   { 0x.e800000000000000p0L, -0x.19335e5d594988aep0L },  /* 232/256 */
>   { 0x.e900000000000000p0L, -0x.18197e2f40e3f01cp0L },  /* 233/256 */
>   { 0x.ea00000000000000p0L, -0x.1700d30aeac0e0f4p0L },  /* 234/256 */
>   { 0x.eb00000000000000p0L, -0x.15e95a4d9791cb7dp0L },  /* 235/256 */
>   { 0x.ec00000000000000p0L, -0x.14d3115d207eac5ep0L },  /* 236/256 */
>   { 0x.ed00000000000000p0L, -0x.13bdf5a7d1ee642fp0L },  /* 237/256 */
>   { 0x.ee00000000000000p0L, -0x.12aa04a44717a48cp0L },  /* 238/256 */
>   { 0x.ef00000000000000p0L, -0x.11973bd1465566d1p0L },  /* 239/256 */
>   { 0x.f000000000000000p0L, -0x.108598b59e3a0689p0L },  /* 240/256 */
>   { 0x.f100000000000000p0L, -0x.f7518e0035c3dd83p-4L },  /* 241/256 */
>   { 0x.f200000000000000p0L, -0x.e65b9e6eed965c37p-4L },  /* 242/256 */
>   { 0x.f300000000000000p0L, -0x.d5779687d887e0d2p-4L },  /* 243/256 */
>   { 0x.f400000000000000p0L, -0x.c4a550a4fd9a19a9p-4L },  /* 244/256 */
>   { 0x.f500000000000000p0L, -0x.b3e4a796a5dac208p-4L },  /* 245/256 */
>   { 0x.f600000000000000p0L, -0x.a33576a16f1f4c64p-4L },  /* 246/256 */
>   { 0x.f700000000000000p0L, -0x.9297997c68c1f4d7p-4L },  /* 247/256 */
>   { 0x.f800000000000000p0L, -0x.820aec4f3a222381p-4L },  /* 248/256 */
>   { 0x.f900000000000000p0L, -0x.718f4bb052abc632p-4L },  /* 249/256 */
>   { 0x.fa00000000000000p0L, -0x.612494a3232afa2ep-4L },  /* 250/256 */
>   { 0x.fb00000000000000p0L, -0x.50caa49660330273p-4L },  /* 251/256 */
>   { 0x.fc00000000000000p0L, -0x.408159624d611d28p-4L },  /* 252/256 */
>   { 0x.fd00000000000000p0L, -0x.3048914711455441p-4L },  /* 253/256 */
>   { 0x.fe00000000000000p0L, -0x.20202aeb11bce252p-4L },  /* 254/256 */
>   { 0x.ff00000000000000p0L, -0x.10080559588b357ep-4L }   /* 255/256 */
>};
>/*
> *  Chebyshev reduced coefficients in power series.
> */
>#define M00  0x.ff805515885e0250p-8L     /*  0.389864041565732301e-2L */
>#define M01  0x.ff00ff00ff0591f1p-8L     /*  0.389105058365765103e-2L */
>#define M02 -0x.7f017e027d037c00p-16L     /* -0.757013732229102636e-5L */
>#define M03  0x.545751e06188aaabp-24L     /*  0.196371909896990662e-7L */
>#define M04 -0x.3f027b08b2000000p-32L     /* -0.573069790783262763e-10L */
>#define M05  0x.327f5b4288888889p-40L     /*  0.179403110200483648e-12L */
>#define M05  0x.327f5b4288888889p-40L     /*  0.179403110200483648e-12L */
>#define M06 -0x.29ae215555555555p-48L     /* -0.578428917628910401e-15L */
>#define M07 -0x.29ae215555555555p-48L     /* -0.578428917628910401e-15L */
>#define M08 -0x.1f00000000000000p-64L     /* -0.656450534122082763e-20L */
>#define M09  0x.1c71c71c71c71c72p-72L     /*  0.235286929792861205e-22L */
>
>static long double zero = 0.0L;
>
>long double logl(long double x) {
>
>   int i, j, n;
>   long double ff, d, t, val;
>
>   /* log(+Inf) = +Inf */
>   if (isinf(x))
>      return x*x+x;
>
>   /* log(NaN) =  NaN with signal */
>   if (isnan(x))
>      return (x+x);
>
>   /* log(x) = sNaN if x < 0. */
>   if (x < zero)
>      return (x - x) / (x - x);
>
>   /* log(+-0) = -Inf with signal */
>   if (x == zero)
>      return (- 1.0L / zero);
>
>   /* Break x into ff * 2**n */
>   ff = frexpl(x, &n);
>
>   /* Find index for the relevant interval */
>   i = 2 * ff * INTERVALS - INTERVALS;
>     
>   d = (ff - f[i][0]) / f[i][0];   /* d in [0,1/INTERVALS] */
>
>   /* Polynomial mapped into in [-1,1] */
>   t = 2 * INTERVALS * d - 1.0L;
>
>   val = M08 + M09 * t;
>   val = M07 + val * t;
>   val = M06 + val * t;
>   val = M05 + val * t;
>   val = M04 + val * t;
>   val = M03 + val * t;
>   val = M02 + val * t;
>   val = M01 + val * t;
>   val = M00 + val * t;
>
>   val += n * LOGL2 + f[i][1];
>
>   return val;
>
>}

Actions: View

Attachments on bug 82654: 54631 | 54632 | 54633 | 54634 | 54635 | 54636