slatec.pd

use strict;
use warnings;
use PDL::Types;

our $VERSION = '2.097';
pp_setversion($VERSION);

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

=head1 NAME

PDL::Slatec - PDL interface to some LINPACK and EISPACK routines - DEPRECATED

=head1 SYNOPSIS

 use PDL::Slatec;

 ($ndeg, $r, $ierr, $c) = polyfit($x, $y, $w, $maxdeg, $eps);

=head1 DESCRIPTION

This module is now deprecated in favour of L<PDL::LinearAlgebra>.

This module serves the dual purpose of providing an interface to
parts of the slatec library and showing how to interface PDL
to an external library.
Using this library requires a Fortran compiler; the source for the routines
is provided for convenience.

Currently available are routines to:
manipulate matrices; calculate FFT's; 
and fit data using polynomials.

=head2 Piecewise cubic Hermite interpolation (PCHIP)

These routines are now in L<PDL::Primitive> as of PDL 2.096.

=cut

EOD

# add function definitions after finishing the first pp_addpm(), since this
# adds a '=head1 FUNCTIONS' line at the end of the text

pp_addpm(<<'END');
=head2 eigsys

=for ref

Eigenvalues and eigenvectors of a real positive definite symmetric matrix.

=for usage

 ($eigvals,$eigvecs) = eigsys($mat)

Note: this function should be extended to calculate only eigenvalues if called 
in scalar context!

This is the EISPACK routine C<rs>.

=head2 matinv

=for ref

Inverse of a square matrix

=for usage

 ($inv) = matinv($mat)

=head2 polyfit

Convenience wrapper routine about the C<polfit> C<slatec> function.
Separates supplied arguments and return values.

=for ref

Fit discrete data in a least squares sense by polynomials
in one variable.  Handles broadcasting correctly--one can pass in a 2D PDL (as C<$y>)
and it will pass back a 2D PDL, the rows of which are the polynomial regression
results (in C<$r>) corresponding to the rows of $y.

=for usage

 ($ndeg, $r, $ierr, $c, $coeffs, $rms) = polyfit($x, $y, $w, $maxdeg, [$eps]);

 $coeffs = polyfit($x,$y,$w,$maxdeg,[$eps]);

where on input:

C<$x> and C<$y> are the values to fit to a polynomial.
C<$w> are weighting factors
C<$maxdeg> is the maximum degree of polynomial to use and 
C<$eps> is the required degree of fit.

and the output switches on list/scalar context.  

In list context: 

C<$ndeg> is the degree of polynomial actually used
C<$r> is the values of the fitted polynomial 
C<$ierr> is a return status code, and
C<$c> is some working array or other (preserved for historical purposes)
C<$coeffs> is the polynomial coefficients of the best fit polynomial.
C<$rms> is the rms error of the fit.

In scalar context, only $coeffs is returned.

Historically, C<$eps> was modified in-place to be a return value of the
rms error.  This usage is deprecated, and C<$eps> is an optional parameter now.
It is still modified if present.

C<$c> is a working array accessible to Slatec - you can feed it to several
other Slatec routines to get nice things out.  It does not broadcast
correctly and should probably be fixed by someone.  If you are 
reading this, that someone might be you.

=for bad

This version of polyfit handles bad values correctly.  Bad values in
$y are ignored for the fit and give computed values on the fitted
curve in the return.  Bad values in $x or $w are ignored for the fit and
result in bad elements in the output.

=head2 polycoef

Convenience wrapper routine around the C<pcoef> C<slatec> function.
Separates supplied arguments and return values.                               

=for ref

Convert the C<polyfit>/C<polfit> coefficients to Taylor series form.

=for usage

 $tc = polycoef($l, $c, $x);

=head2 detslatec

=for ref

compute the determinant of an invertible matrix

=for example

  $mat = identity(5); # unity matrix
  $det = detslatec $mat;

Usage:

=for usage

  $determinant = detslatec $matrix;

=for sig

  Signature: detslatec(mat(n,m); [o] det())

