# Reinch's modification to the Clenshaw recurrence gives no improvement

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:

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

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.

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