Renaming coefficient arrays in cephes

This commit is contained in:
Kent English
2016-11-17 15:41:39 -05:00
parent 5e2f090b11
commit 1dde9ce04a
2 changed files with 14 additions and 37 deletions

View File

@@ -22,7 +22,7 @@ const (
* 1/sqrt(2) <= x < sqrt(2) * 1/sqrt(2) <= x < sqrt(2)
* Theoretical peak relative error = 2.32e-20 * Theoretical peak relative error = 2.32e-20
*/ */
var LP = []float64{ var lP = []float64{
4.5270000862445199635215E-5, 4.5270000862445199635215E-5,
4.9854102823193375972212E-1, 4.9854102823193375972212E-1,
6.5787325942061044846969E0, 6.5787325942061044846969E0,
@@ -32,7 +32,7 @@ var LP = []float64{
2.0039553499201281259648E1, 2.0039553499201281259648E1,
} }
var LQ = []float64{ var lQ = []float64{
1.5062909083469192043167E1, 1.5062909083469192043167E1,
8.3047565967967209469434E1, 8.3047565967967209469434E1,
2.2176239823732856465394E2, 2.2176239823732856465394E2,
@@ -48,7 +48,7 @@ func log1p(x float64) float64 {
return math.Log(z) return math.Log(z)
} }
z = x * x z = x * x
z = -0.5*z + x*(z*polevl(x, LP, 6)/p1evl(x, LQ, 6)) z = -0.5*z + x*(z*polevl(x, lP, 6)/p1evl(x, lQ, 6))
return x + z return x + z
} }
@@ -76,13 +76,13 @@ func log1pmx(x float64) float64 {
* -0.5 <= x <= 0.5 * -0.5 <= x <= 0.5
*/ */
var EP = []float64{ var eP = []float64{
1.2617719307481059087798E-4, 1.2617719307481059087798E-4,
3.0299440770744196129956E-2, 3.0299440770744196129956E-2,
9.9999999999999999991025E-1, 9.9999999999999999991025E-1,
} }
var EQ = []float64{ var eQ = []float64{
3.0019850513866445504159E-6, 3.0019850513866445504159E-6,
2.5244834034968410419224E-3, 2.5244834034968410419224E-3,
2.2726554820815502876593E-1, 2.2726554820815502876593E-1,
@@ -101,8 +101,8 @@ func expm1(x float64) float64 {
return math.Exp(x) - 1 return math.Exp(x) - 1
} }
xx := x * x xx := x * x
r := x * polevl(xx, EP, 2) r := x * polevl(xx, eP, 2)
r = r / (polevl(xx, EQ, 3) - r) r = r / (polevl(xx, eQ, 3) - r)
return r + r return r + r
} }

View File

@@ -13,23 +13,11 @@ package cephes
import "math" import "math"
/* /*
*
* Riemann zeta function of two arguments * Riemann zeta function of two arguments
* *
* *
*
* SYNOPSIS:
*
* double x, q, y, zeta();
*
* y = zeta( x, q );
*
*
*
* DESCRIPTION: * DESCRIPTION:
* *
*
*
* inf. * inf.
* - -x * - -x
* zeta(x,q) = > (k+q) * zeta(x,q) = > (k+q)
@@ -56,11 +44,6 @@ import "math"
* zeta(x,1) = zetac(x) + 1. * zeta(x,1) = zetac(x) + 1.
* *
* *
*
* ACCURACY:
*
*
*
* REFERENCE: * REFERENCE:
* *
* Gradshteyn, I. S., and I. M. Ryzhik, Tables of Integrals, * Gradshteyn, I. S., and I. M. Ryzhik, Tables of Integrals,
@@ -68,18 +51,12 @@ import "math"
* *
*/ */
/*
* Cephes Math Library Release 2.0: April, 1987
* Copyright 1984, 1987 by Stephen L. Moshier
* Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
/* Expansion coefficients /* Expansion coefficients
* for Euler-Maclaurin summation formula * for Euler-Maclaurin summation formula
* (2k)! / B2k * (2k)! / B2k
* where B2k are Bernoulli numbers * where B2k are Bernoulli numbers
*/ */
var A = []float64{ var zetaCoefs = []float64{
12.0, 12.0,
-720.0, -720.0,
30240.0, 30240.0,
@@ -115,18 +92,18 @@ func Zeta(x, q float64) float64 {
} }
/* Asymptotic expansion /* Asymptotic expansion
* http://dlmf.nist.gov/25.11#E43 * http://dlmf.nist.gov/25.11#E43
*/ */
if q > 1e8 { if q > 1e8 {
return (1/(x-1) + 1/(2*q)) * math.Pow(q, 1-x) return (1/(x-1) + 1/(2*q)) * math.Pow(q, 1-x)
} }
/* Euler-Maclaurin summation formula */ // Euler-Maclaurin summation formula
/* Permit negative q but continue sum until n+q > +9 . /* Permit negative q but continue sum until n+q > +9 .
* This case should be handled by a reflection formula. * This case should be handled by a reflection formula.
* If q<0 and x is an integer, there is a relation to * If q<0 and x is an integer, there is a relation to
* the polyGamma function. * the polyGamma function.
*/ */
s := math.Pow(q, -x) s := math.Pow(q, -x)
a := q a := q
@@ -150,7 +127,7 @@ func Zeta(x, q float64) float64 {
for i = 0; i < 12; i++ { for i = 0; i < 12; i++ {
a *= x + k a *= x + k
b /= w b /= w
t := a * b / A[i] t := a * b / zetaCoefs[i]
s = s + t s = s + t
t = math.Abs(t / s) t = math.Abs(t / s)
if t < machEp { if t < machEp {