From a92cbf583f152095b3a78cf47530b23d75be27bd Mon Sep 17 00:00:00 2001
From: Ed J <mohawk2@users.noreply.github.com>
Date: Thu, 8 Feb 2024 20:27:24 +0000
Subject: [PATCH] Math::polyval to use kernel, not polyev

---
 Basic/Math/cpoly.c | 13 ++++++-------
 Basic/Math/cpoly.h |  2 --
 Basic/Math/math.pd | 11 ++++++++---
 3 files changed, 14 insertions(+), 12 deletions(-)

diff --git a/Basic/Math/cpoly.c b/Basic/Math/cpoly.c
index 94a66e566..5a84030ea 100644
--- a/Basic/Math/cpoly.c
+++ b/Basic/Math/cpoly.c
@@ -24,6 +24,7 @@ static int fxshft(int l2, int nn, complex double shc[], complex double qpc[], co
 static int vrshft(int l3, int nn, complex double qpc[], complex double pc[], complex double qhc[], complex double hc[], complex double *tc, complex double *sc, complex double *pvc, complex double *zc);
 static int calct(int nn, complex double sc, complex double pvc, complex double qhc[], complex double hc[], complex double *tc);
 static void nexth(int boolvar, int nn, complex double tc, complex double qhc[], complex double qpc[], complex double hc[]);
+static complex double polyev(int nn, complex double sc, complex double pc[], complex double qc[]);
 static double errev(int nn, double ms, double mp, complex double qc[]);
 static double cauchy(int nn, complex double pc[]);
 static double scale(int nn, complex double pc[]);
@@ -514,19 +515,17 @@ static void nexth(int boolvar, int nn, complex double tc, complex double qhc[],
   }
 }
 
-complex double polyev(int nn, complex double sc, complex double pc[],
+static complex double polyev(int nn, complex double sc, complex double pc[],
 	    complex double qc[])
      /* Evaluates a polynomial  p  at  s  by the Horner recurrence,
-	optionally placing the partial sums in q, returns the computed value
+	placing the partial sums in q, returns the computed value
      */
 {
   int i;
   complex double vc = pc[0];
-  if (qc) qc[0] = vc;
-  for (i=1;i<nn;i++) {
-    vc = vc*sc + pc[i];
-    if (qc) qc[i] = vc;
-  }
+  qc[0] = vc;
+  for (i=1;i<nn;i++)
+    qc[i] = vc = vc*sc + pc[i];
   return vc;
 }
 
diff --git a/Basic/Math/cpoly.h b/Basic/Math/cpoly.h
index cf66309ed..5d1d20768 100644
--- a/Basic/Math/cpoly.h
+++ b/Basic/Math/cpoly.h
@@ -4,8 +4,6 @@ void prtz(int n,double zr[], double zi[]);
 #endif
 char *cpoly(double opr[], double opi[], int degree,
 	   double zeror[], double zeroi[]);
-complex double polyev(int nn, complex double sc, complex double pc[],
-	   complex double qc[]);
 
 #if !defined(FALSE)
 #define FALSE (0)
diff --git a/Basic/Math/math.pd b/Basic/Math/math.pd
index 88a49db6e..364f68276 100644
--- a/Basic/Math/math.pd
+++ b/Basic/Math/math.pd
@@ -427,7 +427,12 @@ native-complex data. Currently C<O(n^2)>.
 pp_def("polyval",
       Pars => 'c(n); x(); [o]y();',
       GenericTypes => ['CD'],
-      Code => '$y() = polyev($SIZE(n), $x(), $P(c), NULL);',
+      Code => <<'EOF',
+$GENERIC(y) vc = $c(n=>0), sc = $x();
+PDL_Indx i;
+for (i=1; i<$SIZE(n); i++) vc = vc*sc + $c(n=>i);
+$y() = vc;
+EOF
       PMCode => pp_line_numbers(__LINE__, <<'EOF'),
 sub PDL::polyval {
   my @args = map PDL->topdl($_), @_;
@@ -457,8 +462,8 @@ EOF
 =for ref
 
 Complex value of a complex polynomial at given point, given coefficients
-in order of decreasing powers. Added in 2.086, works with native-complex
-data.
+in order of decreasing powers. Uses Horner recurrence. Added in 2.086,
+works with native-complex data.
 
 =for usage