add compensated_horner method

Author: jverzani

Project.toml

@@ -2,7 +2,7 @@ name = "Polynomials"
 uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 license = "MIT"
 author = "JuliaMath"
-version = "1.1.6"
+version = "1.1.7"
 
 [deps]
 Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5"

src/polynomials/standard-basis.jl

@@ -297,3 +297,91 @@ function vander(P::Type{<:StandardBasisPolynomial}, x::AbstractVector{T}, n::Int
     end
     return A
 end
+
+
+## --------------------------------------------------
+"""
+    compensated_horner(p::P, x)
+    compensated_horner(ps, x)
+
+Evaluate `p(x)` using a compensation scheme of S. Graillat, Ph. Langlois, N. Louve [Compensated Horner Scheme](https://cadxfem.org/cao/Compensation-horner.pdf). Either a `Polynomial` `p` or it coefficients may be passed in.
+
+The Horner scheme has relative error given by
+
+`|(p(x) - p̂(x))/p(x)| ≤ α ⋅ u ⋅ cond(p, x)`, where `u` is the precision (`2⁻⁵³`).
+
+The compensated Horner scheme has relative error bounded by
+
+`|(p(x) - p̂(x))/p(x)| ≤  u +  β(n) ⋅ u² ⋅ cond(p, x)`.
+
+As noted, this reflects the accuracy of a backward stable computation performed in doubled working precision `u²`. (E.g., polynomial evaluation of a `Polynomial{Float64}` object through `compensated_horner` is as accurate as evaluation of a `Polynomial{Double64}` object (using the `DoubleFloat` package), but significantly faster.
+
+Pointed out in https://discourse.julialang.org/t/more-accurate-evalpoly/45932/5.
+"""
+function compensated_horner(p::P, x) where {T, P <: Polynomials.StandardBasisPolynomial{T}}
+    compensated_horner(coeffs(p), x)
+end
+
+# rexpressed from paper to compute horner_sum in same pass
+# sᵢ -- horner sum
+# c -- compensating term
+@inline function compensated_horner(ps, x) 
+    n, T = length(ps), eltype(ps)
+    aᵢ = ps[end]
+    sᵢ = aᵢ * _one(x)
+    c = zero(T) * _one(x)
+    for i in n-1:-1:1
+	aᵢ = ps[i]
+        pᵢ, πᵢ = two_product_fma(sᵢ, x)
+	sᵢ, σᵢ = two_sum(pᵢ, aᵢ)
+        c = fma(c, x, πᵢ + σᵢ)
+    end
+    sᵢ + c
+end
+
+function compensated_horner(ps::Tuple, x::S) where {S}
+    ps == () && return zero(S)
+    if @generated
+        n = length(ps.parameters)
+        sσᵢ =:(ps[end] * _one(x), zero(S))
+        c = :(zero(S) * _one(x))
+        for i in n-1:-1:1
+            pπᵢ = :(two_product_fma($sσᵢ[1], x))
+	    sσᵢ = :(two_sum($pπᵢ[1], ps[$i]))
+            Σ = :($pπᵢ[2] + $sσᵢ[2])
+            c = :(fma($c, x, $Σ))
+        end
+        s = :($sσᵢ[1] + $c)
+        s
+    else
+        compensated_horner(ps, x)
+    end
+end
+
+# error-free-transformations (EFT) of a∘b = x + y where x, y floating point, ∘ mathematical
+# operator (not ⊕ or ⊗)
+@inline function two_sum(a, b)
+    x = a + b
+    z = x - a
+    y = (a - (x-z)) + (b-z)
+    x, y
+end
+
+# return x, y, floats,  with x + y = a * b
+@inline function two_product_fma(a, b)
+    x = a * b
+    y = fma(a, b, -x)
+    x, y
+end
+
+
+# Condition number of a standard basis polynomial
+# rule of thumb: p̂ a compute value
+# |p(x) - p̃(x)|/|p(x)| ≤ α(n)⋅u ⋅ cond(p,x), where u = finite precision of compuation (2^-p)
+function LinearAlgebra.cond(p::P, x) where {P <: Polynomials.StandardBasisPolynomial}
+    p̃ = Polynomials.:⟒(P)(abs.(coeffs(p)), p.var)
+    p̃(abs(x))/ abs(p(x))
+end
+
+
+

test/StandardBasis.jl

@@ -342,8 +342,10 @@ end
         # issue #209
         ps  = [P([0,1]), P([0,0,1])]
         @test Polynomials.evalpoly.(1/2, ps) ≈ [p(1/2)  for  p  in ps]
+
     end
 
+    
     # constant polynomials and type
     Ts = (Int, Float32, Float64, Complex{Int}, Complex{Float64})
     for P in (Polynomial, ImmutablePolynomial, SparsePolynomial)
@@ -374,7 +376,18 @@ end
         @test eltype(p3) == eltype(p2)
     end
 
-
+    # compensated_horner
+    # polynomial evaluation for polynomials with large condition numbers
+    for P in (Polynomial, ImmutablePolynomial, SparsePolynomial)
+        x = variable(P{Float64})
+        f(x) = (x - 1)^20
+        p = f(x)
+        e₁ = abs( (f(4/3) - p(4/3))/ p(4/3) )
+        e₂ = abs( (f(4/3) - Polynomials.compensated_horner(p, 4/3))/ p(4/3) )
+        @test cond(p, 4/3) > 1/eps()
+        @test e₁ > sqrt(eps())
+        @test e₂ <= 4eps()        
+    end
 end
 
 @testset "Conversion" begin