Make pochhammer more robust #72

putianyi889 · 2019-06-15T20:44:54Z

pochhammer can now deal with arbitrary ::Number without throwing unexpected errors.
I did this by adding a modified version of gamma function named newgamma which returns Inf at negative integers.

Added a test file.

Edit: following @dlfivefifty 's suggestion, I used ogamma which returns 0 at negative integers.

Added function newgamma which is identical to gamma but returns Inf at negative integers.

Update specialfunctions.jl

dlfivefifty · 2019-06-16T11:43:50Z

Make sure your tests pass on Travis. I believe you forgot an import FastTransforms: pochhammer.

Also, it would be more natural to have a special function specifically for 1/gamma(x), say called ogamma, instead of relying on floating point Inf calculations. I'm surprised this doesn't already exist (from what I can tell).

MikaelSlevinsky · 2019-06-16T19:24:49Z

Thanks for taking a look at this! Special functions are always tricky to implement and my early version of the Pochhammer symbol is quite naive. While we're at it, we may as well consider all cases. Let me start mapping out some cases that come to mind for real numbers.

Formally, (x)_n = gamma(x+n)/gamma(x), but this runs into trouble when x or x+n (or both!) are near a negative integer.

For x < 0, we may use the Euler reflection formula, gamma(x)gamma(1-x) = pi/sinpi(x) to relate a Pochhamer evaluation to a positive argument call (n >= 0):

(-x)_n = sinpi(x)/sinpi(x-n) (x-n+1)_n

This formula is a start, but more details arise in the implementation of the ratio of sines (or cardinal sines), especially as n tends to an integer. Fortunately, when n is an integer, this reduces to the more well-known formula (-x)_n = (-1)ⁿ (x-n+1)_n. Otherwise, two options include argument reduction (to bring both arguments as close as possible), and the sine addition formula.

Now, if x+n is a negative integer and x is not, then I think the correct output is arguably analogous to:

julia> log2(-1.0)
ERROR: DomainError with -1.0:
NaN result for non-NaN input.

because gamma(x) is then just a constant finite scaling of gamma(x+n) which is evaluated at a simple pole. Since the pole is simple, the infinite limit depends on the direction of approach, so we cannot simply settle for one Inf or the other. Therefore, NaN seems appropriate, and to conform with Julia's other special functions would suggest returning a DomainError.

Finally, if x+n and x are both negative integers, there is a ratio of residues that we should be able to reason with. For example, if n::Integer >= 0 we may safely use (-x)_n = (-1)ⁿ (x-n+1)_n, but if n<0 then we may want to use (x)_n = 1/(x+n)_-n first.

codecov · 2019-06-16T21:39:05Z

Codecov Report

Merging #72 into master will increase coverage by 0.1%.
The diff coverage is 100%.

@@            Coverage Diff            @@
##           master      #72     +/-   ##
=========================================
+ Coverage    52.3%   52.41%   +0.1%     
=========================================
  Files          33       33             
  Lines        3185     3188      +3     
=========================================
+ Hits         1666     1671      +5     
+ Misses       1519     1517      -2

Impacted Files	Coverage Δ
src/specialfunctions.jl	`74.32% <100%> (+0.83%)`	⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update d2add70...8336ae1. Read the comment docs.

codecov · 2019-06-16T21:39:05Z

Codecov Report

Merging #72 into master will increase coverage by 20.04%.
The diff coverage is 100%.

@@             Coverage Diff             @@
##           master      #72       +/-   ##
===========================================
+ Coverage    52.3%   72.35%   +20.04%     
===========================================
  Files          33       33               
  Lines        3185     2308      -877     
===========================================
+ Hits         1666     1670        +4     
+ Misses       1519      638      -881

Impacted Files	Coverage Δ
src/specialfunctions.jl	`83.22% <100%> (+9.72%)`	⬆️
src/clenshawcurtis.jl	`60% <0%> (-15%)`	⬇️
src/TriangularHarmonics/slowplan.jl	`0% <0%> (ø)`	⬆️
src/TriangularHarmonics/trifunctions.jl	`0% <0%> (ø)`	⬆️
src/chebyshevtransform.jl	`97% <0%> (+4.61%)`	⬆️
src/cjt.jl	`76.71% <0%> (+6.71%)`	⬆️
src/hierarchical.jl	`48.06% <0%> (+6.83%)`	⬆️
src/SphericalHarmonics/sphfunctions.jl	`78.94% <0%> (+13.72%)`	⬆️
src/SphericalHarmonics/slowplan.jl	`46.79% <0%> (+14.35%)`	⬆️
... and 20 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update d2add70...055248f. Read the comment docs.

dlfivefifty · 2019-06-20T09:31:13Z

src/specialfunctions.jl

@@ -61,6 +61,14 @@ function pochhammer(x::Number,n::UnitRange{T}) where T<:Real
    ret
 end

+function ogamma(x::Number)
+    if isinteger(x) && x<0
+        0.0


