lib/libm/noieee_src/n_log.c

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