add sine and cosine integrals #32
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,246 @@ | ||
| const pidiv2 = π/2 | ||
|
|
||
| import Base.Math: @horner | ||
|
There was a problem hiding this comment.
There was a problem hiding this comment. changed to |
||
|
|
||
| # 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 | ||
| absx = abs(x) | ||
| if absx ≤ 6.0 | ||
|
There was a problem hiding this comment. Seems like you could omit the |
||
| 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 absx ≤ 12.0 | ||
| invt = inv(t) | ||
|
There was a problem hiding this comment. i think these are equivalent There was a problem hiding this comment. better then to just use There was a problem hiding this comment. 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. It can, i.e these are all equivalent and get compiled to the same code There was a problem hiding this comment. @MikaelSlevinsky, it seems like that definition would fail (spurious underflow) for There was a problem hiding this comment.
There was a problem hiding this comment. Still fails for There was a problem hiding this comment. Yes, sorry about the dispatch loop. Anyways, surely there are multiple corner cases, but it's a reasonable request that 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)
There was a problem hiding this comment. Sure, that would be a nice property to have. |
||
| return sign(x)*pidiv2 - 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 absx < Inf | ||
| invt = inv(t) | ||
| return sign(x) * pidiv2 - 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 sign(x) * pidiv2 | ||
|
There was a problem hiding this comment.
|
||
| 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.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`` | ||
|
There was a problem hiding this comment. What is gamma here? Also, should put newlines after the function signatures in docstrings. There was a problem hiding this comment. gamma is the Euler-Mascheroni constant, I just tried to replicate the style of current repository https://github.com/JuliaMath/SpecialFunctions.jl/blob/master/src/erf.jl#L40, but i now see that it's inconsistent. There was a problem hiding this comment. Yeah, I'd add something like "where |
||
| for real `x > 0`. | ||
| """ | ||
| 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.
I believe this is unnecessary, i.e just using the constant
π/2is sufficient and will be compiled to a constant in the function anyway... can check comparing@code_llvmor similarThere 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.
removed