lib/libm/noieee_src/n_log.c

   1 /*      $NetBSD: n_log.c,v 1.8 2014/10/10 20:58:09 martin Exp $ */
   2 /*
   3  * Copyright (c) 1992, 1993
   4  *      The Regents of the University of California.  All rights reserved.
   5  *
   6  * Redistribution and use in source and binary forms, with or without
   7  * modification, are permitted provided that the following conditions
   8  * are met:
   9  * 1. Redistributions of source code must retain the above copyright
  10  *    notice, this list of conditions and the following disclaimer.
  11  * 2. Redistributions in binary form must reproduce the above copyright
  12  *    notice, this list of conditions and the following disclaimer in the
  13  *    documentation and/or other materials provided with the distribution.
  14  * 3. Neither the name of the University nor the names of its contributors
  15  *    may be used to endorse or promote products derived from this software
  16  *    without specific prior written permission.
  17  *
  18  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  19  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  20  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  21  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  22  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  23  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  24  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  25  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  26  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  27  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  28  * SUCH DAMAGE.
  29  */
  30
  31 #ifndef lint
  32 #if 0
  33 static char sccsid[] = "@(#)log.c       8.2 (Berkeley) 11/30/93";
  34 #endif
  35 #endif /* not lint */
  36
  37 #include "../src/namespace.h"
  38
  39 #include <math.h>
  40 #include <errno.h>
  41
  42 #include "mathimpl.h"
  43
  44 #ifdef __weak_alias
  45 __weak_alias(log, _log);
  46 __weak_alias(_logl, _log);
  47 __weak_alias(logf, _logf);
  48 #endif
  49
  50 /* Table-driven natural logarithm.
  51  *
  52  * This code was derived, with minor modifications, from:
  53  *      Peter Tang, "Table-Driven Implementation of the
  54  *      Logarithm in IEEE Floating-Point arithmetic." ACM Trans.
  55  *      Math Software, vol 16. no 4, pp 378-400, Dec 1990).
  56  *
  57  * Calculates log(2^m*F*(1+f/F)), |f/j| <= 1/256,
  58  * where F = j/128 for j an integer in [0, 128].
  59  *
  60  * log(2^m) = log2_hi*m + log2_tail*m
  61  * since m is an integer, the dominant term is exact.
  62  * m has at most 10 digits (for subnormal numbers),
  63  * and log2_hi has 11 trailing zero bits.
  64  *
  65  * log(F) = logF_hi[j] + logF_lo[j] is in tabular form in log_table.h
  66  * logF_hi[] + 512 is exact.
  67  *
  68  * log(1+f/F) = 2*f/(2*F + f) + 1/12 * (2*f/(2*F + f))**3 + ...
  69  * the leading term is calculated to extra precision in two
  70  * parts, the larger of which adds exactly to the dominant
  71  * m and F terms.
  72  * There are two cases:
  73  *      1. when m, j are non-zero (m | j), use absolute
  74  *         precision for the leading term.
  75  *      2. when m = j = 0, |1-x| < 1/256, and log(x) ~= (x-1).
  76  *         In this case, use a relative precision of 24 bits.
  77  * (This is done differently in the original paper)
  78  *
  79  * Special cases:
  80  *      0       return signalling -Inf
  81  *      neg     return signalling NaN
  82  *      +Inf    return +Inf
  83 */
  84
  85 #if defined(__vax__) || defined(tahoe)
  86 #define _IEEE           0
  87 #define TRUNC(x)        x = (double) (float) (x)
  88 #else
  89 #define _IEEE           1
  90 #define endian          (((*(int *) &one)) ? 1 : 0)
  91 #define TRUNC(x)        *(((int *) &x) + endian) &= 0xf8000000
  92 #define infnan(x)       0.0
  93 #endif
  94
  95 #define N 128
  96
  97 /* Table of log(Fj) = logF_head[j] + logF_tail[j], for Fj = 1+j/128.
  98  * Used for generation of extend precision logarithms.
  99  * The constant 35184372088832 is 2^45, so the divide is exact.
 100  * It ensures correct reading of logF_head, even for inaccurate
 101  * decimal-to-binary conversion routines.  (Everybody gets the
 102  * right answer for integers less than 2^53.)
 103  * Values for log(F) were generated using error < 10^-57 absolute
 104  * with the bc -l package.
 105 */
 106 static const double     A1 =      .08333333333333178827;
 107 static const double     A2 =      .01250000000377174923;
 108 static const double     A3 =     .002232139987919447809;
 109 static const double     A4 =    .0004348877777076145742;
 110
 111 static const double logF_head[N+1] = {
 112         0.,
 113         .007782140442060381246,
 114         .015504186535963526694,
 115         .023167059281547608406,
 116         .030771658666765233647,
 117         .038318864302141264488,
 118         .045809536031242714670,
 119         .053244514518837604555,
 120         .060624621816486978786,
 121         .067950661908525944454,
 122         .075223421237524235039,
 123         .082443669210988446138,
 124         .089612158689760690322,
 125         .096729626458454731618,
 126         .103796793681567578460,
 127         .110814366340264314203,
 128         .117783035656430001836,
 129         .124703478501032805070,
 130         .131576357788617315236,
 131         .138402322859292326029,
 132         .145182009844575077295,
 133         .151916042025732167530,
 134         .158605030176659056451,
 135         .165249572895390883786,
 136         .171850256926518341060,
 137         .178407657472689606947,
 138         .184922338493834104156,
 139         .191394852999565046047,
 140         .197825743329758552135,
 141         .204215541428766300668,
 142         .210564769107350002741,
 143         .216873938300523150246,
 144         .223143551314024080056,
 145         .229374101064877322642,
 146         .235566071312860003672,
 147         .241719936886966024758,
 148         .247836163904594286577,
 149         .253915209980732470285,
 150         .259957524436686071567,
 151         .265963548496984003577,
 152         .271933715484010463114,
 153         .277868451003087102435,
 154         .283768173130738432519,
 155         .289633292582948342896,
 156         .295464212893421063199,
 157         .301261330578199704177,
 158         .307025035294827830512,
 159         .312755710004239517729,
 160         .318453731118097493890,
 161         .324119468654316733591,
 162         .329753286372579168528,
 163         .335355541920762334484,
 164         .340926586970454081892,
 165         .346466767346100823488,
 166         .351976423156884266063,
 167         .357455888922231679316,
 168         .362905493689140712376,
 169         .368325561158599157352,
 170         .373716409793814818840,
 171         .379078352934811846353,
 172         .384411698910298582632,
 173         .389716751140440464951,
 174         .394993808240542421117,
 175         .400243164127459749579,
 176         .405465108107819105498,
 177         .410659924985338875558,
 178         .415827895143593195825,
 179         .420969294644237379543,
 180         .426084395310681429691,
 181         .431173464818130014464,
 182         .436236766774527495726,
 183         .441274560805140936281,
 184         .446287102628048160113,
 185         .451274644139630254358,
 186         .456237433481874177232,
 187         .461175715122408291790,
 188         .466089729924533457960,
 189         .470979715219073113985,
 190         .475845904869856894947,
 191         .480688529345570714212,
 192         .485507815781602403149,
 193         .490303988045525329653,
 194         .495077266798034543171,
 195         .499827869556611403822,
 196         .504556010751912253908,
 197         .509261901790523552335,
 198         .513945751101346104405,
 199         .518607764208354637958,
 200         .523248143765158602036,
 201         .527867089620485785417,
 202         .532464798869114019908,
 203         .537041465897345915436,
 204         .541597282432121573947,
 205         .546132437597407260909,
 206         .550647117952394182793,
 207         .555141507540611200965,
 208         .559615787935399566777,
 209         .564070138285387656651,
 210         .568504735352689749561,
 211         .572919753562018740922,
 212         .577315365035246941260,
 213         .581691739635061821900,
 214         .586049045003164792433,
 215         .590387446602107957005,
 216         .594707107746216934174,
 217         .599008189645246602594,
 218         .603290851438941899687,
 219         .607555250224322662688,
 220         .611801541106615331955,
 221         .616029877215623855590,
 222         .620240409751204424537,
 223         .624433288012369303032,
 224         .628608659422752680256,
 225         .632766669570628437213,
 226         .636907462236194987781,
 227         .641031179420679109171,
 228         .645137961373620782978,
 229         .649227946625615004450,
 230         .653301272011958644725,
 231         .657358072709030238911,
 232         .661398482245203922502,
 233         .665422632544505177065,
 234         .669430653942981734871,
 235         .673422675212350441142,
 236         .677398823590920073911,
 237         .681359224807238206267,
 238         .685304003098281100392,
 239         .689233281238557538017,
 240         .693147180560117703862
 241 };
 242
 243 static const double logF_tail[N+1] = {
 244         0.,
 245         -.00000000000000543229938420049,
 246          .00000000000000172745674997061,
 247         -.00000000000001323017818229233,
 248         -.00000000000001154527628289872,
 249         -.00000000000000466529469958300,
 250          .00000000000005148849572685810,
 251         -.00000000000002532168943117445,
 252         -.00000000000005213620639136504,
 253         -.00000000000001819506003016881,
 254          .00000000000006329065958724544,
 255          .00000000000008614512936087814,
 256         -.00000000000007355770219435028,
 257          .00000000000009638067658552277,
 258          .00000000000007598636597194141,
 259          .00000000000002579999128306990,
 260         -.00000000000004654729747598444,
 261         -.00000000000007556920687451336,
 262          .00000000000010195735223708472,
 263         -.00000000000017319034406422306,
 264         -.00000000000007718001336828098,
 265          .00000000000010980754099855238,
 266         -.00000000000002047235780046195,
 267         -.00000000000008372091099235912,
 268          .00000000000014088127937111135,
 269          .00000000000012869017157588257,
 270          .00000000000017788850778198106,
 271          .00000000000006440856150696891,
 272          .00000000000016132822667240822,
 273         -.00000000000007540916511956188,
 274         -.00000000000000036507188831790,
 275          .00000000000009120937249914984,
 276          .00000000000018567570959796010,
 277         -.00000000000003149265065191483,
 278         -.00000000000009309459495196889,
 279          .00000000000017914338601329117,
 280         -.00000000000001302979717330866,
 281          .00000000000023097385217586939,
 282          .00000000000023999540484211737,
 283          .00000000000015393776174455408,
 284         -.00000000000036870428315837678,
 285          .00000000000036920375082080089,
 286         -.00000000000009383417223663699,
 287          .00000000000009433398189512690,
 288          .00000000000041481318704258568,
 289         -.00000000000003792316480209314,
 290          .00000000000008403156304792424,
 291         -.00000000000034262934348285429,
 292          .00000000000043712191957429145,
 293         -.00000000000010475750058776541,
 294         -.00000000000011118671389559323,
 295          .00000000000037549577257259853,
 296          .00000000000013912841212197565,
 297          .00000000000010775743037572640,
 298          .00000000000029391859187648000,
 299         -.00000000000042790509060060774,
 300          .00000000000022774076114039555,
 301          .00000000000010849569622967912,
 302         -.00000000000023073801945705758,
 303          .00000000000015761203773969435,
 304          .00000000000003345710269544082,
 305         -.00000000000041525158063436123,
 306          .00000000000032655698896907146,
 307         -.00000000000044704265010452446,
 308          .00000000000034527647952039772,
 309         -.00000000000007048962392109746,
 310          .00000000000011776978751369214,
 311         -.00000000000010774341461609578,
 312          .00000000000021863343293215910,
 313          .00000000000024132639491333131,
 314          .00000000000039057462209830700,
 315         -.00000000000026570679203560751,
 316          .00000000000037135141919592021,
 317         -.00000000000017166921336082431,
 318         -.00000000000028658285157914353,
 319         -.00000000000023812542263446809,
 320          .00000000000006576659768580062,
 321         -.00000000000028210143846181267,
 322          .00000000000010701931762114254,
 323          .00000000000018119346366441110,
 324          .00000000000009840465278232627,
 325         -.00000000000033149150282752542,
 326         -.00000000000018302857356041668,
 327         -.00000000000016207400156744949,
 328          .00000000000048303314949553201,
 329         -.00000000000071560553172382115,
 330          .00000000000088821239518571855,
 331         -.00000000000030900580513238244,
 332         -.00000000000061076551972851496,
 333          .00000000000035659969663347830,
 334          .00000000000035782396591276383,
 335         -.00000000000046226087001544578,
 336          .00000000000062279762917225156,
 337          .00000000000072838947272065741,
 338          .00000000000026809646615211673,
 339         -.00000000000010960825046059278,
 340          .00000000000002311949383800537,
 341         -.00000000000058469058005299247,
 342         -.00000000000002103748251144494,
 343         -.00000000000023323182945587408,
 344         -.00000000000042333694288141916,
 345         -.00000000000043933937969737844,
 346          .00000000000041341647073835565,
 347          .00000000000006841763641591466,
 348          .00000000000047585534004430641,
 349          .00000000000083679678674757695,
 350         -.00000000000085763734646658640,
 351          .00000000000021913281229340092,
 352         -.00000000000062242842536431148,
 353         -.00000000000010983594325438430,
 354          .00000000000065310431377633651,
 355         -.00000000000047580199021710769,
 356         -.00000000000037854251265457040,
 357          .00000000000040939233218678664,
 358          .00000000000087424383914858291,
 359          .00000000000025218188456842882,
 360         -.00000000000003608131360422557,
 361         -.00000000000050518555924280902,
 362          .00000000000078699403323355317,
 363         -.00000000000067020876961949060,
 364          .00000000000016108575753932458,
 365          .00000000000058527188436251509,
 366         -.00000000000035246757297904791,
 367         -.00000000000018372084495629058,
 368          .00000000000088606689813494916,
 369          .00000000000066486268071468700,
 370          .00000000000063831615170646519,
 371          .00000000000025144230728376072,
 372         -.00000000000017239444525614834
 373 };
 374
 375 double
 376 log(double x)
 377 {
 378         int m, j;
 379         double F, f, g, q, u, u2, v, zero = 0.0, one = 1.0;
 380         volatile double u1;
 381
 382         /* Catch special cases */
 383         if (x <= 0) {
 384                 if (_IEEE && x == zero) /* log(0) = -Inf */
 385                         return (-one/zero);
 386                 else if (_IEEE)         /* log(neg) = NaN */
 387                         return (zero/zero);
 388                 else if (x == zero)     /* NOT REACHED IF _IEEE */
 389                         return (infnan(-ERANGE));
 390                 else
 391                         return (infnan(EDOM));
 392         } else if (!finite(x)) {
 393                 if (_IEEE)              /* x = NaN, Inf */
 394                         return (x+x);
 395                 else
 396                         return (infnan(ERANGE));
 397         }
 398
 399         /* Argument reduction: 1 <= g < 2; x/2^m = g;   */
 400         /* y = F*(1 + f/F) for |f| <= 2^-8              */
 401
 402         m = logb(x);
 403         g = ldexp(x, -m);
 404         if (_IEEE && m == -1022) {
 405                 j = logb(g), m += j;
 406                 g = ldexp(g, -j);
 407         }
 408         j = N*(g-1) + .5;
 409         F = (1.0/N) * j + 1;    /* F*128 is an integer in [128, 512] */
 410         f = g - F;
 411
 412         /* Approximate expansion for log(1+f/F) ~= u + q */
 413         g = 1/(2*F+f);
 414         u = 2*f*g;
 415         v = u*u;
 416         q = u*v*(A1 + v*(A2 + v*(A3 + v*A4)));
 417
 418     /* case 1: u1 = u rounded to 2^-43 absolute.  Since u < 2^-8,
 419      *         u1 has at most 35 bits, and F*u1 is exact, as F has < 8 bits.
 420      *         It also adds exactly to |m*log2_hi + log_F_head[j] | < 750
 421     */
 422         if (m | j)
 423                 u1 = u + 513, u1 -= 513;
 424
 425     /* case 2:  |1-x| < 1/256. The m- and j- dependent terms are zero;
 426      *          u1 = u to 24 bits.
 427     */
 428         else
 429                 u1 = u, TRUNC(u1);
 430         u2 = (2.0*(f - F*u1) - u1*f) * g;
 431                         /* u1 + u2 = 2f/(2F+f) to extra precision.      */
 432
 433         /* log(x) = log(2^m*F*(1+f/F)) =                                */
 434         /* (m*log2_hi+logF_head[j]+u1) + (m*log2_lo+logF_tail[j]+q);    */
 435         /* (exact) + (tiny)                                             */
 436
 437         u1 += m*logF_head[N] + logF_head[j];            /* exact */
 438         u2 = (u2 + logF_tail[j]) + q;                   /* tiny */
 439         u2 += logF_tail[N]*m;
 440         return (u1 + u2);
 441 }
 442
 443 /*
 444  * Extra precision variant, returning struct {double a, b;};
 445  * log(x) = a+b to 63 bits, with a is rounded to 26 bits.
 446  */
 447 struct Double
 448 __log__D(double x)
 449 {
 450         int m, j;
 451         double F, f, g, q, u, v, u2;
 452         volatile double u1;
 453         struct Double r;
 454
 455         /* Argument reduction: 1 <= g < 2; x/2^m = g;   */
 456         /* y = F*(1 + f/F) for |f| <= 2^-8              */
 457
 458         m = logb(x);
 459         g = ldexp(x, -m);
 460         if (_IEEE && m == -1022) {
 461                 j = logb(g), m += j;
 462                 g = ldexp(g, -j);
 463         }
 464         j = N*(g-1) + .5;
 465         F = (1.0/N) * j + 1;
 466         f = g - F;
 467
 468         g = 1/(2*F+f);
 469         u = 2*f*g;
 470         v = u*u;
 471         q = u*v*(A1 + v*(A2 + v*(A3 + v*A4)));
 472         if (m | j)
 473                 u1 = u + 513, u1 -= 513;
 474         else
 475                 u1 = u, TRUNC(u1);
 476         u2 = (2.0*(f - F*u1) - u1*f) * g;
 477
 478         u1 += m*logF_head[N] + logF_head[j];
 479
 480         u2 +=  logF_tail[j]; u2 += q;
 481         u2 += logF_tail[N]*m;
 482         r.a = u1 + u2;                  /* Only difference is here */
 483         TRUNC(r.a);
 484         r.b = (u1 - r.a) + u2;
 485         return (r);
 486 }
 487
 488 float
 489 logf(float x)
 490 {
 491         return(log((double)x));
 492 }