use strict;

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.

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.

=cut

### Kludge for backwards compatibility with older scripts
### This should be deleted at some point later than 21-Nov-2003.
BEGIN {use PDL::MatrixOps;}

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 ('ndtri',
"=for ref

The value for which the area under the
Gaussian probability density function (integrated from
minus infinity) is equal to the argument (cf L<erfi|/erfi>).");

doco(qw/bessj0 bessj1 bessy0 bessy1/, 
     "The standard Bessel function." );

doco( qw/bessjn bessyn/, 
'=for ref

The standard Bessel function.
This has a second integer
argument which gives the order of the function required.
');

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.

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 C<infs> and C<nans> in C<$a> to the corresponding value in C<$b>.

badmask can be run with C<$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 C<$mask> true if C<$a> is not a C<NaN> or C<inf> (either positive or negative). Works inplace.',
       BadDoc =>
'Bad values are treated as C<NaN> or C<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 => '

=for ref

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

=for usage

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

',);

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

=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).
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.

=cut

EOD
pp_done();