move csqrt* to Math - #527

PDLPorters · Jan 22, 2025 · 4719490 · 4719490
1 parent 81a6883
commit 4719490
Show file tree

Hide file tree

Showing 5 changed files with 76 additions and 74 deletions.
diff --git a/Changes b/Changes
@@ -21,7 +21,7 @@
 - fix IO::Pic problem with older IO::GD installed (#522) - thanks @shawnlaffan
 - r2C and i2C now handle long double both real and complex (#524) - thanks @wlmb for report
 - add "update_data_from" method as streamlined get_dataref/upd_data
-- add Ops::csqrt{,_up} (#524) - thanks @wlmb
+- add Math::csqrt{,_up} (#524) - thanks @wlmb
 - Core::upd_data now sets DATACHANGED so will work for flowing ops
 - remove CopyBadStatusCode key from PP (#517)
 - badflag now not propagated back to inputs if set on outputs (#517)

diff --git a/lib/PDL/Math.pd b/lib/PDL/Math.pd
@@ -38,6 +38,8 @@ my $F = [map $_->ppsym, grep $_->real && !$_->integer, types()];
 my $C = [ppdefs_complex()];
 my @Rtypes = grep $_->real, types();
 my @Ctypes = grep !$_->real, types();
+my $AF = [map $_->ppsym, grep !$_->integer, types];
+$AF = [(grep $_ ne 'D', @$AF), 'D']; # so defaults to D if non-float given
 
 pp_addpm({At=>'Top'},<<'EOD');
 use strict;
@@ -487,6 +489,39 @@ works with native-complex data.
  ($yr, $yi) = polyval($cr, $ci, $xr, $xi);
 ',);
 
+pp_def('csqrt',
+    GenericTypes => $AF,
+    Pars => 'i(); complex [o] o()',
+    Doc => <<'EOF',
+Take the complex square root of a number choosing that with
+non-negative real part, i.e., a square root with a branch cut
+'infinitesimally' below the negative real axis, the standard choice.
+EOF
+    Code => q{
+        $TFDEGCH(PDL_CFloat,PDL_CDouble,PDL_CLDouble,PDL_CFloat,PDL_CDouble,PDL_CLDouble) tmp=$i();
+        tmp=csqrt(tmp);
+        $o() = tmp;
+    },
+);
+
+pp_def('csqrt_up',
+    GenericTypes => $AF,
+    Pars => 'i(); complex [o] o()',
+    Doc => <<'EOF',
+Take the complex square root of a number choosing that whose imaginary
+part is not negative, i.e., it is a square root with a branch cut
+'infinitesimally' below the positive real axis.
+EOF
+    Code => q{
+        $TFDEGCH(PDL_CFloat,PDL_CDouble,PDL_CLDouble,PDL_CFloat,PDL_CDouble,PDL_CLDouble) tmp=$i();
+        tmp=csqrt(tmp);
+        if(cimag(tmp)<0){
+             tmp = -tmp;
+        }
+        $o() = tmp;
+    },
+);
+
 pp_addpm({At=>'Bot'},<<'EOD');
 
 =head1 BUGS

diff --git a/lib/PDL/Ops.pd b/lib/PDL/Ops.pd
@@ -422,39 +422,6 @@ EOF
   Code => '$c() = $r() + $i() * I;'
 );
 
-pp_def('csqrt',
-    GenericTypes => $AF,
-    Pars => 'i(); complex [o] o()',
-    Doc => <<'EOF',
-Take the complex square root of a number choosing that with
-non-negative real part, i.e., a square root with a branch cut
-'infinitesimally' below the negative real axis, the standard choice.
-EOF
-    Code => q{
-        $TFDEGCH(PDL_CFloat,PDL_CDouble,PDL_CLDouble,PDL_CFloat,PDL_CDouble,PDL_CLDouble) tmp=$i();
-        tmp=csqrt(tmp);
-        $o() = tmp;
-    },
-);
-
-pp_def('csqrt_up',
-    GenericTypes => $AF,
-    Pars => 'i(); complex [o] o()',
-    Doc => <<'EOF',
-Take the complex square root of a number choosing that whose imaginary
-part is not negative, i.e., it is a square root with a branch cut
-'infinitesimally' below the positive real axis.
-EOF
-    Code => q{
-        $TFDEGCH(PDL_CFloat,PDL_CDouble,PDL_CLDouble,PDL_CFloat,PDL_CDouble,PDL_CLDouble) tmp=$i();
-        tmp=csqrt(tmp);
-        if(cimag(tmp)<0){
-             tmp = -tmp;
-        }
-        $o() = tmp;
-    },
-);
-
 pp_def('ipow',
    Inplace => [qw(a ans)],
    Doc => qq{

diff --git a/t/math.t b/t/math.t
@@ -35,6 +35,46 @@ is_pdl erf(0.5), pdl(1.-erfc(0.5)), "erf and erfc";
 is_pdl erf(erfi(0.5)), pdl(0.5), "erfi (both ways)";
 is_pdl erfi(erf(0.5)), pdl(0.5), "erfi (both ways)";
 
+{   # csqrt
+  my $pi=4*atan2(1,1);
+  my $eiO = exp(i()*(sequence(8)-3)*$pi/4);
+  my $eiO2 = exp(i()*(sequence(8)-3)*$pi/8);
+  is_pdl csqrt($eiO), $eiO2, "csqrt of complex";
+  is_pdl csqrt(-1), i(), "csqrt of real -1";
+  my $squares="-9 -4 -1 0 1 4 9";
+  my $roots="3i 2i i 0 1 2 3";
+  is_pdl long($squares)->csqrt,     cdouble($roots), "csqrt of long";
+  is_pdl longlong($squares)->csqrt, cdouble($roots), "csqrt of longlong";
+  is_pdl float($squares)->csqrt,    cfloat($roots), "csqrt of float";
+  is_pdl double($squares)->csqrt,   cdouble($roots), "csqrt of double";
+  is_pdl ldouble($squares)->csqrt,  cldouble($roots), "csqrt of ldouble";
+  is_pdl cfloat($squares)->csqrt,   cfloat($roots), "csqrt of cfloat";
+  is_pdl cdouble($squares)->csqrt,  cdouble($roots), "csqrt of cdouble";
+  is_pdl cldouble($squares)->csqrt, cldouble($roots), "csqrt of cldouble";
+  is_pdl pdl('-2i')->csqrt, pdl('1-i');
+}
+
+{   # csqrt_up
+  my $pi=4*atan2(1,1);
+  my $eiO = exp(i()*sequence(8)*$pi/4);
+  my $eiO2 = exp(i()*sequence(8)*$pi/8);
+  my $sqrt=csqrt_up($eiO);
+  is_pdl($sqrt, $eiO2, "Square of csqrt_up of complex");
+  my $i=csqrt_up(-1);
+  is_pdl($i, i(), "csqrt_up of real -1");
+  my $squares="-9 -4 -1 0 1 4 9";
+  my $roots="3i 2i i 0 1 2 3";
+  is_pdl long($squares)->csqrt_up,     cdouble($roots), "csqrt_up of long";
+  is_pdl longlong($squares)->csqrt_up, cdouble($roots), "csqrt_up of longlong";
+  is_pdl float($squares)->csqrt_up,    cfloat($roots), "csqrt_up of float";
+  is_pdl double($squares)->csqrt_up,   cdouble($roots), "csqrt_up of double";
+  is_pdl ldouble($squares)->csqrt_up,  cldouble($roots), "csqrt_up of ldouble";
+  is_pdl cfloat($squares)->csqrt_up,   cfloat($roots), "csqrt_up of cfloat";
+  is_pdl cdouble($squares)->csqrt_up,  cdouble($roots), "csqrt_up of cdouble";
+  is_pdl cldouble($squares)->csqrt_up, cldouble($roots), "csqrt_up of cldouble";
+  is_pdl pdl('-2i')->csqrt_up, pdl('-1+i');
+}
+
 {
 my $pa = pdl(0.0,30.0);
 $pa->inplace->erf;

diff --git a/t/ops.t b/t/ops.t
@@ -88,46 +88,6 @@ if ($can_complex_power) {
 }
 }
 
-{   # csqrt
-  my $pi=4*atan2(1,1);
-  my $eiO = exp(i()*(sequence(8)-3)*$pi/4);
-  my $eiO2 = exp(i()*(sequence(8)-3)*$pi/8);
-  is_pdl csqrt($eiO), $eiO2, "csqrt of complex";
-  is_pdl csqrt(-1), i(), "csqrt of real -1";
-  my $squares="-9 -4 -1 0 1 4 9";
-  my $roots="3i 2i i 0 1 2 3";
-  is_pdl long($squares)->csqrt,     cdouble($roots), "csqrt of long";
-  is_pdl longlong($squares)->csqrt, cdouble($roots), "csqrt of longlong";
-  is_pdl float($squares)->csqrt,    cfloat($roots), "csqrt of float";
-  is_pdl double($squares)->csqrt,   cdouble($roots), "csqrt of double";
-  is_pdl ldouble($squares)->csqrt,  cldouble($roots), "csqrt of ldouble";
-  is_pdl cfloat($squares)->csqrt,   cfloat($roots), "csqrt of cfloat";
-  is_pdl cdouble($squares)->csqrt,  cdouble($roots), "csqrt of cdouble";
-  is_pdl cldouble($squares)->csqrt, cldouble($roots), "csqrt of cldouble";
-  is_pdl pdl('-2i')->csqrt, pdl('1-i');
-}
-
-{   # csqrt_up
-  my $pi=4*atan2(1,1);
-  my $eiO = exp(i()*sequence(8)*$pi/4);
-  my $eiO2 = exp(i()*sequence(8)*$pi/8);
-  my $sqrt=csqrt_up($eiO);
-  is_pdl($sqrt, $eiO2, "Square of csqrt_up of complex");
-  my $i=csqrt_up(-1);
-  is_pdl($i, i(), "csqrt_up of real -1");
-  my $squares="-9 -4 -1 0 1 4 9";
-  my $roots="3i 2i i 0 1 2 3";
-  is_pdl long($squares)->csqrt_up,     cdouble($roots), "csqrt_up of long";
-  is_pdl longlong($squares)->csqrt_up, cdouble($roots), "csqrt_up of longlong";
-  is_pdl float($squares)->csqrt_up,    cfloat($roots), "csqrt_up of float";
-  is_pdl double($squares)->csqrt_up,   cdouble($roots), "csqrt_up of double";
-  is_pdl ldouble($squares)->csqrt_up,  cldouble($roots), "csqrt_up of ldouble";
-  is_pdl cfloat($squares)->csqrt_up,   cfloat($roots), "csqrt_up of cfloat";
-  is_pdl cdouble($squares)->csqrt_up,  cdouble($roots), "csqrt_up of cdouble";
-  is_pdl cldouble($squares)->csqrt_up, cldouble($roots), "csqrt_up of cldouble";
-  is_pdl pdl('-2i')->csqrt_up, pdl('-1+i');
-}
-
 {
   is_pdl(r2C(long(1)), cdouble(1), "r2C of long");
   is_pdl(r2C(longlong(1)), cdouble(1), "r2C of longlong");