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