BUG: function p1evl in povevl.h not making what's described

[semoi]

Describe your issue.

The description of the function p1evl says if consider the coefficient of x^N is 1.0.
But the description at the beginning of the special/cephes/polevl.h file says it assumes coef[N] = 1.0
And coef[N] is not the coefficient of x^N, but of x^0 (because coeff are stored in reverse order)

I don't know what is the right description, but the code seems to do none of them.

special/cephes/polevl.h

* DESCRIPTION:
 *
 * Evaluates polynomial of degree N:
 *
 *                     2          N
 * y  =  C  + C x + C x  +...+ C x
 *        0    1     2          N
 *
 * Coefficients are stored in reverse order:
 *
 * coef[0] = C  , ..., coef[N] = C  .
 *            N                   0
 *
 *  The function p1evl() assumes that coef[N] = 1.0 and is
 * omitted from the array.  Its calling arguments are
 * otherwise the same as polevl().

/*                                                     p1evl() */
/*                                          
 * Evaluate polynomial when coefficient of x^N  is 1.0.
 * Otherwise same as polevl.
 */

static NPY_INLINE double p1evl(double x, const double coef[], int N)
{
    double ans;
    const double *p;
    int i;

    p = coef;
    ans = x + *p++;
    i = N - 1;

    do
	ans = ans * x + *p++;
    while (--i);

    return (ans);
}

Reproducing Code Example

I've tested with the T[] array and polevl(x, T, 4) and p1evl(x, T, 4) from the special/cephes/ndtr.c file.

coeff = [9, 90, 2232, 7003, 55592] (rounded here for clarity)
Polinomial expected for N=4 : 55592 + 7003x + 2232x2 + 90x3 + 9x4

Function polevl(x, T, 4) :
ans=9
i=N=4
do/while loop :
ans=9*x + 90
i=3
ans=(9*x + 90)*x + 2232
i=2
ans=((9*x + 90)*x + 2232)*x + 7003
i=1
ans=(((9*x + 90)*x + 2232)*x + 7003)*x + 55592

Function p1evl(x, T, 4) :
ans=x + 9
i=N-1=3
do/while loop :
ans=(x + 9)*x + 90
i=2
ans=((x + 9)*x + 90)*x + 2232
i=1
ans=(((x + 9)*x + 90)*x + 2232)*x + 7003

So we see that the coeff of x^4 is 1 (which is ok with 2nd description), but we have a shift of the coeff as we have 7003 instead of 7003x and 55592 not used.

Maybe it's normal because of the 1st description says "p1evl() assumes that coef[N] = 1.0 and is omitted from the array"

Thanks in advance for your feedback on this behaviour / description problem.

Error message

.

SciPy/NumPy/Python version information

1.8.0 1.21.5 sys.version_info(major=3, minor=9, micro=12, releaselevel='final', serial=0)

[WarrenWeckesser]

@semoi, sorry for the lack of response here.

You're not the first person to be confused by p1evl. It appears that the description

 *  The function p1evl() assumes that coef[N] = 1.0 and is
 * omitted from the array. 

is not correct. Presumably the author intended that to be "The function p1evl() assumes that C_N = 1.0", where C_N is the coefficient of $x^N$, as shown in the previous lines:

 *                     2          N
 * y  =  C  + C x + C x  +...+ C x
 *        0    1     2          N
 *
 * Coefficients are stored in reverse order:
 *
 * coef[0] = C  , ..., coef[N] = C  .
 *            N                   0

So for p1evl, the input array coef must have length N (rather than N+1 as in polevl), and the values in coeff must be

    coef[0] = C_{N-1}
    coef[1] = C_{N-2}
    ...
    coef[N-2] = C_1
    coef[N-1] = C_0

I created a pull request to improve the comments about p1evl: #17535