Christian Soeller > PDL-2.4.1 > PDL::Math

Download:
PDL-2.4.1.tar.gz

Annotate this POD

CPAN RT

Open  0
View/Report Bugs
Source   Latest Release: PDL-2.007_04

NAME ^

PDL::Math - extended mathematical operations and special functions

SYNOPSIS ^

 use PDL::Math;

 use PDL::Graphics::TriD;
 imag3d [SURF2D,bessj0(rvals(zeroes(50,50))/2)];

DESCRIPTION ^

This module extends PDL with more advanced mathematical functions than provided by standard Perl.

All the functions have one input pdl, and one output, unless otherwise stated.

Many of the functions are linked from the system maths library or the Cephes maths library (determined when PDL is compiled); a few are implemented entirely in PDL.

log gamma function

This returns 2 piddles -- the first set gives the log(gamma) values, while the second set, of integer values, gives the sign of the gamma function. This is useful for determining factorials, amongst other things.

EOD

} # if: $^O !~ win32

pp_addhdr(' #include <math.h> #include "protos.h" /* Change names when fixing glibc-2.1 bug */ #ifdef MY_FIXY0 #define y0(a) fixy0(a) extern double fixy0(double a); #endif #ifdef MY_FIXYN #define yn(a,b) fixyn(a,b) extern double fixyn(int a, double b); #endif ');

## handle various cases of 'finite' # if ($^O =~ /MSWin/) { # _finite in VC++ 4.0 pp_addhdr(' #define finite _finite #include <float.h> '); }

# patch from Albert Chin if ($^O =~ /hpux/) { pp_addhdr(' #ifdef isfinite #define finite isfinite #endif '); }

# Standard `-lm' my (@ufuncs1) = qw(acos asin atan cosh sinh tan tanh); # F,D only my (@ufuncs1g) = qw(ceil floor rint); # Any type

# Note: # ops.pd has a power() function that does the same thing # (although it has OtherPars => 'int swap;' as well) # - left this in for now. # my (@bifuncs1) = qw(pow); # Any type

# Extended `-lm' my (@ufuncs2) = qw(acosh asinh atanh erf erfc); # F,D only my (@besufuncs) = qw(j0 j1 y0 y1); # " my (@besbifuncs) = qw(jn yn); # " # Need igamma, ibeta, and a fall-back implementation of the above

sub code_ufunc { return '$b() = ' . $_[0] . '($a());'; } sub badcode_ufunc { my $name = $_[0]; return 'if ( $ISBAD(a()) ) { $SETBAD(b()); } else { $b() = ' . $name . '($a()); }'; }

sub code_bifunc { my $name = $_[0]; my $a = $_[1] || 'a'; my $b = $_[2] || 'b'; my $c = $_[3] || 'c'; return "\$$c() = $name(\$$a(),\$$b());"; } sub badcode_bifunc { my $name = $_[0]; my $a = $_[1] || 'a'; my $b = $_[2] || 'b'; my $c = $_[3] || 'c'; return 'if ( $ISBAD('.$a.'()) || $ISBAD('.$b.'()) ) { $SETBAD('.$c.'()); } else { ' . "\$$c() = $name(\$$a(),\$$b()); }"; }

sub inplace_doc { my $func = shift; return "$doco{$func} Works inplace."; }

my $func; foreach $func (@ufuncs1) { pp_def($func, HandleBad => 1, NoBadifNaN => 1, GenericTypes => ['F','D'], Pars => 'a(); [o]b();', Inplace => 1, Doc => inplace_doc( $func ), Code => code_ufunc($func), BadCode => badcode_ufunc($func), ); }

foreach $func (@ufuncs1g) { pp_def($func, HandleBad => 1, NoBadifNaN => 1, Pars => 'a(); [o]b();', Inplace => 1, Doc => inplace_doc( $func ), Code => code_ufunc($func), BadCode => badcode_ufunc($func), ); }

foreach $func (@bifuncs1) { pp_def($func, HandleBad => 1, NoBadifNaN => 1, Pars => 'a(); b(); [o]c();', Inplace => [ 'a' ], Doc => inplace_doc( $func ), Code => code_bifunc($func), BadCode => badcode_bifunc($func), ); }

# Functions provided by extended -lm foreach $func (@ufuncs2) { pp_def($func, HandleBad => 1, NoBadifNaN => 1, GenericTypes => ['F','D'], Pars => 'a(); [o]b();', Inplace => 1, Doc => inplace_doc( $func ), Code => code_ufunc($func), BadCode => badcode_ufunc($func), ); }

foreach $func (@besufuncs) { my $fname = "bess$func"; pp_def($fname, HandleBad => 1, NoBadifNaN => 1, GenericTypes => ['F','D'], Pars => 'a(); [o]b();', Inplace => 1, Doc => inplace_doc( $fname ), Code => code_ufunc($func), BadCode => badcode_ufunc($func), ); }

