Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add upper_sqrt #524

Merged
merged 2 commits into from
Jan 16, 2025
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 33 additions & 0 deletions lib/PDL/Ops.pd
Original file line number Diff line number Diff line change
Expand Up @@ -422,6 +422,39 @@ 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);
mohawk2 marked this conversation as resolved.
Show resolved Hide resolved
$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{
Expand Down
56 changes: 56 additions & 0 deletions t/ops.t
Original file line number Diff line number Diff line change
Expand Up @@ -112,6 +112,62 @@ 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);
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");
Expand Down
Loading