From 671bc20167e40f687b8a9613e0c0ec8ab56e1cbe Mon Sep 17 00:00:00 2001
From: Ed J <mohawk2@users.noreply.github.com>
Date: Wed, 7 Feb 2024 21:33:37 +0000
Subject: [PATCH] add Math::polyval

---
 Basic/Math/cpoly.c | 15 ++++++++-------
 Basic/Math/cpoly.h |  2 ++
 Basic/Math/math.pd | 46 ++++++++++++++++++++++++++++++++++++++++++++--
 Changes            |  1 +
 t/math.t           |  8 ++++++++
 5 files changed, 63 insertions(+), 9 deletions(-)

diff --git a/Basic/Math/cpoly.c b/Basic/Math/cpoly.c
index 697630522..bcb3c04a2 100644
--- a/Basic/Math/cpoly.c
+++ b/Basic/Math/cpoly.c
@@ -24,8 +24,6 @@ static int fxshft(int l2, complex double *zc, int nn, complex double *shc, compl
 static int vrshft(int l3, complex double *zc, int nn, complex double *qpc, complex double *pc, complex double *qhc, complex double *hc, complex double *tc, complex double *sc, complex double *pvc);
 static int calct(int nn, complex double *qhc, complex double *hc, complex double *tc, complex double *sc, complex double *pvc);
 static void nexth(int boolvar, int nn, complex double *qhc, complex double *qpc, complex double *hc, complex double tc);
-static complex double polyev(int nn, complex double sc, complex double pc[],
-	    complex double qc[]);
 static double errev(int nn, complex double qc[], double ms, double mp);
 static double cauchy(int nn, complex double pc[]);
 static double scale(int nn, complex double pc[]);
@@ -513,16 +511,19 @@ static void nexth(int boolvar, int nn, complex double *qhc, complex double *qpc,
   }
 }
 
-static complex double polyev(int nn, complex double sc, complex double pc[],
+complex double polyev(int nn, complex double sc, complex double pc[],
 	    complex double qc[])
      /* Evaluates a polynomial  p  at  s  by the Horner recurrence,
-	placing the partial sums in q and the computed value in pv
+	optionally placing the partial sums in q, returns the computed value
      */
 {
   int i;
-  complex double vc = qc[0] = pc[0];
-  for (i=1;i<nn;i++)
-    qc[i] = vc = vc*sc + pc[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;
+  }
   return vc;
 }
 
diff --git a/Basic/Math/cpoly.h b/Basic/Math/cpoly.h
index 5d1d20768..cf66309ed 100644
--- a/Basic/Math/cpoly.h
+++ b/Basic/Math/cpoly.h
@@ -4,6 +4,8 @@ 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 f2cf7ef9d..c244492d0 100644
--- a/Basic/Math/math.pd
+++ b/Basic/Math/math.pd
@@ -347,7 +347,7 @@ pp_def("polyroots",
   if (fail)
      $CROAK("cpoly: %s", fail);
 EOF
-      PMCode => <<'EOF',
+      PMCode => pp_line_numbers(__LINE__, <<'EOF'),
 sub PDL::polyroots {
   my @args = map PDL->topdl($_), @_;
   my $natcplx = !$args[0]->type->real;
@@ -366,8 +366,50 @@ of decreasing powers. As of 2.086, works with native-complex data.
 
 =for usage
 
- ($rr, $ri) = polyroots($cr, $ci);
  $roots = polyroots($coeffs); # native complex
+ ($rr, $ri) = polyroots($cr, $ci);
+',);
+
+pp_def("polyval",
+      Pars => 'c(n); x(); [o]y();',
+      GenericTypes => ['CD'],
+      Code => '$y() = polyev($SIZE(n), $x(), $P(c), NULL);',
+      PMCode => pp_line_numbers(__LINE__, <<'EOF'),
+sub PDL::polyval {
+  my @args = map PDL->topdl($_), @_;
+  my $natcplx = !$args[0]->type->real;
+  barf "need array context" if !$natcplx and !(wantarray//1);
+  if (!$natcplx) {
+    splice @args, 0, 2, $args[0]->czip($args[1]); # c
+    splice @args, 1, 2, $args[1]->czip($args[2]); # x
+  }
+  my @ins = splice @args, 0, 2;
+  my $explicit_out = my @outs = @args;
+  if ($natcplx) {
+    $_ //= PDL->null for $outs[0];
+  } else {
+    $_ //= PDL->null for @outs[0,1];
+  }
+  my @args_out = $natcplx ? @outs : PDL->null;
+  PDL::_polyval_int(@ins, @args_out);
+  if (!$natcplx) {
+    $outs[0] .= $args_out[0]->re;
+    $outs[1] .= $args_out[0]->im;
+  }
+  $natcplx ? $outs[0] : @outs;
+}
+EOF
+      Doc => '
+=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.
+
+=for usage
+
+ $y = polyval($coeffs, $x); # native complex
+ ($yr, $yi) = polyval($cr, $ci, $xr, $xi);
 ',);
 
 pp_addpm({At=>'Bot'},<<'EOD');
diff --git a/Changes b/Changes
index bfeb358be..25588dfef 100644
--- a/Changes
+++ b/Changes
@@ -5,6 +5,7 @@
 - add gv command to shells
 - add PDL::Trans::VTable
 - Math::polyroots now handles native-complex
+- add Math::polyval
 
 2.085 2024-01-30
 - switch FFT code to use heap, not VLA (#436) - thanks @HaraldJoerg for report
diff --git a/t/math.t b/t/math.t
index b7441e106..c99435150 100644
--- a/t/math.t
+++ b/t/math.t
@@ -63,6 +63,14 @@ ok all(approx $got=qsort(polyroots($coeffs)->re), $roots), 'polyroots native com
 polyroots $coeffs, $got=null; $got=$got->re->qsort;
 ok all(approx $got, $roots), 'polyroots native complex explicit output args' or diag $got;
 
+my ($coeffs2, $x, $exp_val) = (cdouble(3,2,1), cdouble(5,7,9), cdouble(86,162,262));
+ok all(approx $got=polyval($coeffs2, $x), $exp_val), 'polyval natcom no output' or diag $got;
+polyval($coeffs2, $x, $got=null);
+ok all(approx $got, $exp_val), 'polyval natcom explicit output' or diag $got;
+ok all(approx $got=(polyval($coeffs2->re, zeroes(3), $x->re, zeroes(3)))[0], $exp_val->re), 'polyval legacy no output' or diag $got;
+polyval($coeffs2->re, zeroes(3), $x->re, zeroes(3), $got=null);
+ok all(approx $got, $exp_val->re), 'polyval legacy explicit output' or diag $got;
+
 {
 my $pa = sequence(41) - 20;
 $pa /= 4;