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 1 commit
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
106 changes: 60 additions & 46 deletions src/sincosint.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,8 @@ import Base.Math: @horner
# This program computes the sine integral ∫_0^x sin(t)/t dt using the rational approximations of A.J. MacLeod, Numer. Algor., 12:259--272, 1996.
function sinint(x::Float64)
t = x*x
if abs(x) ≤ 6.0
absx = abs(x)
if absx ≤ 6.0
Copy link
Member

Choose a reason for hiding this comment

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

Seems like you could omit the abs(x) call and just compare t ≤ 36 etcetera?

return x * @horner(t, 1.00000000000000000000E0,
-0.44663998931312457298E-1,
0.11209146443112369449E-2,
Expand All @@ -22,24 +23,24 @@ function sinint(x::Float64)
0.10378561511331814674E-11,
0.13754880327250272679E-14,
0.10223981202236205703E-17)
elseif abs(x) ≤ 12.0
elseif absx ≤ 12.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 sign(x)*pidiv2 - cos(x)/x*@horner(invt, 0.99999999962173909991E0,
return sign(x)*pidiv2 - cos(x) * @horner(invt, 0.99999999962173909991E0,
0.36451060338631902917E3,
0.44218548041288440874E5,
0.22467569405961151887E7,
0.49315316723035561922E8,
0.43186795279670283193E9,
0.11847992519956804350E10,
0.45573267593795103181E9) /
@horner(invt, 1.00000000000000000000E0,
(x * @horner(invt, 1.00000000000000000000E0,
0.36651060273229347594E3,
0.44927569814970692777E5,
0.23285354882204041700E7,
0.53117852017228262911E8,
0.50335310667241870372E9,
0.16575285015623175410E10,
0.11746532837038341076E10) -
0.11746532837038341076E10)) -
sin(x)*invt * @horner(invt, 0.99999999920484901956E0,
0.51385504875307321394E3,
0.92293483452013810811E5,
Expand All @@ -58,9 +59,10 @@ function sinint(x::Float64)
0.54570971054996441467E11,
0.18241750166645704670E12,
0.15407148148861454434E12)
else
elseif absx < Inf
invt = inv(t)
return sign(x)*pidiv2 - cos(x)/x*(1.0 - @horner(invt, 0.19999999999999978257E1,
return sign(x) * pidiv2 - cos(x) / x * (1.0 -
@horner(invt, 0.19999999999999978257E1,
0.22206119380434958727E4,
0.84749007623988236808E6,
0.13959267954823943232E9,
Expand All @@ -76,7 +78,7 @@ function sinint(x::Float64)
0.20157980379272098841E12,
0.26229141857684496445E13,
0.87852907334918467516E13)*invt) -
sin(x)*invt * (1.0 - @horner(invt, 0.59999999999999993089E1,
sin(x) * invt * (1.0 - @horner(invt, 0.59999999999999993089E1,
0.96527746044997139158E4,
0.56077626996568834185E7,
0.15022667718927317198E10,
Expand All @@ -94,17 +96,24 @@ function sinint(x::Float64)
0.85134283716950697226E14,
0.11304079361627952930E16,
0.42519841479489798424E16)*invt)
elseif isnan(x)
return NaN
else
return sign(x) * pidiv2
Copy link
Member

@stevengj stevengj Jun 13, 2017 •
edited

Choose a reason for hiding this comment

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

copysign(pidiv2, x), and similarly elsewhere.

end
end


# This program computes the cosine integral γ + log x + ∫_0^x (cos(t)-1)/t dt using the rational approximations of A.J. MacLeod, Numer. Algor., 12:259--272, 1996.
function cosint(x::Float64)
t, r0, r1 = x*x, 0.616505485620716233797110404100, 3.384180422851186426397851146402
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) + (t-r0^2) * @horner(t, -0.24607411378767540707E0,
return log(x/r0) + ((x - r01) - r02) * (x + r0) *
@horner(t, -0.24607411378767540707E0,
0.72113492241301534559E-2,
-0.11867127836204767056E-3,
0.90542655466969866243E-6,
Expand All @@ -117,7 +126,8 @@ function cosint(x::Float64)
0.73191677761328838216E-9,
0.94351174530907529061E-12)
elseif x ≤ 6.0
return log(x/r1) + (t-r1^2) * @horner(t, -0.15684781827145408780E0,
return log(x/r1) + ((x - r11) - r12) * (x + r1) *
@horner(t, -0.15684781827145408780E0,
0.66253165609605468916E-2,
-0.12822297297864512864E-3,
0.12360964097729408891E-5,
Expand All @@ -134,41 +144,41 @@ function cosint(x::Float64)
0.27706844497155995398E-15)
elseif x ≤ 12.0
invt = inv(t)
return sin(x)/x * @horner(invt, 0.99999999962173909991E0,
0.36451060338631902917E3,
0.44218548041288440874E5,
0.22467569405961151887E7,
0.49315316723035561922E8,
0.43186795279670283193E9,
0.11847992519956804350E10,
0.45573267593795103181E9) /
@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)
else
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,
Expand Down Expand Up @@ -204,6 +214,10 @@ function cosint(x::Float64)
0.85134283716950697226E14,
0.11304079361627952930E16,
0.42519841479489798424E16)*invt)
elseif isnan(x)
return NaN
else
return 0.0
end
end

Expand Down
28 changes: 19 additions & 9 deletions test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,12 @@ using Base.Test

const SF = SpecialFunctions

# useful test functions for relative error, which differ from isapprox
# in that relerrc separately looks at the real and imaginary parts
relerr(z, x) = z == x ? 0.0 : abs(z - x) / abs(x)
relerrc(z, x) = max(relerr(real(z),real(x)), relerr(imag(z),imag(x)))
≅(a,b) = relerrc(a,b) ≤ 1e-13

@testset "error functions" begin
@test SF.erf(Float16(1)) ≈ 0.84270079294971486934
@test SF.erf(1) ≈ 0.84270079294971486934
Expand Down Expand Up @@ -53,9 +59,9 @@ end
sinintvals = [0.9460830703671830149, 1.605412976802694849, 1.848652527999468256, 1.75820313894905306, 1.54993124494467414, 1.4246875512805065, 1.4545966142480936, 1.5741868217069421, 1.665040075829602, 1.658347594218874, 1.578306806945727416, 1.504971241526373371, 1.499361722862824564, 1.556211050077665054, 1.618194443708368739, 1.631302268270032886, 1.590136415870701122, 1.536608096861185462, 1.518630031769363932, 1.548241701043439840]
cosintvals = [0.3374039229009681347, 0.4229808287748649957, 0.119629786008000328, -0.14098169788693041, -0.19002974965664388, -0.06805724389324713, 0.07669527848218452, 0.122433882532010, 0.0553475313331336, -0.045456433004455, -0.08956313549547997948, -0.04978000688411367560, 0.02676412556403455504, 0.06939635592758454727, 0.04627867767436043960, -0.01420019012019002240, -0.05524268226081385053, -0.04347510299950100478, 0.00515037100842612857, 0.04441982084535331654]
for x in 1:20
@test SF.sinint(x) ≈ sinintvals[x]
@test SF.sinint(-x) ≈ -sinintvals[x]
@test SF.cosint(x) ≈ cosintvals[x]
@test SF.sinint(x) ≅ sinintvals[x]
@test SF.sinint(-x) ≅ -sinintvals[x]
@test SF.cosint(x) ≅ cosintvals[x]
end

@test SF.sinint(1.f0) == Float32(SF.sinint(1.0))
Expand All @@ -67,6 +73,16 @@ end
@test SF.sinint(1//2) == SF.sinint(0.5)
@test SF.cosint(1//2) == SF.cosint(0.5)

@test SF.sinint(1e300) ≅ SF.pidiv2
@test SF.cosint(1e300) ≅ -8.17881912115908554103E-301
@test SF.sinint(1e-300) ≅ 1.0E-300
@test SF.cosint(1e-300) ≅ -690.1983122333121

@test SF.sinint(Inf) == SF.pidiv2
@test SF.cosint(Inf) == 0.0
@test isnan(SF.sinint(NaN))
@test isnan(SF.cosint(NaN))

@test_throws ErrorException SF.sinint(big(1.0))
@test_throws ErrorException SF.cosint(big(1.0))
@test_throws DomainError SF.cosint(-1.0)
Copy link
Member

Choose a reason for hiding this comment

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

Other things to test: a very small value, like x = 1e-300, and a very large value like x = 1e300. I'm a bit concerned that the latter will suffer from spurious overflow.

Copy link
Member

Choose a reason for hiding this comment

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

Also test Inf (should give π/2 or 0) and NaN (should give NaN)

Expand Down Expand Up @@ -283,12 +299,6 @@ end
end
end

# useful test functions for relative error, which differ from isapprox
# in that relerrc separately looks at the real and imaginary parts
relerr(z, x) = z == x ? 0.0 : abs(z - x) / abs(x)
relerrc(z, x) = max(relerr(real(z),real(x)), relerr(imag(z),imag(x)))
≅(a,b) = relerrc(a,b) ≤ 1e-13

@testset "gamma and friends" begin
@testset "digamma" begin
@testset "$elty" for elty in (Float32, Float64)
Expand Down