simplify factorials; change scalar_product_function

src/chebyshev.jl

@@ -21,18 +21,15 @@ const ChebyshevBasis{P} = ChebyshevBasisFirstKind{P}
 
 degree_one_univariate_polynomial(::Type{<:ChebyshevBasisFirstKind}, variable::MP.AbstractVariable) = MA.@rewrite(variable + 0)
 
-function scalar_product_function(::Type{<:ChebyshevBasisFirstKind})
-    function sp(i::Int)
-        if i == 0
-            return π
-        elseif mod(i, 2) == 1
-            return 0
-        else
-            n = Int(i/2)
-            return π*factorial(i)/(2^i * factorial(n)^2)
-        end
+function scalar_product_function(::Type{<:ChebyshevBasisFirstKind}, i::Int)
+    if i == 0
+        return π
+    elseif isodd(i)
+        return 0
+    else
+        n = div(i, 2)
+        return (π/2^i)*prod(n+1:i)/factorial(n)
     end
-    return sp
 end
 
 """
@@ -48,18 +45,13 @@ end
 
 degree_one_univariate_polynomial(::Type{<:ChebyshevBasisSecondKind}, variable::MP.AbstractVariable) = MA.@rewrite(2variable + 0)
 
-
-function scalar_product_function(::Type{<:ChebyshevBasisSecondKind})
-    function sp(i::Int)
-        if i == 0
-            return π/2
-        elseif mod(i, 2) == 1
-            return 0
-        else
-            n = Int(i/2)
-            return π*factorial(i)/(2^(i+1) * factorial(n) * factorial(n +1))
-        end
+function scalar_product_function(::Type{<:ChebyshevBasisSecondKind}, i::Int)
+    if i == 0
+        return π/2
+    elseif isodd(i)
+        return 0
+    else
+        n = div(i, 2)
+        return π/(2^(i+1))*prod(n+2:i)/factorial(n)
     end
-    return sp
 end
-

src/hermite.jl

@@ -21,18 +21,15 @@ reccurence_first_coef(::Type{<:ProbabilistsHermiteBasis}, degree) = 1
 reccurence_third_coef(::Type{<:ProbabilistsHermiteBasis}, degree) = -(degree - 1)
 degree_one_univariate_polynomial(::Type{<:ProbabilistsHermiteBasis}, variable::MP.AbstractVariable) = MA.@rewrite(1variable)
 
-function scalar_product_function(::Type{<:ProbabilistsHermiteBasis})
-    function sp(i::Int)
-       if i == 0
-           return √(2*π)
-       elseif mod(i, 2) == 1
-           return 0
-       else
-           n = Int(i/2)
-           return factorial(2*n)/((2^n)*factorial(n))*√(2*π)
-       end
-   end
-   return sp
+function scalar_product_function(::Type{<:ProbabilistsHermiteBasis}, i::Int)
+    if i == 0
+        return √(2*π)
+    elseif isodd(i)
+        return 0
+    else
+        n = div(i, 2)
+        return (√(2*π)/(2^n))*prod(n+1:2*n)
+    end
 end
 
 """
@@ -49,16 +46,13 @@ reccurence_first_coef(::Type{<:PhysicistsHermiteBasis}, degree) = 2
 reccurence_third_coef(::Type{<:PhysicistsHermiteBasis}, degree) = -2(degree - 1)
 degree_one_univariate_polynomial(::Type{<:PhysicistsHermiteBasis}, variable::MP.AbstractVariable) = MA.@rewrite(2variable)
 
-function scalar_product_function(::Type{<:PhysicistsHermiteBasis})
-    function sp(i::Int)
-       if i == 0
-            return √(π)
-        elseif mod(i, 2) == 1
-            return 0
-        else
-            n = Int(i/2)
-            return √(π)*factorial(i)/(factorial(n)*2^i)
-        end
+function scalar_product_function(::Type{<:PhysicistsHermiteBasis}, i::Int)
+    if i == 0
+        return √(π)
+    elseif isodd(i)
+        return 0
+    else
+        n = div(i, 2)
+        return (√(π)/(2^i))*prod(n+1:i)
     end
