Math::polyval to use kernel, not polyev

PDLPorters · Feb 8, 2024 · a92cbf5 · a92cbf5
1 parent c327e5b
commit a92cbf5
Show file tree

Hide file tree

Showing 3 changed files with 14 additions and 12 deletions.
diff --git 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
@@ -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
@@ -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