Skip to content

Commit 0e81579

Browse files
MikaelSlevinskyararslan
MikaelSlevinsky
authored and
ararslan
committed
add sine and cosine integrals (#32)
1 parent 533a549 commit 0e81579

File tree

6 files changed

+310
-10
lines changed

6 files changed

+310
-10
lines changed

‎README.md

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,8 @@
11
# SpecialFunctions.jl
22

3-
Special mathematical functions in Julia, including Bessel, Hankel, Airy, error, Dawson, eta, zeta,
4-
digamma, inverse digamma, trigamma, and polygamma functions.
5-
These functions were formerly part of Base.
3+
Special mathematical functions in Julia, including Bessel, Hankel, Airy, error, Dawson, sine and cosine integrals,
4+
eta, zeta, digamma, inverse digamma, trigamma, and polygamma functions.
5+
Most of these functions were formerly part of Base.
66

77
[![Travis](https://travis-ci.org/JuliaMath/SpecialFunctions.jl.svg?branch=master)](https://travis-ci.org/JuliaMath/SpecialFunctions.jl)
88
[![AppVeyor](https://ci.appveyor.com/api/projects/status/ccfgkm2cjcggu158/branch/master?svg=true)](https://ci.appveyor.com/project/ararslan/specialfunctions-jl/branch/master)

‎docs/src/index.md

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,8 @@ libraries.
1313
| [`erfi(x)`](@ref SpecialFunctions.erfi) | imaginary error function defined as `-im * erf(x * im)`, where `im` is the imaginary unit |
1414
| [`erfcx(x)`](@ref SpecialFunctions.erfcx) | scaled complementary error function, i.e. accurate `exp(x^2) * erfc(x)` for large `x` |
1515
| [`dawson(x)`](@ref SpecialFunctions.dawson) | scaled imaginary error function, a.k.a. Dawson function, i.e. accurate `exp(-x^2) * erfi(x) * sqrt(pi) / 2` for large `x` |
16+
| [`sinint(x)`](@ref SpecialFunctions.sinint) | [sine integral](https://en.wikipedia.org/wiki/Trigonometric_integral) at `x` |
17+
| [`cosint(x)`](@ref SpecialFunctions.cosint) | [cosine integral](https://en.wikipedia.org/wiki/Trigonometric_integral) at `x` |
1618
| [`digamma(x)`](@ref SpecialFunctions.digamma) | [digamma function](https://en.wikipedia.org/wiki/Digamma_function) (i.e. the derivative of `lgamma` at `x`) |
1719
| [`eta(x)`](@ref SpecialFunctions.eta) | [Dirichlet eta function](https://en.wikipedia.org/wiki/Dirichlet_eta_function) at `x` |
1820
| [`zeta(x)`](@ref SpecialFunctions.zeta) | [Riemann zeta function](https://en.wikipedia.org/wiki/Riemann_zeta_function) at `x` |
@@ -51,7 +53,7 @@ from the Julia REPL.
5153

5254
## Note
5355

54-
Prior to Julia 0.6, these functions were available in Julia's Base module.
56+
Prior to Julia 0.6, most of these functions were available in Julia's Base module.
5557
Because of this, the symbols from this package are not exported on Julia 0.5
5658
to avoid name conflicts.
5759
In this case, the symbols will need to be explicitly imported or called

‎docs/src/special.md

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,8 @@ SpecialFunctions.erfi
1212
SpecialFunctions.dawson
1313
SpecialFunctions.erfinv
1414
SpecialFunctions.erfcinv
15+
SpecialFunctions.sinint
16+
SpecialFunctions.cosint
1517
SpecialFunctions.digamma
1618
SpecialFunctions.invdigamma
1719
SpecialFunctions.trigamma

‎src/SpecialFunctions.jl

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -59,6 +59,9 @@ if VERSION >= v"0.6.0-dev.2767"
5959
end
6060
end
6161

62+
export sinint,
63+
cosint
64+
6265
if isdefined(Base.Math, :openspecfun)
6366
const openspecfun = Base.Math.openspecfun
6467
else
@@ -67,6 +70,7 @@ end
6770

6871
include("bessel.jl")
6972
include("erf.jl")
73+
include("sincosint.jl")
7074
include("gamma.jl")
7175
include("deprecated.jl")
7276

‎src/sincosint.jl

Lines changed: 258 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,258 @@
1+
using Base.Math.@horner
2+
3+
# Compute the sine integral: ∫_0^x sin(t)/t dt,
4+
# and the cosine integral: γ + log x + ∫_0^x (cos(t)-1)/t dt,
5+
# using the rational approximants tabulated in:
6+
# A.J. MacLeod, "Rational approximations, software and test methods for
7+
# sine and cosine integrals," Numer. Algor. 12, pp. 259--272 (1996).
8+
# http://dx.doi.org/10.1007/BF02142806
9+
# https://link.springer.com/article/10.1007/BF02142806
10+
#
11+
# Note: the second zero of Ci(x) has a typo that is fixed:
12+
#
13+
# r1 = 3.38418 04228 51186 42639 78511 46402 in the article, but is in fact:
14+
#
15+
# r1 = 3.38418 04225 51186 42639 78511 46402.
16+
#
17+
18+
function sinint(x::Float64)
19+
t = x*x
20+
if t ≤ 36.0
21+
return x * @horner(t, 1.00000000000000000000E0,
22+
-0.44663998931312457298E-1,
23+
0.11209146443112369449E-2,
24+
-0.13276124407928422367E-4,
25+
0.85118014179823463879E-7,
26+
-0.29989314303147656479E-9,
27+
0.55401971660186204711E-12,
28+
-0.42406353433133212926E-15) /
29+
@horner(t, 1.00000000000000000000E0,
30+
0.10891556624243098264E-1,
31+
0.59334456769186835896E-4,
32+
0.21231112954641805908E-6,
33+
0.54747121846510390750E-9,
34+
0.10378561511331814674E-11,
35+
0.13754880327250272679E-14,
36+
0.10223981202236205703E-17)
37+
elseif t ≤ 144.0
38+
invt = inv(t)
39+
return copysign(π/2, x) - cos(x) *
40+
@horner(invt, 0.99999999962173909991E0,
41+
0.36451060338631902917E3,
42+
0.44218548041288440874E5,
43+
0.22467569405961151887E7,
44+
0.49315316723035561922E8,
45+
0.43186795279670283193E9,
46+
0.11847992519956804350E10,
47+
0.45573267593795103181E9) /
48+
(x * @horner(invt, 1.00000000000000000000E0,
49+
0.36651060273229347594E3,
50+
0.44927569814970692777E5,
51+
0.23285354882204041700E7,
52+
0.53117852017228262911E8,
53+
0.50335310667241870372E9,
54+
0.16575285015623175410E10,
55+
0.11746532837038341076E10)) -
56+
sin(x)*invt * @horner(invt, 0.99999999920484901956E0,
57+
0.51385504875307321394E3,
58+
0.92293483452013810811E5,
59+
0.74071341863359841727E7,
60+
0.28142356162841356551E9,
61+
0.49280890357734623984E10,
62+
0.35524762685554302472E11,
63+
0.79194271662085049376E11,
64+
0.17942522624413898907E11) /
65+
@horner(invt, 1.00000000000000000000E0,
66+
0.51985504708814870209E3,
67+
0.95292615508125947321E5,
68+
0.79215459679762667578E7,
69+
0.31977567790733781460E9,
70+
0.62273134702439012114E10,
71+
0.54570971054996441467E11,
72+
0.18241750166645704670E12,
73+
0.15407148148861454434E12)
74+
elseif t < Inf
75+
invt = inv(t)
76+
return copysign(π/2, x) - cos(x) / x * (1.0 -
77+
@horner(invt, 0.19999999999999978257E1,
78+
0.22206119380434958727E4,
79+
0.84749007623988236808E6,
80+
0.13959267954823943232E9,
81+
0.10197205463267975592E11,
82+
0.30229865264524075951E12,
83+
0.27504053804288471142E13,
84+
0.21818989704686874983E13) /
85+
@horner(invt, 1.00000000000000000000E0,
86+
0.11223059690217167788E4,
87+
0.43685270974851313242E6,
88+
0.74654702140658116258E8,
89+
0.58580034751805687471E10,
90+
0.20157980379272098841E12,
91+
0.26229141857684496445E13,
92+
0.87852907334918467516E13)*invt) -
93+
sin(x) * invt * (1.0 - @horner(invt, 0.59999999999999993089E1,
94+
0.96527746044997139158E4,
95+
0.56077626996568834185E7,
96+
0.15022667718927317198E10,
97+
0.19644271064733088465E12,
98+
0.12191368281163225043E14,
99+
0.31924389898645609533E15,
100+
0.25876053010027485934E16,
101+
0.12754978896268878403E16) /
102+
@horner(invt, 1.00000000000000000000E0,
103+
0.16287957674166143196E4,
104+
0.96636303195787870963E6,
105+
0.26839734750950667021E9,
106+
0.37388510548029219241E11,
107+
0.26028585666152144496E13,
108+
0.85134283716950697226E14,
109+
0.11304079361627952930E16,
110+
0.42519841479489798424E16)*invt)
111+
elseif isnan(x)
112+
return NaN
113+
else
114+
return copysign(π/2, x)
115+
end
116+
end
117+
118+
function cosint(x::Float64)
119+
t, r0, r1 = x*x, 0.616505485620716233797110404100, 3.384180422551186426397851146402
120+
r01, r02 = 0.6162109375, 0.29454812071623379711E-3
121+
r11, r12 = 3.3837890625, 0.39136005118642639785E-3
122+
if x < 0.0
123+
return throw(DomainError())
124+
elseif x ≤ 3.0
125+
return log(x/r0) + ((x - r01) - r02) * (x + r0) *
126+
@horner(t, -0.24607411378767540707E0,
127+
0.72113492241301534559E-2,
128+
-0.11867127836204767056E-3,
129+
0.90542655466969866243E-6,
130+
-0.34322242412444409037E-8,
131+
0.51950683460656886834E-11) /
132+
@horner(t, 1.00000000000000000000E0,
133+
0.12670095552700637845E-1,
134+
0.78168450570724148921E-4,
135+
0.29959200177005821677E-6,
136+
0.73191677761328838216E-9,
137+
0.94351174530907529061E-12)
138+
elseif x ≤ 6.0
139+
return log(x/r1) + ((x - r11) - r12) * (x + r1) *
140+
@horner(t, -0.15684781827145408780E0,
141+
0.66253165609605468916E-2,
142+
-0.12822297297864512864E-3,
143+
0.12360964097729408891E-5,
144+
-0.66450975112876224532E-8,
145+
0.20326936466803159446E-10,
146+
-0.33590883135343844613E-13,
147+
0.23686934961435015119E-16) /
148+
@horner(t, 1.00000000000000000000E0,
149+
0.96166044388828741188E-2,
150+
0.45257514591257035006E-4,
151+
0.13544922659627723233E-6,
152+
0.27715365686570002081E-9,
153+
0.37718676301688932926E-12,
154+
0.27706844497155995398E-15)
155+
elseif x ≤ 12.0
156+
invt = inv(t)
157+
return sin(x) * @horner(invt, 0.99999999962173909991E0,
158+
0.36451060338631902917E3,
159+
0.44218548041288440874E5,
160+
0.22467569405961151887E7,
161+
0.49315316723035561922E8,
162+
0.43186795279670283193E9,
163+
0.11847992519956804350E10,
164+
0.45573267593795103181E9) /
165+
(x * @horner(invt, 1.00000000000000000000E0,
166+
0.36651060273229347594E3,
167+
0.44927569814970692777E5,
168+
0.23285354882204041700E7,
169+
0.53117852017228262911E8,
170+
0.50335310667241870372E9,
171+
0.16575285015623175410E10,
172+
0.11746532837038341076E10)) -
173+
cos(x) * invt * @horner(invt, 0.99999999920484901956E0,
174+
0.51385504875307321394E3,
175+
0.92293483452013810811E5,
176+
0.74071341863359841727E7,
177+
0.28142356162841356551E9,
178+
0.49280890357734623984E10,
179+
0.35524762685554302472E11,
180+
0.79194271662085049376E11,
181+
0.17942522624413898907E11) /
182+
@horner(invt, 1.00000000000000000000E0,
183+
0.51985504708814870209E3,
184+
0.95292615508125947321E5,
185+
0.79215459679762667578E7,
186+
0.31977567790733781460E9,
187+
0.62273134702439012114E10,
188+
0.54570971054996441467E11,
189+
0.18241750166645704670E12,
190+
0.15407148148861454434E12)
191+
elseif x < Inf
192+
invt = inv(t)
193+
return sin(x)/x * (1.0 - @horner(invt, 0.19999999999999978257E1,
194+
0.22206119380434958727E4,
195+
0.84749007623988236808E6,
196+
0.13959267954823943232E9,
197+
0.10197205463267975592E11,
198+
0.30229865264524075951E12,
199+
0.27504053804288471142E13,
200+
0.21818989704686874983E13) /
201+
@horner(invt, 1.00000000000000000000E0,
202+
0.11223059690217167788E4,
203+
0.43685270974851313242E6,
204+
0.74654702140658116258E8,
205+
0.58580034751805687471E10,
206+
0.20157980379272098841E12,
207+
0.26229141857684496445E13,
208+
0.87852907334918467516E13)*invt) -
209+
cos(x)*invt * (1.0 - @horner(invt, 0.59999999999999993089E1,
210+
0.96527746044997139158E4,
211+
0.56077626996568834185E7,
212+
0.15022667718927317198E10,
213+
0.19644271064733088465E12,
214+
0.12191368281163225043E14,
215+
0.31924389898645609533E15,
216+
0.25876053010027485934E16,
217+
0.12754978896268878403E16) /
218+
@horner(invt, 1.00000000000000000000E0,
219+
0.16287957674166143196E4,
220+
0.96636303195787870963E6,
221+
0.26839734750950667021E9,
222+
0.37388510548029219241E11,
223+
0.26028585666152144496E13,
224+
0.85134283716950697226E14,
225+
0.11304079361627952930E16,
226+
0.42519841479489798424E16)*invt)
227+
elseif isnan(x)
228+
return NaN
229+
else
230+
return 0.0
231+
end
232+
end
233+
234+
for f in (:sinint, :cosint)
235+
@eval begin
236+
($f)(x::Float32) = Float32(($f)(Float64(x)))
237+
($f)(x::Float16) = Float16(($f)(Float64(x)))
238+
($f)(x::Real) = ($f)(float(x))
239+
($f)(x::AbstractFloat) = error("not implemented for ", typeof(x))
240+
end
241+
end
242+
243+
244+
"""
245+
sinint(x)
246+
247+
Compute the sine integral function of `x`, defined by ``\\operatorname{Si}(x) := \\int_0^x\\frac{\\sin t}{t} dt``
248+
for real `x`.
249+
"""
250+
sinint
251+
252+
"""
253+
cosint(x)
254+
255+
Compute the cosine integral function of `x`, defined by ``\\operatorname{Ci}(x) := \\gamma + \\log x + \\int_0^x \\frac{\\cos t - 1}{t} dt``
256+
for real `x > 0`, where ``\\gamma`` is the Euler-Mascheroni constant.
257+
"""
258+
cosint

0 commit comments

Comments
 (0)