-    return sp
 end

src/laguerre.jl

@@ -20,9 +20,7 @@ reccurence_deno_coef(::Type{<:LaguerreBasis}, degree) = degree
 
 degree_one_univariate_polynomial(::Type{<:LaguerreBasis}, variable::MP.AbstractVariable) = MA.@rewrite(1 - variable)
 
-function scalar_product_function(::Type{<:LaguerreBasis})
-    function sp(i::Int)
-        return factorial(i) 
-    end
+function scalar_product_function(::Type{<:LaguerreBasis}, i::Int)
+    return factorial(i) 
 end
 

src/legendre.jl

@@ -29,14 +29,10 @@ reccurence_deno_coef(::Type{<:LegendreBasis}, degree) = degree
 
 degree_one_univariate_polynomial(::Type{<:LegendreBasis}, variable::MP.AbstractVariable) = MA.@rewrite(variable + 0)
 
-function scalar_product_function(::Type{<:LegendreBasis})
-    function sp(i::Int)
-        if mod(i, 2) == 1
-            return 0
-        else
-            return 2/(i+1)
-        end
+function scalar_product_function(::Type{<:LegendreBasis}, i::Int)
+    if isodd(i)
+        return 0
+    else
+        return 2/(i+1)
     end
-    return sp
 end
-

src/orthogonal.jl

@@ -127,20 +127,20 @@ function basis_covering_monomials(B::Type{<:AbstractMultipleOrthogonalBasis}, mo
     return _basis_from_monomials(B, variables(monos), MP.monovec(collect(m)))
 end
 
-function scalar_product_function(::Type{<:AbstractMultipleOrthogonalBasis}) end
+function scalar_product_function(::Type{<:AbstractMultipleOrthogonalBasis}, i::Int) end
 
 LinearAlgebra.dot(p, q, basis_type::Type{<:AbstractMultipleOrthogonalBasis}) = _integral(p*q, basis_type)
 
 function _integral(p::Number, basis_type::Type{<:AbstractMultipleOrthogonalBasis})
-    return p*scalar_product_function(basis_type)(0)
+    return p*scalar_product_function(basis_type, 0)
 end
 
 function _integral(p::MP.AbstractVariable, basis_type::Type{<:AbstractMultipleOrthogonalBasis})
-    return scalar_product_function(basis_type)(1)
+    return scalar_product_function(basis_type, 1)
 end
 
 function _integral(p::MP.AbstractMonomial, basis_type::Type{<:AbstractMultipleOrthogonalBasis})
-    return prod([scalar_product_function(basis_type)(i) for i in exponents(p)])
+    return prod([scalar_product_function(basis_type, i) for i in exponents(p)])
 end
 
 function _integral(p::MP.AbstractTerm, basis_type::Type{<:AbstractMultipleOrthogonalBasis})

test/chebyshev.jl

@@ -9,7 +9,7 @@ using MultivariateBases
         4x^3 - 3x,
         8x^4 - 8x^2 + 1
     ], true)
-    univ_orthogonal_test(ChebyshevBasis, i -> π*(1/2 + (0^i)/2))
+    univ_orthogonal_test(ChebyshevBasis, i -> π*(1/2 + (0^i)/2), atol = 1e-12)
     orthogonal_test(ChebyshevBasisSecondKind, x -> [
         1,
         2x,
@@ -18,7 +18,7 @@ using MultivariateBases
         16x^4 - 12x^2 + 1
     ], true)
 
-    univ_orthogonal_test(ChebyshevBasisSecondKind, i -> π/2)
+    univ_orthogonal_test(ChebyshevBasisSecondKind, i -> π/2, atol = 1e-12)
 end
 
 @testset "API degree = $degree" for degree in 0:3