File: spence.c

package info (click to toggle)
python-scipy 0.18.1-2
  • links: PTS, VCS
  • area: main
  • in suites: stretch
  • size: 75,464 kB
  • ctags: 79,406
  • sloc: python: 143,495; cpp: 89,357; fortran: 81,650; ansic: 79,778; makefile: 364; sh: 265
file content (126 lines) | stat: -rw-r--r-- 2,243 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
/*                                                     spence.c
 *
 *     Dilogarithm
 *
 *
 *
 * SYNOPSIS:
 *
 * double x, y, spence();
 *
 * y = spence( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Computes the integral
 *
 *                    x
 *                    -
 *                   | | log t
 * spence(x)  =  -   |   ----- dt
 *                 | |   t - 1
 *                  -
 *                  1
 *
 * for x >= 0.  A rational approximation gives the integral in
 * the interval (0.5, 1.5).  Transformation formulas for 1/x
 * and 1-x are employed outside the basic expansion range.
 *
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      0,4         30000       3.9e-15     5.4e-16
 *
 *
 */

/*                                                     spence.c */


/*
 * Cephes Math Library Release 2.1:  January, 1989
 * Copyright 1985, 1987, 1989 by Stephen L. Moshier
 * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
 */

#include "mconf.h"

static double A[8] = {
    4.65128586073990045278E-5,
    7.31589045238094711071E-3,
    1.33847639578309018650E-1,
    8.79691311754530315341E-1,
    2.71149851196553469920E0,
    4.25697156008121755724E0,
    3.29771340985225106936E0,
    1.00000000000000000126E0,
};

static double B[8] = {
    6.90990488912553276999E-4,
    2.54043763932544379113E-2,
    2.82974860602568089943E-1,
    1.41172597751831069617E0,
    3.63800533345137075418E0,
    5.03278880143316990390E0,
    3.54771340985225096217E0,
    9.99999999999999998740E-1,
};

extern double MACHEP;

double spence(x)
double x;
{
    double w, y, z;
    int flag;

    if (x < 0.0) {
	mtherr("spence", DOMAIN);
	return (NPY_NAN);
    }

    if (x == 1.0)
	return (0.0);

    if (x == 0.0)
	return (NPY_PI * NPY_PI / 6.0);

    flag = 0;

    if (x > 2.0) {
	x = 1.0 / x;
	flag |= 2;
    }

    if (x > 1.5) {
	w = (1.0 / x) - 1.0;
	flag |= 2;
    }

    else if (x < 0.5) {
	w = -x;
	flag |= 1;
    }

    else
	w = x - 1.0;


    y = -w * polevl(w, A, 7) / polevl(w, B, 7);

    if (flag & 1)
	y = (NPY_PI * NPY_PI) / 6.0 - log(x) * log(1.0 - x) - y;

    if (flag & 2) {
	z = log(x);
	y = -0.5 * z * z - y;
    }

    return (y);
}