From 34616d4c2f2ea3aa0579953346af4ad129ea1e81 Mon Sep 17 00:00:00 2001
From: Luis Mochan <mochan@fis.unam.mx>
Date: Mon, 13 Jan 2025 19:31:16 -0600
Subject: [PATCH 1/2] Add csqrt_up and tests

---
 lib/PDL/Ops.pd | 18 ++++++++++++++++++
 t/ops.t        | 28 ++++++++++++++++++++++++++++
 2 files changed, 46 insertions(+)

diff --git a/lib/PDL/Ops.pd b/lib/PDL/Ops.pd
index cf3eaffd1..9fe0ff9ce 100644
--- a/lib/PDL/Ops.pd
+++ b/lib/PDL/Ops.pd
@@ -422,6 +422,24 @@ EOF
   Code => '$c() = $r() + $i() * I;'
 );
 
+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/ops.t b/t/ops.t
index b22db2afc..34a9bb55f 100644
--- a/t/ops.t
+++ b/t/ops.t
@@ -112,6 +112,34 @@ if ($can_complex_power) {
 is $pa->at(0), '16', 'sqrt orig value ok';
 }
 
+{   # 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";
+    my $lsqrt=long($squares)->csqrt_up;
+    is_pdl($lsqrt, cdouble($roots), "csqrt_up of long");
+    my $llsqrt=longlong($squares)->csqrt_up;
+    is_pdl($llsqrt, cdouble($roots), "csqrt_up of longlong");
+    my $fsqrt=float($squares)->csqrt_up;
+    is_pdl($fsqrt, cfloat($roots), "csqrt_up of float");
+    my $dsqrt=double($squares)->csqrt_up;
+    is_pdl($dsqrt,cdouble($roots), "csqrt_up of double");
+    my $ldsqrt=ldouble($squares)->csqrt_up;
+    is_pdl($ldsqrt, cldouble($roots), "csqrt_up of ldouble");
+    my $cfsqrt=cfloat($squares)->csqrt_up;
+    is_pdl($cfsqrt, cfloat($roots), "csqrt_up of cfloat");
+    my $cdsqrt=cdouble($squares)->csqrt_up;
+    is_pdl($cdsqrt,cdouble($roots), "csqrt_up of cdouble");
+    my $cldsqrt=cldouble($squares)->csqrt_up;
+    is_pdl($cldsqrt,cldouble($roots), "csqrt_up of cldouble");
+}
+
 {
     is_pdl(r2C(long(1)), cdouble(1), "r2C of long");
     is_pdl(r2C(longlong(1)), cdouble(1), "r2C of longlong");

From 59b535f4ee68d5c69d9de239767f5a7df88dd754 Mon Sep 17 00:00:00 2001
From: Luis Mochan <mochan@fis.unam.mx>
Date: Wed, 15 Jan 2025 12:05:01 -0600
Subject: [PATCH 2/2] Added csqrt and tests

---
 lib/PDL/Ops.pd | 15 +++++++++++++++
 t/ops.t        | 28 ++++++++++++++++++++++++++++
 2 files changed, 43 insertions(+)

diff --git a/lib/PDL/Ops.pd b/lib/PDL/Ops.pd
index 9fe0ff9ce..3f101beba 100644
--- a/lib/PDL/Ops.pd
+++ b/lib/PDL/Ops.pd
@@ -422,6 +422,21 @@ 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()',
diff --git a/t/ops.t b/t/ops.t
index 34a9bb55f..c878b6603 100644
--- a/t/ops.t
+++ b/t/ops.t
@@ -112,6 +112,34 @@ if ($can_complex_power) {
 is $pa->at(0), '16', 'sqrt orig value ok';
 }
 
+{   # 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);
+    my $sqrt=csqrt($eiO);
+    is_pdl($sqrt, $eiO2, "csqrt of complex");
+    my $i=csqrt(-1);
+    is_pdl($i, i(), "csqrt of real -1");
+    my $squares="-9 -4 -1 0 1 4 9";
+    my $roots="3i 2i i 0 1 2 3";
+    my $lsqrt=long($squares)->csqrt;
+    is_pdl($lsqrt, cdouble($roots), "csqrt of long");
+    my $llsqrt=longlong($squares)->csqrt;
+    is_pdl($llsqrt, cdouble($roots), "csqrt of longlong");
+    my $fsqrt=float($squares)->csqrt;
+    is_pdl($fsqrt, cfloat($roots), "csqrt of float");
+    my $dsqrt=double($squares)->csqrt;
+    is_pdl($dsqrt,cdouble($roots), "csqrt of double");
+    my $ldsqrt=ldouble($squares)->csqrt;
+    is_pdl($ldsqrt, cldouble($roots), "csqrt of ldouble");
+    my $cfsqrt=cfloat($squares)->csqrt;
+    is_pdl($cfsqrt, cfloat($roots), "csqrt of cfloat");
+    my $cdsqrt=cdouble($squares)->csqrt;
+    is_pdl($cdsqrt,cdouble($roots), "csqrt of cdouble");
+    my $cldsqrt=cldouble($squares)->csqrt;
+    is_pdl($cldsqrt,cldouble($roots), "csqrt of cldouble");
+}
+
 {   # csqrt_up
     my $pi=4*atan2(1,1);
     my $eiO = exp(i()*sequence(8)*$pi/4);