0
$\begingroup$

The classical Clenshaw recurrence (see Algorithm 3.1 here) is less accurate as $x\to \pm 1$. So Reich proposed a modification to it, which is discussed as Algorithm 3.2 as well as by Oliver.

While trying to derive a point of transition between the two algorithms I generated ULP plots for these two algorithms, and was surprised to find I could not see any difference at all, anywhere.

First I tried 64 Unif(0,1) coefficients, with perturbed abscissas and perturbed coefficients:

modified_clenshaw_unif01_perturbed_abscissas

classic_clenshaw

Next, 64 Unif(0,1) exact coefficients, perturbed abscissas:

modified_exact_u01

enter image description here

I also tried Unif(-1,1) coefficients, perturbed and non-perturbed, as well as decaying coefficients more representative of smooth functions. But each time there was no difference in the appearance of the ULP plot.

What data will actually demonstrate the superiority of Reinch's modification to the Clenshaw recurrence?

My code, for reference:

template<class Real>
inline Real chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const Real& x)
{
    if (length < 2)
    {
        if (length == 0)
        {
            return 0;
        }
        return c[0]/2;
    }
    Real b2 = 0;
    Real b1 = c[length -1];
    for(size_t j = length - 2; j >= 1; --j)
    {
        Real tmp = 2*x*b1 - b2 + c[j];
        b2 = b1;
        b1 = tmp;
    }
    return x*b1 - b2 + c[0]/2;
}

template<class Real>
inline Real modified_chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const Real& x)
{
    using std::abs;
    if (length < 2)
    {
        if (length == 0)
        {
            return 0;
        }
        return c[0]/2;
    }
    Real cutoff = 0.0;
    if (abs(x) <= cutoff)
    {
        Real b2 = 0;
        Real b1 = c[length -1];
        for(size_t j = length - 2; j >= 1; --j)
        {
            Real tmp = 2*x*b1 - b2 + c[j];
            b2 = b1;
            b1 = tmp;
        }
        return x*b1 - b2 + c[0]/2;
    }
    else if (x > cutoff)
    {
        Real b = c[length -1];
        Real d = b;
        Real b2 = b;
        for (size_t r = length - 2; r >= 1; --r)
        {
            d = 2*(x-1)*b + d + c[r];
            b2 = b;
            b = d + b;
        }
        return x*b - b2 + c[0]/2;
    }
    else
    {
        Real b = c[length -1];
        Real d = b;
        Real b2 = b;
        for (size_t r = length - 2; r >= 1; --r)
        {
            d = 2*(x+1)*b - d + c[r];
            b2 = b;
            b = d - b;
        }
        return x*b - b2 + c[0]/2;
    }
}

Also, I perturbed coefficients and abscissas by generating them in long double and rounding them to float, and calculated the "exact" number with respect to long double precision.

$\endgroup$
1
$\begingroup$

I have found that in fact the modified version is better about 10% of the time under a broad set of perturbations (abscissas, coefficients of various distributions, so on), but it's quite difficult to see. I had to simply brute force over all 32 bit representables and just calculate how many times one version was better than the other, as well as compute total ULP distance.

However, if you need to translate your Chebyshev series from $[a,b]\to [-1,1]$, then the superiority becomes apparent:

enter image description here

(On the lhs, the orange abscissas cover the blue ones; this is an artifact of the rendering order; they accuracy is roughly the same on the lhs.)

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.