view release on metacpan or  search on metacpan

complex.pd  view on Meta::CPAN

           $c(m=>1) = 0.5 * atan2 (2*ai, 1 - ar*ar - i2);
        ^
;

pp_def 'Cproj',
	Pars => 'a(m=2); [o]c(m=2)',
	Inplace => 1,
        GenericTypes => $F,
        Doc => 'compute the projection of a complex number to the riemann sphere. Works inplace',
        Code => q^
           $GENERIC() ar = $a(m=>0), ai = $a(m=>1);

           double den = ar*ar + ai*ai + 1;

           $c(m=>0) = 2*ar / den;
           $c(m=>1) = 2*ai / den;
        ^
;

pp_def 'Croots',
	Pars => 'a(m=2); [o]c(m=2,n)',
        OtherPars => 'int n => n',
        GenericTypes => $F,
        Doc => 'Compute the C<n> roots of C<a>. C<n> must be a positive integer. The result will always be a complex type!',
        PMCode => q^sub PDL::Complex::Croots($$) {
           my ($pdl, $n) = @_;
           my $r = PDL->null;
           &PDL::Complex::_Croots_int($pdl, $r, $n);
           bless $r;
        }^,
        Code => q^
           double s, c;
           double ar = $a(m=>0), ai = $a(m=>1),
                  n1 = 1 / (double)$COMP(n),
                  rr = pow (CABS (ar, ai), n1), /* do not optimize the sqrt out of this expr! */
                  at = atan2 (ai, ar) * n1,
                  ti = M_2PI * n1;

           loop(n) %{
               SINCOS (at, s, c);

               $c(m=>0) = rr * c;
               $c(m=>1) = rr * s;

               at += ti;
            %}
        ^
;

pp_addpm <<'EOD';

=head2 re, im

Return the real or imaginary part of the complex number(s) given.

These are slicing operators, so data flow works. The real and
imaginary parts are returned as ndarrays (ref eq PDL).

=cut

sub re($) { $_[0]->slice("(0)") }
sub im($) { $_[0]->slice("(1)") }

{
no warnings 'redefine';
# if the argument does anything other than pass through 0-th dim, re-bless
sub slice :lvalue {
  my $first = ref $_[1] ? $_[1][0] : (split ',', $_[1])[0];
  my $class = ($first//'') =~ /^[:x]?$/i ? ref($_[0]) : 'PDL';
  my $ret = bless $_[0]->SUPER::slice(@_[1..$#_]), $class;
  $ret;
}
}

EOD

pp_def 'rCpolynomial',
       Pars => 'coeffs(n); x(c=2,m); [o]out(c=2,m)',
       Doc => 'evaluate the polynomial with (real) coefficients C<coeffs> at the (complex) position(s) C<x>. C<coeffs[0]> is the constant term.',
       GenericTypes => $F,
       PMCode=> q!
sub rCpolynomial {
    my $coeffs = shift;
    my $x = shift;
    my $out = $x->copy;
    _rCpolynomial_int($coeffs,$x,$out);
    return PDL::complex($out);
    }
!,
       Code => q!
          loop(m) %{
             double xr = 1;
             double xi = 0;
             double or = 0;
             double oi = 0;
             double Xr;
             loop(n) %{
                or += $coeffs() * xr;
                oi += $coeffs() * xi;
                Xr = xr;
                xr = Xr * $x(c=>0) - xi * $x(c=>1);
                xi = xi * $x(c=>0) + Xr * $x(c=>1);
             %}
             $out(c=>0) = or;
             $out(c=>1) = oi;
          %}
       !
;

pp_def('Ctricpy',
  Pars => 'A(c=2,m,n);[o] C(c=2,m,n)',
  OtherPars => 'int uplo',
  OtherParsDefaults => {uplo => 0},
  ArgOrder => [qw(A uplo C)],
  Code => '
if ($COMP(uplo)) { /* lower */
  broadcastloop %{ loop(n,m=:n+1,c) %{ $C() = $A(); %} %}
} else {
  broadcastloop %{ loop(n,m=n,c) %{ $C() = $A(); %} %}
}
  ',
  Doc => <<'EOT',