Christian Soeller > PDL-2.4.3 > PDL::Math

PDL-2.4.3.tar.gz

Annotate this POD

# CPAN RT

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

# 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.

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

# 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);

',);