pp_addpm({At=>Top},<<'EOD');
=head1 NAME
PDL::Math - extended mathematical operations and special functions
=head1 SYNOPSIS
use PDL::Math;
use PDL::Graphics::TriD;
imag3d [SURF2D,bessj0(rvals(zeroes(50,50))/2)];
=head1 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.
The functions are usually available from the system
maths library, however if they are not (determined
when PDL is compiled) a version from the Cephes math
library is used.
=cut
EOD
# Internal doc util
my %doco;
sub doco {
my @funcs = @_;
my $doc = pop @funcs;
for (@funcs) { $doco{$_} = $doc }
}
doco (qw/acos asin atan tan/,
'The usual trigonometric function.');
doco (qw/cosh sinh tanh acosh asinh atanh/,
'The standard hyperbolic function.');
doco (qw/ceil floor rint/,
'Round to integral values in floating-point format');
doco( 'pow',"Synonym for `**'");
doco ('erf',"The error function");
doco ('erfc',"The complement of the error function");
doco ('erfi',"The inverse of the error function");
doco(qw/bessj0 bessj1 bessy0 bessy1 bessjn bessyn/, <<'EOD');
=for ref
The standard Bessel Functions
=cut
EOD
$doco{'bessjn'} .= <<'EOD';
=pod
This has a second integer
argument which gives the order of the function required.
=cut
EOD
$doco{'bessyn'} .= <<'EOD';
=pod
This has a second integer
argument which gives the order of the function required.
=cut
EOD
if ($^O !~ /win32/i) { # doesn't seem to be in the MS VC lib
doco( 'lgamma' ,<<'EOD');
=for ref
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.
=cut
EOD
}
pp_addhdr('
#include <math.h>
#include "protos.h"
');
if ($^O =~ /MSWin/) {
# _finite in VC++ 4.0
pp_addhdr('
#define finite _finite
#include <float.h>
');
}
# Standard `-lm'
my (@ufuncs1) = qw(acos asin atan cosh sinh tan tanh); # F,D only
my (@ufuncs1g) = qw(ceil floor rint); # Any type
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
foreach $func (@ufuncs1) {
pp_def($func,
GenericTypes => [F,D],
Pars => 'a(); [o]b();',
Doc => $doco{$func},
Code => '
$b() = '.$func.'($a());
');
}
foreach $func (@ufuncs1g) {
pp_def($func,
Pars => 'a(); [o]b();',
Doc => $doco{$func},
Code => '
$b() = '.$func.'($a());
');
}
foreach $func (@bifuncs1) {
pp_def($func,
Pars => 'a(); b(); [o]c();',
Doc => $doco{$func},
Code => '
$c() = '.$func.'($a(),$b());
');
}
# Functions provided by extended -lm
foreach $func (@ufuncs2) {
pp_def($func,
GenericTypes => [F,D],
Pars => 'a(); [o]b();',
Doc => $doco{$func},
Code => '
$b() = '.$func.'($a());
');
}
foreach $func (@besufuncs) {
pp_def("bess$func",
GenericTypes => [F,D],
Pars => 'a(); [o]b();',
Doc => $doco{"bess$func"},
Code => '
$b() = '.$func.'($a());
');
}
foreach $func (@besbifuncs) {
pp_def("bess$func",
GenericTypes => [F,D],
Pars => 'a(); int n(); [o]b();',
Doc => $doco{"bess$func"},
Code => '
$b() = '.$func.'($n(),$a());
');
}
if ($^O !~ /win32/i) {
pp_def("lgamma",
Pars => 'a(); [o]b(); int[o]s()',
Doc => $doco{"lgamma"},
Code => '
extern int signgam;
$b() = lgamma($a());
$s() = signgam;
');
}
pp_def('badmask',
Pars => 'a(); b(); [o]c();',
Doc => 'Clears all infs and nans in a to the corresponding value in b',
Code => '
$c() = finite($a()) ? $a() : $b();
');
# Extra functions from cephes
pp_def("erfi",
GenericTypes => [F,D],
Pars => 'a(); [o]b()',
Doc => $doco{"erfi"},
Code => '
extern double ndtri(double);
$b() = ndtri((double)$a());
');
pp_def("svd",
Pars => 'a(n,m); [o]u(n,m); [o,phys]z(n); [o]v(n,n);',
GenericTypes => [D],
Code => '
extern void SVD( double *W, double *Z, int nRow, int nCol );
int sm = $SIZE(m), sn = $SIZE(n), i;
double *w, *t, zv;
t = w = (double *) malloc(sn*(sm+sn)*sizeof(double));
loop (m) %{
loop(n) %{
*t++ = $a();
%}
%}
SVD(w, $P(z), sm, sn);
t = w;
loop (n) %{
zv = sqrt($z());
$z() = zv;
%}
loop (m) %{
loop (n) %{
$u() = *t++/$z();
%}
%}
loop (n) %{
for (i=0;i<sn;i++) {
$v(n0=>i, n1=>n) = *t++;
}
%}
free(w);
',
, Doc => '
=for ref
Singular value decomposition of array.
=for usage
($u, $s, $v) = svd($a);
',);
# XXX The next two really need driver routines...
pp_def("eigens",
Pars => '[phys]a(m); [o,phys]ev(n,n); [o,phys]e(n)',
GenericTypes => [D],
Code => '
extern void eigens( double *A, double *RR, double *E, int N );
register int sn = $SIZE(n);
if($SIZE(m) != (sn * (sn + 1))/2) {
barf("Wrong sized args for eigens");
}
eigens($P(a), $P(ev), $P(e), sn);
',
PMCode =>'
sub PDL::eigens {
my ($a) = @_;
my (@d) = $a->dims;
barf "Need real square matrix for eigens"
if $#d != 1 or $d[0] != $d[1];
my ($n) = $d[0];
my ($sym) = 0.5*($a + $a->mv(0,1));
my ($err) = PDL::max(abs($sym));
barf "Need symmetric component non-zero for eigens"
if $err == 0;
$err = PDL::max(abs($a-$sym))/$err;
warn "Using symmetrized version of the matrix in eigens"
if $err > 1e-5;
# Get lower triangular form --
my ($lt) = PDL::where($sym,$sym->xvals <= $sym->yvals)->copy;
my ($ev) = PDL->zeroes($n,$n);
my ($e) = PDL->zeroes($n);
&PDL::_eigens_int($lt, $ev, $e);
return $ev, $e;
}
'
, Doc => '
=for ref
Eigenvalues and -vectors of a symmetric square matrix. If passed
an asymmetric matrix, the routine will warn and symmetrize it.
=for usage
($e, $ev) = eigens($a);
',);
# XXX Destroys a!!!
# To use the new a again, must store both a and ips.
pp_def("simq",
Pars => '[phys]a(n,n); [phys]b(n); [o,phys]x(n); int [o,phys]ips(n)',
OtherPars => 'int flag;',
GenericTypes => [D],
Code => '
extern int simq( double *A, double *B, double *X,
int n, int flag, int *IPS );
simq($P(a),$P(b),$P(x),$SIZE(n),$COMP(flag),$P(ips));
',
Doc => '
=for ref
Solution of simultaneous linear equations, a x = b.
a is an n x n matrix (i.e., a vector of length n*n), stored row-wise:
that is, a(i,j) = a[ij], where ij = i*n + j. While this is the
transpose of the normal column-wise storage, this corresponds to
normal PDL usage. The contents of matrix a may be altered (but may be
required for subsequent calls with flag = -1).
b, x, ips are vectors of length n.
Set flag=0 to solve. Set flag=-1 to do a new back substitution for
different b vector using the same a matrix previously reduced when
flag=0 (the ips vector generated in the previous solution is also
required).
');
pp_def("squaretotri",
Pars => 'a(n,n); b(m)',
Code => '
register int mna=0, nb=0, ns = $SIZE(n);
#if (PERL_VERSION >= 5) && (PERL_SUBVERSION >= 57)
dXSARGS;
#endif
threadloop %{
if($SIZE(m) != (ns * (ns+1))/2) {
barf("Wrong sized args for squaretotri");
}
loop(m) %{
$b() = $a(n0 => mna, n1 => nb);
mna++; if(mna > nb) {mna = 0; nb ++;}
%}
%}
',
Doc => '
=for ref
Convert a symmetric square matrix to triangular vector storage
=cut
',
);
pp_addpm({At=>Bot},<<'EOD');
sub eigen_c {
my($mat) = @_;
my $s = $mat->getdim(0);
my $z = zeroes($s * ($s+1) / 2);
my $ev = zeroes($s);
squaretotri($mat,$z);
my $k = 0 * $mat;
PDL::eigens($z, $k, $ev);
return ($ev, $k);
}
=head1 BUGS
Hasn't been tested on all platforms to ensure Cephes
versions are picked up automatically and used correctly.
=head1 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).
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 copyright notice should be included in the file.
=cut
EOD
pp_done();