From cd90acb05a2c9da562e7fbeb10f9f386b05e9798 Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Mon, 21 Oct 2024 07:37:13 +0100
Subject: [PATCH 1/8] fix an inexact error

---
 src/classical/jacobi.jl | 17 +++++++++++++++--
 test/test_jacobi.jl     |  3 +++
 2 files changed, 18 insertions(+), 2 deletions(-)

diff --git a/src/classical/jacobi.jl b/src/classical/jacobi.jl
index 0bdc730..45008d3 100644
--- a/src/classical/jacobi.jl
+++ b/src/classical/jacobi.jl
@@ -312,12 +312,25 @@ function _jacobi_convert_b(a, b, k, T) # Jacobi(a, b+k) \ Jacobi(a, b)
     end
 end
 
+"""
+    _approx_integer(x; kw...)
+
+Returns `y::Integer` such that `isapprox(x,y; kw...)` is `true`. Throws `InexactError` if `y` does not exist.
+"""
+function _approx_integer(x; kw...)
+    y = round(x)
+    if !isapprox(x,y; kw...)
+        throw(InexactError(_approx_integer, x))
+    end
+    Integer(y)
+end
+
 function \(A::Jacobi, B::Jacobi)
     T = promote_type(eltype(A), eltype(B))
     aa, ab = A.a, A.b
     ba, bb = B.a, B.b
-    ka = Integer(aa-ba)
-    kb = Integer(ab-bb)
+    ka = _approx_integer(aa-ba)
+    kb = _approx_integer(ab-bb)
     if ka >= 0
         C1 = _jacobi_convert_a(ba, ab, ka, T)
         if kb >= 0
diff --git a/test/test_jacobi.jl b/test/test_jacobi.jl
index 7ad1866..fb3ab13 100644
--- a/test/test_jacobi.jl
+++ b/test/test_jacobi.jl
@@ -147,6 +147,9 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMa
             # special weighted conversions
             W = (JacobiWeight(-1,-1) .* Jacobi(0,0)) \ Jacobi(0,0)
             @test ((JacobiWeight(-1,-1) .* Jacobi(0,0)) * W)[0.1,1:10] ≈ Jacobi(0,0)[0.1,1:10]
+
+            # potential inexact errors
+            Jacobi(1.9,1.9) \ Jacobi(0.9,0.9) # isinteger(1.9-0.9) = false
         end
 
         @testset "Derivative" begin

From 36a346aefd520a7aec75ac3501c1f2be2c147f86 Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Mon, 21 Oct 2024 07:46:44 +0100
Subject: [PATCH 2/8] include doc

---
 docs/src/index.md | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/docs/src/index.md b/docs/src/index.md
index 3cbb8bd..b00bd6b 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -233,6 +233,9 @@ ClassicalOrthogonalPolynomials.WeightedOPLayout
 ClassicalOrthogonalPolynomials.interlace!
 ```
 ```@docs
+ClassicalOrthogonalPolynomials._approx_integer
+```
+```@docs
 ClassicalOrthogonalPolynomials._tritrunc
 ```
 ```@docs

From fd5d0657bcfdc9fbcdb1b7025d93fbd81c3ebeca Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Mon, 21 Oct 2024 18:02:53 +0100
Subject: [PATCH 3/8] Update jacobi.jl

---
 src/classical/jacobi.jl | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/src/classical/jacobi.jl b/src/classical/jacobi.jl
index 45008d3..f957aba 100644
--- a/src/classical/jacobi.jl
+++ b/src/classical/jacobi.jl
@@ -292,6 +292,9 @@ end
 # the type specification is also because LazyArrays.jl is slow otherwise
 # getindex and print could be very slow. This needs to be fixed by LazyArrays.jl
 function _jacobi_convert_a(a, b, k, T) # Jacobi(a+k, b) \ Jacobi(a, b)
+    j = round(k)
+    @assert j ≈ k
+    k = Integer(j)
     if iszero(k)
         Eye{T}(∞)
     elseif isone(k)
@@ -302,6 +305,9 @@ function _jacobi_convert_a(a, b, k, T) # Jacobi(a+k, b) \ Jacobi(a, b)
     end
 end
 function _jacobi_convert_b(a, b, k, T) # Jacobi(a, b+k) \ Jacobi(a, b)
+    j = round(k)
+    @assert j ≈ k
+    k = Integer(j)
     if iszero(k)
         Eye{T}(∞)
     elseif isone(k)
@@ -329,8 +335,8 @@ function \(A::Jacobi, B::Jacobi)
     T = promote_type(eltype(A), eltype(B))
     aa, ab = A.a, A.b
     ba, bb = B.a, B.b
