-
Notifications
You must be signed in to change notification settings - Fork 99
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.
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
add sine and cosine integrals #32
Merged
Merged
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 fbe070a
fix sincosint
MikaelSlevinsky c4e0d19
absx => t, and sign(x) => copysign
MikaelSlevinsky 8ad92b7
reformulate comments for clearer citation, including a note regarding…
MikaelSlevinsky c6d7089
add colon
MikaelSlevinsky 9a8b0e5
add commas
MikaelSlevinsky 3685fcf
remove pidiv2, change import to using
MikaelSlevinsky e14d654
int(t) => 1.0/t
MikaelSlevinsky 41a7a91
Revert "int(t) => 1.0/t"
MikaelSlevinsky 1f2a786
add newlines and add ``\\gamma`` in the cosine integral docstring
MikaelSlevinsky File filter
Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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) | ||
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 |
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
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
There was a problem hiding this comment.
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
There was a problem hiding this comment.
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
sinceinv
is more for matrix inversesThere was a problem hiding this comment.
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.
There was a problem hiding this comment.
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)
There was a problem hiding this comment.
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)
.There was a problem hiding this comment.
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))
There was a problem hiding this comment.
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)
(givesNaN + NaN*im
instead ofInf - Inf*im
). (And I assume that the nestedx/y
calls the existing implementation, since otherwise you have a dispatch loop.)There was a problem hiding this comment.
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):There was a problem hiding this comment.
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.