C<detslatec> computes the determinant of an invertible matrix and barfs if
the matrix argument provided is non-invertible. The matrix broadcasts as usual.

This routine was previously known as C<det> which clashes now with
L<det|PDL::MatrixOps/det> which is provided by L<PDL::MatrixOps>.

=head2 fft

=for ref

Fast Fourier Transform

=for example

  $v_in = pdl(1,0,1,0);
  ($azero,$x,$y) = PDL::Slatec::fft($v_in);

C<PDL::Slatec::fft> is a convenience wrapper for L</ezfftf>, and
performs a Fast Fourier Transform on an input vector C<$v_in>. The
return values are the same as for L</ezfftf>.

=head2 rfft

=for ref

reverse Fast Fourier Transform

=for example

  $v_out = PDL::Slatec::rfft($azero,$x,$y);
  print $v_in, $vout
  [1 0 1 0] [1 0 1 0]

C<PDL::Slatec::rfft> is a convenience wrapper for L</ezfftb>,
and performs a reverse Fast Fourier Transform. The input is the same
as the output of L</PDL::Slatec::fft>, and the output
of C<rfft> is a data vector, similar to what is input into
L</PDL::Slatec::fft>.

=cut

END

use strict;

sub intcopy {
  my ($cvarname, $rvalue) = @_; # closures
  (sub {"&".$cvarname->(@_)}, sub {"PDL_LongLong *"},
    sub {"PDL_LongLong ".$cvarname->(@_)." = ".$rvalue->(@_).";"});
}
# defslatec( $pdlname, $funcnames, $argstr, $docstring, $funcref )
#
# $pdlname is the name of the PDL function to be created
# $funcnames is a reference to a hash array, whose keys define
# the single (S), and double precision (D) names of the
# SLATEC routines to be linked to.
#
# $argstr is a list of arguments expected by the SLATEC routine
# - some of the allowed type names are:
#   FuncRet
#     - specifies that this is a function, not a subroutine, and
#       that the output of the function should be stored in this
#       variable
#   IntVal
#     - used for passing in a value to func, either calculated or a constant,
#       but not a caller-supplied arg
#   Logical
#     - A Fortran Boolean, defined as same size as REAL/INTEGER so 'long'
#
# $docstring gives the text to be used as the function documentation
#
# $funcref gets placed within a '=for ref' pod statement at the
# start of the documentation - ie it is placed before the
# text within $docstring. This string gets printed out
# in the perldl shell after a '?? string' command
#
# [Par_forcetype(undef=none,ref=OtherPar), C_argin, C_paramtype, C_decl_value]
my %KIND2I = (
  Mat => ['', sub{"\$P($_[0][2])"}, sub {PDL::Type->new($_[0])->ctype." *"}],
  FuncRet => ['', sub {()}, sub {()}],
  Integer => ['longlong ', sub{"\$P($_[0][2])"}, sub {"PDL_LongLong *"}],
  # scalar for calling C functions
  IntScalar => ['longlong ', sub{"\$$_[0][2]()"}, sub {"PDL_LongLong "}],
  MatScalar => ['', sub{"\$$_[0][2]()"}, sub {PDL::Type->new($_[0])->ctype." "}],
  Logical => ['int ', sub{"\$P($_[0][2])"}, sub {"PDL_Long *"}],
  OtherParDim => [\'PDL_LongLong ', intcopy(sub {"dim_$_[0][2]"}, sub {"\$SIZE($_[0][2])"})],
  IntVal => [undef, intcopy(sub {$_[0][2]}, sub {$_[0][4]})],
);
my %ignore_ppar = map +($_=>1), qw(IntVal OtherParDim);
my %prototype = ( F => "float", D => "double" );
my %ftypes = (S => 'F', D => 'D');
our $debug = 0;  # print out calls to pp_def; use "local $debug=1" before func
sub defslatec {
    my ($pname,$fnames,$argstr,$docstring,$funcref) = @_;
    my @args = grep /\S/, split ';', $argstr;
    my @args2 = map {
		/^\s*([a-zA-Z]+)\s+	# "Type name"
		  ((?:\[[^]]*\])?)\s*	# Options
		  ([a-zA-Z]+)\s*	# Par name
		  (?:
		    ((?:\(.*\))?)\s*$	# Dims
		    |
		    =\s*((?:.*?)?)\s*$	# value for IntVal
		  )
		 /x or die("Invalid slatec par $_");
		[$1,$2,$3,$4,$5]} @args;
    # is this for a function (Type name eq "FuncRet")
    # or a subroutine?
    die "Only one FuncRet allowed in pars list.\n"
      if (my @funcret = grep $_->[0] eq "FuncRet", @args2) > 1;
    my $fpar = @funcret ? "\$$funcret[0][2]()" : undef;
    my $pars = join '; ', map
      +($KIND2I{$_->[0]}[0] // die "Invalid ppars ",join(',',@$_),"\n") .
        join(' ', grep length, $_->[1], "$_->[2]$_->[3]"),
      grep !$ignore_ppar{$_->[0]}, @args2;
    my $otherpars = join ';', map
      ${$KIND2I{$_->[0]}[0]} . join(' => ', @$_[2,2]),
      grep ref $KIND2I{$_->[0]}[0], @args2;
    my $argorder;
    if ($otherpars) {
      my @mand = map $_->[2], grep $_->[1] !~ /\bo\b/, @args2;
      my @opt = map $_->[2], grep $_->[1] =~ /\bo\b/, @args2;
      $argorder = [@mand, @opt];
    }
    my @talts = map
      [$ftypes{$_} // die("FTYPE $_ NOT THERE\n"),$fnames->{$_}],
      sort keys %$fnames;
    my $func = "PDL_FORTRAN(\$T".join('',map {$_->[0]} @talts) . "(" .
      join(',',map qq{$_->[1]}, @talts)."))";
    if ( defined $fpar ) { $func = "$fpar = $func"; }
    my $funcargs = join ',', map $KIND2I{$_->[0]}[1]->($_), @args2;
    my @protos = map $KIND2I{$_->[0]}[2], @args2;
    my @cheaders;
    foreach my $t ( @talts ) {
        my @decl_args;
        for my $ind (0..$#protos) {
          my $arg = $protos[$ind];
          next unless my @ret = $arg->($t->[0]);
          push @decl_args, "$ret[0]$args2[$ind][2]";
        }
        push @cheaders, (defined $fpar ? $prototype{$t->[0]} : "void")
          . " PDL_FORTRAN($t->[1])(" .
          join(',', @decl_args) . ");";
    }
    # add on the function reference, if supplied, to the start of
    # the doc string
    if ( defined $docstring ) {
      $docstring = "=for ref\n\n$funcref\n\n$docstring" if defined $funcref;
    } else {
      $docstring = '';
    }
    # IntVals
    my @intcopy_code = map +($KIND2I{$_->[0]}[3]||sub{})->($_), @args2;
    my $code = join("\n", @intcopy_code, "$func($funcargs);");
    print <<"ENDDBG" if $debug;
pp_def('$pname',
  Pars => '$pars',@{[ !$otherpars ? '' : "
  OtherPars => '$otherpars',"]}@{[ !$argorder ? '' : "
  ArgOrder => [qw(@$argorder)],"]}
  Code => pp_line_numbers(__LINE__, <<'EOF'),
$code
EOF
  CHeader => '@cheaders',
  GenericTypes => [qw(@{[join ", ", map $_->[0], @talts]})],
  Doc => <<'EOF',
$docstring
EOF
);
ENDDBG
    pp_def($pname,
      Pars => $pars,
      ($otherpars ? (OtherPars => $otherpars) : ()),
      $argorder ? (ArgOrder => $argorder) : (),
      Code => $code,
      CHeader => join("\n", @cheaders),
      GenericTypes => [map $_->[0], @talts],
      Doc => $docstring,
    );
}

pp_addhdr(qq|
static void MAIN__ () {
   /* Cheat to define MAIN__ symbol */                                          
   croak("This should never happen");                                           
}

/* these two need to be done manually because we don't use defslatec for them */
void PDL_FORTRAN(polfit)(
  PDL_LongLong *N, PDL_Float *X, PDL_Float *Y, PDL_Float *W, PDL_LongLong *MAXDEG,
  PDL_LongLong *NDEG, PDL_Float *EPS, PDL_Float *R, PDL_LongLong *IERR, PDL_Float *A
);
void PDL_FORTRAN(dpolft)(
  PDL_LongLong *N, PDL_Double *X, PDL_Double *Y, PDL_Double *W, PDL_LongLong *MAXDEG,
  PDL_LongLong *NDEG, PDL_Double *EPS, PDL_Double *R, PDL_LongLong *IERR, PDL_Double *A
);
void PDL_FORTRAN(slatecbarf)();

|);

pp_add_exported('',"eigsys matinv polyfit polycoef");

pp_addpm(<<'END');

use PDL::Core;
use PDL::Basic;
use PDL::Primitive;
use PDL::Ufunc;
use strict;

*eigsys = \&PDL::eigsys;

sub PDL::eigsys {
	my($h) = @_;
	$h = float($h);
	rs($h,
		my $eigval=PDL->null,
		1,my $eigmat=PDL->null,
		my $errflag=PDL->null
	);
#	print $covar,$eigval,$eigmat,$fvone,$fvtwo,$errflag;
	if(sum($errflag) > 0) {
		barf("Non-positive-definite matrix given to eigsys: $h\n");
	}
	return ($eigval,$eigmat);
}

*matinv = \&PDL::matinv;

sub PDL::matinv {
	my($m) = @_;
	my(@dims) = $m->dims;

	# Keep from dumping core (FORTRAN does no error checking)
	barf("matinv requires a 2-D square matrix")
		unless( @dims >= 2 && $dims[0] == $dims[1] );
  
	$m = $m->copy(); # Make sure we don't overwrite :(
	my ($ipvt,$info) = gefa($m);
	if(sum($info) > 0) {
		barf("Uninvertible matrix given to inv: $m\n");
	}
	gedi($m,$ipvt,1);
	$m;
}

*detslatec = \&PDL::detslatec;
sub PDL::detslatec {
	my($m) = @_;
	$m = $m->copy(); # Make sure we don't overwrite :(
	gefa($m,(my $ipvt=null),(my $info=null));
	if(sum($info) > 0) {
		barf("Uninvertible matrix given to inv: $m\n");
	}
	my ($det) = gedi($m,$ipvt,10);
	return $det->slice('(0)')*10**$det->slice('(1)');
}

sub prepfft {
	my($n) = @_;
	my $ifac = ezffti($n,my $wsave = PDL->zeroes(float(),$n*3));
	return ($wsave, $ifac);
}

sub fft (;@) {
	my ($v) = @_;
	ezfftf($v, prepfft($v->getdim(0)));
}

sub rfft {
	my ($az,$x,$y) = @_;
	ezfftb($az,$x,$y,prepfft($x->getdim(0)));
}

# polynomial fitting routines
# simple wrappers around the SLATEC implementations

*polyfit = \&PDL::polyfit;
sub PDL::polyfit {
  barf 'Usage: polyfit($x, $y, $w, $maxdeg, [$eps]);'
    unless (@_ == 5 || @_==4 );

  my ($x_in, $y_in, $w_in, $maxdeg_in, $eps_io) = map PDL->topdl($_), @_;

  my $template_ind = maximum_n_ind([(map $_->ndims-1, $x_in, $y_in, $w_in), $maxdeg_in->ndims, defined $eps_io ? $eps_io->ndims : -1], 1)->sclr;
  my $template = $_[$template_ind];
  # if $w_in does not match the data vectors ($x_in, $y_in), then we can resize
  # it to match the size of $y_in :
  $w_in = $w_in + $template->zeroes;
  $eps_io = $eps_io + $template->slice('(0)')->zeroes; # also needs to match but with one less dim
  my $max_maxdeg = $maxdeg_in->max->sclr;

  # Now call polfit
  my $rms = pdl($eps_io);
  my ($ndeg, $r, $ierr, $a1, $coeffs) = polfit($x_in, $y_in, $w_in, $maxdeg_in, $rms, $max_maxdeg+1);
  # Preserve historic compatibility by flowing rms error back into the argument
  $eps_io .= $rms if UNIVERSAL::isa($_[4],'PDL');

  # Return the arrays
  wantarray ? ($ndeg, $r, $ierr, $a1, $coeffs, $rms) : $coeffs;
}


*polycoef = \&PDL::polycoef;
sub PDL::polycoef {
  barf 'Usage: polycoef($l, $c, $x);'
    unless @_ == 3;

  # Allocate memory for return PDL
  # Simply l + 1 but cant see how to get PP to do this - TJ
  # Not sure the return type since I do not know
  # where PP will get the information from
  my $tc = PDL->zeroes( abs($_[0]) + 1 );                                     

  # Run the slatec routine
  pcoef($_[0], $_[1], $tc, $_[2]);

  # Return results
  return $tc;

}
END

defslatec('svdc',
	{S => 'ssvdc'},
	'Mat		x	(n,p);
	 IntVal		ldx=$SIZE(n);
	 IntVal		n=$SIZE(n);
	 IntVal		p=$SIZE(p);
	 Mat	[o]	s	(p);
	 Mat	[o]	e	(p);
	 Mat	[o]	u	(n,p);
	 IntVal		ldu=$SIZE(n);
	 Mat	[o]	v	(p,p);
	 IntVal		ldv=$SIZE(p);
	 Mat	[t]	work	(n);
	 Integer	job	();
	 Integer [o]	info	();
	',
'singular value decomposition of a matrix'
);

defslatec('poco',
	{S => 'spoco', D => 'dpoco'},
	'Mat [io]	a	(n,n);
	 IntVal		lda=$SIZE(n);
	 IntVal		n=$SIZE(n);
	 Mat [o]	rcond	();
	 Mat [o]	z	(n);
	 Integer [o]	info	();
	',
'Factor a real symmetric positive definite matrix
and estimate the condition number of the matrix.'
);

defslatec('geco',
	{S => 'sgeco', D => 'dgeco'},
	'Mat		a	(n,n);
	 IntVal		lda=$SIZE(n);
	 IntVal		n=$SIZE(n);
	 Integer [o]	ipvt	(n);
	 Mat	 [o]	rcond	();
	 Mat	 [o]	z	(n);
	',
'Factor a matrix using Gaussian elimination and estimate
the condition number of the matrix.'
);

defslatec('gefa',
	{S => 'sgefa', D => 'dgefa'},
	'Mat [io]	a	(n,n);
	 IntVal		lda=$SIZE(n);
	 IntVal		n=$SIZE(n);
	 Integer [o]	ipvt	(n);
	 Integer [o]	info	();
	',
'Factor a matrix using Gaussian elimination.'
);

# pofa and sqrdc aren't (yet?) implemented
#
defslatec('podi',
	{S => 'spodi', D => 'dpodi'},
	'Mat [io]	a	(n,n);
	 IntVal		lda=$SIZE(n);
	 IntVal		n=$SIZE(n);
	 Mat	[o]	det	(two=2);
	 Integer	job	();
	',
'Compute the determinant and inverse of a certain real
symmetric positive definite matrix using the factors
computed by L</poco>.'
);

defslatec('gedi',
	{S => 'sgedi', D => 'dgedi'},
	'Mat [io]	a	(n,n);
	 IntVal		lda=$SIZE(n);
	 IntVal		n=$SIZE(n);
	 Integer	ipvt	(n);
	 Mat	 [o]	det	(two=2);
	 Mat	 [t]	work	(n);
	 Integer	job	();
	',
'Compute the determinant and inverse of a matrix using the
factors computed by L</geco> or L</gefa>.'
);

defslatec('gesl',
	{S => 'sgesl', D => 'dgesl'},
	'Mat		a	(lda,n);
	 IntVal		lda=$SIZE(lda);
	 IntVal		n=$SIZE(n);
	 Integer	ipvt	(n);
	 Mat [io]	b	(n);
	 Integer	job	();
	',
'Solve the real system C<A*X=B> or C<TRANS(A)*X=B> using the
factors computed by L</geco> or L</gefa>.'
);

defslatec('rs',
	 {S => 'rsfoo'},
	'IntVal		lda=$SIZE(n);
	 IntVal		n=$SIZE(n);
	 Mat		a	(n,n);
	 Mat	[o]	w	(n);
	 Integer	matz	();
	 Mat	[o]	z	(n,n);
	 Mat	[t]	fvone	(n);
	 Mat	[t]	fvtwo	(n);
	 Integer [o]	ierr	();
	',
'This subroutine calls the recommended sequence of
subroutines from the eigensystem subroutine package (EISPACK)
to find the eigenvalues and eigenvectors (if desired)
of a REAL SYMMETRIC matrix.'

);

# XXX wsave : at least 3n
defslatec('ezffti',
	{S => 'ezffti'},
	'Integer	n	();
	 Mat [io]	wsave(foo);
	 Integer [o]	ifac(ni=15);
	',
'Subroutine ezffti initializes the work array C<wsave(3n or more)>
and C<ifac()>
which is used in both L</ezfftf> and L</ezfftb>.
The prime factorization
of C<n> together with a tabulation of the trigonometric functions
are computed and stored in C<wsave()>.'
);

# XXX Correct for azero, a and b
defslatec('ezfftf',
	{S => 'ezfftf'},
	'IntVal		n=$SIZE(n);
	 Mat		r(n);
	 Mat [o]	azero();
	 Mat [o]	a(n);
	 Mat [o]	b(n);
	 Mat		wsave(foo);
	 Integer	ifac(ni=15);
	'
);

defslatec('ezfftb',
	{S => 'ezfftb'},
	'IntVal		n=$SIZE(n);
	 Mat  [o]	r(n);
	 Mat		azero();
	 Mat		a(n);
	 Mat		b(n);
	 Mat		wsave(foo);
	 Integer	ifac(ni=15);
	'
);

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

defslatec('pcoef',
      {S => 'pcoef', D => 'dpcoef'},
      '
      Integer   l ();
      Mat      c ();
      Mat [o]  tc (bar);
      Mat      a (foo);
      ',
'Convert the C<polfit> coefficients to Taylor series form.
C<c> and C<a()> must be of the same type.'
);                                                                            

defslatec('polyvalue',
      {S => 'pvalue', D => 'dp1vlu'},
      '
      Integer     l ();
      OtherParDim nder;
      Mat         x    ();
      Mat [o]     yfit ();
      Mat [o]     yp   (nder);
      Mat         a    (foo);
      ', <<'EOF',
=for ref

Use the coefficients C<c> generated by L</polyfit> (or L</polfit>) to evaluate
the polynomial fit of degree C<l>, along with the first C<nder> of its
derivatives, at a specified point C<x>.
Broadcasts correctly over multiple C<x> positions.

=for usage

 ($yfit, $yp) = polyvalue($l, $nder, $x, $c);
EOF
);                                                                            

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

# This version of polfit accepts bad values and also allows broadcasting

# indices: 
#  n   runs across input points; 
#  foo runs across wacky Slatec buffer size;
#  bar runs across polynomial coefficients.
pp_def('polfit',
  Pars => 'x(n); y(n); w(n); longlong maxdeg(); longlong [o]ndeg(); [io]eps(); [o]r(n); longlong [o]ierr(); [o]a(foo=CALC(3*($SIZE(n) + $SIZE(bar)))); [o]coeffs(bar);[t]xtmp(n);[t]ytmp(n);[t]wtmp(n);[t]rtmp(n)',
  OtherPars => 'IV max_maxdeg_plus1 => bar',
  GenericTypes => ['F','D'],
  HandleBad => 1,
  NoBadifNaN => 1,
  CHeader => <<'EOF',
void PDL_FORTRAN(dpcoef)(PDL_LongLong *l,PDL_Double *c,PDL_Double *tc,PDL_Double *a);
void PDL_FORTRAN(pcoef)(PDL_LongLong *l,PDL_Float *c,PDL_Float *tc,PDL_Float *a);
EOF
  Code => 'PDL_LongLong howmany = PDL_IF_BAD(0,$SIZE(n)), i;
              if($SIZE(n)<$maxdeg())
                $CROAK("Must have at least <n> points to fit <n> coefficients");
              PDL_IF_BAD(
              loop(n) %{   /* get rid of bad values.  Call polfit with [xyw]tmp instead of [xyz]. */
                if ($ISGOOD(y()) && $ISGOOD(x()) && $ISGOOD(w())) {
                  $xtmp(n=>howmany) = $x();
                  $ytmp(n=>howmany) = $y();
                  $wtmp(n=>howmany) = $w();
                  howmany++;
                }
              %}
	      if (howmany <= $maxdeg()) {
		/* Not enough good datapoints -- set this whole row to BAD. */
                loop(n) %{ $SETBAD(r()); %}
                $ierr() = 2;
              } else,) {
                  /* Found enough good datapoints for a fit -- do the fit */
                  $GENERIC() zero = 0;
                /* Do the fit */
                PDL_FORTRAN($TFD(polfit,dpolft))(
                  &howmany,
                  PDL_IF_BAD($P(xtmp),$P(x)),
                  PDL_IF_BAD($P(ytmp),$P(y)),
                  PDL_IF_BAD($P(wtmp),$P(w)),
                  $P(maxdeg),$P(ndeg),$P(eps),
                  PDL_IF_BAD($P(rtmp),$P(r)),
                  $P(ierr),$P(a)
                );
		PDL_LongLong maxord = $a(foo=>0)+0.5;
		PDL_LongLong ord = PDLMIN($maxdeg(),$a(foo=>(maxord*3)+2));
		/* Extract the polynomial coefficients into coeffs -- used for bad values */
                PDL_FORTRAN($TFD(,d)pcoef)(&(ord), &(zero), $P(coeffs), $P(a));
                for(i=ord+1; i<=$maxdeg(); i++)
                   $coeffs(bar=>i) = 0;
                PDL_IF_BAD(
                howmany=0;
                loop(n) %{  /* replace bad values */
                  if ($ISGOOD(y())) {
                    $r() = $rtmp(n=>howmany);
                    howmany++;
                  } else if($ISGOOD(x())) {
		     /* Bad values are omitted from the call to polfit, so we must reconstitute them on return */
	             /* (just because a value is bad in y, does not mean the fit itself is bad there) */
                     /* */
                     /* The number in ord is not the number of coefficients in the polynomial, it is the highest */
                     /* order coefficient -- so 3 for a cubic, which has 4 coefficients. */
		     PDL_LongLong ii;
                     $GENERIC() acc = 0;
                     for( ii=ord; ii>0; ii-- ) {
                        acc += $coeffs(bar=>ii);
                        acc *= $x();
                     }
                     acc += $coeffs(bar=>0);
                     $r() = acc;
                  } else {
                    /* x and y are bad here... */
		    $SETBAD(r());
                  }
                %},)
              }',
  Doc => 'Fit discrete data in a least squares sense by polynomials
          in one variable. C<x()>, C<y()> and C<w()> must be of the same type.
	  This version handles bad values appropriately',
);

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

pp_addpm(<<'EOD');
=head1 AUTHOR

Copyright (C) 1997 Tuomas J. Lukka. 
Copyright (C) 2000 Tim Jenness, Doug Burke.            
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();