-    ka = _approx_integer(aa-ba)
-    kb = _approx_integer(ab-bb)
+    ka = aa - ba
+    kb = ab - bb
     if ka >= 0
         C1 = _jacobi_convert_a(ba, ab, ka, T)
         if kb >= 0
@@ -345,6 +351,7 @@ function \(A::Jacobi, B::Jacobi)
         else
             list = (C1, C2)
             ApplyArray{T,2,typeof(*),typeof(list)}(*, list)
+            C1 * C2
         end
     else
         inv(B \ A)

From f42af977f6e906c371fa4f8e915d778b0ca26216 Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Mon, 21 Oct 2024 18:03:22 +0100
Subject: [PATCH 4/8] Update index.md

---
 docs/src/index.md | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/docs/src/index.md b/docs/src/index.md
index b00bd6b..3cbb8bd 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -233,9 +233,6 @@ ClassicalOrthogonalPolynomials.WeightedOPLayout
 ClassicalOrthogonalPolynomials.interlace!
 ```
 ```@docs
-ClassicalOrthogonalPolynomials._approx_integer
-```
-```@docs
 ClassicalOrthogonalPolynomials._tritrunc
 ```
 ```@docs

From 7777e19d6939797892aa6207c1770ca2f3bf7ceb Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Mon, 21 Oct 2024 18:52:33 +0100
Subject: [PATCH 5/8] Update jacobi.jl

---
 src/classical/jacobi.jl | 13 -------------
 1 file changed, 13 deletions(-)

diff --git a/src/classical/jacobi.jl b/src/classical/jacobi.jl
index f957aba..462d0d4 100644
--- a/src/classical/jacobi.jl
+++ b/src/classical/jacobi.jl
@@ -318,19 +318,6 @@ function _jacobi_convert_b(a, b, k, T) # Jacobi(a, b+k) \ Jacobi(a, b)
     end
 end
 
-"""
-    _approx_integer(x; kw...)
-
-Returns `y::Integer` such that `isapprox(x,y; kw...)` is `true`. Throws `InexactError` if `y` does not exist.
-"""
-function _approx_integer(x; kw...)
-    y = round(x)
-    if !isapprox(x,y; kw...)
-        throw(InexactError(_approx_integer, x))
-    end
-    Integer(y)
-end
-
 function \(A::Jacobi, B::Jacobi)
     T = promote_type(eltype(A), eltype(B))
     aa, ab = A.a, A.b

From d49412bf30ed075c1670a44a9f605fd18f7dd359 Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Mon, 21 Oct 2024 23:03:39 +0100
Subject: [PATCH 6/8] Update jacobi.jl

---
 src/classical/jacobi.jl | 41 +++++------------------------------------
 1 file changed, 5 insertions(+), 36 deletions(-)

diff --git a/src/classical/jacobi.jl b/src/classical/jacobi.jl
index 462d0d4..4cf4faa 100644
--- a/src/classical/jacobi.jl
+++ b/src/classical/jacobi.jl
@@ -288,34 +288,17 @@ function _jacobi_convert_b(a, b) # Jacobi(a, b+1) \ Jacobi(a, b)
     end
 end
 
-# TODO: implement ApplyArray(*, A) in LazyArrays.jl, then these can be simplified
-# the type specification is also because LazyArrays.jl is slow otherwise
-# getindex and print could be very slow. This needs to be fixed by LazyArrays.jl
 function _jacobi_convert_a(a, b, k, T) # Jacobi(a+k, b) \ Jacobi(a, b)
     j = round(k)
     @assert j ≈ k
     k = Integer(j)
-    if iszero(k)
-        Eye{T}(∞)
-    elseif isone(k)
-        _jacobi_convert_a(a, b)
-    else
-        list = tuple([_jacobi_convert_a(a+j, b) for j in k-1:-1:0]...)
-        ApplyArray{T,2,typeof(*),typeof(list)}(*, list)
-    end
+    reduce(*, [_jacobi_convert_a(a+j, b) for j in k-1:-1:0], init=Eye{T}(∞))
 end
 function _jacobi_convert_b(a, b, k, T) # Jacobi(a, b+k) \ Jacobi(a, b)
     j = round(k)
     @assert j ≈ k
     k = Integer(j)
-    if iszero(k)
-        Eye{T}(∞)
-    elseif isone(k)
-        _jacobi_convert_b(a, b)
-    else
-        list = tuple([_jacobi_convert_b(a, b+j) for j in k-1:-1:0]...)
-        ApplyArray{T,2,typeof(*),typeof(list)}(*, list)
-    end
+    reduce(*, [_jacobi_convert_b(a, b+j) for j in k-1:-1:0], init=Eye{T}(∞))
 end
 
 function \(A::Jacobi, B::Jacobi)