foreach $func (@besbifuncs) { my $fname = "bess$func"; pp_def($fname, HandleBad => 1, NoBadifNaN => 1, GenericTypes => ['F','D'], Pars => 'a(); int n(); [o]b();', Inplace => [ 'a' ], Doc => inplace_doc( $fname ), Code => code_bifunc($func,'n','a','b'), BadCode => badcode_bifunc($func,'n','a','b'), ); }

if ($^O !~ /win32/i) { pp_def("lgamma", HandleBad => 1, Pars => 'a(); [o]b(); int[o]s()', Doc => $doco{"lgamma"}, Code => 'extern int signgam; $b() = lgamma($a()); $s() = signgam;', # what happens to signgam if $a() is bad? BadCode => 'extern int signgam; if ( $ISBAD(a()) ) { $SETBAD(b()); $SETBAD(s()); } else { $b() = lgamma($a()); $s() = signgam; }', ); } # if: os !~ win32

pp_def( 'badmask', Pars => 'a(); b(); [o]c();', Inplace => [ 'a' ], HandleBad => 1, Code => '$c() = finite($a()) ? $a() : $b();', BadCode => '$c() = ( finite($a()) && $ISGOOD(a()) ) ? $a() : $b();', CopyBadStatusCode => 'if ( a == c && $ISPDLSTATEBAD(a) ) PDL->propogate_badflag( c, 0 ); /* propogate badflag if inplace AND its changed */ $SETPDLSTATEGOOD(c); /* always make sure the output is "good" */ ', Doc => '=for ref

Clears all infs and nans in $a to the corresponding value in $b.

badmask can be run with $a inplace:

  badmask($a->inplace,0);
  $a->inplace->badmask(0);

', BadDoc => 'If bad values are present, these are also cleared.', );

pp_def( 'isfinite', Pars => 'a(); int [o]mask();', Inplace => 1, HandleBad => 1, Code => '$mask() = finite((double) $a()) != 0;', BadCode => '$mask() = finite((double) $a()) != 0 && $ISGOOD($a());', CopyBadStatusCode => 'if ( a == mask && $ISPDLSTATEBAD(a) ) PDL->propogate_badflag( mask, 0 ); /* propogate badflag if inplace AND its changed */ $SETPDLSTATEGOOD(mask); /* always make sure the output is "good" */ ', Doc => 'Sets $mask true if $a is not a NaN or inf (either positive or negative). Works inplace.', BadDoc => 'Bad values are treated as NaN or inf.', );

# Extra functions from cephes pp_def( "erfi", HandleBad => 1, NoBadifNaN => 1, GenericTypes => ['F','D'], Pars => 'a(); [o]b()', Inplace => 1, Doc => inplace_doc( "erfi" ), Code => 'extern double ndtri(double), SQRTH; $b() = SQRTH*ndtri((1+(double)$a())/2);', BadCode => 'extern double ndtri(double), SQRTH; if ( $ISBAD(a()) ) { $SETBAD(b()); } else { $b() = SQRTH*ndtri((1+(double)$a())/2); }', );

pp_def( "ndtri", HandleBad => 1, NoBadifNaN => 1, GenericTypes => ['F','D'], Pars => 'a(); [o]b()', Inplace => 1, Doc => inplace_doc( "ndtri" ), Code => 'extern double ndtri(double); $b() = ndtri((double)$a());', BadCode => 'extern double ndtri(double); if ( $ISBAD(a()) ) { $SETBAD(b()); } else { $b() = ndtri((double)$a()); }', );

pp_def("polyroots", Pars => 'cr(n); ci(n); [o]rr(m); [o]ri(m);', RedoDimsCode => 'int sn = $PDL(cr)->dims[0]; $SIZE(m) = sn-1;', GenericTypes => ['D'], Code => ' extern int cpoly( double *cr, double *ci, int deg, double *rr, double *ri ); int deg = $SIZE(n)-1, i; if (cpoly($P(cr), $P(ci), deg, $P(rr), $P(ri))) barf("PDL::Math::polyroots failed"); ', , Doc => '

Complex roots of a complex polynomial, given coefficients in order of decreasing powers.

 ($rr, $ri) = polyroots($cr, $ci);

',);

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

BUGS ^

Hasn't been tested on all platforms to ensure Cephes versions are picked up automatically and used correctly.

AUTHOR ^

Copyright (C) R.J.R. Williams 1997 (rjrw@ast.leeds.ac.uk), Karl Glazebrook (kgb@aaoepp.aao.gov.au) and Tuomas J. Lukka (Tuomas.Lukka@helsinki.fi). Portions (C) Craig DeForest 2002 (deforest@boulder.swri.edu).

All rights reserved. There is no warranty. You are allowed to redistribute this software / documentation under certain conditions. For details, see the file COPYING in the PDL distribution. If this file is separated from the PDL distribution, the PDL copyright notice should be included in the file.

syntax highlighting: