usr/src/lib/libast/common/uwin/log.c

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