Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

add sine and cosine integrals #32

Merged
merged 10 commits into from
Jun 27, 2017
Merged

add sine and cosine integrals #32

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
2cfc869
add sine and cosine integrals
MikaelSlevinsky Jun 12, 2017
fbe070a
fix sincosint
MikaelSlevinsky Jun 13, 2017
c4e0d19
absx => t, and sign(x) => copysign
MikaelSlevinsky Jun 13, 2017
8ad92b7
reformulate comments for clearer citation, including a note regarding…
MikaelSlevinsky Jun 13, 2017
c6d7089
add colon
MikaelSlevinsky Jun 13, 2017
9a8b0e5
add commas
MikaelSlevinsky Jun 13, 2017
3685fcf
remove pidiv2, change import to using
MikaelSlevinsky Jun 14, 2017
e14d654
int(t) => 1.0/t
MikaelSlevinsky Jun 14, 2017
41a7a91
Revert "int(t) => 1.0/t"
MikaelSlevinsky Jun 15, 2017
1f2a786
add newlines and add ``\\gamma`` in the cosine integral docstring
MikaelSlevinsky Jun 27, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
6 changes: 3 additions & 3 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
# SpecialFunctions.jl

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

[![Travis](https://travis-ci.org/JuliaMath/SpecialFunctions.jl.svg?branch=master)](https://travis-ci.org/JuliaMath/SpecialFunctions.jl)
[![AppVeyor](https://ci.appveyor.com/api/projects/status/ccfgkm2cjcggu158/branch/master?svg=true)](https://ci.appveyor.com/project/ararslan/specialfunctions-jl/branch/master)
Expand Down
4 changes: 3 additions & 1 deletion docs/src/index.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,8 @@ libraries.
| [`erfi(x)`](@ref SpecialFunctions.erfi) | imaginary error function defined as `-im * erf(x * im)`, where `im` is the imaginary unit |
| [`erfcx(x)`](@ref SpecialFunctions.erfcx) | scaled complementary error function, i.e. accurate `exp(x^2) * erfc(x)` for large `x` |
| [`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` |
| [`sinint(x)`](@ref SpecialFunctions.sinint) | [sine integral](https://en.wikipedia.org/wiki/Trigonometric_integral) at `x` |
| [`cosint(x)`](@ref SpecialFunctions.cosint) | [cosine integral](https://en.wikipedia.org/wiki/Trigonometric_integral) at `x` |
| [`digamma(x)`](@ref SpecialFunctions.digamma) | [digamma function](https://en.wikipedia.org/wiki/Digamma_function) (i.e. the derivative of `lgamma` at `x`) |
| [`eta(x)`](@ref SpecialFunctions.eta) | [Dirichlet eta function](https://en.wikipedia.org/wiki/Dirichlet_eta_function) at `x` |
| [`zeta(x)`](@ref SpecialFunctions.zeta) | [Riemann zeta function](https://en.wikipedia.org/wiki/Riemann_zeta_function) at `x` |
Expand Down Expand Up @@ -51,7 +53,7 @@ from the Julia REPL.

## Note

Prior to Julia 0.6, these functions were available in Julia's Base module.
Prior to Julia 0.6, most of these functions were available in Julia's Base module.
Because of this, the symbols from this package are not exported on Julia 0.5
to avoid name conflicts.
In this case, the symbols will need to be explicitly imported or called
Expand Down
2 changes: 2 additions & 0 deletions docs/src/special.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,8 @@ SpecialFunctions.erfi
SpecialFunctions.dawson
SpecialFunctions.erfinv
SpecialFunctions.erfcinv
SpecialFunctions.sinint
SpecialFunctions.cosint
SpecialFunctions.digamma
SpecialFunctions.invdigamma
SpecialFunctions.trigamma
Expand Down
4 changes: 4 additions & 0 deletions src/SpecialFunctions.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,9 @@ if VERSION >= v"0.6.0-dev.2767"
end
end

export sinint,
cosint

if isdefined(Base.Math, :openspecfun)
const openspecfun = Base.Math.openspecfun
else
Expand All @@ -67,6 +70,7 @@ end

include("bessel.jl")
include("erf.jl")
include("sincosint.jl")
include("gamma.jl")
include("deprecated.jl")

Expand Down
258 changes: 258 additions & 0 deletions src/sincosint.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,258 @@
using Base.Math.@horner

# Compute the sine integral: ∫_0^x sin(t)/t dt,
# and the cosine integral: γ + log x + ∫_0^x (cos(t)-1)/t dt,
# using the rational approximants tabulated in:
# A.J. MacLeod, "Rational approximations, software and test methods for
# sine and cosine integrals," Numer. Algor. 12, pp. 259--272 (1996).
# http://dx.doi.org/10.1007/BF02142806
# https://link.springer.com/article/10.1007/BF02142806
#
# Note: the second zero of Ci(x) has a typo that is fixed:
#
# r1 = 3.38418 04228 51186 42639 78511 46402 in the article, but is in fact:
#
# r1 = 3.38418 04225 51186 42639 78511 46402.
#

function sinint(x::Float64)
t = x*x
if t ≤ 36.0
return x * @horner(t, 1.00000000000000000000E0,
-0.44663998931312457298E-1,
0.11209146443112369449E-2,
-0.13276124407928422367E-4,
0.85118014179823463879E-7,
-0.29989314303147656479E-9,
0.55401971660186204711E-12,
-0.42406353433133212926E-15) /
@horner(t, 1.00000000000000000000E0,
0.10891556624243098264E-1,
0.59334456769186835896E-4,
0.21231112954641805908E-6,
0.54747121846510390750E-9,
0.10378561511331814674E-11,
0.13754880327250272679E-14,
0.10223981202236205703E-17)
elseif t ≤ 144.0
invt = inv(t)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why not just 1/t

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i think these are equivalent

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

better then to just use 1/t since inv is more for matrix inverses

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i don't know why it's defined if it can't be used in a completely equivalent case, but i don't care enough to discuss it further: it's changed.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It can, i.e these are all equivalent and get compiled to the same code
1/2.0
1.0/2.0
inv(2.0)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@MikaelSlevinsky, it seems like that definition would fail (spurious underflow) for 0.1/(1e308+0im).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

/(x::Real, y::Complex) = (z = y/x; isinf(z) ? x/y : inv(z))

Copy link
Member

@stevengj stevengj Jun 30, 2017 •
edited

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Still fails for Inf / (1 + 1im) (gives NaN + NaN*im instead of Inf - Inf*im). (And I assume that the nested x/y calls the existing implementation, since otherwise you have a dispatch loop.)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, sorry about the dispatch loop. Anyways, surely there are multiple corner cases, but it's a reasonable request that x/(y+0.0im) == x/y. Many base functions behave like this (at least for input that isn't extreme in some sense):

julia> function realcomplextest(x::AbstractArray, f::Function)
           for y in x
               test = f(y) - f(y+0im)
               if test != 0
                   println(y,"  ",test)
               end
           end
       end
realcomplextest (generic function with 1 method)

julia> realcomplextest(linspace(0,50,10001), sin)

julia> realcomplextest(linspace(0,50,10001), cos)

julia> realcomplextest(linspace(0,50,10001), exp)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sure, that would be a nice property to have.

return copysign(π/2, x) - cos(x) *
@horner(invt, 0.99999999962173909991E0,
0.36451060338631902917E3,
0.44218548041288440874E5,
0.22467569405961151887E7,
0.49315316723035561922E8,
0.43186795279670283193E9,
0.11847992519956804350E10,
0.45573267593795103181E9) /
(x * @horner(invt, 1.00000000000000000000E0,
0.36651060273229347594E3,
0.44927569814970692777E5,
0.23285354882204041700E7,
0.53117852017228262911E8,
0.50335310667241870372E9,
0.16575285015623175410E10,
0.11746532837038341076E10)) -
sin(x)*invt * @horner(invt, 0.99999999920484901956E0,
0.51385504875307321394E3,
0.92293483452013810811E5,
0.74071341863359841727E7,
0.28142356162841356551E9,
0.49280890357734623984E10,
0.35524762685554302472E11,
0.79194271662085049376E11,
0.17942522624413898907E11) /
@horner(invt, 1.00000000000000000000E0,
0.51985504708814870209E3,
0.95292615508125947321E5,
0.79215459679762667578E7,
0.31977567790733781460E9,
0.62273134702439012114E10,
0.54570971054996441467E11,
0.18241750166645704670E12,
0.15407148148861454434E12)
elseif t < Inf
invt = inv(t)
return copysign(π/2, x) - cos(x) / x * (1.0 -
@horner(invt, 0.19999999999999978257E1,
0.22206119380434958727E4,
0.84749007623988236808E6,
0.13959267954823943232E9,
0.10197205463267975592E11,
0.30229865264524075951E12,
0.27504053804288471142E13,
0.21818989704686874983E13) /
@horner(invt, 1.00000000000000000000E0,
0.11223059690217167788E4,
0.43685270974851313242E6,
0.74654702140658116258E8,
0.58580034751805687471E10,
0.20157980379272098841E12,
0.26229141857684496445E13,
0.87852907334918467516E13)*invt) -
sin(x) * invt * (1.0 - @horner(invt, 0.59999999999999993089E1,
0.96527746044997139158E4,
0.56077626996568834185E7,
0.15022667718927317198E10,
0.19644271064733088465E12,
0.12191368281163225043E14,
0.31924389898645609533E15,
0.25876053010027485934E16,
0.12754978896268878403E16) /
@horner(invt, 1.00000000000000000000E0,
0.16287957674166143196E4,
0.96636303195787870963E6,
0.26839734750950667021E9,
0.37388510548029219241E11,
0.26028585666152144496E13,
0.85134283716950697226E14,
0.11304079361627952930E16,
0.42519841479489798424E16)*invt)
elseif isnan(x)
return NaN
else
return copysign(π/2, x)
end
end

function cosint(x::Float64)
t, r0, r1 = x*x, 0.616505485620716233797110404100, 3.384180422551186426397851146402
r01, r02 = 0.6162109375, 0.29454812071623379711E-3
r11, r12 = 3.3837890625, 0.39136005118642639785E-3
if x < 0.0
return throw(DomainError())
elseif x ≤ 3.0
return log(x/r0) + ((x - r01) - r02) * (x + r0) *
@horner(t, -0.24607411378767540707E0,
0.72113492241301534559E-2,
-0.11867127836204767056E-3,
0.90542655466969866243E-6,
-0.34322242412444409037E-8,
0.51950683460656886834E-11) /
@horner(t, 1.00000000000000000000E0,
0.12670095552700637845E-1,
0.78168450570724148921E-4,
0.29959200177005821677E-6,
0.73191677761328838216E-9,
0.94351174530907529061E-12)
elseif x ≤ 6.0
return log(x/r1) + ((x - r11) - r12) * (x + r1) *
@horner(t, -0.15684781827145408780E0,
0.66253165609605468916E-2,
-0.12822297297864512864E-3,
0.12360964097729408891E-5,
-0.66450975112876224532E-8,
0.20326936466803159446E-10,
-0.33590883135343844613E-13,
0.23686934961435015119E-16) /
@horner(t, 1.00000000000000000000E0,
0.96166044388828741188E-2,
0.45257514591257035006E-4,
0.13544922659627723233E-6,
0.27715365686570002081E-9,
0.37718676301688932926E-12,
0.27706844497155995398E-15)
elseif x ≤ 12.0
invt = inv(t)
return sin(x) * @horner(invt, 0.99999999962173909991E0,
0.36451060338631902917E3,
0.44218548041288440874E5,
0.22467569405961151887E7,
0.49315316723035561922E8,
0.43186795279670283193E9,
0.11847992519956804350E10,
0.45573267593795103181E9) /
(x * @horner(invt, 1.00000000000000000000E0,
0.36651060273229347594E3,
0.44927569814970692777E5,
0.23285354882204041700E7,
0.53117852017228262911E8,
0.50335310667241870372E9,
0.16575285015623175410E10,
0.11746532837038341076E10)) -
cos(x) * invt * @horner(invt, 0.99999999920484901956E0,
0.51385504875307321394E3,
0.92293483452013810811E5,
0.74071341863359841727E7,
0.28142356162841356551E9,
0.49280890357734623984E10,
0.35524762685554302472E11,
0.79194271662085049376E11,
0.17942522624413898907E11) /
@horner(invt, 1.00000000000000000000E0,
0.51985504708814870209E3,
0.95292615508125947321E5,
0.79215459679762667578E7,
0.31977567790733781460E9,
0.62273134702439012114E10,
0.54570971054996441467E11,
0.18241750166645704670E12,
0.15407148148861454434E12)
elseif x < Inf
invt = inv(t)
return sin(x)/x * (1.0 - @horner(invt, 0.19999999999999978257E1,
0.22206119380434958727E4,
0.84749007623988236808E6,
0.13959267954823943232E9,
0.10197205463267975592E11,
0.30229865264524075951E12,
0.27504053804288471142E13,
0.21818989704686874983E13) /
@horner(invt, 1.00000000000000000000E0,
0.11223059690217167788E4,
0.43685270974851313242E6,
0.74654702140658116258E8,
0.58580034751805687471E10,
0.20157980379272098841E12,
0.26229141857684496445E13,
0.87852907334918467516E13)*invt) -
cos(x)*invt * (1.0 - @horner(invt, 0.59999999999999993089E1,
0.96527746044997139158E4,
0.56077626996568834185E7,
0.15022667718927317198E10,
0.19644271064733088465E12,
0.12191368281163225043E14,
0.31924389898645609533E15,
0.25876053010027485934E16,
0.12754978896268878403E16) /
@horner(invt, 1.00000000000000000000E0,
0.16287957674166143196E4,
0.96636303195787870963E6,
0.26839734750950667021E9,
0.37388510548029219241E11,
0.26028585666152144496E13,
0.85134283716950697226E14,
0.11304079361627952930E16,
0.42519841479489798424E16)*invt)
elseif isnan(x)
return NaN
else
return 0.0
end
end

for f in (:sinint, :cosint)
@eval begin
($f)(x::Float32) = Float32(($f)(Float64(x)))
($f)(x::Float16) = Float16(($f)(Float64(x)))
($f)(x::Real) = ($f)(float(x))
($f)(x::AbstractFloat) = error("not implemented for ", typeof(x))
end
end


"""
sinint(x)

Compute the sine integral function of `x`, defined by ``\\operatorname{Si}(x) := \\int_0^x\\frac{\\sin t}{t} dt``
for real `x`.
"""
sinint

"""
cosint(x)

Compute the cosine integral function of `x`, defined by ``\\operatorname{Ci}(x) := \\gamma + \\log x + \\int_0^x \\frac{\\cos t - 1}{t} dt``
for real `x > 0`, where ``\\gamma`` is the Euler-Mascheroni constant.
"""
cosint