File Coverage

lib/PDL/Ops.pd

Criterion	Covered	Total	%
statement	16	16	100.0
branch	17	64	26.5
condition			n/a
subroutine			n/a
pod			n/a
total	33	80	41.2

use strict;
use warnings;
use PDL::Types qw(types ppdefs ppdefs_all ppdefs_complex);
require PDL::Core::Dev;

my $A = [ppdefs_all];
my $C = [ppdefs_complex];
my $F = [map $_->ppsym, grep $_->real && !$_->integer, types];
$F = [(grep $_ ne 'D', @$F), 'D']; # so defaults to D if non-float given
my $AF = [map $_->ppsym, grep !$_->integer, types];
$AF = [(grep $_ ne 'D', @$AF), 'D']; # so defaults to D if non-float given
my $T = [map $_->ppsym, grep $_->integer, types];
my $U = [map $_->ppsym, grep $_->unsigned, types];
my $S = [map $_->ppsym, grep $_->real && !$_->unsigned, types];
my %is_real; @is_real{ppdefs()} = ();
my @Rtypes = grep $_->real, types();
my @Ctypes = grep !$_->real, types();
my @Ftypes = grep !$_->integer, types();

pp_addpm({At=>'Top'},<<'EOD');

use strict;
use warnings;

=head1 NAME

PDL::Ops - Fundamental mathematical operators

=head1 DESCRIPTION

This module provides the functions used by PDL to
overload the basic mathematical operators (C<+ - / *>
etc.) and functions (C etc.)

It also includes the function C, which should
be a perl function so that we can overload it!

Matrix multiplication (the operator C) is handled
by the module L.

=head1 SYNOPSIS

none

=cut

EOD

pp_addpm({At=>'Bot'},<<'EOPM');

=head1 AUTHOR

Tuomas J. Lukka ([e-mail]),
Karl Glazebrook ([e-mail]),
Doug Hunt ([e-mail]),
Christian Soeller ([e-mail]),
Doug Burke ([e-mail]),
and Craig DeForest ([e-mail]).

=cut

EOPM