This should be type stable: zero(float(x))

Also, why not x ≤ 0?

Because the original gamma already returns Inf
gamma(0.0)=Inf
gamma(-0.0)=-Inf

dlfivefifty · 2019-06-20T09:31:45Z

src/specialfunctions.jl

+    if isinteger(x) && x<0
+        0.0
+    else
+        1.0/gamma(x)


inv(gamma(x)) is cleaner.

dlfivefifty · 2019-06-20T09:32:27Z

test/specialfunctionstests.jl

+    @test pochhammer(-1,2) == 0
+    @test pochhammer(-5,3) == -60
+    @test pochhammer(-1,-0.5) == 0
+    @test 1.0/pochhammer(-0.5,-0.5) == 0


Add complex tests?

gamma is nonzero and holomorphic except the singularities, so complex tests don't make any difference, I suppose. I skipped array tests, though. Maybe add some array tests

I mean complex types but at integers: ogamma(-1+0im)

dlfivefifty · 2019-06-20T14:28:57Z

test/specialfunctionstests.jl

+    @test pochhammer(-1,-0.5) == 0
+    @test 1.0/pochhammer(-0.5,-0.5) == 0
+    @test pochhammer(-1+0im,-1) == -2
+    @test pochhammer([2,-5,-1],3) == [24,-60,0]


vector-valued arguments should be deprecated in favour of pochhammer.([2,-5,-1],3), so perhaps remove this test.

So can pochhammer(x::AbstractArray{T},n::Integer) be removed? Do abstractarrays work with dot operations?

It can be deprecated, though that can be done separately. Yes abstract arrays support broadcastibg

Merge pull request JuliaApproximation#72 from putianyi889/master

putianyi889 added 8 commits June 15, 2019 17:14

Update specialfunctions.jl

5814f60

Added function newgamma which is identical to gamma but returns Inf at negative integers.

Merge pull request #1 from putianyi889/patch-1

99d1d6b

Update specialfunctions.jl

change pochhammer

e2825ea

Create specialfunctionstests.jl

422a992

Update specialfunctionstests.jl

cdb83dd

Update runtests.jl

2c0dc6e

Update specialfunctionstests.jl

ed9f733

Update specialfunctionstests.jl

c140788

putianyi889 added 4 commits June 16, 2019 19:27

Update specialfunctions.jl

b56c1a5

Update specialfunctionstests.jl

d98a398

Update specialfunctions.jl

8d1839b

Update specialfunctionstests.jl

8336ae1

dlfivefifty requested changes Jun 20, 2019

View reviewed changes

putianyi889 added 3 commits June 20, 2019 12:32

Update specialfunctions.jl

54c0bd0

Update specialfunctionstests.jl

4ab2f19

Update specialfunctionstests.jl

6acfdee

dlfivefifty reviewed Jun 20, 2019

View reviewed changes

putianyi889 added 3 commits June 20, 2019 16:41

Update specialfunctionstests.jl

31d5218

Update specialfunctionstests.jl

2d371f2

Update specialfunctions.jl

055248f

dlfivefifty merged commit 3e9c187 into JuliaApproximation:master Jun 21, 2019

putianyi889 added a commit to putianyi889/FastTransforms.jl that referenced this pull request Jun 22, 2019

Merge pull request #2 from JuliaApproximation/master

4917162

Merge pull request JuliaApproximation#72 from putianyi889/master

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Make pochhammer more robust #72

Make pochhammer more robust #72

Uh oh!

putianyi889 commented Jun 15, 2019 •

edited

Loading

Uh oh!

dlfivefifty commented Jun 16, 2019

Uh oh!

MikaelSlevinsky commented Jun 16, 2019

Uh oh!

codecov bot commented Jun 16, 2019

Uh oh!

codecov bot commented Jun 16, 2019 •

edited

Loading

Uh oh!

dlfivefifty Jun 20, 2019

Uh oh!

dlfivefifty Jun 20, 2019

Uh oh!

putianyi889 Jun 20, 2019

Uh oh!

dlfivefifty Jun 20, 2019

Uh oh!

dlfivefifty Jun 20, 2019

Uh oh!

putianyi889 Jun 20, 2019

Uh oh!

dlfivefifty Jun 20, 2019

Uh oh!

dlfivefifty Jun 20, 2019

Uh oh!

putianyi889 Jun 20, 2019

Uh oh!

dlfivefifty Jun 20, 2019

Uh oh!

Uh oh!

Make pochhammer more robust #72

Make pochhammer more robust #72

Uh oh!

Conversation

putianyi889 commented Jun 15, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

dlfivefifty commented Jun 16, 2019

Uh oh!

MikaelSlevinsky commented Jun 16, 2019

Uh oh!

codecov bot commented Jun 16, 2019

Codecov Report

Uh oh!

codecov bot commented Jun 16, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

putianyi889 commented Jun 15, 2019 •

edited

Loading

codecov bot commented Jun 16, 2019 •

edited

Loading