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

   1 #include "FEATURE/uwin"
   2
   3 #if !_UWIN || _lib_lgamma
   4
   5 void _STUB_lgamma(){}
   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[] = "@(#)lgamma.c    8.2 (Berkeley) 11/30/93";
  39
  40 /*
  41  * Coded by Peter McIlroy, Nov 1992;
  42  *
  43  * The financial support of UUNET Communications Services is greatfully
  44  * acknowledged.
  45  */
  46
  47 #define gamma   ______gamma
  48 #define lgamma  ______lgamma
  49
  50 #include <math.h>
  51 #include <errno.h>
  52 #include "mathimpl.h"
  53
  54 #undef  gamma
  55 #undef  lgamma
  56
  57 /* Log gamma function.
  58  * Error:  x > 0 error < 1.3ulp.
  59  *         x > 4, error < 1ulp.
  60  *         x > 9, error < .6ulp.
  61  *         x < 0, all bets are off. (When G(x) ~ 1, log(G(x)) ~ 0)
  62  * Method:
  63  *      x > 6:
  64  *              Use the asymptotic expansion (Stirling's Formula)
  65  *      0 < x < 6:
  66  *              Use gamma(x+1) = x*gamma(x) for argument reduction.
  67  *              Use rational approximation in
  68  *              the range 1.2, 2.5
  69  *              Two approximations are used, one centered at the
  70  *              minimum to ensure monotonicity; one centered at 2
  71  *              to maintain small relative error.
  72  *      x < 0:
  73  *              Use the reflection formula,
  74  *              G(1-x)G(x) = PI/sin(PI*x)
  75  * Special values:
  76  *      non-positive integer    returns +Inf.
  77  *      NaN                     returns NaN
  78 */
  79 static int endian;
  80 #if defined(vax) || defined(tahoe)
  81 #define _IEEE           0
  82 /* double and float have same size exponent field */
  83 #define TRUNC(x)        x = (double) (float) (x)
  84 #else
  85 #define _IEEE           1
  86 #define TRUNC(x)        *(((int *) &x) + endian) &= 0xf8000000
  87 #define infnan(x)       0.0
  88 #endif
  89
  90 static double small_lgam(double);
  91 static double large_lgam(double);
  92 static double neg_lgam(double);
  93 static double zero = 0.0, one = 1.0;
  94 int signgam;
  95
  96 #define UNDERFL (1e-1020 * 1e-1020)
  97
  98 #define LEFT    (1.0 - (x0 + .25))
  99 #define RIGHT   (x0 - .218)
 100 /*
 101  * Constants for approximation in [1.244,1.712]
 102 */
 103 #define x0      0.461632144968362356785
 104 #define x0_lo   -.000000000000000015522348162858676890521
 105 #define a0_hi   -0.12148629128932952880859
 106 #define a0_lo   .0000000007534799204229502
 107 #define r0      -2.771227512955130520e-002
 108 #define r1      -2.980729795228150847e-001
 109 #define r2      -3.257411333183093394e-001
 110 #define r3      -1.126814387531706041e-001
 111 #define r4      -1.129130057170225562e-002
 112 #define r5      -2.259650588213369095e-005
 113 #define s0       1.714457160001714442e+000
 114 #define s1       2.786469504618194648e+000
 115 #define s2       1.564546365519179805e+000
 116 #define s3       3.485846389981109850e-001
 117 #define s4       2.467759345363656348e-002
 118 /*
 119  * Constants for approximation in [1.71, 2.5]
 120 */
 121 #define a1_hi   4.227843350984671344505727574870e-01
 122 #define a1_lo   4.670126436531227189e-18
 123 #define p0      3.224670334241133695662995251041e-01
 124 #define p1      3.569659696950364669021382724168e-01
 125 #define p2      1.342918716072560025853732668111e-01
 126 #define p3      1.950702176409779831089963408886e-02
 127 #define p4      8.546740251667538090796227834289e-04
 128 #define q0      1.000000000000000444089209850062e+00
 129 #define q1      1.315850076960161985084596381057e+00
 130 #define q2      6.274644311862156431658377186977e-01
 131 #define q3      1.304706631926259297049597307705e-01
 132 #define q4      1.102815279606722369265536798366e-02
 133 #define q5      2.512690594856678929537585620579e-04
 134 #define q6      -1.003597548112371003358107325598e-06
 135 /*
 136  * Stirling's Formula, adjusted for equal-ripple. x in [6,Inf].
 137 */
 138 #define lns2pi  .418938533204672741780329736405
 139 #define pb0      8.33333333333333148296162562474e-02
 140 #define pb1     -2.77777777774548123579378966497e-03
 141 #define pb2      7.93650778754435631476282786423e-04
 142 #define pb3     -5.95235082566672847950717262222e-04
 143 #define pb4      8.41428560346653702135821806252e-04
 144 #define pb5     -1.89773526463879200348872089421e-03
 145 #define pb6      5.69394463439411649408050664078e-03
 146 #define pb7     -1.44705562421428915453880392761e-02
 147
 148 extern __pure double lgamma(double x)
 149 {
 150         double r;
 151
 152         signgam = 1;
 153         endian = ((*(int *) &one)) ? 1 : 0;
 154
 155         if (!finite(x))
 156                 if (_IEEE)
 157                         return (x+x);
 158                 else return (infnan(EDOM));
 159
 160         if (x > 6 + RIGHT) {
 161                 r = large_lgam(x);
 162                 return (r);
 163         } else if (x > 1e-16)
 164                 return (small_lgam(x));
 165         else if (x > -1e-16) {
 166                 if (x < 0)
 167                         signgam = -1, x = -x;
 168                 return (-log(x));
 169         } else
 170                 return (neg_lgam(x));
 171 }
 172
 173 static double
 174 large_lgam(double x)
 175 {
 176         double z, p, x1;
 177         struct Double t, u, v;
 178         u = __log__D(x);
 179         u.a -= 1.0;
 180         if (x > 1e15) {
 181                 v.a = x - 0.5;
 182                 TRUNC(v.a);
 183                 v.b = (x - v.a) - 0.5;
 184                 t.a = u.a*v.a;
 185                 t.b = x*u.b + v.b*u.a;
 186                 if (_IEEE == 0 && !finite(t.a))
 187                         return(infnan(ERANGE));
 188                 return(t.a + t.b);
 189         }
 190         x1 = 1./x;
 191         z = x1*x1;
 192         p = pb0+z*(pb1+z*(pb2+z*(pb3+z*(pb4+z*(pb5+z*(pb6+z*pb7))))));
 193                                         /* error in approximation = 2.8e-19 */
 194
 195         p = p*x1;                       /* error < 2.3e-18 absolute */
 196                                         /* 0 < p < 1/64 (at x = 5.5) */
 197         v.a = x = x - 0.5;
 198         TRUNC(v.a);                     /* truncate v.a to 26 bits. */
 199         v.b = x - v.a;
 200         t.a = v.a*u.a;                  /* t = (x-.5)*(log(x)-1) */
 201         t.b = v.b*u.a + x*u.b;
 202         t.b += p; t.b += lns2pi;        /* return t + lns2pi + p */
 203         return (t.a + t.b);
 204 }
 205
 206 static double
 207 small_lgam(double x)
 208 {
 209         int x_int;
 210         double y, z, t, r = 0, p, q, hi, lo;
 211         struct Double rr;
 212         x_int = (int)(x + .5);
 213         y = x - x_int;
 214         if (x_int <= 2 && y > RIGHT) {
 215                 t = y - x0;
 216                 y--; x_int++;
 217                 goto CONTINUE;
 218         } else if (y < -LEFT) {
 219                 t = y +(1.0-x0);
 220 CONTINUE:
 221                 z = t - x0_lo;
 222                 p = r0+z*(r1+z*(r2+z*(r3+z*(r4+z*r5))));
 223                 q = s0+z*(s1+z*(s2+z*(s3+z*s4)));
 224                 r = t*(z*(p/q) - x0_lo);
 225                 t = .5*t*t;
 226                 z = 1.0;
 227                 switch (x_int) {
 228                 case 6: z  = (y + 5);
 229                 case 5: z *= (y + 4);
 230                 case 4: z *= (y + 3);
 231                 case 3: z *= (y + 2);
 232                         rr = __log__D(z);
 233                         rr.b += a0_lo; rr.a += a0_hi;
 234                         return(((r+rr.b)+t+rr.a));
 235                 case 2: return(((r+a0_lo)+t)+a0_hi);
 236                 case 0: r -= log1p(x);
 237                 default: rr = __log__D(x);
 238                         rr.a -= a0_hi; rr.b -= a0_lo;
 239                         return(((r - rr.b) + t) - rr.a);
 240                 }
 241         } else {
 242                 p = p0+y*(p1+y*(p2+y*(p3+y*p4)));
 243                 q = q0+y*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*q6)))));
 244                 p = p*(y/q);
 245                 t = (double)(float) y;
 246                 z = y-t;
 247                 hi = (double)(float) (p+a1_hi);
 248                 lo = a1_hi - hi; lo += p; lo += a1_lo;
 249                 r = lo*y + z*hi;        /* q + r = y*(a0+p/q) */
 250                 q = hi*t;
 251                 z = 1.0;
 252                 switch (x_int) {
 253                 case 6: z  = (y + 5);
 254                 case 5: z *= (y + 4);
 255                 case 4: z *= (y + 3);
 256                 case 3: z *= (y + 2);
 257                         rr = __log__D(z);
 258                         r += rr.b; r += q;
 259                         return(rr.a + r);
 260                 case 2: return (q+ r);
 261                 case 0: rr = __log__D(x);
 262                         r -= rr.b; r -= log1p(x);
 263                         r += q; r-= rr.a;
 264                         return(r);
 265                 default: rr = __log__D(x);
 266                         r -= rr.b;
 267                         q -= rr.a;
 268                         return (r+q);
 269                 }
 270         }
 271 }
 272
 273 static double
 274 neg_lgam(double x)
 275 {
 276         int xi;
 277         double y, z, one = 1.0, zero = 0.0;
 278         extern double gamma();
 279
 280         /* avoid destructive cancellation as much as possible */
 281         if (x > -170) {
 282                 xi = (int)x;
 283                 if (xi == x)
 284                         if (_IEEE)
 285                                 return(one/zero);
 286                         else
 287                                 return(infnan(ERANGE));
 288                 y = gamma(x);
 289                 if (y < 0)
 290                         y = -y, signgam = -1;
 291                 return (log(y));
 292         }
 293         z = floor(x + .5);
 294         if (z == x) {           /* convention: G(-(integer)) -> +Inf */
 295                 if (_IEEE)
 296                         return (one/zero);
 297                 else
 298                         return (infnan(ERANGE));
 299         }
 300         y = .5*ceil(x);
 301         if (y == ceil(y))
 302                 signgam = -1;
 303         x = -x;
 304         z = fabs(x + z);        /* 0 < z <= .5 */
 305         if (z < .25)
 306                 z = sin(M_PI*z);
 307         else
 308                 z = cos(M_PI*(0.5-z));
 309         z = log(M_PI/(z*x));
 310         y = large_lgam(x);
 311         return (z - y);
 312 }
 313
 314 #endif