Math::polyroots now handles native-complex

Basic/Math/math.pd

@@ -342,23 +342,32 @@ pp_def("polyroots",
       Pars => 'cr(n); ci(n); [o]rr(m); [o]ri(m);',
       RedoDimsCode => '$SIZE(m) = $SIZE(n)-1;',
       GenericTypes => ['D'],
-      Code => '
-              int deg = (int)$SIZE(n)-1;
-              char *fail = cpoly($P(cr), $P(ci), deg, $P(rr), $P(ri));
-              if (fail)
-                 $CROAK("cpoly: %s", fail);
-',
-      , Doc => '
-
+      Code => <<'EOF',
+  char *fail = cpoly($P(cr), $P(ci), $SIZE(m), $P(rr), $P(ri));
+  if (fail)
+     $CROAK("cpoly: %s", fail);
+EOF
+      PMCode => <<'EOF',
+sub PDL::polyroots {
+  my @args = map PDL->topdl($_), @_;
+  my $natcplx = !$args[0]->type->real;
+  barf "need array context" if !$natcplx and !(wantarray//1);
+  splice @args, 0, 1, map $args[0]->$_, qw(re im) if $natcplx;
+  $_ //= PDL->null for @args[2,3];
+  PDL::_polyroots_int(@args);
+  $natcplx ? $args[2]->czip($args[3]) : @args[2,3];
+}
+EOF
+      Doc => '
 =for ref
 
 Complex roots of a complex polynomial, given coefficients in order
-of decreasing powers.
+of decreasing powers. As of 2.086, works with native-complex data.
 
 =for usage
 
  ($rr, $ri) = polyroots($cr, $ci);
-
+ $roots = polyroots($coeffs); # native complex
 ',);
 
 pp_addpm({At=>'Bot'},<<'EOD');

Changes

@@ -4,6 +4,7 @@
 - add MatrixOps::tritosquare
 - add gv command to shells
 - add PDL::Trans::VTable
+- Math::polyroots now handles native-complex
 
 2.085 2024-01-30
 - switch FFT code to use heap, not VLA (#436) - thanks @HaraldJoerg for report

t/math.t

@@ -53,12 +53,15 @@ $pa->inplace->erfi;
 ok( all( approx( $pa, pdl(0.00886,0.0) )), "erfi inplace" );
 }
 
-ok( all approx( qsort(
-	(polyroots(
-			pdl( 1,-55,1320,-18150,157773,-902055, 3416930,-8409500,12753576,-10628640,3628800 ),
-			zeroes(11)
-		))[0]
-), 1+sequence(10) ) );
+my $coeffs = pdl(cdouble, 1,-55,1320,-18150,157773,-902055, 3416930,-8409500,12753576,-10628640,3628800);
+my $roots = 1+sequence(10);
+my $got;
+ok all(approx $got=qsort((polyroots $coeffs->re, $coeffs->im)[0]), $roots), 'polyroots' or diag $got;
+polyroots $coeffs->re, $coeffs->im, $got=null; $got->inplace->qsort;
+ok all(approx $got, $roots), 'polyroots with explicit output args' or diag $got;
+ok all(approx $got=qsort(polyroots($coeffs)->re), $roots), 'polyroots native complex no output args' or diag $got;
+polyroots $coeffs, $got=null; $got=$got->re->qsort;
+ok all(approx $got, $roots), 'polyroots native complex explicit output args' or diag $got;
 
 {
 my $pa = sequence(41) - 20;