pp_addhdr('
#include

#define MOD(X,N) (((N) == 0)   ?    0   :   (   (X) - (PDL_ABS(N))  *  ((long long)((X)/(PDL_ABS(N))) + (   ( ((N) * ((long long)((X)/(N)))) != (X) )   ?   ( ( ((N)<0) ? 1 : 0 )  +  ( (((X)<0) ? -1 : 0)))  :  0 ))))
#define BU_MOD(X,N)(((N) == 0)   ?    0   :   ( (X)-(N)*((uint64_t)((X)/(N))) ))
#define SPACE(A,B)   ( ((A)<(B)) ? -1 : ((A)!=(B)) )
');

my %char2escape = ('>'=>'E','<'=>'E');
my $chars = '(['.join('', map quotemeta, sort keys %char2escape).'])';
sub protect_chars {
  my ($txt) = @_;
  $txt =~ s/$chars/$char2escape{$1}/g;
  return $txt;
}

# simple binary operators

pp_addhdr(pp_line_numbers(__LINE__, <<'EOF'));
#define PDL_BADVAL_WARN_X(datatype, ctype, ppsym, ...) \
  bad_anyval.type = datatype; bad_anyval.value.ppsym = PDL->bvals.ppsym;
#define PDL_BADVAL_WARN(var) \
  { \
    PDL_Anyval bad_anyval = { PDL_INVALID, {0} }; \
    if (!(var->has_badvalue && var->badvalue.type != var->datatype)) { \
      if (var->has_badvalue) \
        bad_anyval = var->badvalue; \
      else { \
        PDL_GENERICSWITCH(PDL_TYPELIST_ALL, var->datatype, PDL_BADVAL_WARN_X, ) \
      } \
    } \
    if (bad_anyval.type < 0) \
      barf("Error getting badvalue, type=%d", bad_anyval.type); \
    complex double bad_c; \
    ANYVAL_TO_CTYPE(bad_c, complex double, bad_anyval); \
    if( bad_c == 0 || bad_c == 1 ) \
      warn(#var " badvalue is set to 0 or 1. This will cause data loss when using badvalues for comparison operators."); \
  }
EOF
sub biop {
    my ($name,$op,$mutator,$doc,%extra) = @_;
    my $optxt = protect_chars ref $op eq 'ARRAY' ? $op->[1] : $op;
    $op = $op->[0] if ref $op eq 'ARRAY';
    $extra{HdrCode} = << 'EOH';
  if (swap) {
    pdl *tmp = a;
    a = b;
    b = tmp;
  }
EOH
    # handle exceptions
    if ( exists $extra{Exception} ) {
# NOTE This option is unused.
	#      See also `ufunc()`.
	delete $extra{Exception};
    }
    if ($extra{Comparison}) {
        my $first_complex = $Ctypes[0]->sym;
	$extra{HdrCode} .= < 1;
  if ((a->datatype >= $first_complex) || (b->datatype >= $first_complex))
    barf("Can't compare complex numbers");
EOF
	$extra{HdrCode} .= "  PDL_BADVAL_WARN(a)\n  PDL_BADVAL_WARN(b)\n";
        delete $extra{Comparison};
    }

    my $bitwise = delete $extra{Bitwise};
    pp_def($name,
	   Pars => 'a(); b(); [o]c();',
	   OtherPars => 'int $swap',
	   OtherParsDefaults => { swap => 0 },
	   HandleBad => 1,
	   NoBadifNaN => 1,
	   Inplace => [ 'a' ],
	   Overload => [$op, $mutator, $bitwise],
	   NoExport => 1,
	   Code => pp_line_numbers(__LINE__, <
PDL_IF_BAD(char anybad = 0;,)
broadcastloop %{
  PDL_IF_BAD(if ( ( \$PDLSTATEISBAD(a) && \$ISBAD(a()) )
               || ( \$PDLSTATEISBAD(b) && \$ISBAD(b()) )) { \$SETBAD(c()); anybad = 1; } else,)
     \$c() = \$a() $op \$b();
%}
PDL_IF_BAD(if (anybad) \$PDLSTATESETBAD(c);,)
EOF
	   %extra,
	   Doc => $doc,
    );
}

#simple binary functions
sub bifunc {
    my ($name,$func,$mutator,$doc,%extra) = @_;
    my $funcov = ref $func eq 'ARRAY' ? $func->[1] : $func;
    my $isop = $funcov =~ s/^op//;
    my $funcovp = protect_chars $funcov;
    $func = $func->[0] if ref $func eq 'ARRAY';
    my $got_complex = PDL::Core::Dev::got_complex_version($func, 2);
    $extra{GenericTypes} = [ grep exists $is_real{$_}, @{$extra{GenericTypes}} ]
	if !$got_complex and $extra{GenericTypes};
    $extra{HdrCode} .= << 'EOH';
  if (swap) {
    pdl *tmp = a;
    a = b;
    b = tmp;
  }
EOH
    # is this one to be used as a function or operator ?

    my $codestr;
    if ($extra{unsigned}){
#a little dance to avoid the MOD macro warnings for byte & ushort datatypes
      my $t = join '', map $_->ppsym, grep $_->real, types();
      my $v = join ',', map
        $_->unsigned ? 'BU_' : '',
        grep $_->real, types();
      $codestr = << "ENDCODE";
\$c() = (\$GENERIC(c))\$T$t($v)$func(\$a(),\$b());
ENDCODE
#end dance
    } else {
      $codestr = '$c() = ($GENERIC(c))'.$func.'($a(),$b());';
    }
    delete $extra{unsigned}; #remove the key so it doesn't get added in pp_def.

    pp_def($name,
	   HandleBad => 1,
	   NoBadifNaN => 1,
	   Pars => 'a(); b(); [o]c();',
	   OtherPars => 'int $swap',
	   OtherParsDefaults => { swap => 0 },
	   Inplace => [ 'a' ],
	   Overload => [$funcov, $mutator],
	   NoExport => 1,
	   Code => pp_line_numbers(__LINE__, <
PDL_IF_BAD(char anybad = 0;,)
broadcastloop %{
  PDL_IF_BAD(if ( \$ISBAD(a()) || \$ISBAD(b()) ) { anybad = 1; \$SETBAD(c()); } else ,) {
     $codestr
  }
%}
PDL_IF_BAD(if (anybad) \$PDLSTATESETBAD(c);,)
EOF
	   %extra,
	   Doc => $doc,
    );
}

# simple unary functions and operators
sub ufunc {
    my ($name,$func,$overload,$doc,%extra) = @_;
    my $funcov = ref $func eq 'ARRAY' ? $func->[1] : $func;
    my $funcovp = protect_chars $funcov;
    $func = $func->[0] if ref $func eq 'ARRAY';
    my $got_complex = PDL::Core::Dev::got_complex_version($func, 1);
    $extra{GenericTypes} = [ grep exists $is_real{$_}, @{$extra{GenericTypes}} ]
	if !$got_complex and $extra{GenericTypes};

    # handle exceptions
    if ( exists $extra{Exception} ) {
#	print "Warning: ignored exception for $name\n";
# NOTE This option is unused.
	#      See also `biop()`.
	delete $extra{Exception};
    }
    my $codestr = '$b() = ($GENERIC(b))'.$func.'($a());';
    if (delete $extra{NoTgmath} and $got_complex) {
        # don't bother if not got complex version
        $codestr = join "\n",
            'types('.join('', map $_->ppsym, @Rtypes).') %{'.$codestr.'%}',
            (map 'types('.$_->ppsym.') %{$b() = c'.$func.$_->floatsuffix.'($a());%}', @Ctypes),
            ;
    }
    # do not have to worry about propagation of the badflag when
    # inplace since only input ndarray is a, hence its badflag
    # won't change
    # UNLESS an exception occurs...
    pp_def($name,
	   Pars => 'a(); [o]b()',
	   HandleBad => 1,
	   NoBadifNaN => 1,
	   Inplace => 1,
	   !$overload ? () : (Overload => $funcov),
	   NoExport => 1,
	   Code => pp_line_numbers(__LINE__, <
PDL_IF_BAD(if ( \$ISBAD(a()) ) \$SETBAD(b()); else {,)
  $codestr
PDL_IF_BAD(},)
EOF
	   %extra,
	   Doc => $doc,
    );
}

######################################################################

# we trap some illegal operations here -- see the Exception option
# note, for the ufunc()'s, the checks do not work too well
#    for unsigned integer types (ie < 0)
#
# XXX needs thinking about
#    - have to integrate into Code section as well (so
#      12/pdl(2,4,0,3) is trapped and flagged bad)
#      --> complicated
#    - perhaps could use type %{ %} ?
#
# ==> currently have commented out the exception code, since
#     want to see if can use NaN/Inf for bad values
#     (would solve many problems for F,D types)
#
# there is an issue over how we handle comparison operators
# - see Primitive/primitive.pd/zcover() for more discussion
#

## arithmetic ops
biop('plus','+',1,'add two ndarrays',GenericTypes => $A);
biop('mult','*',1,'multiply two ndarrays',GenericTypes => $A);
biop('minus','-',1,'subtract two ndarrays',GenericTypes => $A);
biop('divide','/',1,'divide two ndarrays', Exception => '$b() == 0', GenericTypes => $A);

## note: divide should perhaps trap division by zero as well

## comparison ops
# not defined for complex numbers
biop('gt','>',0,'the binary E (greater than) operation', Comparison => 2);
biop('lt','<',0,'the binary E (less than) operation', Comparison => 2);
biop('le','<=',0,'the binary E= (less equal) operation', Comparison => 2);
biop('ge','>=',0,'the binary E= (greater equal) operation', Comparison => 2);
biop('eq','==',0,'binary I operation (C<==>)', Comparison => 1, GenericTypes => $A);
biop('ne','!=',0,'binary I operation (C)', Comparison => 1, GenericTypes => $A);

## bit ops
# those need to be limited to the right types
biop('shiftleft','<<',1,'bitwise leftshift C<$a> by C<$b>',GenericTypes => $T);
biop('shiftright','>>',1,'bitwise rightshift C<$a> by C<$b>',GenericTypes => $T);
biop('or2','|',1,'bitwise I of two ndarrays',GenericTypes => $T,
      Bitwise => 1);
biop('and2','&',1,'bitwise I of two ndarrays',GenericTypes => $T,
      Bitwise => 1);
biop('xor','^',1,'bitwise I of two ndarrays',GenericTypes => $T,
      Bitwise => 1);

pp_addpm(
"=head2 xor2\n\n=for ref\n\nSynonym for L.\n\n=cut\n
*PDL::xor2 = *xor2 = \\&PDL::xor;"
    );

# some standard binary functions
bifunc('power',['pow','op**'],1,'raise ndarray C<$a> to the power C<$b>',GenericTypes => [@$C, @$F]);
bifunc('atan2','atan2',0,'elementwise C of two ndarrays',GenericTypes => $F);
bifunc('modulo',['MOD','op%'],1,'elementwise C operation',unsigned=>1);
bifunc('spaceship',['SPACE','op<=>'],0,'elementwise "<=>" operation');

# some standard unary functions
ufunc('bitnot','~',1,'unary bitwise negation',GenericTypes => $T);
ufunc('sqrt','sqrt',1,'elementwise square root', GenericTypes => $A); # Exception => '$a() < 0');
ufunc('sin','sin',1,'the sin function', GenericTypes => $A);
ufunc('cos','cos',1,'the cos function', GenericTypes => $A);
ufunc('not','!',1,'the elementwise I operation');
ufunc('exp','exp',1,'the exponential function',GenericTypes => [@$C, @$F]);
ufunc('log','log',1,'the natural logarithm',GenericTypes => [@$C, @$F], Exception => '$a() <= 0');

# no export these because clash with Test::Deep (re) or internal (_*abs)
cfunc('re', 'creal', 1, 0, 'Returns the real part of a complex number.',
  '$complexv() = $b() + I * cimag($complexv());'
);
cfunc('im', 'cimag', 1, 0, 'Returns the imaginary part of a complex number.',
  '$complexv() = creal($complexv()) + I * $b();'
);
cfunc('_cabs', 'fabs', 1, 0, 'Returns the absolute (length) of a complex number.', undef,
    PMFunc=>'',
);
my $rabs_code = '
  types('.join('', @$U).') %{ $b()=$a(); %}
  types('.join('', @$S).') %{ $b()=PDL_ABS($a()); %}
';
pp_def ( '_rabs',
	Pars=>'a(); [o]b()',
	HandleBad => 1,
	NoBadifNaN => 1,
	    Inplace => 1,
	NoExport => 1,
	Code => pp_line_numbers(__LINE__-1, qq{
PDL_IF_BAD(if ( \$ISBAD(a()) ) \$SETBAD(b()); else,)
  $rabs_code
	}),
	Doc=>undef,
	PMFunc=>'',
);

# make log10() work on scalars (returning scalars)
# as well as ndarrays
ufunc('log10','log10',0,'the base 10 logarithm', GenericTypes => $A,
      Exception => '$a() <= 0',
      NoTgmath => 1, # glibc for at least GCC 8.3.0 won't tgmath log10 though 7.1.0 did
      NoExport => 0,
      PMCode => <<'EOF',
sub PDL::log10 {
    my ($x, $y) = @_;
    return log($x) / log(10) if !UNIVERSAL::isa($x,"PDL");
    barf "inplace but output given" if $x->is_inplace and defined $y;
    if ($x->is_inplace) { $x->set_inplace(0); $y = $x; }
    elsif (!defined $y) { $y = $x->initialize; }
    &PDL::_log10_int( $x, $y );
    $y;
};
EOF
);

pp_def(
       'assgn',
       HandleBad => 1,
       GenericTypes => $A,
       Pars => 'a(); [o]b();',
       Code => pp_line_numbers(__LINE__-1, q{
PDL_IF_BAD(char anybad = 0;,)
broadcastloop %{
  PDL_IF_BAD(if ($ISBAD(a())) { anybad = 1; $SETBAD(b()); continue; },)
  $b() = $a();
%}
PDL_IF_BAD(if (anybad) $PDLSTATESETBAD(b);,)
       }),
       Doc =>
'Plain numerical assignment. This is used to implement the ".=" operator',
);

# special functions for complex data types that don't work well with
# the ufunc/bifunc logic
sub cfunc {
    my ($name, $func, $make_real, $force_complex, $doc, $backcode, %extra) = @_;
    my $codestr = pp_line_numbers(__LINE__-1,"\$b() = $func(\$complexv());");
    pp_def($name,
	   GenericTypes=>$C,
	   Pars => ($force_complex ? '!real ' : '').'complexv(); '.($make_real ? 'real' : '').' [o]b()',
	   HandleBad => 1,
	   NoBadifNaN => 1,
	   (($make_real || $force_complex) ? () : (Inplace => 1)),
	   NoExport => 1,
	   Code => pp_line_numbers(__LINE__-1, qq{
PDL_IF_BAD(if ( \$ISBAD(complexv()) ) \$SETBAD(b()); else,)
  $codestr
	   }),
	   !$backcode ? () : (
	     DefaultFlow => 1,
	     TwoWay => 1,
	     BackCode => pp_line_numbers(__LINE__-1, qq{
		PDL_IF_BAD(if ( \$ISBAD(b()) ) \$SETBAD(complexv()); else {,)
		   $backcode
		PDL_IF_BAD(},)
	     }),
	   ),
	   %extra,
	   Doc => $doc . (!$backcode ? '' : ' Flows data back & forth.'),
    );
}

cfunc('carg', 'carg', 1, 1, 'Returns the polar angle of a complex number.', undef, NoExport => 0);
cfunc('conj', 'conj', 0, 0, 'complex conjugate.', undef, NoExport => 0);

pp_def('czip',
  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 ndarray to "i" times the C ndarray. Only
takes real ndarrays as input.
EOF
  Code => '$c() = $r() + $i() * I;'
);

pp_def('ipow',
   Inplace => [qw(a ans)],
   Doc => qq{
=for ref

raise ndarray C<\$a> to integer power C<\$b>

Algorithm from L
},
   Pars => 'a(); longlong b(); [o] ans()',
   GenericTypes => [qw(P Q), @$AF],
   Code => pp_line_numbers(__LINE__-1, q{
$GENERIC(b) n = $b();
if (n == 0) {
  $ans() = 1;
  continue;
}
$GENERIC() y = 1;
$GENERIC() x = $a();
if (n < 0) {
  x = 1 / x;
  n = -n;
}
while (n > 1) {
  if (n % 2) {
    y *= x;
    n -= 1;
  }
  x *= x;
  n /= 2;
}
$ans() = x * y;
})
);

pp_addpm(<<'EOPM');

=head2 abs

=for ref

Returns the absolute value of a number.

=cut

sub PDL::abs { $_[0]->type->real ? goto &PDL::_rabs : goto &PDL::_cabs }
EOPM

pp_def('abs2',
  GenericTypes=>$A,
  HandleBad => 1,
  Pars => 'a(); real [o]b()',
  Doc => 'Returns the square of the absolute value of a number.',
  Code => <<'EOF',
PDL_IF_BAD(if ($ISBAD(a())) { $SETBAD(b()); continue; },)
$b() = PDL_IF_GENTYPE_REAL(
  $a()*$a(),
  creall($a())*creall($a()) + cimagl($a())*cimagl($a())
);
EOF
);

pp_def('r2C',
  GenericTypes=>$AF,
  Pars => 'r(); complex [o]c()',
  Doc => 'convert real to native complex, with an imaginary part of zero',
  PMCode => << 'EOF',
sub PDL::r2C ($) {
  return $_[0] if UNIVERSAL::isa($_[0], 'PDL') and !$_[0]->type->real;
  my $r = $_[1] // PDL->nullcreate($_[0]);
  PDL::_r2C_int($_[0], $r);
  $r;
}
EOF
  Code => '$c() = $r();'
);

pp_def('i2C',
  GenericTypes=>$AF,
  Pars => 'i(); complex [o]c()',
  Doc => 'convert imaginary to native complex, with a real part of zero',
  PMCode => << 'EOF',
sub PDL::i2C ($) {
  return $_[0] if UNIVERSAL::isa($_[0], 'PDL') and !$_[0]->type->real;
  my $r = $_[1] // PDL->nullcreate($_[0]);
  PDL::_i2C_int($_[0], $r);
  $r;
}
EOF
  Code => '$c() = $i() * I;'
);

pp_addpm(<<'EOF');
# This is to used warn if an operand is non-numeric or non-PDL.
sub warn_non_numeric_op_wrapper {
  require Scalar::Util;
  my ($cb, $op_name) = @_;
  return sub {
    my ($op1, $op2) = @_;
    warn "'$op2' is not numeric nor a PDL in operator $op_name"
      unless Scalar::Util::looks_like_number($op2)
            || ( Scalar::Util::blessed($op2) && $op2->isa('PDL') );
    $cb->(@_);
  }
}

{ package # hide from MetaCPAN
    PDL;
  use overload
    "eq"    => PDL::Ops::warn_non_numeric_op_wrapper(\&PDL::eq, 'eq'),
    ".="    => sub {
      my @args = !$_[2] ? @_[1,0] : @_[0,1];
      PDL::Ops::assgn(@args);
      return $args[1];
    },
    'abs' => sub { PDL::abs($_[0]) },
    '++' => sub { $_[0] += ($PDL::Core::pdl_ones[$_[0]->get_datatype]//barf "Couldn't find 'one' for type ", $_[0]->get_datatype) },
    '--' => sub { $_[0] -= ($PDL::Core::pdl_ones[$_[0]->get_datatype]//barf "Couldn't find 'one' for type ", $_[0]->get_datatype) },
    ;
}
EOF

pp_done();