Support Jacobi(-1,0)\Ultraspherical(-1/2) and Jacobi(0,-1)

[DanielVandH]

On the current version,

julia> Jacobi(-1,0)\Ultraspherical(-1/2)
(ℵ₀×ℵ₀ LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}} with indices OneToInf()×OneToInf()) * (ℵ₀×ℵ₀ LinearAlgebra.Diagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{SubArray{Float64, 1, Ultraspherical{Float64, Float64}, Tuple{Float64, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}, SubArray{Float64, 1, Jacobi{Float64, Float64}, Tuple{Float64, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}}}} with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf():
 1.0  NaN       ⋅      ⋅      ⋅      ⋅      ⋅      ⋅   …
  ⋅   NaN    NaN       ⋅      ⋅      ⋅      ⋅      ⋅
  ⋅      ⋅   NaN    NaN       ⋅      ⋅      ⋅      ⋅
  ⋅      ⋅      ⋅   NaN    NaN       ⋅      ⋅      ⋅
  ⋅      ⋅      ⋅      ⋅   NaN    NaN       ⋅      ⋅
  ⋅      ⋅      ⋅      ⋅      ⋅   NaN    NaN       ⋅   …
  ⋅      ⋅      ⋅      ⋅      ⋅      ⋅   NaN    NaN
  ⋅      ⋅      ⋅      ⋅      ⋅      ⋅      ⋅   NaN
  ⋅      ⋅      ⋅      ⋅      ⋅      ⋅      ⋅      ⋅
  ⋅      ⋅      ⋅      ⋅      ⋅      ⋅      ⋅      ⋅
  ⋅      ⋅      ⋅      ⋅      ⋅      ⋅      ⋅      ⋅   …
  ⋅      ⋅      ⋅      ⋅      ⋅      ⋅      ⋅      ⋅
  ⋅      ⋅      ⋅      ⋅      ⋅      ⋅      ⋅      ⋅
 ⋮                                  ⋮                  ⋱

This adds a fix, noting that

$$\begin{align*} U_0^{(-1/2)}(x) &= P_0^{(-1,0)}(x), \\\ U_1^{(-1/2)}(x) &= -P_0^{(-1,0)}(x) - 2P_1^{(-1,0)}(x), \\\ U_n^{(-1/2)}(x) &= -\frac{2}{2n-1}\left(P_{n-1}^{(-1,0)}(x) + P_n^{(-1,0)}(x)\right), \quad n \geq 2. \end{align*}$$

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 91.03%. Comparing base (15d5f64) to head (cbc7424).
⚠️ Report is 1 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #245      +/-   ##
==========================================
+ Coverage   90.99%   91.03%   +0.03%     
==========================================
  Files          19       19              
  Lines        1744     1751       +7     
==========================================
+ Hits         1587     1594       +7     
  Misses        157      157              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[DanielVandH]

Also added Jacobi(0,-1)

[dlfivefifty]

Add a comment why this special case is needed.

Also, where in (A\B̃)*Diagonal(B[1,:]./B̃[1,:]) are the NaNs coming in?

[DanielVandH]

The NaNs come from B̃[1,:], since

julia> Jacobi(-1, -1)[1, :]
ℵ₀-element view(::Jacobi{Float64, Int64}, 1.0, :) with eltype Float64 with indices OneToInf():
   1.0
   0.0
 NaN
 NaN
 NaN
 NaN
 NaN
 NaN
 NaN
 NaN
   ⋮

(since we haven't yet given a definition to Jacobi(-1, -1) in the code that works, although it'd be all zero anyway -except for the first two entries- I think?)

[dlfivefifty]

P_k^(-1,-1)(x) are well-defined:

I think using the series or Hypergeometric formula. They are just degenerate since P_1^(-1,-1)(x) = 0, so we cannot use the 3-term recurrence to build them.

So really we should be fixing evaluation for degenerate Jacobi polynomials.

[DanielVandH]

Won't it still lead to a NaN anyway since $P_1^{(-1,-1)}(x) = 0$? We are dividing by it

[dlfivefifty]

Touché

[dlfivefifty]

update the comment to mention even if it is defined the above won't work because of division by zero.