/*
* CDDL HEADER START
*
* The contents of this file are subject to the terms of the
* Common Development and Distribution License (the "License").
* You may not use this file except in compliance with the License.
*
* You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
* See the License for the specific language governing permissions
* and limitations under the License.
*
* When distributing Covered Code, include this CDDL HEADER in each
* file and include the License file at usr/src/OPENSOLARIS.LICENSE.
* If applicable, add the following below this CDDL HEADER, with the
* fields enclosed by brackets "[]" replaced with your own identifying
* information: Portions Copyright [yyyy] [name of copyright owner]
*
* CDDL HEADER END
*/
/*
* Copyright 2011 Nexenta Systems, Inc. All rights reserved.
*/
/*
* Copyright 2005 Sun Microsystems, Inc. All rights reserved.
* Use is subject to license terms.
*/
/*
* double __k_lgamma(double x, int *signgamp);
*
* K.C. Ng, March, 1989.
*
* Part of the algorithm is based on W. Cody's lgamma function.
*/
#include "libm.h"
static const double
/*
* Numerator and denominator coefficients for rational minimax Approximation
* P/Q over (0.5,1.5).
*/
/*
* Numerator and denominator coefficients for rational minimax Approximation
* G/H over (1.5,4.0).
*/
/*
* Numerator and denominator coefficients for rational minimax Approximation
* U/V over (4.0,12.0).
*/
/*
* Coefficients for minimax approximation over (12, INF).
*/
/*
* Return sin(pi*x). We assume x is finite and negative, and if it
* is an integer, then the sign of the zero returned doesn't matter.
*/
static double
sin_pi(double x) {
double y, z;
int n;
y = -x;
if (y <= 0.25)
if (y >= two52)
return (zero);
z = floor(y);
if (y == z)
return (zero);
/* argument reduction: set y = |x| mod 2 */
y *= 0.5;
y = 2.0 * (y - floor(y));
/* now floor(y * 4) tells which octant y is in */
n = (int)(y * 4.0);
switch (n) {
case 0:
break;
case 1:
case 2:
break;
case 3:
case 4:
break;
case 5:
case 6:
break;
default:
break;
}
return (-y);
}
static double
double t, p;
/*
* written by K.C. Ng, Feb 2, 1989.
*
* Since
* we have
* G(-z) = -pi/(sin(pi*z)*G(z)*z)
* = pi/(sin(pi*(-z))*G(z)*z)
* Algorithm
* z = |z|
* t = sin_pi(z); ...note that when z>2**52, z is an int
* and hence t=0.
*
* if (t == 0.0) return 1.0/0.0;
* if (t< 0.0) *signgamp = -1; else t= -t;
* if (z+1.0 == 1.0) ...tiny z
* return -log(z);
* else
* return log(pi/(t*z))-__k_lgamma(z, signgamp);
*/
t = sin_pi(z); /* t := sin(pi*z) */
if (t == zero) /* return 1.0/0.0 = +INF */
z = -z;
p = z + one;
if (p == one)
p = -log(z);
else
if (t < zero)
*signgamp = -1;
return (p);
}
double
double t, p, q, cr, y;
/* purge off +-inf, NaN and negative arguments */
if (!finite(x))
return (x * x);
*signgamp = 1;
if (signbit(x))
/* lgamma(x) ~ log(1/x) for really tiny x */
t = one + x;
if (t == one) {
if (x == zero)
return (one / x);
return (-log(x));
}
/* for tiny < x < inf */
if (x <= 1.5) {
if (x < 0.6796875) {
y = x;
} else {
y = x - one;
}
if (x <= 0.5 || x >= 0.6796875) {
if (x == one)
return (zero);
(q7+y)))))));
} else {
y = x - one;
(h7+y)))))));
}
} else if (x <= 4.0) {
if (x == 2.0)
return (zero);
y = x - 2.0;
return (y*(D2+y*(p/q)));
} else if (x <= 12.0) {
y = x - 4.0;
return (D4+y*(p/q));
} else if (x <= 1.0e17) { /* x ~< 2**(prec+3) */
t = one / x;
y = t * t;
q = log(x);
return (x*(q-one)-(0.5*q-p));
} else { /* may overflow */
return (x * (log(x) - 1.0));
}
}