add Pars type-spec "!real", fix Math::polyroots - #511

Changes

@@ -1,6 +1,8 @@
 - fix test that assumed acosh(0)->byte, i.e. nan()->byte, was always 0 (#514) - thanks @eserte for report
 - separate PDL::Type POD documentation
 - partly restoring pre-2.096 xform type-selection: if xform given no typed outputs, and non-available (greater than last-given type) typed inputs, use last-given (#511, https://github.com/moocow-the-bovine/PDL-CCS/issues/18)
+- add Pars type-spec "!real" which makes it an error to supply real values (#511)
+- fix Math::polyroots with native-complex input and supplied null output
 
 2.098 2025-01-03
 - fix Windows build problems

lib/PDL/Core/pdl.h.PL

@@ -293,6 +293,8 @@ PDL_TYPELIST_ALL(X)
 #define  PDL_PARENTDATACHANGED   (1 << 1)
   /* Parent dims or incs has been altered without changing this pdl.        */
 #define  PDL_PARENTDIMSCHANGED   (1 << 2)
+  /* Marked read-only by user; throw error if given as output to xform.     */
+#define  PDL_READONLY            (1 << 3)
   /* Physical data representation of the parent has changed (e.g.           */
   /* physical transposition), so incs etc. need to be recalculated.         */
 #define  PDL_ANYCHANGED          (PDL_PARENTDATACHANGED|PDL_PARENTDIMSCHANGED)
@@ -327,6 +329,7 @@ PDL_TYPELIST_ALL(X)
  X(PDL_ALLOCATED) \
  X(PDL_PARENTDATACHANGED) \
  X(PDL_PARENTDIMSCHANGED) \
+ X(PDL_READONLY) \
  X(PDL_DATAFLOW_F) \
  X(PDL_NOMYDIMS) \
  X(PDL_MYDIMS_TRANS) \
@@ -419,6 +422,7 @@ typedef struct pdl_transvtable {
 #define PDL_PARAM_ISIGNORE           (1 << 10)
 #define PDL_PARAM_ISNOTCOMPLEX       (1 << 11)
 #define PDL_PARAM_ALLOW_NULL         (1 << 12)
+#define PDL_PARAM_ISNOTREAL          (1 << 13)
 
 #define PDL_LIST_FLAGS_PARAMS(X) \
  X(PDL_PARAM_ISREAL) \
@@ -433,7 +437,8 @@ typedef struct pdl_transvtable {
  X(PDL_PARAM_ISPHYS) \
  X(PDL_PARAM_ISIGNORE) \
  X(PDL_PARAM_ISNOTCOMPLEX) \
- X(PDL_PARAM_ALLOW_NULL)
+ X(PDL_PARAM_ALLOW_NULL) \
+ X(PDL_PARAM_ISNOTREAL)
 
 /* All trans must start with this */
 

lib/PDL/Core/pdlapi.c

@@ -1126,6 +1126,10 @@ static inline pdl_error pdl__transtype_select(
       continue;
     short flags = vtable->par_flags[i];
     pdl_datatypes dtype = pdl->datatype;
+    if (flags & PDL_PARAM_ISNOTREAL && dtype < PDL_CF)
+      return pdl_make_error(PDL_EUSERERROR,
+        "%s: ndarray %s must be complex, but is type %s",
+        vtable->name, vtable->par_names[i], PDL_TYPENAME(dtype));
     if (flags & PDL_PARAM_ISNOTCOMPLEX && dtype >= PDL_CF)
       return pdl_make_error(PDL_EUSERERROR,
         "%s: ndarray %s must be real, but is type %s",

lib/PDL/Math.pd

@@ -360,11 +360,20 @@ EOF
 sub PDL::polyroots {
   my @args = map PDL->topdl($_), @_;
   my $natcplx = !$args[0]->type->real;
-  barf "need array context" if !$natcplx and !(wantarray//1);
+  barf "need array context if give real data and no outputs"
+    if !$natcplx and @_ < 3 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];
+  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 ? (map PDL->null, 1..2) : @outs; # opposite from polyfromroots
+  PDL::_polyroots_int(@ins, @args_out);
+  return @args_out if !$natcplx;
+  $outs[0] .= PDL::czip(@args_out[0,1]);
 }
 EOF
       Doc => '
@@ -379,17 +388,15 @@ As of 2.086, works with native-complex data.
 =for usage
 
  $roots = polyroots($coeffs); # native complex
+ polyroots($coeffs, $roots=null); # native complex
  ($rr, $ri) = polyroots($cr, $ci);
+ polyroots($cr, $ci, $rr, $ri);
 ',);
 
 pp_def("polyfromroots",
       Pars => 'r(m); [o]c(n=CALC($SIZE(m)+1));',
       GenericTypes => ['CD'],
       Code => <<'EOF',
-/*
-algorithm from octave poly.m, O(n^2)
-per https://cs.stackexchange.com/questions/116643/what-is-the-most-efficient-algorithm-to-compute-polynomial-coefficients-from-its, FFT allows O(n*log(n)^2)
-*/
 $c(n=>0) = 1.0;
 loop(m) %{ $c(n=>m+1) = 0.0; %}
 PDL_Indx k;
@@ -427,7 +434,11 @@ EOF
 
 Calculates the complex coefficients of a polynomial from its complex
 roots, in order of decreasing powers. Added in 2.086, works with
-native-complex data. Currently C<O(n^2)>.
+native-complex data.
+
+Algorithm is from Octave poly.m, O(n^2), per
+L<https://cs.stackexchange.com/questions/116643/what-is-the-most-efficient-algorithm-to-compute-polynomial-coefficients-from-its>;
+using an FFT would allow O(n*log(n)^2).
 
 =for usage
 

lib/PDL/Ops.pd

@@ -480,7 +480,7 @@ sub cfunc {
     my $codestr = pp_line_numbers(__LINE__-1,"\$b() = $func(\$complexv());");
     pp_def($name,
 	   GenericTypes=>$C,
-	   Pars => 'complexv(); '.($make_real ? 'real' : '').' [o]b()',
+	   Pars => ($make_real ? '!real ' : '').'complexv(); '.($make_real ? 'real' : '').' [o]b()',
 	   HandleBad => 1,
 	   NoBadifNaN => 1,
 	   ($make_real ? () : (Inplace => 1)),
@@ -506,7 +506,7 @@ cfunc('carg', 'carg', 1, 'Returns the polar angle of a complex number.', undef);
 cfunc('conj', 'conj', 0, 'complex conjugate.', undef);
 
 pp_def('czip',
-  Pars => 'r(); i(); complex [o]c()',
+  Pars => '!complex r(); !complex i(); complex [o]c()',
   Doc => <<'EOF',
 convert real, imaginary to native complex, (sort of) like LISP zip
 function. Will add the C<r> ndarray to "i" times the C<i> ndarray. Only

lib/PDL/PP.pod

@@ -912,6 +912,9 @@ and C<Doc> omitted):
 As of 2.096, there is also a C<!complex> type, which means if a
 complex-valued ndarray is supplied to the operation, an error will
 be thrown. In the normal case, its type will be as if C<real> were given.
+As of 2.099, there is further a C<!real> type, which is the opposite
+of the above; the idea is to protect a user who e.g. through a bug
+supplied a real ndarray to L<PDL::Ops/im>.
 
 Finally, there are the C<type+> qualifiers. Let's illustrate the C<int+>
 qualifier with the actual definition of sumover:

lib/PDL/PP/PdlParObj.pm

@@ -11,7 +11,7 @@ our %INVALID_PAR = map +($_=>1), qw(
 );
 
 my $typeregex = join '|', map $_->ppforcetype, types;
-my $complex_regex = join '|', qw(real complex !complex);
+my $complex_regex = join '|', qw(real complex !real !complex);
 our $sqbr_re = qr/\[([^]]*)\]/x;
 our $pars_re = qr/^
 	\s*(?:($complex_regex|$typeregex)\b([+]*)|)\s*	# $1,2: first option then plus
@@ -28,12 +28,14 @@ my %flag2info = (
   phys => [[qw(FlagPhys)]],
   real => [[qw(FlagTypeOverride FlagReal)]],
   complex => [[qw(FlagTypeOverride FlagComplex)]],
+  '!real' => [[qw(FlagTypeOverride FlagNotReal)]],
   '!complex' => [[qw(FlagTypeOverride FlagNotComplex)]],
   (map +($_->ppforcetype => [[qw(FlagTypeOverride FlagTyped)], 'Type']), types),
 );
 my %flag2c = qw(
   FlagReal PDL_PARAM_ISREAL
   FlagComplex PDL_PARAM_ISCOMPLEX
+  FlagNotReal PDL_PARAM_ISNOTREAL
   FlagNotComplex PDL_PARAM_ISNOTCOMPLEX
   FlagTyped PDL_PARAM_ISTYPED
   FlagTplus PDL_PARAM_ISTPLUS

t/nat_complex.t

@@ -21,6 +21,8 @@ ok !$x->type->real, 'complex type not real';
 ok double->real, 'real type is real';
 ok !$x->sumover->type->real, 'sumover type=complex';
 
+eval {double(5)->re};
+like $@, qr/must be complex/, 'error if given real data to re';
 $x = cdouble(2,3);
 $x-=i2C(3);
 is type($x), 'cdouble', 'type promotion ndarray - i';