@@ -331,15 +314,7 @@ function \(A::Jacobi, B::Jacobi)
         else
             C2 = inv(_jacobi_convert_b(ba, ab, -kb, T))
         end
-        if iszero(ka)
-            C2
-        elseif iszero(kb)
-            C1
-        else
-            list = (C1, C2)
-            ApplyArray{T,2,typeof(*),typeof(list)}(*, list)
-            C1 * C2
-        end
+        C1 * C2
     else
         inv(B \ A)
     end
@@ -355,14 +330,8 @@ function \(w_A::WeightedJacobi, B::Jacobi)
     w_A \ (JacobiWeight(zero(a),zero(b)) .* B)
 end
 
-function \(A::AbstractJacobi, w_B::WeightedJacobi)
-    Ã = Jacobi(A)
-    (A \ Ã) * (Ã \ w_B)
-end
-function \(w_A::WeightedJacobi, B::AbstractJacobi)
-    B̃ = Jacobi(B)
-    (w_A \ B̃) * (B̃ \ B)
-end
+\(A::AbstractJacobi, w_B::WeightedJacobi) = Jacobi(A) \ w_B
+\(w_A::WeightedJacobi, B::AbstractJacobi) = w_A \ Jacobi(B)
 
 
 function broadcastbasis(::typeof(+), w_A::WeightedJacobi, w_B::WeightedJacobi)

From 8e6d7500ca1f8ec990ca3cb3f19b48d159d6942a Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Tue, 22 Oct 2024 08:42:33 +0100
Subject: [PATCH 7/8] Update jacobi.jl

---
 src/classical/jacobi.jl | 13 +++----------
 1 file changed, 3 insertions(+), 10 deletions(-)

diff --git a/src/classical/jacobi.jl b/src/classical/jacobi.jl
index 4cf4faa..6ab97ec 100644
--- a/src/classical/jacobi.jl
+++ b/src/classical/jacobi.jl
@@ -122,6 +122,7 @@ OrthogonalPolynomial(w::JacobiWeight) = Jacobi(w.a, w.b)
 orthogonalityweight(P::Jacobi) = JacobiWeight(P.a, P.b)
 
 const WeightedJacobi{T} = WeightedBasis{T,<:JacobiWeight,<:Jacobi}
+WeightedJacobi(P::Jacobi{T}) where T = JacobiWeight(zero(T),zero(T)) .* P
 
 """
     HalfWeighted{lr}(Jacobi(a,b))
@@ -320,16 +321,8 @@ function \(A::Jacobi, B::Jacobi)
     end
 end
 
-function \(A::Jacobi, w_B::WeightedJacobi)
-    a,b = A.a,A.b
-    (JacobiWeight(zero(a),zero(b)) .* A) \ w_B
-end
-
-function \(w_A::WeightedJacobi, B::Jacobi)
-    a,b = B.a,B.b
-    w_A \ (JacobiWeight(zero(a),zero(b)) .* B)
-end
-
+\(A::Jacobi, w_B::WeightedJacobi) = WeightedJacobi(A) \ w_B
+\(w_A::WeightedJacobi, B::Jacobi) = w_A \ WeightedJacobi(B)
 \(A::AbstractJacobi, w_B::WeightedJacobi) = Jacobi(A) \ w_B
 \(w_A::WeightedJacobi, B::AbstractJacobi) = w_A \ Jacobi(B)
 

From cf2ec38dd43e18ad09b03f6ba8723022fe53c9b2 Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Fri, 15 Nov 2024 17:08:35 +0000
Subject: [PATCH 8/8] fix

---
 src/classical/jacobi.jl | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/src/classical/jacobi.jl b/src/classical/jacobi.jl
index 4da3d15..188cec0 100644
--- a/src/classical/jacobi.jl
+++ b/src/classical/jacobi.jl
@@ -295,13 +295,17 @@ end
 
 function _jacobi_convert_a(a, b, k, T) # Jacobi(a+k, b) \ Jacobi(a, b)
     j = round(k)
-    @assert j ≈ k
+    if j ≉ k
+        throw(ArgumentError("non-integer conversions not supported"))
+    end
     k = Integer(j)
     reduce(*, [_jacobi_convert_a(a+j, b) for j in k-1:-1:0], init=Eye{T}(∞))
 end
 function _jacobi_convert_b(a, b, k, T) # Jacobi(a, b+k) \ Jacobi(a, b)
     j = round(k)
-    @assert j ≈ k
+    if j ≉ k
+        throw(ArgumentError("non-integer conversions not supported"))
+    end
     k = Integer(j)
     reduce(*, [_jacobi_convert_b(a, b+j) for j in k-1:-1:0], init=Eye{T}(∞))
 end