diff --git a/Changes b/Changes index 402a2e9c0..144582e88 100644 --- a/Changes +++ b/Changes @@ -28,6 +28,7 @@ - finish, then split PDL::IO::ENVI out to separate distro - split PDL::Graphics::TriD out to separate distro - fix approx_artol for complex values (#507) - thanks @wlmb +- split PDL::GSL out to separate distro 2.095 2024-11-03 - add PDL_GENTYPE_IS_{REAL,FLOATREAL,COMPLEX,SIGNED,UNSIGNED}_##ppsym (#502) diff --git a/Libtmp/GSL/.gitignore b/Libtmp/GSL/.gitignore deleted file mode 100644 index ca34462d6..000000000 --- a/Libtmp/GSL/.gitignore +++ /dev/null @@ -1,27 +0,0 @@ -lib/PDL/GSL/CDF.c -lib/PDL/GSL/CDF.pm -lib/PDL/GSL/CDF.xs -lib/PDL/GSL/DIFF.c -lib/PDL/GSL/DIFF.pm -lib/PDL/GSL/DIFF.xs -lib/PDL/GSL/INTEG.c -lib/PDL/GSL/INTEG.pm -lib/PDL/GSL/INTEG.xs -lib/PDL/GSL/INTERP.c -lib/PDL/GSL/INTERP.pm -lib/PDL/GSL/INTERP.xs -lib/PDL/GSL/LINALG.c -lib/PDL/GSL/LINALG.pm -lib/PDL/GSL/LINALG.xs -lib/PDL/GSL/MROOT.c -lib/PDL/GSL/MROOT.pm -lib/PDL/GSL/MROOT.xs -lib/PDL/GSL/RNG.c -lib/PDL/GSL/RNG.pm -lib/PDL/GSL/RNG.xs -lib/PDL/GSL/SF.c -lib/PDL/GSL/SF.pm -lib/PDL/GSL/SF.xs -lib/PDL/Stats/Distr.c -lib/PDL/Stats/Distr.pm -lib/PDL/Stats/Distr.xs diff --git a/Libtmp/GSL/Makefile.PL b/Libtmp/GSL/Makefile.PL deleted file mode 100644 index 004ac1509..000000000 --- a/Libtmp/GSL/Makefile.PL +++ /dev/null @@ -1,62 +0,0 @@ -use strict; -use warnings; -use ExtUtils::MakeMaker; -use PDL::Core::Dev; -use File::Spec::Functions; - -sub get_gsl_config { - my ($flags) = @_; - no warnings 'exec'; - `gsl-config $flags`; -} - -# Version check -chomp (my $version = get_gsl_config('--version') // '0'); -my $new_enough = 0; -if (!$version) { - $version = 'UNKNOWN VERSION'; -} else { - my $major = (split /\./,$version)[0]; - $new_enough= $major >= 2; -} -if (!$new_enough) { - write_dummy_make("Not building GSL modules: GSL version $version found, but need at least 2.0"); - return; -} - -# the real stuff happens in the subdirs -sub get_gsl_libs { - my $lib = ($ENV{GSL_LIBS} || - get_gsl_config('--libs') || - warn "\tno GSL link info (libgsl probably not available)\n"); - my $inc = ($ENV{GSL_INC} || - get_gsl_config('--cflags') || - warn "\tno GSL include info (libgsl probably not available)\n\n"); - chomp $lib; chomp $inc; - ($inc,$lib); -} - -# these will be used in the subdirs -my ($GSL_includes, $GSL_libs) = get_gsl_libs(); - -my @pd_srcs; -undef &MY::init_PM; # suppress warning -*MY::init_PM = sub { - package MY; # so that "SUPER" works right - my ($self) = @_; - $self->SUPER::init_PM; - @pd_srcs = ::pdlpp_eumm_update_deep($self); -}; - -undef &MY::postamble; # suppress warning -*MY::postamble = sub { pdlpp_postamble_int(@pd_srcs); }; - -WriteMakefile( - NAME => 'PDL::GSL', - VERSION_FROM => 'lib/PDL/GSL/CDF.pd', - MIN_PERL_VERSION => '5.014', - INC => join(' ', "-I".curdir(), $GSL_includes), - LIBS => [$GSL_libs], - clean => { FILES => join ' ', qw(MANIFEST.bak) }, - NO_MYMETA => 1, -); diff --git a/Libtmp/GSL/README b/Libtmp/GSL/README deleted file mode 100644 index 1c23bad1f..000000000 --- a/Libtmp/GSL/README +++ /dev/null @@ -1,12 +0,0 @@ -PDL interface to GSL by Christian Pellegrin , -Copyleft, 1999 see COPYING in the root of the PDL build tree. -The documentation is in the .pd file and is automaticaly generated. - -We need: - - - an example for every function (C<=for example> section in docs). - - a test for every SF function (to go into t/gsl_sf.t) - -Patches by GSL users are most welcome! - -2002 Christian Soeller diff --git a/Libtmp/GSL/gslerr.h b/Libtmp/GSL/gslerr.h deleted file mode 100644 index c55a05eec..000000000 --- a/Libtmp/GSL/gslerr.h +++ /dev/null @@ -1,8 +0,0 @@ -#include - -#define GSLERR(x,y) \ - { \ - int status = x y; \ - if (status) \ - return PDL->make_error(PDL_EUSERERROR, "Error in %s: %s", #x, gsl_strerror(status)); \ - } diff --git a/Libtmp/GSL/lib/PDL/Demos/GSL_CDF.pm b/Libtmp/GSL/lib/PDL/Demos/GSL_CDF.pm deleted file mode 100644 index cd203eff4..000000000 --- a/Libtmp/GSL/lib/PDL/Demos/GSL_CDF.pm +++ /dev/null @@ -1,28 +0,0 @@ -package PDL::Demos::GSL_CDF; - -sub info {('gsl_cdf', 'GSL cumulative distribution functions')} - -my @demo = ( -[act => q| -# This demo illustrates the PDL::GSL::CDF module. -# It shows the power of PDL with a concise way to generate a table of -# more extreme small p-values, and the associated Z scores. -use PDL::GSL::CDF; -$pvalue = ipow(pdl(10),-(sequence(32) + 1)); -$z = gsl_cdf_ugaussian_Qinv($pvalue); -$pdl = $pvalue->cat($z)->transpose; -print $pdl->splitdim(1,8)->mv(2,1)->clump(-2)->string("%4.4g"); -|], - -[act => q| -# And more extreme high Z scores, and the associated p-values. -$z = sequence(24) + 1; -$pvalue = gsl_cdf_ugaussian_Q($z); -$pdl = $z->cat($pvalue)->transpose; -print $pdl->splitdim(1,8)->mv(2,1)->clump(-2)->string("%4.4g"); -|], -); - -sub demo { @demo } - -1; diff --git a/Libtmp/GSL/lib/PDL/Demos/GSL_RNG.pm b/Libtmp/GSL/lib/PDL/Demos/GSL_RNG.pm deleted file mode 100644 index 2db5dede5..000000000 --- a/Libtmp/GSL/lib/PDL/Demos/GSL_RNG.pm +++ /dev/null @@ -1,112 +0,0 @@ -package PDL::Demos::GSL_RNG; - -sub info {('gsl_rng', 'GSL randomness functions (Req.: PDL::Graphics::Simple)')} - -my @demo = ( -[act => q| -# This demo illustrates the PDL::GSL::RNG module. -# It shows the power of PDL with a concise way to generate graphs of -# different random number distributions. -use PDL::Graphics::Simple; -use PDL::GSL::RNG; -$x = zeroes(100)->xlinvals(-5,5); -$w = pgswin(); -# Exponential Power Distribution -$w->plot( - with=>'lines', key=>'a=1 b=2.5', $x, ran_exppow_pdf($x, 1, 2.5), - with=>'lines', key=>'a=1 b=0.5', $x, ran_exppow_pdf($x, 1, 0.5), - {le=>'tr', yrange=>[0,0.8], title=>'Exponential Power Distribution', - xlabel=>'x', ylabel=>'p(x)'} -); -|], - -[act => q| -# Cauchy Distribution -$w->plot( - with=>'lines', key=>'a=1', $x, ran_cauchy_pdf($x, 1), - with=>'lines', key=>'a=2', $x, ran_cauchy_pdf($x, 2), - {le=>'tr', yrange=>[0,0.4], title=>'Cauchy Distribution', - xlabel=>'x', ylabel=>'p(x)'} -); -|], - -[act => q| -# Rayleigh Tail Distribution -$x = zeroes(100)->xlinvals(0,5); -$w->plot( - with=>'lines', key=>'a=1 sigma=1', $x, ran_rayleigh_tail_pdf($x, 1, 1), - with=>'lines', key=>'a=0.5 sigma=2', $x, ran_rayleigh_tail_pdf($x, 0.5, 2), - {le=>'tr', yrange=>[0,1.1], title=>'Rayleigh Tail Distribution', - xlabel=>'x', ylabel=>'p(x)'} -); -|], - -[act => q| -# Gamma Distribution -$x = zeroes(100)->xlinvals(0,5); -$w->plot( - with=>'lines', key=>'a=1 b=1', $x, ran_gamma_pdf($x, 1, 1), - with=>'lines', key=>'a=2 b=1', $x, ran_gamma_pdf($x, 2, 1), - with=>'lines', key=>'a=3 b=1', $x, ran_gamma_pdf($x, 3, 1), - {le=>'tr', yrange=>[0,1], title=>'Gamma Distribution', - xlabel=>'x', ylabel=>'p(x)'} -); -|], - -[act => q| -# Bivariate Gaussian Distribution -# inspired by https://www.perlmonks.org/?node_id=11104262 -$points = pdl '[219 88 2.7; 38 95 1.7; 45 268 0.8]'; -($XSIZE, $YSIZE) = (300, 300); -($xcoord, $ycoord, $weight) = $points # xyw nweights - ->slice(",*$XSIZE,*$YSIZE,") # xyw nx ny nweights - ->using(0..2); # nx ny nweights -$xbase = xvals($XSIZE)->slice(",*$YSIZE"); # nx ny -$ybase = xvals($YSIZE)->slice("*$XSIZE,"); # nx ny -for (1..90) { - $h = ( - $weight * ran_bivariate_gaussian_pdf( - $xcoord-$xbase, $ycoord-$ybase, $_, $_, 0 - ) # nx ny nweights - )->mv(-1,0)->sumover; # nx ny - $w->plot(with=>'image', $h, {title=>'Bivariate Gaussian Distribution',j=>1}); -} -|], - -[act => q| -# Same, but with a colourful heatmap (if you have the right libraries) -sub as_heatmap { - my ($d) = @_; - my $max = $d->max; - die "as_heatmap: can't work if max == 0" if $max == 0; - $d /= $max; # negative OK - my $hue = (1 - $d)*240; - $d = cat($hue, pdl(1), pdl(1)); - (hsv_to_rgb($d->mv(-1,0)) * 255)->byte->mv(0,-1); -} -if (eval 'use PDL::Graphics::ColorSpace; 1') { - for (1..90) { - $h = ( - $weight * ran_bivariate_gaussian_pdf( - $xcoord-$xbase, $ycoord-$ybase, $_, $_, 0 - ) # nx ny nweights - )->mv(-1,0)->sumover; # nx ny - $w->plot( - with=>'image', as_heatmap($h), - {title=>'Bivariate Gaussian Distribution (heatmap)',j=>1} - ); - } -} -|], - -[comment => q| -See https://www.gnu.org/software/gsl/doc/html/randist.html for more. -|], -); - -sub demo { @demo } -sub done {' - undef $w; -'} - -1; diff --git a/Libtmp/GSL/lib/PDL/GSL/CDF.pd b/Libtmp/GSL/lib/PDL/GSL/CDF.pd deleted file mode 100644 index 52c85fc3c..000000000 --- a/Libtmp/GSL/lib/PDL/GSL/CDF.pd +++ /dev/null @@ -1,181 +0,0 @@ -use strict; -use warnings; - -our $VERSION = '2.096'; - -pp_addpm({At=>'Top'},<<'EOD'); - -use strict; -use warnings; - -=head1 NAME - -PDL::GSL::CDF - PDL interface to GSL Cumulative Distribution Functions - -=head1 DESCRIPTION - -This is an interface to the Cumulative Distribution Function package present in the GNU Scientific Library. - -Let us have a continuous random number distributions are defined by a probability density function C. - -The cumulative distribution function for the lower tail C is defined by the integral of C, and -gives the probability of a variate taking a value less than C. These functions are named B. - -The cumulative distribution function for the upper tail C is defined by the integral of C, and -gives the probability of a variate taking a value greater than C. These functions are named B. - -The upper and lower cumulative distribution functions are related by C and -satisfy C<0 E= P(x) E= 1> and C<0 E= Q(x) E= 1>. - -The inverse cumulative distributions, C and C give the values of C which correspond -to a specific value of C

or C. They can be used to find confidence limits from probability values. -These functions are named B and B. - -For discrete distributions the probability of sampling the integer value C is given by C, where -C. The cumulative distribution for the lower tail C of a discrete distribution is -defined as, where the sum is over the allowed range of the distribution less than or equal to C. - -The cumulative distribution for the upper tail of a discrete distribution C is defined as giving the sum -of probabilities for all values greater than C. These two definitions satisfy the identity C. - -If the range of the distribution is C<1> to C inclusive then C, C -while C, C. - -=head1 SYNOPSIS - - use PDL; - use PDL::GSL::CDF; - - my $p = gsl_cdf_tdist_P( $t, $df ); - - my $t = gsl_cdf_tdist_Pinv( $p, $df ); - -=cut - -EOD - -pp_addhdr(' -#include -#include - -static char *funcname; - -static void cdf_error_handler(const char *reason, const char *file, int line, int status) { - char buf[200]; - sprintf(buf,"Error in %s: %s", funcname, gsl_strerror(status)); - barf(buf); -} - -'); - -use Text::ParseWords qw(shellwords); -chomp(my $header = `gsl-config --cflags`); -$header =~ s#\\#/#g; # win32 -($header) = map {s/^-I//;$_} grep /^-I/, shellwords $header; -$header .= '/gsl/gsl_cdf.h'; -my $h = do { open my $fh, "<", $header or die "$header: $!"; local $/ = undef; <$fh> }; -my @functions = $h =~ m/(double\s+gsl_cdf_.+?\)\s*;\s*)\n/xmsg; -my @func_defs; - -my %p_type = ( - 'double' => 'double', - 'unsigned int' => 'ulonglong', -); - -my %func_desc = ( - gsl_cdf_gaussian => ['The Gaussian Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Gaussian distribution with standard deviation I.'], - gsl_cdf_ugaussian => ['The Unit Gaussian Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the unit Gaussian distribution.'], - gsl_cdf_exponential => ['The Exponential Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the exponential distribution with mean I.'], - gsl_cdf_laplace => ['The Laplace Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Laplace distribution with width I.'], - gsl_cdf_exppow => ['The Exponential Power Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) for the exponential power distribution with parameters I and I.'], - gsl_cdf_cauchy => ['The Cauchy Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Cauchy distribution with scale parameter I.'], - gsl_cdf_rayleigh => ['The Rayleigh Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Rayleigh distribution with scale parameter I.'], - gsl_cdf_gamma => ['The Gamma Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the gamma distribution with parameters I and I.'], - gsl_cdf_flat => ['The Flat (Uniform) Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for a uniform distribution from I to I.'], - gsl_cdf_lognormal => ['The Lognormal Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the lognormal distribution with parameters I and I.'], - gsl_cdf_chisq => ['The Chi-squared Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the chi-squared distribution with I degrees of freedom.'], - gsl_cdf_fdist => ['The F-distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the F-distribution with I and I degrees of freedom.'], - gsl_cdf_tdist => ['The t-distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the t-distribution with I degrees of freedom.'], - gsl_cdf_beta => ['The Beta Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the beta distribution with parameters I and I.'], - gsl_cdf_logistic => ['The Logistic Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the logistic distribution with scale parameter I.'], - gsl_cdf_pareto => ['The Pareto Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Pareto distribution with exponent I and scale I.'], - gsl_cdf_weibull => ['The Weibull Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Weibull distribution with scale I and exponent I.'], - gsl_cdf_gumbel1 => ['The Type-1 Gumbel Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Type-1 Gumbel distribution with parameters I and I.'], - gsl_cdf_gumbel2 => ['The Type-2 Gumbel Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Type-2 Gumbel distribution with parameters I and I.'], - gsl_cdf_poisson => ['The Poisson Distribution', 'These functions compute the cumulative distribution functions P(k), Q(k) for the Poisson distribution with parameter I.'], - gsl_cdf_binomial => ['The Binomial Distribution', 'These functions compute the cumulative distribution functions P(k), Q(k) for the binomial distribution with parameters I

and I.'], - gsl_cdf_negative_binomial => ['The Negative Binomial Distribution', 'These functions compute the cumulative distribution functions P(k), Q(k) for the negative binomial distribution with parameters I

and I.'], - gsl_cdf_pascal => ['The Pascal Distribution', 'These functions compute the cumulative distribution functions P(k), Q(k) for the Pascal distribution with parameters I

and I.'], - gsl_cdf_geometric => ['The Geometric Distribution', 'These functions compute the cumulative distribution functions P(k), Q(k) for the geometric distribution with parameter I

.'], - gsl_cdf_hypergeometric => ['The Hypergeometric Distribution', 'These functions compute the cumulative distribution functions P(k), Q(k) for the hypergeometric distribution with parameters I, I and I.'], -); - -for (@functions) { - s/\n\s+/ /xmsg; - s/const //g; - if (m/^(\w+)\s+(\w+)\s*\(\s*(.+)\s*\)\s*\;$/s) { - my ($out_type, $function, $pars) = ($1, $2, $3); - my @pars = split /,/, $pars; - for (@pars) { - if (m/^(.+)( \w+)$/) { - my ($type, $par) = ($1, $2); - s/^ | $//g for ($type, $par); - $par = lc $par; - $_ = [$p_type{$type}, $par]; - } - } - push @func_defs, [ $out_type, $function, \@pars ]; - } -} - -@func_defs = sort { $a->[1] cmp $b->[1] } @func_defs; # sort by function name - -for my $f (@func_defs) { - my ($out_type, $function, $pars) = @$f; - my $func_short = join '_', (split '_', $function)[0..2]; - my ($p, $code) = print_ppdef($out_type, $function, @$pars); - my $desc = delete $func_desc{$func_short}; - pp_addpm(qq{\n=head2 $desc->[0] (${func_short}_*)\n\n$desc->[1]\n\n=cut\n\n}) if $desc; - pp_def($function, - HandleBad => 1, - GenericTypes => ['D'], - Pars => $p, - Code => $code, - Doc => '', - ); -} - -sub print_ppdef { - my ($out_type, $function, @pars) = @_; - @pars = map $_->[0] eq 'double' ? [$_->[1]] : $_, @pars; - my $pars = join ' ', map "@$_();", @pars, ['[o]out']; - my $code = pp_line_numbers(__LINE__, <[-1]()", @pars]}); -gsl_set_error_handler(current_handler); -EOF - $code = - "PDL_IF_BAD(if ( !(" . join(' && ', map { "\$ISGOOD($_->[-1]())" } @pars ) . ") ) {" . - "\$SETBAD(out()); }" . - 'else,) {' . $code . "}\n"; - return ($pars, $code); -} - -pp_addpm({At=>'Bot'},<<'EOD'); - -=head1 AUTHOR - -Copyright (C) 2009 Maggie J. Xiong - -The GSL CDF module was written by J. Stover. - -All rights reserved. There is no warranty. You are allowed to redistribute this software / documentation as described in the file COPYING in the PDL distribution. - -=cut - -EOD - -pp_done(); - diff --git a/Libtmp/GSL/lib/PDL/GSL/DIFF-FUNC.c b/Libtmp/GSL/lib/PDL/GSL/DIFF-FUNC.c deleted file mode 100644 index 895d63977..000000000 --- a/Libtmp/GSL/lib/PDL/GSL/DIFF-FUNC.c +++ /dev/null @@ -1,44 +0,0 @@ -#include "EXTERN.h" -#include "perl.h" -#include "pdlcore.h" -#include - -SV* ext_funname; - -void set_funname(SV *fn) { - ext_funname = fn; -} - -double FUNC(double x,void * p){ - - double res; - int count; - - dSP; - SV* funname = ext_funname; /* get function name on the perl side */ - - ENTER; - SAVETMPS; - - PUSHMARK(SP); - - XPUSHs(sv_2mortal(newSVnv(x))); - - PUTBACK; - - count=call_sv(funname,G_SCALAR); - - SPAGAIN; - - if (count!=1) - croak("error calling perl function\n"); - - /* recover output value */ - res = POPn; - - PUTBACK; - FREETMPS; - LEAVE; - - return res; -} diff --git a/Libtmp/GSL/lib/PDL/GSL/DIFF.pd b/Libtmp/GSL/lib/PDL/GSL/DIFF.pd deleted file mode 100644 index 011945db4..000000000 --- a/Libtmp/GSL/lib/PDL/GSL/DIFF.pd +++ /dev/null @@ -1,174 +0,0 @@ -use strict; -use warnings; - -{ no warnings 'once'; # pass info back to Makefile.PL -$PDL::Core::Dev::EXTRAS{$::PDLMOD}{OBJECT} .= join '', map " $::PDLBASE-$_\$(OBJ_EXT)", qw(FUNC); -} - -pp_addpm({At=>'Top'},<<'EOD'); - -use strict; -use warnings; - -=head1 NAME - -PDL::GSL::DIFF - PDL interface to numerical differentiation routines in GSL - -=head1 DESCRIPTION - -This is an interface to the numerical differentiation package present in the -GNU Scientific Library. - -=head1 SYNOPSIS - - use PDL; - use PDL::GSL::DIFF; - - my $x0 = 3.3; - - my @res = gsldiff(\&myfunction,$x0); - # same as above: - @res = gsldiff(\&myfunction,$x0,{Method => 'central'}); - - # use only values greater than $x0 to get the derivative - @res = gsldiff(\&myfunction,$x0,{Method => 'forward'}); - - # use only values smaller than $x0 to get the derivative - @res = gsldiff(\&myfunction,$x0,{Method => 'backward'}); - - sub myfunction{ - my ($x) = @_; - return $x**2; - } - -=cut -EOD - -pp_addpm({At=>'Bot'},<<'EOD'); # the rest of FUNCTIONS section -=head2 gsldiff - -=for ref - -This functions serves as an interface to the three differentiation -functions present in GSL: gsl_diff_central, gsl_diff_backward and -gsl_diff_forward. To compute the derivative, the central method uses -values greater and smaller than the point at which the derivative is -to be evaluated, while backward and forward use only values smaller -and greater respectively. gsldiff() returns both the derivative -and an absolute error estimate. The default method is 'central', -others can be specified by passing an option. - -Please check the GSL documentation for more information. - -=for usage - -Usage: - - ($d,$abserr) = gsldiff($function_ref,$x,{Method => $method}); - -=for example - -Example: - - #derivative using default method ('central') - ($d,$abserr) = gsldiff(\&myf,3.3); - - #same as above with method set explicitly - ($d,$abserr) = gsldiff(\&myf,3.3,{Method => 'central'}); - - #using backward & forward methods - ($d,$abserr) = gsldiff(\&myf,3.3,{Method => 'backward'}); - ($d,$abserr) = gsldiff(\&myf,3.3,{Method => 'forward'}); - - sub myf{ - my ($x) = @_; - return exp($x); - } - -=head1 BUGS - -Feedback is welcome. Log bugs in the PDL bug database (the -database is always linked from L). - -=head1 SEE ALSO - -L - -The GSL documentation is online at - - http://www.gnu.org/software/gsl/manual/ - -=head1 AUTHOR - -This file copyright (C) 2003 Andres Jordan -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. - -The GSL differentiation routines were written by David Morrison. - -=cut - -EOD - -pp_addhdr(' -#include -#include -#include -#include - -double FUNC(double x,void * p); -void set_funname(SV *fn); -'); - -pp_addpm(' -sub gsldiff{ - my $opt; - if (ref($_[$#_]) eq \'HASH\'){ $opt = pop @_; } - else{ $opt = {Method => \'central\'}; } - die \'Usage: gsldiff(function_ref, x, {Options} )\' - if $#_<1 || $#_>2; - my ($f,$x) = @_; - my ($res,$abserr); - if($$opt{Method}=~/cent/i){ - ($res,$abserr) = PDL::GSL::DIFF::diff_central($x,$f); - } - elsif($$opt{Method}=~/back/i){ - ($res,$abserr) = PDL::GSL::DIFF::diff_backward($x,$f); - } - elsif($$opt{Method}=~/forw/i){ - ($res,$abserr) = PDL::GSL::DIFF::diff_forward($x,$f); - } - else{ - barf("Unknown differentiation method $$opt{Method} in gsldiff\n"); - } - return ($res,$abserr); -} - -'); - -pp_add_exported('gsldiff'); - -sub diff_func { - my ($which) = @_; - pp_def("diff_$which", - Pars => 'double x(); double [o] res(); double [o] abserr();', - OtherPars => 'SV* function;', - Doc => undef, - Code => < - -#define PDL PDL_GSL_INTEG -extern Core *PDL; - -#define max_nested_integrals 20 - -static SV* ext_funname[max_nested_integrals]; -static int current_fun = -1; - -void set_funname(SV *fn) { - current_fun++; - if (current_fun > (max_nested_integrals)) barf("Too many nested integrals, sorry!\n"); - ext_funname[current_fun] = fn; -} -void dec_func() { - current_fun--; -} - -void my_handler (const char * reason, - const char * file, - int line, - int gsl_errno){ - printf("Warning: %s at line %d of GSL file %s\n",reason,line,file); -} - -double FUNC(double x,void * p){ - - SV* funname; - - int count; - - I32 ax ; - - double res; - double* resp; - - dSP; - - resp = &res; - - - ENTER; - SAVETMPS; - - /* get function name on the perl side */ - funname = ext_funname[current_fun]; - - PUSHMARK(SP); - - XPUSHs(sv_2mortal(newSVnv(x))); - - PUTBACK; - - count=call_sv(funname,G_SCALAR); - - SPAGAIN; - SP -= count ; - ax = (SP - PL_stack_base) + 1 ; - - if (count!=1) - croak("error calling perl function\n"); - - /* recover output value */ - /*res = POPn;*/ - *resp = SvNV(ST(0)); - - PUTBACK; - FREETMPS; - LEAVE; - - return res; -} diff --git a/Libtmp/GSL/lib/PDL/GSL/INTEG.pd b/Libtmp/GSL/lib/PDL/GSL/INTEG.pd deleted file mode 100644 index e2a01f5b8..000000000 --- a/Libtmp/GSL/lib/PDL/GSL/INTEG.pd +++ /dev/null @@ -1,953 +0,0 @@ -use strict; -use warnings; - -{ no warnings 'once'; # pass info back to Makefile.PL -$PDL::Core::Dev::EXTRAS{$::PDLMOD}{OBJECT} .= join '', map " $::PDLBASE-$_\$(OBJ_EXT)", qw(FUNC); -} - -pp_bless('PDL::GSL::INTEG'); - -pp_addpm({At=>'Top'},<<'EOD'); -use strict; -use warnings; - -=head1 NAME - -PDL::GSL::INTEG - PDL interface to numerical integration routines in GSL - -=head1 DESCRIPTION - -This is an interface to the numerical integration package present in the -GNU Scientific Library, which is an implementation of QUADPACK. - -Functions are named B where {algorithm} -is the QUADPACK naming convention. The available functions are: - -=over 3 - -=item gslinteg_qng: Non-adaptive Gauss-Kronrod integration - -=item gslinteg_qag: Adaptive integration - -=item gslinteg_qags: Adaptive integration with singularities - -=item gslinteg_qagp: Adaptive integration with known singular points - -=item gslinteg_qagi: Adaptive integration on infinite interval of the form (-\infty,\infty) - -=item gslinteg_qagiu: Adaptive integration on infinite interval of the form (la,\infty) - -=item gslinteg_qagil: Adaptive integration on infinite interval of the form (-\infty,lb) - -=item gslinteg_qawc: Adaptive integration for Cauchy principal values - -=item gslinteg_qaws: Adaptive integration for singular functions - -=item gslinteg_qawo: Adaptive integration for oscillatory functions - -=item gslinteg_qawf: Adaptive integration for Fourier integrals - -=back - -Each algorithm computes an approximation to the integral, I, -of the function f(x)w(x), where w(x) is a weight function -(for general integrands w(x)=1). The user provides absolute -and relative error bounds (epsabs,epsrel) which specify -the following accuracy requirement: - -|RESULT - I| <= max(epsabs, epsrel |I|) - - -The routines will fail to converge if the -error bounds are too stringent, but always return the best -approximation obtained up to that stage - -All functions return the result, and estimate of the absolute -error and an error flag (which is zero if there were no problems). -You are responsible for checking for any errors, no warnings are issued -unless the option {Warn => 'y'} is specified in which case -the reason of failure will be printed. - -You can nest integrals up to 20 levels. If you find yourself in -the unlikely situation that you need more, you can change the value -of 'max_nested_integrals' in the first line of the file 'FUNC.c' -and recompile. - -=head1 NOMENCLATURE - -Throughout this documentation we strive to use the same variables that -are present in the original GSL documentation (see L). Oftentimes those variables are called C and -C. Since good Perl coding practices discourage the use of Perl -variables C<$a> and C<$b>, here we refer to Parameters C and C -as C<$pa> and C<$pb>, respectively, and Limits (of domain or -integration) as C<$la> and C<$lb>. - -=head1 SYNOPSIS - - use PDL; - use PDL::GSL::INTEG; - - my $la = 1.2; - my $lb = 3.7; - my $epsrel = 0; - my $epsabs = 1e-6; - - # Non adaptive integration - my ($res,$abserr,$ierr,$neval) = gslinteg_qng(\&myf,$la,$lb,$epsrel,$epsabs); - # Warnings on - my ($res,$abserr,$ierr,$neval) = gslinteg_qng(\&myf,$la,$lb,$epsrel,$epsabs,{Warn=>'y'}); - - # Adaptive integration with warnings on - my $limit = 1000; - my $key = 5; - my ($res,$abserr,$ierr) = gslinteg_qag(\&myf,$la,$lb,$epsrel, - $epsabs,$limit,$key,{Warn=>'y'}); - - sub myf{ - my ($x) = @_; - return exp(-$x**2); - } - -EOD - -pp_addpm({At=>'Bot'},<<'EOD'); -=head1 BUGS - -Feedback is welcome. Log bugs in the PDL bug database (the -database is always linked from L). - -=head1 SEE ALSO - -L - -The GSL documentation for numerical integration is online at -L - -=head1 AUTHOR - -This file copyright (C) 2003,2005 Andres Jordan -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. - -The GSL integration routines were written by Brian Gough. QUADPACK -was written by Piessens, Doncker-Kapenga, Uberhuber and Kahaner. - -=cut - - -EOD - -pp_addhdr(' -#include -#include -#include -#include - -double FUNC(double x,void * p); -void set_funname(SV *fn); -void dec_func(); -void my_handler (const char * reason, const char * file, int line, int gsl_errno); - -'); - -pp_def('gslinteg_qng', - Pars => 'a(); b(); epsabs(); epsrel(); int gslwarn(); - [o] result(); [o] abserr(); int [o] neval(); int [o] ierr();', - OtherPars => 'SV* function;', - GenericTypes => ['D'], - PMCode => <<'EOF', -sub gslinteg_qng{ - my $opt = ref($_[-1]) eq 'HASH' ? pop @_ : {Warn => 'n'}; - barf 'Usage: gslinteg_qng($function_ref,$la,$lb,$epsabs,$epsrel,[opt])' - unless (@_ == 5); - my $warn = $$opt{Warn}=~/y/i ? 1 : 0; - my ($f,$la,$lb,$epsabs,$epsrel) = @_; - $_=PDL->null for my ($res,$abserr,$neval,$ierr); - _gslinteg_qng_int($la,$lb,$epsabs,$epsrel,$warn,$res,$abserr,$neval,$ierr,$f); - return ($res,$abserr,$ierr,$neval); -} -EOF - Code => <<'EOF', -gsl_error_handler_t * old_handler; -if ($gslwarn() == 1) { old_handler = gsl_set_error_handler(&my_handler); } -else { gsl_set_error_handler_off ();} -set_funname($COMP(function)); -gsl_function F; -F.function = &FUNC; -F.params = 0; -$ierr() = gsl_integration_qng(&F,$a(),$b(),$epsabs(),$epsrel(),$P(result),$P(abserr),(size_t *) $P(neval)); -dec_func(); -EOF - Doc => <<'EOF', -=for ref - -Non-adaptive Gauss-Kronrod integration - -This function applies the Gauss-Kronrod 10-point, 21-point, 43-point -and 87-point integration rules in succession until an estimate of the -integral of f over ($la,$lb) is achieved within the desired absolute -and relative error limits, $epsabs and $epsrel. It is meant for fast -integration of smooth functions. It returns an array with the result, -an estimate of the absolute error, an error flag and the number of -function evaluations performed. - -=for usage - -Usage: - - ($res,$abserr,$ierr,$neval) = gslinteg_qng($function_ref,$la,$lb, - $epsrel,$epsabs,[{Warn => $warn}]); - -=for example - -Example: - - my ($res,$abserr,$ierr,$neval) = gslinteg_qng(\&f,0,1,0,1e-9); - # with warnings on - my ($res,$abserr,$ierr,$neval) = gslinteg_qng(\&f,0,1,0,1e-9,{Warn => 'y'}); - - sub f{ - my ($x) = @_; - return ($x**2.6)*log(1.0/$x); - } -EOF -); - -pp_def('gslinteg_qag', - Pars => 'a(); b(); epsabs();epsrel(); - int limit(); int key(); int n(); int gslwarn(); - [o] result(); [o] abserr(); int [o] ierr();', - GenericTypes => ['D'], - OtherPars => 'SV* function;', - Code => <<'EOF', -gsl_error_handler_t * old_handler; -if ($gslwarn() == 1) { old_handler = gsl_set_error_handler(&my_handler); } -else { gsl_set_error_handler_off ();} -{gsl_integration_workspace *w; -set_funname($COMP(function)); -gsl_function F; -F.function = &FUNC; -F.params = 0; -w = gsl_integration_workspace_alloc((size_t) $n()); -$ierr() = gsl_integration_qag(&F,$a(),$b(),$epsabs(),$epsrel(),(size_t) $limit(),$key(),w,$P(result),$P(abserr)); -gsl_integration_workspace_free(w); -dec_func(); -} -EOF - PMCode => <<'EOF', -sub gslinteg_qag { - my $opt = ref($_[-1]) eq 'HASH' ? pop @_ : {Warn => 'n'}; - my $warn = $$opt{Warn}=~/y/i ? 1 : 0; - my ($f,$la,$lb,$epsabs,$epsrel,$limit,$key) = @_; - barf 'Usage: gslinteg_qag($function_ref,$la,$lb,$epsabs,$epsrel,$limit,$key,[opt]) ' - unless ($#_ == 6); - $_ = PDL->null for my ($res,$abserr,$ierr); - _gslinteg_qag_int($la,$lb,$epsabs,$epsrel,$limit,$key,$limit,$warn,$res,$abserr,$ierr,$f); - return ($res,$abserr,$ierr); -} -EOF - Doc => <<'EOF', -=for ref - -Adaptive integration - -This function applies an integration rule adaptively until an estimate of -the integral of f over ($la,$lb) is achieved within the desired absolute and -relative error limits, $epsabs and $epsrel. On each iteration the adaptive -integration strategy bisects the interval with the largest error estimate; -the maximum number of allowed subdivisions is given by the parameter $limit. -The integration rule is determined by the -value of $key, which has to be one of (1,2,3,4,5,6) and correspond to -the 15, 21, 31, 41, 51 and 61 point Gauss-Kronrod rules respectively. -It returns an array with the result, an estimate of the absolute error -and an error flag. - -=for usage - -Usage: - - ($res,$abserr,$ierr) = gslinteg_qag($function_ref,$la,$lb,$epsrel, - $epsabs,$limit,$key,[{Warn => $warn}]); - -=for example - -Example: - - my ($res,$abserr,$ierr) = gslinteg_qag(\&f,0,1,0,1e-10,1000,1); - # with warnings on - my ($res,$abserr,$ierr) = gslinteg_qag(\&f,0,1,0,1e-10,1000,1,{Warn => 'y'}); - - sub f{ - my ($x) = @_; - return ($x**2.6)*log(1.0/$x); - } -EOF -); - -pp_def('gslinteg_qags', - Pars => 'a(); b(); epsabs(); epsrel(); - int limit(); int n(); int gslwarn(); - [o] result(); [o] abserr(); int [o] ierr();', - GenericTypes => ['D'], - OtherPars => 'SV* function;', - PMCode => <<'EOF', -sub gslinteg_qags{ - my $opt = ref($_[-1]) eq 'HASH' ? pop @_ : {Warn => 'n'}; - my $warn = $$opt{Warn}=~/y/i ? 1 : 0; - my ($f,$la,$lb,$epsabs,$epsrel,$limit) = @_; - barf 'Usage: gslinteg_qags($function_ref,$la,$lb,$epsabs,$epsrel,$limit,[opt]) ' - unless ($#_ == 5); - $_ = PDL->null for my ($res,$abserr,$ierr); - _gslinteg_qags_int($la,$lb,$epsabs,$epsrel,$limit,$limit,$warn,$res,$abserr,$ierr,$f); - return ($res,$abserr,$ierr); -} -EOF - Code => <<'EOF', -gsl_error_handler_t * old_handler; -if ($gslwarn() == 1) { old_handler = gsl_set_error_handler(&my_handler); } -else { gsl_set_error_handler_off ();} -{gsl_integration_workspace *w; -set_funname($COMP(function)); -gsl_function F; -F.function = &FUNC; -F.params = 0; -w = gsl_integration_workspace_alloc((size_t) $n()); -$ierr() = gsl_integration_qags(&F,$a(),$b(),$epsabs(),$epsrel(),(size_t) $limit(),w,$P(result),$P(abserr)); -gsl_integration_workspace_free(w); -dec_func(); -} -EOF - Doc => <<'EOF', -=for ref - -Adaptive integration with singularities - -This function applies the Gauss-Kronrod 21-point integration rule -adaptively until an estimate of the integral of f over ($la,$lb) is -achieved within the desired absolute and relative error limits, -$epsabs and $epsrel. The algorithm is such that it -accelerates the convergence of the integral in the presence of -discontinuities and integrable singularities. -The maximum number of allowed subdivisions done by the adaptive -algorithm must be supplied in the parameter $limit. - -=for usage - -Usage: - - ($res,$abserr,$ierr) = gslinteg_qags($function_ref,$la,$lb,$epsrel, - $epsabs,$limit,[{Warn => $warn}]); - -=for example - -Example: - - my ($res,$abserr,$ierr) = gslinteg_qags(\&f,0,1,0,1e-10,1000); - # with warnings on - ($res,$abserr,$ierr) = gslinteg_qags(\&f,0,1,0,1e-10,1000,{Warn => 'y'}); - - sub f{ - my ($x) = @_; - return ($x)*log(1.0/$x); - } -EOF -); - -pp_def('gslinteg_qagp', - Pars => 'pts(l); epsabs(); epsrel();int limit();int n();int gslwarn(); - [o] result(); [o] abserr();int [o] ierr();', - GenericTypes => ['D'], - OtherPars => 'SV* function;', - Code => <<'EOF', -gsl_error_handler_t * old_handler; -if ($gslwarn() == 1) { old_handler = gsl_set_error_handler(&my_handler); } -else { gsl_set_error_handler_off ();} -{gsl_integration_workspace *w; -set_funname($COMP(function)); -gsl_function F; -F.function = &FUNC; -F.params = 0; -w = gsl_integration_workspace_alloc((size_t) $n()); -$ierr() = gsl_integration_qagp(&F,$P(pts),(size_t) $SIZE(l),$epsabs(),$epsrel(),(size_t) $limit(),w,$P(result),$P(abserr)); -gsl_integration_workspace_free(w); -dec_func(); -} -EOF - PMCode => <<'EOF', -sub gslinteg_qagp{ - my $opt = ref($_[-1]) eq 'HASH' ? pop @_ : {Warn => 'n'}; - my $warn = $$opt{Warn}=~/y/i ? 1 : 0; - my ($f,$points,$epsabs,$epsrel,$limit) = @_; - barf 'Usage: gslinteg_qagp($function_ref,$points,$epsabs,$epsrel,$limit,[opt]) ' - unless ($#_ == 4); - $_ = PDL->null for my ($res,$abserr,$ierr); - _gslinteg_qagp_int($points,$epsabs,$epsrel,$limit,$limit,$warn,$res,$abserr,$ierr,$f); - return ($res,$abserr,$ierr); -} -EOF - Doc => <<'EOF', -=for ref - -Adaptive integration with known singular points - -This function applies the adaptive integration algorithm used by -gslinteg_qags taking into account the location of singular points -until an estimate of -the integral of f over ($la,$lb) is achieved within the desired absolute and -relative error limits, $epsabs and $epsrel. -Singular points are supplied in the ndarray $points, whose endpoints -determine the integration range. -So, for example, if the function has singular points at x_1 and x_2 and the -integral is desired from a to b (a < x_1 < x_2 < b), $points = pdl(a,x_1,x_2,b). -The maximum number of allowed subdivisions done by the adaptive -algorithm must be supplied in the parameter $limit. - -=for usage - -Usage: - - ($res,$abserr,$ierr) = gslinteg_qagp($function_ref,$points,$epsabs, - $epsrel,$limit,[{Warn => $warn}]) - -=for example - -Example: - - my $points = pdl(0,1,sqrt(2),3); - my ($res,$abserr,$ierr) = gslinteg_qagp(\&f,$points,0,1e-3,1000); - # with warnings on - ($res,$abserr,$ierr) = gslinteg_qagp(\&f,$points,0,1e-3,1000,{Warn => 'y'}); - - sub f{ - my ($x) = @_; - my $x2 = $x**2; - my $x3 = $x**3; - return $x3 * log(abs(($x2-1.0)*($x2-2.0))); - } - -EOF -); - -pp_def('gslinteg_qagi', - Pars => 'epsabs(); epsrel(); int limit(); int n();int gslwarn(); - [o] result(); [o] abserr(); int [o] ierr();', - GenericTypes => ['D'], - OtherPars => 'SV* function;', - Code => <<'EOF', -gsl_error_handler_t * old_handler; -if ($gslwarn() == 1) { old_handler = gsl_set_error_handler(&my_handler); } -else { gsl_set_error_handler_off ();} -{gsl_integration_workspace *w; -set_funname($COMP(function)); -gsl_function F; -F.function = &FUNC; -F.params = 0; -w = gsl_integration_workspace_alloc((size_t) $n()); -$ierr() = gsl_integration_qagi(&F,$epsabs(),$epsrel(),(size_t) $limit(),w,$P(result),$P(abserr)); -gsl_integration_workspace_free(w); -dec_func(); -} -EOF - PMCode => <<'EOF', -sub gslinteg_qagi { - my $opt = ref($_[-1]) eq 'HASH' ? pop @_ : {Warn => 'n'}; - my $warn = $$opt{Warn}=~/y/i ? 1 : 0; - my ($f,$epsabs,$epsrel,$limit) = @_; - barf 'Usage: gslinteg_qagi($function_ref,$epsabs,$epsrel,$limit,[opt])' - unless ($#_ == 3); - $_ = PDL->null for my ($res,$abserr,$ierr); - _gslinteg_qagi_int($epsabs,$epsrel,$limit,$limit,$warn,$res,$abserr,$ierr,$f); - return ($res,$abserr,$ierr); -} -EOF - Doc => <<'EOF', -=for ref - -Adaptive integration on infinite interval - -This function estimates the integral of the function f over the -infinite interval (-\infty,+\infty) within the desired absolute and -relative error limits, $epsabs and $epsrel. -After a transformation, the algorithm -of gslinteg_qags with a 15-point Gauss-Kronrod rule is used. -The maximum number of allowed subdivisions done by the adaptive -algorithm must be supplied in the parameter $limit. - -=for usage - -Usage: - - ($res,$abserr,$ierr) = gslinteg_qagi($function_ref,$epsabs, - $epsrel,$limit,[{Warn => $warn}]); - -=for example - -Example: - - my ($res,$abserr,$ierr) = gslinteg_qagi(\&myfn,1e-7,0,1000); - # with warnings on - ($res,$abserr,$ierr) = gslinteg_qagi(\&myfn,1e-7,0,1000,{Warn => 'y'}); - - sub myfn{ - my ($x) = @_; - return exp(-$x - $x*$x) ; - } -EOF -); - -pp_def('gslinteg_qagiu', - Pars => 'a(); epsabs(); epsrel();int limit();int n();int gslwarn(); - [o] result(); [o] abserr();int [o] ierr();', - GenericTypes => ['D'], - OtherPars => 'SV* function;', - Doc => <<'EOF', -=for ref - -Adaptive integration on infinite interval - -This function estimates the integral of the function f over the -infinite interval (la,+\infty) within the desired absolute and -relative error limits, $epsabs and $epsrel. -After a transformation, the algorithm -of gslinteg_qags with a 15-point Gauss-Kronrod rule is used. -The maximum number of allowed subdivisions done by the adaptive -algorithm must be supplied in the parameter $limit. - -=for usage - -Usage: - - ($res,$abserr,$ierr) = gslinteg_qagiu($function_ref,$la,$epsabs, - $epsrel,$limit,[{Warn => $warn}]); - -=for example - -Example: - - my $alfa = 1; - my ($res,$abserr,$ierr) = gslinteg_qagiu(\&f,99.9,1e-7,0,1000); - # with warnings on - ($res,$abserr,$ierr) = gslinteg_qagiu(\&f,99.9,1e-7,0,1000,{Warn => 'y'}); - - sub f{ - my ($x) = @_; - if (($x==0) && ($alfa == 1)) {return 1;} - if (($x==0) && ($alfa > 1)) {return 0;} - return ($x**($alfa-1))/((1+10*$x)**2); - } -EOF - PMCode => <<'EOF', -sub gslinteg_qagiu{ - my $opt = ref($_[-1]) eq 'HASH' ? pop @_ : {Warn => 'n'}; - my $warn = $$opt{Warn}=~/y/i ? 1 : 0; - my ($f,$la,$epsabs,$epsrel,$limit) = @_; - barf 'Usage: gslinteg_qagiu($function_ref,$la,$epsabs,$epsrel,$limit,[opt])' - unless ($#_ == 4); - $_ = PDL->null for my ($res,$abserr,$ierr); - _gslinteg_qagiu_int($la,$epsabs,$epsrel,$limit,$limit,$warn,$res,$abserr,$ierr,$f); - return ($res,$abserr,$ierr); -} -EOF - Code =>' -gsl_error_handler_t * old_handler; -if ($gslwarn() == 1) { old_handler = gsl_set_error_handler(&my_handler); } -else { gsl_set_error_handler_off ();} -{gsl_integration_workspace *w; -set_funname($COMP(function)); -gsl_function F; -F.function = &FUNC; -F.params = 0; -w = gsl_integration_workspace_alloc((size_t) $n()); -$ierr() = gsl_integration_qagiu(&F,$a(),$epsabs(),$epsrel(),(size_t) $limit(),w,$P(result),$P(abserr)); -gsl_integration_workspace_free(w); -dec_func(); -} -'); - -pp_def('gslinteg_qagil', - Pars => 'b(); epsabs(); epsrel();int limit();int n();int gslwarn(); - [o] result(); [o] abserr();int [o] ierr();', - GenericTypes => ['D'], - OtherPars => 'SV* function;', - Doc => <<'EOF', -=for ref - -Adaptive integration on infinite interval - -This function estimates the integral of the function f over the -infinite interval (-\infty,lb) within the desired absolute and -relative error limits, $epsabs and $epsrel. -After a transformation, the algorithm -of gslinteg_qags with a 15-point Gauss-Kronrod rule is used. -The maximum number of allowed subdivisions done by the adaptive -algorithm must be supplied in the parameter $limit. - -=for usage - -Usage: - - ($res,$abserr,$ierr) = gslinteg_qagl($function_ref,$lb,$epsabs, - $epsrel,$limit,[{Warn => $warn}]); - -=for example - -Example: - - my ($res,$abserr,$ierr) = gslinteg_qagil(\&myfn,1.0,1e-7,0,1000); - # with warnings on - ($res,$abserr,$ierr) = gslinteg_qagil(\&myfn,1.0,1e-7,0,1000,{Warn => 'y'}); - - sub myfn{ - my ($x) = @_; - return exp($x); - } -EOF - PMCode => <<'EOF', -sub gslinteg_qagil{ - my $opt = ref($_[-1]) eq 'HASH' ? pop @_ : {Warn => 'n'}; - my $warn = $$opt{Warn}=~/y/i ? 1 : 0; - my ($f,$lb,$epsabs,$epsrel,$limit) = @_; - barf 'Usage: gslinteg_qagil($function_ref,$lb,$epsabs,$epsrel,$limit,[opt])' - unless ($#_ == 4); - $_ = PDL->null for my ($res,$abserr,$ierr); - _gslinteg_qagil_int($lb,$epsabs,$epsrel,$limit,$limit,$warn,$res,$abserr,$ierr,$f); - return ($res,$abserr,$ierr); -} -EOF - Code =>' -gsl_error_handler_t * old_handler; -if ($gslwarn() == 1) { old_handler = gsl_set_error_handler(&my_handler); } -else { gsl_set_error_handler_off ();} -{gsl_integration_workspace *w; -set_funname($COMP(function)); -gsl_function F; -F.function = &FUNC; -F.params = 0; -w = gsl_integration_workspace_alloc((size_t) $n()); -$ierr() = gsl_integration_qagil(&F,$b(),$epsabs(),$epsrel(),(size_t) $limit(),w,$P(result),$P(abserr)); -gsl_integration_workspace_free(w); -dec_func(); -} -'); - -pp_def('gslinteg_qawc', - Pars => 'a(); b(); c(); epsabs(); epsrel();int limit();int n();int gslwarn(); - [o] result(); [o] abserr();int [o] ierr();', - GenericTypes => ['D'], - OtherPars => 'SV* function;', - PMCode => <<'EOF', -sub gslinteg_qawc{ - my $opt = ref($_[-1]) eq 'HASH' ? pop @_ : {Warn => 'n'}; - my $warn = $$opt{Warn}=~/y/i ? 1 : 0; - my ($f,$la,$lb,$c,$epsabs,$epsrel,$limit) = @_; - barf 'Usage: gslinteg_qawc($function_ref,$la,$lb,$c,$epsabs,$epsrel,$limit,[opt])' - unless ($#_ == 6); - $_ = PDL->null for my ($res,$abserr,$ierr); - _gslinteg_qawc_int($la,$lb,$c,$epsabs,$epsrel,$limit,$limit,$warn,$res,$abserr,$ierr,$f); - return ($res,$abserr,$ierr); -} -EOF - Doc => <<'EOF', -=for ref - -Adaptive integration for Cauchy principal values - -This function computes the Cauchy principal value of the integral of f over (la,lb), -with a singularity at c, I = \int_{la}^{lb} dx f(x)/(x - c). The integral is -estimated within the desired absolute and relative error limits, $epsabs and $epsrel. -The maximum number of allowed subdivisions done by the adaptive -algorithm must be supplied in the parameter $limit. - -=for usage - -Usage: - - ($res,$abserr,$ierr) = gslinteg_qawc($function_ref,$la,$lb,$c,$epsabs,$epsrel,$limit) - -=for example - -Example: - - my ($res,$abserr,$ierr) = gslinteg_qawc(\&f,-1,5,0,0,1e-3,1000); - # with warnings on - ($res,$abserr,$ierr) = gslinteg_qawc(\&f,-1,5,0,0,1e-3,1000,{Warn => 'y'}); - - sub f{ - my ($x) = @_; - return 1.0 / (5.0 * $x * $x * $x + 6.0) ; - } -EOF - Code =>' -gsl_error_handler_t * old_handler; -if ($gslwarn() == 1) { old_handler = gsl_set_error_handler(&my_handler); } -else { gsl_set_error_handler_off ();} -{gsl_integration_workspace *w; -set_funname($COMP(function)); -gsl_function F; -F.function = &FUNC; -F.params = 0; -w = gsl_integration_workspace_alloc((size_t) $n()); -$ierr() = gsl_integration_qawc(&F,$a(),$b(),$c(),$epsabs(),$epsrel(),(size_t) $limit(),w,$P(result),$P(abserr)); -gsl_integration_workspace_free(w); -dec_func(); -} -'); - -pp_def('gslinteg_qaws', - Pars => 'a(); b(); epsabs(); epsrel();int limit(); - int n(); alpha(); beta(); int mu(); int nu(); int gslwarn(); - [o] result(); [o] abserr(); int [o] ierr();', - GenericTypes => ['D'], - OtherPars => 'SV* function;', - PMCode => <<'EOF', -sub gslinteg_qaws{ - my $opt = ref($_[-1]) eq 'HASH' ? pop @_ : {Warn => 'n'}; - my $warn = $$opt{Warn}=~/y/i ? 1 : 0; - my ($f,$alpha,$beta,$mu,$nu,$la,$lb,$epsabs,$epsrel,$limit) = @_; - barf 'Usage: gslinteg_qaws($function_ref,$alpha,$beta,$mu,$nu,$la,$lb,$epsabs,$epsrel,$limit,[opt])' - unless ($#_ == 9); - $_ = PDL->null for my ($res,$abserr,$ierr); - _gslinteg_qaws_int($la,$lb,$epsabs,$epsrel,$limit,$limit,$alpha,$beta,$mu,$nu,$warn,$res,$abserr,$ierr,$f); - return ($res,$abserr,$ierr); -} -EOF - Doc => <<'EOF', -=for ref - -Adaptive integration for singular functions - -The algorithm in gslinteg_qaws is designed for integrands with algebraic-logarithmic -singularities at the end-points of an integration region. -Specifically, this function computes the integral given by -I = \int_{la}^{lb} dx f(x) (x-la)^alpha (lb-x)^beta log^mu (x-la) log^nu (lb-x). -The integral is -estimated within the desired absolute and relative error limits, $epsabs and $epsrel. -The maximum number of allowed subdivisions done by the adaptive -algorithm must be supplied in the parameter $limit. - -=for usage - -Usage: - - ($res,$abserr,$ierr) = - gslinteg_qaws($function_ref,$alpha,$beta,$mu,$nu,$la,$lb, - $epsabs,$epsrel,$limit,[{Warn => $warn}]); - -=for example - -Example: - - my ($res,$abserr,$ierr) = gslinteg_qaws(\&f,0,0,1,0,0,1,0,1e-7,1000); - # with warnings on - ($res,$abserr,$ierr) = gslinteg_qaws(\&f,0,0,1,0,0,1,0,1e-7,1000,{Warn => 'y'}); - - sub f{ - my ($x) = @_; - if($x==0){return 0;} - else{ - my $u = log($x); - my $v = 1 + $u*$u; - return 1.0/($v*$v); - } - } -EOF - Code => ' -gsl_error_handler_t * old_handler; -if ($gslwarn() == 1) { old_handler = gsl_set_error_handler(&my_handler); } -else { gsl_set_error_handler_off ();} -{gsl_integration_qaws_table * qtab; -gsl_integration_workspace *w; - -qtab = gsl_integration_qaws_table_alloc($alpha(),$beta(),$mu(),$nu()); - -set_funname($COMP(function)); -gsl_function F; -F.function = &FUNC; -F.params = 0; -w = gsl_integration_workspace_alloc((size_t) $n()); -$ierr() = gsl_integration_qaws(&F,$a(),$b(),qtab,$epsabs(),$epsrel(),(size_t) $limit(),w,$P(result),$P(abserr)); -gsl_integration_workspace_free(w); -gsl_integration_qaws_table_free(qtab); -dec_func(); -} -'); - -pp_def('gslinteg_qawo', - Pars => 'a(); b(); epsabs(); epsrel();int limit();int n(); - int sincosopt(); omega(); L(); int nlevels();int gslwarn(); - [o] result(); [o] abserr();int [o] ierr();', - GenericTypes => ['D'], - OtherPars => 'SV* function;', - PMCode => <<'EOF', -sub gslinteg_qawo{ - my $opt = ref($_[-1]) eq 'HASH' ? pop @_ : {Warn => 'n'}; - my $warn = $$opt{Warn}=~/y/i ? 1 : 0; - my ($f,$omega,$sincosopt,$la,$lb,$epsabs,$epsrel,$limit) = @_; - barf 'Usage: gslinteg_qawo($function_ref,$omega,$sin_or_cos,$la,$lb,$epsabs,$epsrel,$limit,[opt])' - unless ($#_ == 7); - my $OPTION_SIN_COS; - if($sincosopt=~/cos/i){ $OPTION_SIN_COS = 0;} - elsif($sincosopt=~/sin/i){ $OPTION_SIN_COS = 1;} - else { barf("Error in argument 3 of function gslinteg_qawo: specify 'cos' or 'sin'\n");} - - my $L = $lb - $la; - my $nlevels = $limit; - $_ = PDL->null for my ($res,$abserr,$ierr); - _gslinteg_qawo_int($la,$lb,$epsabs,$epsrel,$limit,$limit,$OPTION_SIN_COS,$omega,$L,$nlevels,$warn,$res,$abserr,$ierr,$f); - return ($res,$abserr,$ierr); -} -EOF - Doc => <<'EOF', -=for ref - -Adaptive integration for oscillatory functions - -This function uses an adaptive algorithm to compute the integral of f over -(la,lb) with the weight function sin(omega*x) or cos(omega*x) -- which of -sine or cosine is used is determined by the parameter $opt ('cos' or 'sin'). -The integral is -estimated within the desired absolute and relative error limits, $epsabs and $epsrel. -The maximum number of allowed subdivisions done by the adaptive -algorithm must be supplied in the parameter $limit. - -=for usage - -Usage: - - ($res,$abserr,$ierr) = gslinteg_qawo($function_ref,$omega,$sin_or_cos, - $la,$lb,$epsabs,$epsrel,$limit,[opt]) - -=for example - -Example: - - my $PI = 3.14159265358979323846264338328; - my ($res,$abserr,$ierr) = PDL::GSL::INTEG::gslinteg_qawo(\&f,10*$PI,'sin',0,1,0,1e-7,1000); - # with warnings on - ($res,$abserr,$ierr) = PDL::GSL::INTEG::gslinteg_qawo(\&f,10*$PI,'sin',0,1,0,1e-7,1000,{Warn => 'y'}); - - sub f{ - my ($x) = @_; - if($x==0){return 0;} - else{ return log($x);} - } -EOF - Code => ' -gsl_error_handler_t * old_handler; -if ($gslwarn() == 1) { old_handler = gsl_set_error_handler(&my_handler); } -else { gsl_set_error_handler_off ();} -{gsl_integration_qawo_table * qtab; -gsl_integration_workspace *w; -enum gsl_integration_qawo_enum T; - -T = GSL_INTEG_SINE; -if ($sincosopt() == 0){ T = GSL_INTEG_COSINE ;} - -qtab = gsl_integration_qawo_table_alloc($omega(),$L(),T,(size_t) $nlevels()); - -set_funname($COMP(function)); -gsl_function F; -F.function = &FUNC; -F.params = 0; -w = gsl_integration_workspace_alloc((size_t) $n()); -$ierr() = gsl_integration_qawo(&F,$a(),$epsabs(),$epsrel(),(size_t) $limit(),w,qtab,$P(result),$P(abserr)); -gsl_integration_workspace_free(w); -gsl_integration_qawo_table_free(qtab); -dec_func(); -} -'); - -pp_def('gslinteg_qawf', - Pars => 'a(); epsabs();int limit();int n(); - int sincosopt(); omega(); int nlevels();int gslwarn(); - [o] result(); [o] abserr();int [o] ierr();', - GenericTypes => ['D'], - OtherPars => 'SV* function;', - PMCode => <<'EOF', -sub gslinteg_qawf{ - my $opt = ref($_[-1]) eq 'HASH' ? pop @_ : {Warn => 'n'}; - my $warn = $$opt{Warn}=~/y/i ? 1 : 0; - my ($f,$omega,$sincosopt,$la,$epsabs,$limit) = @_; - barf 'Usage: gslinteg_qawf($function_ref,$omega,$sin_or_cos,$la,$epsabs,$limit,[opt])' - unless ($#_ == 5); - my $OPTION_SIN_COS; - if($sincosopt=~/cos/i){ $OPTION_SIN_COS = 0;} - elsif($sincosopt=~/sin/i){ $OPTION_SIN_COS = 1;} - else { barf("Error in argument 3 of function gslinteg_qawf: specify 'cos' or 'sin'\n");} - my $nlevels = $limit; - $_ = PDL->null for my ($res,$abserr,$ierr); - _gslinteg_qawf_int($la,$epsabs,$limit,$limit,$OPTION_SIN_COS,$omega,$nlevels,$warn,$res,$abserr,$ierr,$f); - return ($res,$abserr,$ierr); -} -EOF - Doc => <<'EOF', -=for ref - -Adaptive integration for Fourier integrals - -This function attempts to compute a Fourier integral of the function -f over the semi-infinite interval [la,+\infty). Specifically, it attempts -tp compute I = \int_{la}^{+\infty} dx f(x)w(x), where w(x) is sin(omega*x) -or cos(omega*x) -- which of sine or cosine is used is determined by -the parameter $opt ('cos' or 'sin'). -The integral is -estimated within the desired absolute error limit $epsabs. -The maximum number of allowed subdivisions done by the adaptive -algorithm must be supplied in the parameter $limit. - -=for usage - -Usage: - - gslinteg_qawf($function_ref,$omega,$sin_or_cos,$la,$epsabs,$limit,[opt]) - -=for example - -Example: - - my ($res,$abserr,$ierr) = gslinteg_qawf(\&f,$PI/2.0,'cos',0,1e-7,1000); - # with warnings on - ($res,$abserr,$ierr) = gslinteg_qawf(\&f,$PI/2.0,'cos',0,1e-7,1000,{Warn => 'y'}); - - sub f{ - my ($x) = @_; - if ($x == 0){return 0;} - return 1.0/sqrt($x) - } -EOF - Code => ' -gsl_error_handler_t * old_handler; -if ($gslwarn() == 1) { old_handler = gsl_set_error_handler(&my_handler); } -else { gsl_set_error_handler_off ();} -{gsl_integration_qawo_table * qtab; -gsl_integration_workspace *w; -gsl_integration_workspace *cw; -enum gsl_integration_qawo_enum T; - -T = GSL_INTEG_SINE; -if ($sincosopt() == 0){ T = GSL_INTEG_COSINE ;} - -qtab = gsl_integration_qawo_table_alloc($omega(),1.,T,(size_t) $nlevels()); - -set_funname($COMP(function)); -gsl_function F; -F.function = &FUNC; -F.params = 0; - -w = gsl_integration_workspace_alloc((size_t) $n()); -cw = gsl_integration_workspace_alloc((size_t) $n()); - -$ierr() = gsl_integration_qawf(&F,$a(),$epsabs(),(size_t) $limit(),w,cw,qtab,$P(result),$P(abserr)); - -gsl_integration_workspace_free(w); -gsl_integration_workspace_free(cw); -gsl_integration_qawo_table_free(qtab); -dec_func(); -} -'); - - -pp_done(); diff --git a/Libtmp/GSL/lib/PDL/GSL/INTERP.pd b/Libtmp/GSL/lib/PDL/GSL/INTERP.pd deleted file mode 100644 index f3bda5abd..000000000 --- a/Libtmp/GSL/lib/PDL/GSL/INTERP.pd +++ /dev/null @@ -1,378 +0,0 @@ -use strict; -use warnings; - -pp_bless('PDL::GSL::INTERP'); # make the functions generated go into our namespace - -pp_addpm({At=>'Top'},<<'EOD'); -use strict; -use warnings; - -=head1 NAME - -PDL::GSL::INTERP - PDL interface to Interpolation routines in GSL - -=head1 DESCRIPTION - -This is an interface to the interpolation package present in the -GNU Scientific Library. - -=head1 SYNOPSIS - - use PDL; - use PDL::GSL::INTERP; - - my $x = sequence(10); - my $y = exp($x); - - my $spl = PDL::GSL::INTERP->init('cspline',$x,$y); - - my $res = $spl->eval(4.35); - $res = $spl->deriv(4.35); - $res = $spl->deriv2(4.35); - $res = $spl->integ(2.1,7.4); - -=head1 NOMENCLATURE - -Throughout this documentation we strive to use the same variables that -are present in the original GSL documentation (see L). Oftentimes those variables are called C and -C. Since good Perl coding practices discourage the use of Perl -variables C<$a> and C<$b>, here we refer to Parameters C and C -as C<$pa> and C<$pb>, respectively, and Limits (of domain or -integration) as C<$la> and C<$lb>. -EOD - -pp_addpm({At=>'Bot'},<<'EOD'); -=head1 BUGS - -Feedback is welcome. - -=head1 SEE ALSO - -L - -The GSL documentation for interpolation is online at -L - -=head1 AUTHOR - -This file copyright (C) 2003 Andres Jordan -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. - -The GSL interpolation module was written by Gerard Jungman. - -=cut -EOD - -pp_addhdr(' -#include -#include -#include -#include -#include "gslerr.h" - -typedef gsl_spline GslSpline; -typedef gsl_interp_accel GslAccel; - -'); - -pp_def('init', - Pars => 'double x(n); double y(n);', - OtherPars => 'gsl_spline *spl', - Doc => <<'EOF', -=for ref - -The init method initializes a new instance of INTERP. It needs as -input an interpolation type and two ndarrays holding the x and y -values to be interpolated. The GSL routines require that x be -monotonically increasing and a quicksort is performed by default to -ensure that. You can skip the quicksort by passing the option -{Sort => 0}. - -The available interpolation types are : - -=over 2 - -=item linear - -=item polynomial - -=item cspline (natural cubic spline) - -=item cspline_periodic (periodic cubic spline) - -=item akima (natural akima spline) - -=item akima_periodic (periodic akima spline) - -=back - -Please check the GSL documentation for more information. - -=for usage - -Usage: - - $blessed_ref = PDL::GSL::INTERP->init($interp_method,$x,$y,$opt); - -=for example - -Example: - - $x = sequence(10); - $y = exp($x); - - $spl = PDL::GSL::INTERP->init('cspline',$x,$y) - $spl = PDL::GSL::INTERP->init('cspline',$x,$y,{Sort => 1}) #same as above - - # no sorting done on x, user is certain that x is monotonically increasing - $spl = PDL::GSL::INTERP->init('cspline',$x,$y,{Sort => 0}); -=cut -EOF - Code =>'GSLERR(gsl_spline_init,($COMP(spl),$P(x),$P(y),$SIZE(n)));', - PMCode => <<'EOF', -sub init { - my $opt; - if (ref($_[$#_]) eq 'HASH'){ $opt = pop @_; } - else{ $opt = {Sort => 1}; } - my ($class,$type,$x,$y) = @_; - if( (ref($x) ne 'PDL') || (ref($y) ne 'PDL') ){ - barf("Have to pass ndarrays as arguments to init method\n"); - } - if($$opt{Sort} != 0){ - my $idx = PDL::Ufunc::qsorti($x); - $x = $x->index($idx); - $y = $y->index($idx); - } - my $ene = nelem($x); - my $obj1 = new_spline($type,$ene); - my $obj2 = new_accel(); - _init_int($x,$y,$$obj1); - my @ret_a = ($obj1,$obj2); - return bless(\@ret_a, $class); -} -EOF -); - -pp_def('eval', - Pars => 'double x(); double [o] out();', - OtherPars => 'gsl_spline *spl;gsl_interp_accel *acc;', - Doc => <<'EOF', -=for ref - -The function eval returns the interpolating function at a given point. -It will barf with an "input domain error" if you try to extrapolate. - -=for usage - -Usage: - - $result = $spl->eval($points); - -=for example - -Example: - - my $res = $spl->eval($x) - -=cut -EOF - PMCode => <<'EOF', -sub eval { - my $opt; - my ($obj,$x) = @_; - my $s_obj = $$obj[0]; - my $a_obj = $$obj[1]; - _eval_int($x,my $o=PDL->null,$$s_obj,$$a_obj); - $o; -} -EOF - HandleBad => 1, - Code => ' - PDL_IF_BAD(if ($ISBAD($x())) { - $out() = $x(); - } else,) { - GSLERR(gsl_spline_eval_e,($COMP(spl), $x(), $COMP(acc), $P(out))); - } - ', -); - -pp_def('deriv', - Pars => 'double x(); double [o] out();', - OtherPars => 'gsl_spline *spl;gsl_interp_accel *acc;', - PMCode => <<'EOF', -sub deriv { - my ($obj,$x) = @_; - my $s_obj = $$obj[0]; - my $a_obj = $$obj[1]; - _deriv_int($x,my $o=PDL->null,$$s_obj,$$a_obj); - $o; -} -EOF - Doc => <<'EOF', -=for ref - -The deriv function returns the derivative of the -interpolating function at a given point. -It will barf with an "input domain error" if you try to extrapolate. - -=for usage - -Usage: - - $result = $spl->deriv($points); - -=for example - -Example: - - my $res = $spl->deriv($x) - -=cut -EOF - Code =>' -GSLERR(gsl_spline_eval_deriv_e,($COMP(spl), $x(), $COMP(acc), $P(out))); -'); - -pp_def('deriv2', - Pars => 'double x(); double [o] out();', - OtherPars => 'gsl_spline *spl;gsl_interp_accel *acc;', - Doc => <<'EOF', -=for ref - -The deriv2 function returns the second derivative -of the interpolating function at a given point. -It will barf with an "input domain error" if you try to extrapolate. - -=for usage - -Usage: - - $result = $spl->deriv2($points); - -=for example - -Example: - - my $res = $spl->deriv2($x) - -=cut -EOF - PMCode => <<'EOF', -sub deriv2 { - my ($obj,$x) = @_; - my $s_obj = $$obj[0]; - my $a_obj = $$obj[1]; - _deriv2_int($x,my $o=PDL->null,$$s_obj,$$a_obj); - $o; -} -EOF - Code =>' -GSLERR(gsl_spline_eval_deriv2_e,($COMP(spl), $x(), $COMP(acc), $P(out))); -'); - -pp_def('integ', - Pars => 'double a(); double b(); double [o] out();', - OtherPars => 'gsl_spline *spl;gsl_interp_accel *acc;', - Doc => <<'EOF', -=for ref - -The integ function returns the integral -of the interpolating function between two points. -It will barf with an "input domain error" if you try to extrapolate. - -=for usage - -Usage: - - $result = $spl->integ($la,$lb); - -=for example - -Example: - - my $res = $spl->integ($la,$lb) - -=cut -EOF - PMCode => <<'EOF', -sub integ { - my ($obj,$la,$lb) = @_; - my $s_obj = $$obj[0]; - my $a_obj = $$obj[1]; - _integ_int($la,$lb,my $o=PDL->null,$$s_obj,$$a_obj); - $o; -} -EOF - Code =>' -GSLERR(gsl_spline_eval_integ_e,($COMP(spl), $a(), $b(), $COMP(acc),$P(out))); -'); - -# XS functions for the INTERP objects - -pp_addxs('',' -MODULE = PDL::GSL::INTERP PACKAGE = PDL::GSL::INTERP - -#define DEF_INTERP(X) if (!strcmp(TYPE,#X)) spline=gsl_spline_alloc( gsl_interp_ ## X , ene); strcat(ula,#X ", "); - -GslSpline * -new_spline (TYPE,ene) - char *TYPE - int ene - CODE: - GslSpline * spline = NULL; - char ula[100]; -strcpy(ula,""); -DEF_INTERP(linear); -DEF_INTERP(polynomial); -DEF_INTERP(cspline); -DEF_INTERP(cspline_periodic); -DEF_INTERP(akima); -DEF_INTERP(akima_periodic); - if (spline==NULL) { - barf("Unknown interpolation type, please use one of the following: %s", ula); - } - else - RETVAL = spline; - OUTPUT: - RETVAL - -GslAccel * -new_accel () - CODE: - GslAccel * accel = NULL; - accel = gsl_interp_accel_alloc(); - if (accel == NULL){ - barf("Problem allocating accelerator object\n"); - } - RETVAL = accel; - OUTPUT: - RETVAL - -MODULE = PDL::GSL::INTERP PACKAGE = GslSplinePtr PREFIX = spl_ - -void -spl_DESTROY(spline) - GslSpline * spline - CODE: - gsl_spline_free(spline); - -MODULE = PDL::GSL::INTERP PACKAGE = GslAccelPtr PREFIX = acc_ - -void -acc_DESTROY(accel) - GslAccel * accel - CODE: - gsl_interp_accel_free(accel); - -'); - -pp_export_nothing; - -pp_add_boot('gsl_set_error_handler_off(); -'); - -pp_done(); diff --git a/Libtmp/GSL/lib/PDL/GSL/LINALG.pd b/Libtmp/GSL/lib/PDL/GSL/LINALG.pd deleted file mode 100644 index de819e710..000000000 --- a/Libtmp/GSL/lib/PDL/GSL/LINALG.pd +++ /dev/null @@ -1,186 +0,0 @@ -use strict; -use warnings; - -pp_addpm({At=>'Top'},<<'EOD'); -use strict; -use warnings; - -=head1 NAME - -PDL::GSL::LINALG - PDL interface to linear algebra routines in GSL - -=head1 SYNOPSIS - - use PDL::LiteF; - use PDL::MatrixOps; # for 'x' - use PDL::GSL::LINALG; - my $A = pdl [ - [0.18, 0.60, 0.57, 0.96], - [0.41, 0.24, 0.99, 0.58], - [0.14, 0.30, 0.97, 0.66], - [0.51, 0.13, 0.19, 0.85], - ]; - my $B = sequence(2,4); # column vectors - LU_decomp(my $lu=$A->copy, my $p=null, my $signum=null); - # transpose so first dim means is vector, higher dims broadcast - LU_solve($lu, $p, $B->transpose, my $x=null); - $x = $x->inplace->transpose; # now can be matrix-multiplied - -=head1 DESCRIPTION - -This is an interface to the linear algebra package present in the -GNU Scientific Library. Functions are named as in GSL, but with the -initial C removed. They are provided in both real and -complex double precision. - -Currently only LU decomposition interfaces here. Pull requests welcome! - -EOD - -pp_addpm({At=>'Bot'},<<'EOD'); -=head1 SEE ALSO - -L - -The GSL documentation for linear algebra is online at -L - -=cut - -EOD - -pp_addhdr(' -#include -#include "gslerr.h" - -#define MATRIX_SETUP(m, rows, cols, lda, datap) \\ - m.size1 = rows; \\ - m.size2 = cols; \\ - m.tda = lda; \\ - m.data = (double *)datap; \\ - m.owner = 0; - -/* rely on optimiser for this to produce sensible code */ -#define PERM_SETUP(p, psize, datap) \\ - gsl_permutation *p ## _local = gsl_permutation_alloc(psize); \\ - size_t *p ## _local_data = p ## _local->data; \\ - size_t p ## _i_local = 0; \\ - for (p ## _i_local = 0; p ## _i_local < (size_t)psize; p ## _i_local++) { \\ - p ## _local_data[p ## _i_local] = datap[p ## _i_local]; \\ - } \\ - p.data = p ## _local_data; \\ - p.size = psize; -#define PERM_RETURN(p, psize, datap) \\ - for (p ## _i_local = 0; p ## _i_local < (size_t)psize; p ## _i_local++) { \\ - datap[p ## _i_local] = p ## _local_data[p ## _i_local]; \\ - } \\ - gsl_permutation_free(p ## _local); - -#define VECTOR_SETUP(v, vsize, datap) \\ - v.size = vsize; \\ - v.stride = 1; \\ - v.data = (double *)datap; \\ - v.owner = 0; -'); - -pp_def('LU_decomp', - HandleBad => 0, - Pars => '[io,phys]A(n,m); indx [o,phys]ipiv(p=CALC($PDL(A)->ndims > 1 ? PDLMIN($PDL(A)->dims[0], $PDL(A)->dims[1]) : 1)); int [o,phys]signum()', - GenericTypes => [qw(C D)], - Code => <<'EOF', -/* make sure the PDL data types match */ -gsl_matrix$TDC(,_complex) m; -gsl_permutation p; -int s; -MATRIX_SETUP(m, $SIZE(m), $SIZE(n), $SIZE(n), $P(A)) -PERM_SETUP(p, $SIZE(p), $P(ipiv)) -GSLERR(gsl_linalg$TDC(,_complex)_LU_decomp, (&m, &p, &s)) -PERM_RETURN(p, $SIZE(p), $P(ipiv)) -$signum() = s; -EOF - Doc => <<'EOF', -=for ref - -LU decomposition of the given (real or complex) matrix. - -EOF -); - -pp_def('LU_solve', - HandleBad => 0, - Pars => '[phys]LU(n,m); indx [phys]ipiv(p); [phys]B(n); [o,phys]x(n)', - GenericTypes => [qw(C D)], - Code => <<'EOF', -gsl_matrix$TDC(,_complex) m; -gsl_permutation p; -gsl_vector$TDC(,_complex) b, x; -int s; -MATRIX_SETUP(m, $SIZE(m), $SIZE(n), $SIZE(n), $P(LU)) -PERM_SETUP(p, $SIZE(p), $P(ipiv)) -VECTOR_SETUP(b, $SIZE(n), $P(B)) -VECTOR_SETUP(x, $SIZE(n), $P(x)) -GSLERR(gsl_linalg$TDC(,_complex)_LU_solve, (&m, &p, &b, &x)) -PERM_RETURN(p, $SIZE(p), $P(ipiv)) -EOF - Doc => <<'EOF', -=for ref - -Solve C using the LU and permutation from L, real -or complex. - -EOF -); - -pp_def('LU_det', - HandleBad => 0, - Pars => '[phys]LU(n,m); int [phys]signum(); [o]det()', - GenericTypes => [qw(C D)], - Code => <<'EOF', -gsl_matrix$TDC(,_complex) m; -MATRIX_SETUP(m, $SIZE(m), $SIZE(n), $SIZE(n), $P(LU)) -types (D) %{ - $det() = gsl_linalg_LU_det(&m, $signum()); -%} -types (C) %{ - gsl_complex z = gsl_linalg_complex_LU_det(&m, $signum()); - $det() = GSL_REAL(z) + I*GSL_IMAG(z); -%} -EOF - Doc => <<'EOF', -=for ref - -Find the determinant from the LU decomp. - -EOF -); - -pp_def('solve_tridiag', - HandleBad => 0, - Pars => '[phys]diag(n); [phys]superdiag(n); [phys]subdiag(n); [phys]B(n); [o,phys]x(n)', - GenericTypes => [qw(D)], - Code => <<'EOF', -gsl_vector d, sup, sub, b, x; -VECTOR_SETUP(d, $SIZE(n), $P(diag)) -VECTOR_SETUP(sup, $SIZE(n)-1, $P(superdiag)) -VECTOR_SETUP(sub, $SIZE(n)-1, $P(subdiag)) -VECTOR_SETUP(b, $SIZE(n), $P(B)) -VECTOR_SETUP(x, $SIZE(n), $P(x)) -#define CONST_VEC (const gsl_vector *) -GSLERR(gsl_linalg_solve_tridiag, ( - CONST_VEC &d, CONST_VEC &sup, CONST_VEC &sub, CONST_VEC &b, &x -)) -#undef CONST_VEC -EOF - Doc => <<'EOF', -=for ref - -Solve C where A is a tridiagonal system. Real only, because -GSL does not have a complex function. - -EOF -); - -pp_add_boot('gsl_set_error_handler_off(); -'); - -pp_done(); diff --git a/Libtmp/GSL/lib/PDL/GSL/MROOT-FUNC.c b/Libtmp/GSL/lib/PDL/GSL/MROOT-FUNC.c deleted file mode 100644 index 18703b891..000000000 --- a/Libtmp/GSL/lib/PDL/GSL/MROOT-FUNC.c +++ /dev/null @@ -1,145 +0,0 @@ -// This file copyright (C) 2006 Andres Jordan -// and Simon Casassus -// 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. - -#include -#include -#include "EXTERN.h" -#include "perl.h" -#include "pdl.h" -#include "pdlcore.h" - -#define PDL PDL_GSL_MROOT -extern Core *PDL; - -static SV* ext_funname1; -static PDL_Indx ene; -void set_funname(SV *fn, PDL_Indx n) { - ext_funname1 = fn; - ene = n; -} - -void DFF(double* xval, double* vector){ - //this version tries just to get the output - dSP; - ENTER; - SAVETMPS; - - pdl* px = PDL->pdlnew(); - if (!px) PDL->pdl_barf("Failed to create pdl"); - SV* pxsv = sv_newmortal(); - PDL->SetSV_PDL(pxsv, px); - int ndims = 1; - PDL_Indx pdims[] = { (PDL_Indx) ene }; - PDL->barf_if_error(PDL->setdims(px,pdims,ndims)); - px->datatype = PDL_D; - px->data = (void *) xval; - px->state |= PDL_ALLOCATED | PDL_DONTTOUCHDATA; - - /* get function name on the perl side */ - PUSHMARK(SP); - XPUSHs(pxsv); - PUTBACK; - int count=call_sv(ext_funname1,G_SCALAR); - SPAGAIN; - SP -= count ; - I32 ax = (SP - PL_stack_base) + 1 ; - if (count!=1) - croak("error calling perl function\n"); - - /* recover output value */ - pdl* pvector = PDL->SvPDLV(ST(0)); - PDL->barf_if_error(PDL->make_physical(pvector)); - - double *xpass = (double *) pvector->data; - PDL_Indx i; - for(i=0;idata = NULL; - - PUTBACK; - FREETMPS; - LEAVE; -} - - -int my_f (const gsl_vector * v, void * params, gsl_vector * df) -{ - double dp = *((double *)params); - int nelem = (int) dp; - double xfree[nelem], vector[nelem]; - int iloop; - for (iloop=0;iloop< nelem;iloop++) { - xfree[iloop] = gsl_vector_get(v, iloop); - vector[iloop] = gsl_vector_get(v, iloop) * gsl_vector_get(v, iloop); - } - DFF(xfree, vector); - for (iloop=0;iloop< nelem;iloop++) { - gsl_vector_set(df, iloop, vector[iloop]); - } - return GSL_SUCCESS; -} - -int -print_state (size_t iter, gsl_multiroot_fsolver * s) -{ - printf ("iter = %3zu x = % .3f % .3f " - "f(x) = % .3e % .3e\n", - iter, - gsl_vector_get (s->x, 0), - gsl_vector_get (s->x, 1), - gsl_vector_get (s->f, 0), - gsl_vector_get (s->f, 1)); - return 1; -} - -int fsolver (double *xfree, int nelem, double epsabs, int method) -{ - gsl_multiroot_fsolver_type *T; - gsl_multiroot_fsolver *s; - int status; - size_t i, iter = 0; - size_t n = nelem; - double p[1] = { nelem }; - int iloop; - // struct func_params p = {1.0, 10.0}; - gsl_multiroot_function func = {&my_f, n, p}; - gsl_vector *x = gsl_vector_alloc (n); - for (iloop=0;iloopf, epsabs); - } - while (status == GSL_CONTINUE && iter < 1000); - if (status) - warn ("Final status = %s\n", gsl_strerror (status)); - for (iloop=0;iloopx, iloop); - } - gsl_multiroot_fsolver_free (s); - gsl_vector_free (x); - return GSL_SUCCESS; -} diff --git a/Libtmp/GSL/lib/PDL/GSL/MROOT.pd b/Libtmp/GSL/lib/PDL/GSL/MROOT.pd deleted file mode 100644 index 77ec6a86f..000000000 --- a/Libtmp/GSL/lib/PDL/GSL/MROOT.pd +++ /dev/null @@ -1,143 +0,0 @@ -use strict; -use warnings; - -{ no warnings 'once'; # pass info back to Makefile.PL -$PDL::Core::Dev::EXTRAS{$::PDLMOD}{OBJECT} .= join '', map " $::PDLBASE-$_\$(OBJ_EXT)", qw(FUNC); -} - -pp_bless('PDL::GSL::MROOT'); - -pp_add_exported('','gslmroot_fsolver'); - -pp_addhdr(' -#include -#include -#include -#include -#include - -void set_funname(SV *fn, PDL_Indx n); -int my_f (const gsl_vector * v, void * params, gsl_vector * df); -int fsolver (double *xfree, int nelem, double epsabs, int method); -'); - -pp_def('gslmroot_fsolver', - Pars => 'double [io]xfree(n); double epsabs(); int method();', - OtherPars => 'SV* function1;', - Doc => <<'EOF', -=for ref - -Multidimensional root finder without using derivatives - -This function provides an interface to the multidimensional root finding algorithms -in the GSL library. It takes a minimum of two arguments: an ndarray $init with an -initial guess for the roots of the system and a reference to a function. The latter -function must return an ndarray whose i-th element is the i-th equation evaluated at -the vector x (an ndarray which is the sole input to this function). See the example in -the Synopsis above for an illustration. The function returns an ndarray with the roots -for the system of equations. - -Two optional arguments can be specified as shown below. One is B, which can -take the values 0,1,2,3. They correspond to the 'hybrids', 'hybrid', 'dnewton' and -'broyden' algorithms respectively (see GSL documentation for details). The other -optional argument is B, which sets the absolute accuracy to which the roots -of the system of equations are required. The default value for Method is 0 ('hybrids' -algorithm) and the default for Epsabs is 1e-3. - -=for usage - - $res = gslmroot_fsolver($init, $function_ref, - [{Method => $method, Epsabs => $epsabs}]); - -EOF - Code =>' -set_funname($COMP(function1), $SIZE(n)); -if (fsolver($P(xfree), $SIZE(n), $epsabs(), $method()) != GSL_SUCCESS) - $CROAK("Something is wrong: could not assign fsolver type...\n"); -', - PMCode => <<'EOF', -sub gslmroot_fsolver { - my ($x, $f_vect) = @_; - my $opt = ref($_[-1]) eq 'HASH' ? pop @_ : {Method => 0, EpsAbs => 1e-3}; - if( (ref($x) ne 'PDL')){ - barf("Have to pass ndarray as first argument to fsolver\n"); - } - my $res = $x->copy; - _gslmroot_fsolver_int($res, $$opt{'EpsAbs'}, $$opt{'Method'}, $f_vect); - return $res; -} -EOF -); - -pp_addpm({At=>'Top'},<<'EOD'); -use strict; -use warnings; - -=head1 NAME - -PDL::GSL::MROOT - PDL interface to multidimensional root-finding routines in GSL - -=head1 DESCRIPTION - -This is an interface to the multidimensional root-finding package present in the -GNU Scientific Library. - -At the moment there is a single function B which provides an interface -to the algorithms in the GSL library that do not use derivatives. - - -=head1 SYNOPSIS - - use PDL; - use PDL::GSL::MROOT; - - my $init = pdl (-10.00, -5.0); - my $epsabs = 1e-7; - - - $res = gslmroot_fsolver($init, \&rosenbrock, - {Method => 0, EpsAbs => $epsabs}); - - - sub rosenbrock{ - my ($x) = @_; - my $c = 1; - my $d = 10; - my $y = zeroes($x); - - my $y0 = $y->slice(0); - $y0 .= $c * (1 - $x->slice(0)); - - my $y1 = $y->slice(1); - $y1 .= $d * ($x->slice(1) - $x->slice(0)**2); - - return $y; - } -EOD - -pp_addpm({At=>'Bot'},<<'EOD'); -=head1 SEE ALSO - -L - -The GSL documentation is online at - - http://www.gnu.org/software/gsl/manual/ - -=head1 AUTHOR - -This file copyright (C) 2006 Andres Jordan -and Simon Casassus -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_add_boot('gsl_set_error_handler_off(); -'); - -pp_done(); diff --git a/Libtmp/GSL/lib/PDL/GSL/RNG.pd b/Libtmp/GSL/lib/PDL/GSL/RNG.pd deleted file mode 100644 index a15ef5986..000000000 --- a/Libtmp/GSL/lib/PDL/GSL/RNG.pd +++ /dev/null @@ -1,1282 +0,0 @@ -use strict; -use warnings; - -use PDL::Types qw(types ppdefs_all); - -pp_bless('PDL::GSL::RNG'); # make the functions generated go into our namespace, and - # not PDL's namespace - -pp_addpm({At=>'Top'},<<'EOD'); -use strict; -use warnings; - -=head1 NAME - -PDL::GSL::RNG - PDL interface to RNG and randist routines in GSL - -=head1 DESCRIPTION - -This is an interface to the rng and randist packages present -in the GNU Scientific Library. - -=head1 SYNOPSIS - - use PDL; - use PDL::GSL::RNG; - - $rng = PDL::GSL::RNG->new('taus'); - - $rng->set_seed(time()); - - $x=zeroes(5,5,5) - - $rng->get_uniform($x); # inplace - - $y=$rng->get_uniform(3,4,5); # creates new pdl - -=head1 NOMENCLATURE - -Throughout this documentation we strive to use the same variables that -are present in the original GSL documentation (see L). Oftentimes those variables are called C and -C. Since good Perl coding practices discourage the use of Perl -variables C<$a> and C<$b>, here we refer to Parameters C and C -as C<$pa> and C<$pb>, respectively, and Limits (of domain or -integration) as C<$la> and C<$lb>. - -=cut -EOD - -pp_addpm({At=>'Middle'},<<'EOD'); -=head2 new - -=for ref - -The new method initializes a new instance of the RNG. - -The available RNGs are: - -=over - -=item coveyou - -=item cmrg - -=item fishman18 - -=item fishman20 - -=item fishman2x - -=item gfsr4 - -=item knuthran - -=item knuthran2 - -=item knuthran2002 - -=item lecuyer21 - -=item minstd - -=item mrg - -=item mt19937 - -=item mt19937_1999 - -=item mt19937_1998 - -=item r250 - -=item ran0 - -=item ran1 - -=item ran2 - -=item ran3 - -=item rand - -=item rand48 - -=item random128_bsd - -=item random128_glibc2 - -=item random128_libc5 - -=item random256_bsd - -=item random256_glibc2 - -=item random256_libc5 - -=item random32_bsd - -=item random32_glibc2 - -=item random32_libc5 - -=item random64_bsd - -=item random64_glibc2 - -=item random64_libc5 - -=item random8_bsd - -=item random8_glibc2 - -=item random8_libc5 - -=item random_bsd - -=item random_glibc2 - -=item random_libc5 - -=item randu - -=item ranf - -=item ranlux - -=item ranlux389 - -=item ranlxd1 - -=item ranlxd2 - -=item ranlxs0 - -=item ranlxs1 - -=item ranlxs2 - -=item ranmar - -=item slatec - -=item taus - -=item taus2 - -=item taus113 - -=item transputer - -=item tt800 - -=item uni - -=item uni32 - -=item vax - -=item waterman14 - -=item zuf - -=item default - -=back - -The last one (default) uses the environment variable GSL_RNG_TYPE. - -Note that only a few of these rngs are recommended for general -use. Please check the GSL documentation for more information. - -=for usage - -Usage: - - $blessed_ref = PDL::GSL::RNG->new($RNG_name); - -Example: - -=for example - - $rng = PDL::GSL::RNG->new('taus'); - -=head2 set_seed - -=for ref - -Sets the RNG seed. - -Usage: - -=for usage - - $rng->set_seed($integer); - # or - $rng = PDL::GSL::RNG->new('taus')->set_seed($integer); - -Example: - -=for example - - $rng->set_seed(666); - -=head2 min - -=for ref - -Return the minimum value generable by this RNG. - -Usage: - -=for usage - - $integer = $rng->min(); - -Example: - -=for example - - $min = $rng->min(); $max = $rng->max(); - -=head2 max - -=for ref - -Return the maximum value generable by the RNG. - -Usage: - -=for usage - - $integer = $rng->max(); - -Example: - -=for example - - $min = $rng->min(); $max = $rng->max(); - -=head2 name - -=for ref - -Returns the name of the RNG. - -Usage: - -=for usage - - $string = $rng->name(); - -Example: - -=for example - - $name = $rng->name(); - -=head2 ran_shuffle - -=for ref - -Shuffles values in ndarray, treating it as flat. - -Usage: - -=for usage - - $rng->ran_shuffle($ndarray); - -=head2 ran_shuffle_vec - -=for ref - -Returns values in Perl list, shuffled. - -Usage: - -=for usage - - @shuffled = $rng->ran_shuffle_vec(@vec); - -=head2 ran_choose - -=for ref - -Chooses values from C<$inndarray> to C<$outndarray>, treating both as flat. - -Usage: - -=for usage - - $rng->ran_choose($inndarray,$outndarray); - -=head2 ran_choose_vec - -=for ref - -Chooses C<$n> values from C<@vec>. - -Usage: - -=for usage - - @chosen = $rng->ran_choose_vec($n,@vec); - -=head2 ran_dir - -=for ref - -Returns C<$n> random vectors in C<$ndim> dimensions. - -Usage: - -=for usage - - $ndarray = $rng->ran_dir($ndim,$n); - -Example: - -=for example - - $o = $rng->ran_dir($ndim,$n); - -=head2 ran_discrete_preproc - -=for ref - -This method returns a handle that must be used when calling -L. You specify the probability of the integer number -that are returned by L. - -Usage: - -=for usage - - $discrete_dist_handle = $rng->ran_discrete_preproc($double_ndarray_prob); - -Example: - -=for example - - $prob = pdl [0.1,0.3,0.6]; - $ddh = $rng->ran_discrete_preproc($prob); - $o = $rng->ran_discrete($discrete_dist_handle,100); - -=cut -EOD - -pp_addpm({At=>'Bot'},<<'EOD'); -=head1 BUGS - -Feedback is welcome. Log bugs in the PDL bug database (the -database is always linked from L). - -=head1 SEE ALSO - -L - -The GSL documentation for random number distributions is online at -L - -=head1 AUTHOR - -This file copyright (C) 1999 Christian Pellegrin -Docs mangled by C. Soeller. 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. - -The GSL RNG and randist modules were written by James Theiler. - -=cut -EOD - -# PP interface to RNG - - -############################## -# -# make_get_sub generates a wrapper PDL subroutine that handles the -# fill-a-PDL and create-a-PDL cases for each of the GSL functions. -# --CED -# -sub make_get_sub { - my ($fname,$par) =@_; - my $s = ' -sub ' . $fname . ' { -my ($obj,' . $par . '@var) = @_;'; - - if ($par ne '') { - my $ss=$par; - - $ss =~ s/,//; - $s .= 'if (!(' . $ss . '>0)) {barf("first parameter must be an int >0")};'; - } - -$s .= 'if (ref($var[0]) eq \'PDL\') { - _' . $fname . '_int($var[0],' . $par . '$obj); - return $var[0]; -} -else { - my $p; - - $p = zeroes @var; - _' . $fname . '_int($p,' . $par . '$obj); - return $p; -} -} -' -} - -pp_addhdr(' -#include - -#include "gsl/gsl_rng.h" -#include "gsl/gsl_randist.h" -'); - -use Config; -my %p_type = ( - 'double' => 'double', - 'unsigned int' => $Config{intsize} == 4 ? 'uint' : - $Config{intsize} == 8 ? 'ulonglong' : die("Unknown intsize, not [48]"), - 'size_t' => 'SIZE', -); -use Text::ParseWords qw(shellwords); -chomp(my $header = `gsl-config --cflags`); -$header =~ s#\\#/#g; # win32 -($header) = map {s/^-I//;$_} grep /^-I/, shellwords $header; -$header .= '/gsl/gsl_randist.h'; -my $h = do { open my $fh, "<", $header or die "$header: $!"; local $/ = undef; <$fh> }; -my @defs; -DEF: for my $def ($h =~ m/(double\s+gsl_ran_\S*_pdf\s*\(.+?\))/sg) { - local $_ = $def; - s/\n\s+/ /xmsg; - s/const //g; - if (m/^(\w+)\s+(\w+)\s*\(\s*(.+)\s*\)/s) { - my ($out_type, $function, $pars) = ($1, $2, $3); - my @pars = split /,/, $pars; - for (@pars) { - die "Couldn't parse '$_' in '$def'" if !m/^(.+?)(\w+)\s*(\[\s*\])?\s*$/; - my ($type, $par, $vec) = ($1, $2, $3); - next DEF if $type =~ /gsl_ran_discrete_t/; - s/^ | $//g for ($type, $par); - $par = lc $par; - # numbers interfere with $ISGOOD in BadCode - $par =~ s/1/a/g; - $par =~ s/2/b/g; - die("unknown type '$type' in '$def'") if !$p_type{$type}; - $_ = [$p_type{$type}, $par, $vec ? 1 : ()]; - } - push @defs, [ $out_type, $function, \@pars ]; - } -} -@defs = sort { $a->[1] cmp $b->[1] } @defs; # sort by function name -my @export_names; -for my $f (@defs) { - my ($out_type, $function, $pars) = @$f; - my ($p, $code) = gen_def($out_type, $function, @$pars); - (my $pp_name = $function) =~ s/^gsl_//; - push @export_names, $pp_name; - pp_def($pp_name, - HandleBad => 1, - GenericTypes => ['D'], - Pars => $p, - Code => $code, - Doc => '', - ); -} -sub gen_def { - my ($out_type, $function, @pars) = @_; - my $SIZEVAR; - my @formal_pars = map { - if ($_->[0] eq 'SIZE') { - die "$function got second sizevar (@$_) but already got '$SIZEVAR'" - if $SIZEVAR; - $SIZEVAR = $_->[1]; - () - } else { - my $prefix = $_->[0] eq 'double' ? '' : "$_->[0] "; - my $suffix = !$_->[2] ? '()' : - $SIZEVAR ? "($SIZEVAR)" : - die "$function got vector par (@$_) but no sizevar"; - $prefix.$_->[1].$suffix - } - } @pars; - my $pars = join '; ', @formal_pars, '[o]out()'; - my @call_pars = map $_->[0] eq 'SIZE' ? "\$SIZE($_->[1])" : - $_->[2] ? "\$P($_->[1])" : "\$$_->[1]()", - @pars; - my $code = join "\n", "PDL_IF_BAD(char anybad=0;,)", - (map $_->[2] ? - "PDL_IF_BAD(loop ($SIZEVAR) %{ if (\$ISBAD($_->[1]())) { anybad=1; break; } %},)" : - "PDL_IF_BAD(if ( \$ISBAD($_->[1]())) anybad=1;,)", - grep $_->[0] ne 'SIZE', @pars), - "PDL_IF_BAD(if (anybad) { \$SETBAD(out()); continue; },)", - "\$out() = $function(@{[join ', ', @call_pars]});\n"; - return ($pars, $code); -} - -sub pp_defnd { # hide the docs - my ($name, %hash) = @_; - pp_def($name,%hash,Doc=>undef); -} - -pp_def('ran_shuffle_1d', - Pars => '[io]a(n)', - GenericTypes => [ppdefs_all()], - OtherPars => 'gsl_rng *rng', - Code => ' -gsl_ran_shuffle($COMP(rng), $P(a), $SIZE(n), sizeof($GENERIC())); -', - PMCode => <<'EOF', -sub ran_shuffle_1d { _ran_shuffle_1d_int(@_[1,0]) } -EOF - Doc => <<'EOF', -=for ref - -Takes n-dimensional ndarray, and shuffles it along its zero-th dimension. - -Usage: - -=for usage - - $vec2d = sequence(10,10); - $rng->ran_shuffle_1d($vec2d); -EOF -); - -pp_def('get_uniform', - Pars => '[o]a()', - GenericTypes => ['F','D'], - OtherPars => 'gsl_rng *rng', - Code => ' -$a() = gsl_rng_uniform($COMP(rng));', - PMCode => make_get_sub('get_uniform',''), - Doc => <<'EOF', -=for ref - -This function creates an ndarray with given dimensions or accepts an -existing ndarray and fills it. get_uniform() returns values 0<=x<1, - -Usage: - -=for usage - - $ndarray = $rng->get_uniform($list_of_integers) - $rng->get_uniform($ndarray); - -Example: - -=for example - - $x = zeroes 5,6; $max=100; - $o = $rng->get_uniform(10,10); $rng->get_uniform($x); -EOF -); - -pp_def('get_uniform_pos', - Pars => '[o]a()', - GenericTypes => ['F','D'], - OtherPars => 'gsl_rng *rng', - Code => ' -$a() = gsl_rng_uniform_pos($COMP(rng));', - PMCode => make_get_sub('get_uniform_pos',''), - Doc => <<'EOF', -=for ref - -This function creates an ndarray with given dimensions or accepts an -existing ndarray and fills it. get_uniform_pos() returns values 0get_uniform_pos($list_of_integers) - $rng->get_uniform_pos($ndarray); - -Example: - -=for example - - $x = zeroes 5,6; - $o = $rng->get_uniform_pos(10,10); $rng->get_uniform_pos($x); -EOF -); - -pp_def('get', - Pars => '[o]a()', - GenericTypes => ['F','D'], - OtherPars => 'gsl_rng *rng', - Code => ' -$a() = gsl_rng_get($COMP(rng));', - PMCode => make_get_sub('get',''), - Doc => <<'EOF', -=for ref - -This function creates an ndarray with given dimensions or accepts an -existing ndarray and fills it. get() returns integer values -between a minimum and a maximum specific to every RNG. - -Usage: - -=for usage - - $ndarray = $rng->get($list_of_integers) - $rng->get($ndarray); - -Example: - -=for example - - $x = zeroes 5,6; - $o = $rng->get(10,10); $rng->get($x); -EOF -); - -pp_def('get_int', - Pars => '[o]a()', - GenericTypes => ['F','D'], - OtherPars => 'IV n; gsl_rng *rng', - Code => ' -$a() = gsl_rng_uniform_int($COMP(rng),$COMP(n));', - PMCode => make_get_sub('get_int','$n,'), - Doc => <<'EOF', -=for ref - -This function creates an ndarray with given dimensions or accepts an -existing ndarray and fills it. get_int() returns integer values -between 0 and $max. - -Usage: - -=for usage - - $ndarray = $rng->get($max, $list_of_integers) - $rng->get($max, $ndarray); - -Example: - -=for example - - $x = zeroes 5,6; $max=100; - $o = $rng->get(10,10); $rng->get($x); -EOF -); - -# randist stuff - -sub add_randist { - my ($name,$doc,$params) = @_; - my $npar = @$params; - my $pars1=join '; ', (map "double $_", @$params), 'gsl_rng *rng'; - my $fcall1=join ',', map "\$COMP($_)", @$params; - my $arglist=join '', map "\$$_,", @$params; - my $pars2=join ';', map "$_()", @$params; - my $fcall2=join ',', map "\$$_()", @$params; - my $arglist2=join ',', map "\$${_}_ndarray", @$params; - - pp_def( - 'ran_' . $name, - Pars => '[o]output()', - GenericTypes => ['F','D'], - OtherPars => $pars1, - Code =>' -$output() = gsl_ran_' . $name . '($COMP(rng),' . $fcall1 . ');', - PMCode =>' -sub ran_' . $name . ' { -my ($obj,' . $arglist . '@var) = @_; -if (ref($var[0]) eq \'PDL\') { - _ran_' . $name . '_int($var[0],' . $arglist . '$obj); - return $var[0]; -} -else { - my $p; - $p = zeroes @var; - _ran_' . $name . '_int($p,' . $arglist . '$obj); - return $p; -} -} -', - Doc => <ran_$name(@{[join ', ', map qq{\$$_}, @$params]},[list of integers = output ndarray dims]); - \$rng->ran_$name(@{[join ', ', map qq{\$$_}, @$params]}, \$output_ndarray); - -Example: - -=for example - - \$o = \$rng->ran_$name(@{[join ', ', map qq{\$$_}, @$params]},10,10); - \$rng->ran_$name(@{[join ', ', map qq{\$$_}, @$params]},\$o); -EOF -); - - pp_def( - 'ran_' . $name . '_var', - Pars => $pars2 . ';[o]output()', - GenericTypes => ['F','D'], - OtherPars => 'gsl_rng *rng', - Code =>' -$output() = gsl_ran_' . $name . '($COMP(rng),' . $fcall2 . ');', - PMCode =>' -sub ran_' . $name . '_var { -my ($obj,@var) = @_; - if (scalar(@var) != ' . $npar . ') {barf("Bad number of parameters!");} - _ran_' . $name . '_var_int(@var,my $x=PDL->null,$obj); - return $x; -} -', - Doc => < except that it takes the distribution -parameters as an ndarray and returns an ndarray of equal dimensions. - -Usage: - -=for usage - - \$ndarray = \$rng->ran_${name}_var($arglist2); -EOF - ); -} - -add_randist('gaussian','values from Gaussian distribution with mean zero and standard deviation C<$sigma>.',[qw(sigma)]); -add_randist('ugaussian_tail', 'variates from the upper tail of a Gaussian distribution with C (AKA unit Gaussian distribution).',[qw(tail)]); -add_randist('exponential', 'variates from the exponential distribution with mean C<$mu>.',[qw(mu)]); -add_randist('laplace', 'variates from the Laplace distribution with width C<$pa>.',[qw(pa)]); -add_randist('exppow', 'variates from the exponential power distribution with scale parameter C<$pa> and exponent C<$pb>.',[qw(pa pb)]); -add_randist('cauchy', 'variates from the Cauchy distribution with scale parameter C<$pa>.',[qw(pa)]); -add_randist('rayleigh', 'variates from the Rayleigh distribution with scale parameter C<$sigma>.',[qw(sigma)]); -add_randist('rayleigh_tail', 'variates from the tail of the Rayleigh distribution -with scale parameter C<$sigma> and a lower limit of C<$la>.',[qw(x sigma)]); -add_randist('levy', 'variates from the Levy symmetric stable distribution with scale C<$c> and exponent C<$alpha>.',[qw(mu x)]); -add_randist('gamma', 'variates from the gamma distribution.',[qw(pa pb)]); -add_randist('flat', 'variates from the flat (uniform) distribution from C<$la> to C<$lb>.',[qw(la lb)]); -add_randist('lognormal', 'variates from the lognormal distribution with parameters C<$mu> (location) and C<$sigma> (scale).',[qw(mu sigma)]); -add_randist('chisq', 'variates from the chi-squared distribution with C<$nu> degrees of freedom.',[qw(nu)]); -add_randist('fdist', 'variates from the F-distribution with degrees of freedom C<$nu1> and C<$nu2>.',[qw(nu1 nu2)]); -add_randist('tdist', q{variates from the t-distribution (AKA Student's -t-distribution) with C<$nu> degrees of freedom.},[qw(nu)]); -add_randist('beta', 'variates from the beta distribution with parameters C<$pa> and C<$pb>.',[qw(pa pb)]); -add_randist('logistic', 'random variates from the logistic distribution.',[qw(m)]); -add_randist('pareto', 'variates from the Pareto distribution of order C<$pa> and scale C<$lb>.',[qw(pa lb)]); -add_randist('weibull', 'variates from the Weibull distribution with scale C<$pa> and exponent C<$pb>. (Some literature uses C for C<$pa> and C for C<$pb>.)',[qw(pa pb)]); -add_randist('gumbel1', 'variates from the Type-1 Gumbel distribution.',[qw(pa pb)]); -add_randist('gumbel2', 'variates from the Type-2 Gumbel distribution.',[qw(pa pb)]); -add_randist('poisson','integer values from the Poisson distribution with mean C<$mu>.',[qw(mu)]); -add_randist('bernoulli', 'values 0 or 1, the result of a Bernoulli trial with probability C<$p>.',[qw(p)]); -add_randist('binomial', 'integer values from the binomial distribution, the number of successes in C<$n> independent trials with probability C<$p>.',[qw(p n)]); -add_randist('negative_binomial', 'integer values from the negative binomial -distribution, the number of failures occurring before C<$n> successes in -independent trials with probability C<$p> of success. Note that C<$n> is -not required to be an integer.',[qw(p n)]); -add_randist('pascal', 'integer values from the Pascal distribution. -The Pascal distribution is simply a negative binomial distribution -(see L) with an integer value of C<$n>.',[qw(p n)]); -add_randist('geometric', 'integer values from the geometric distribution, the number of independent trials with probability C<$p> until the first success.',[qw(p)]); -add_randist('hypergeometric', 'integer values from the hypergeometric distribution. -If a population contains C<$n1> elements of type 1 and C<$n2> elements of -type 2 then the hypergeometric distribution gives the probability of obtaining -C<$x> elements of type 1 in C<$t> samples from the population without replacement.',[qw(n1 n2 t)]); -add_randist('logarithmic', 'integer values from the logarithmic distribution.',[qw(p)]); - -# specific randist - -pp_def('ran_additive_gaussian', - Pars => '[o]x()', - GenericTypes => ['F','D'], - OtherPars => 'double sigma; gsl_rng *rng', - Code =>'$x() += gsl_ran_gaussian($COMP(rng), $COMP(sigma));', - PMCode =>' - sub ran_additive_gaussian { - my ($obj,$sigma,$var) = @_; - barf("In additive gaussian mode you must specify an ndarray!") - if ref($var) ne \'PDL\'; - _ran_additive_gaussian_int($var,$sigma,$obj); - return $var; - } - ', - Doc => <<'EOF', -=for ref - -Add Gaussian noise of given sigma to an ndarray. - -Usage: - -=for usage - - $rng->ran_additive_gaussian($sigma,$ndarray); - -Example: - -=for example - - $rng->ran_additive_gaussian(1,$image); -EOF -); - -pp_def('ran_additive_poisson', - Pars => '[o]x()', - GenericTypes => ['F','D'], - OtherPars => 'double sigma; gsl_rng *rng', - Code =>'$x() += gsl_ran_poisson($COMP(rng), $COMP(sigma));', - PMCode =>' - sub ran_additive_poisson { - my ($obj,$sigma,$var) = @_; - barf("In additive poisson mode you must specify an ndarray!") - if ref($var) ne \'PDL\'; - _ran_additive_poisson_int($var,$sigma,$obj); - return $var; - } - ', - Doc => <<'EOF', -=for ref - -Add Poisson noise of given C<$mu> to a C<$ndarray>. - -Usage: - -=for usage - - $rng->ran_additive_poisson($mu,$ndarray); - -Example: - -=for example - - $rng->ran_additive_poisson(1,$image); -EOF -); - -pp_def('ran_feed_poisson', - Pars => '[o]x()', - GenericTypes => ['F','D'], - OtherPars => 'gsl_rng *rng', - Code =>'$x() = gsl_ran_poisson($COMP(rng), $x());', - PMCode =>' - sub ran_feed_poisson { - my ($obj,$var) = @_; - barf("In poisson mode you must specify an ndarray!") - if ref($var) ne \'PDL\'; - _ran_feed_poisson_int($var,$obj); - return $var; - } - ', - Doc => <<'EOF', -=for ref - -This method simulates shot noise, taking the values of ndarray as -values for C<$mu> to be fed in the poissonian RNG. - -Usage: - -=for usage - - $rng->ran_feed_poisson($ndarray); - -Example: - -=for example - - $rng->ran_feed_poisson($image); -EOF -); - -pp_def('ran_bivariate_gaussian', - Pars => '[o]x(n)', - GenericTypes => ['F','D'], - OtherPars => 'double sigma_x; double sigma_y; double rho; gsl_rng *rng', - Code => <<'EOF', -double xx,yy; -gsl_ran_bivariate_gaussian($COMP(rng), $COMP(sigma_x), $COMP(sigma_y),$COMP(rho), &xx, &yy); -$x(n=>0)=xx; -$x(n=>1)=yy; -EOF - PMCode => <<'EOF', -sub ran_bivariate_gaussian { - my ($obj,$sigma_x,$sigma_y,$rho,$n) = @_; - barf("Not enough parameters for gaussian bivariate!") if $n<=0; - my $p = zeroes(2,$n); - _ran_bivariate_gaussian_int($p,$sigma_x,$sigma_y,$rho,$obj); - return $p; -} -EOF - Doc => <<'EOF', -=for ref - -Generates C<$n> bivariate gaussian random deviates. - -Usage: - -=for usage - - $ndarray = $rng->ran_bivariate_gaussian($sigma_x,$sigma_y,$rho,$n); - -Example: - -=for example - - $o = $rng->ran_bivariate_gaussian(1,2,0.5,1000); -EOF -); - -pp_defnd('ran_dir_2d', - Pars => '[o]x(n)', - GenericTypes => ['F','D'], - OtherPars => 'gsl_rng *rng', - Code =>' -double xx,yy; -gsl_ran_dir_2d($COMP(rng), &xx, &yy); -$x(n=>0)=xx; -$x(n=>1)=yy; -'); - -pp_defnd('ran_dir_3d', - Pars => '[o]x(n)', - GenericTypes => ['F','D'], - OtherPars => 'gsl_rng *rng', - Code =>' -double xx,yy,zz; -gsl_ran_dir_3d($COMP(rng), &xx, &yy, &zz); -$x(n=>0)=xx; -$x(n=>1)=yy; -$x(n=>2)=zz; -'); - -my $MAX_DIMENSIONS = 100; - -pp_defnd('ran_dir_nd', - Pars => '[o]x(n)', - GenericTypes => ['F','D'], - OtherPars => 'IV ns => n; gsl_rng *rng', - Code =>' -double xxx[' . $MAX_DIMENSIONS .']; -gsl_ran_dir_nd($COMP(rng), $COMP(ns), xxx); -loop (n) %{ $x() = xxx[n]; %}'); - -pp_addpm(' - sub ran_dir { - my ($obj,$ndim,$n) = @_; - barf("Not enough parameters for random vectors!") if $n<=0; - my $p = zeroes($ndim,$n); - if ($ndim==2) { ran_dir_2d($p,$obj); } - elsif ($ndim==3) { ran_dir_3d($p,$obj); } - elsif ($ndim>=4 && $ndim<=' . $MAX_DIMENSIONS . ') { ran_dir_nd($p,$ndim,$obj); } - else { barf("Bad number of dimensions!"); } - return $p; - } - '); - -pp_def('ran_discrete', - Pars => '[o]x()', - GenericTypes => ['F','D'], - OtherPars => 'gsl_ran_discrete_t *rng_discrete; gsl_rng *rng', - Code =>' -$x()=gsl_ran_discrete($COMP(rng), $COMP(rng_discrete)); ', - PMCode =>' -sub ran_discrete { -my ($obj, $rdt, @var) = @_; -if (ref($var[0]) eq \'PDL\') { - _ran_discrete_int($var[0], $rdt, $obj); - return $var[0]; -} -else { - my $p; - $p = zeroes @var; - _ran_discrete_int($p, $rdt, $obj); - return $p; -} -} -', - Doc => <<'EOF', -=for ref - -Is used to get the desired samples once a proper handle has been -enstablished (see ran_discrete_preproc()). - -Usage: - -=for usage - - $ndarray = $rng->ran_discrete($discrete_dist_handle,$num); - -Example: - -=for example - - $prob = pdl [0.1,0.3,0.6]; - $ddh = $rng->ran_discrete_preproc($prob); - $o = $rng->ran_discrete($discrete_dist_handle,100); -EOF -); - - pp_addpm(' -sub ran_shuffle_vec { -my ($obj,@in) = @_; -$obj->ran_shuffle(my $p = PDL->sequence(PDL::indx(), 0+@in)); -@in[$p->list]; -} -'); - - pp_addpm(' -sub ran_choose_vec { -my ($obj,$nout,@in) = @_; -$obj->ran_choose(PDL->sequence(PDL::indx(), 0+@in),my $pout = PDL->zeroes(PDL::indx(), $nout)); -@in[$pout->list]; -} -'); - - -pp_def('ran_ver', - Pars => '[o]x(n)', - GenericTypes => ['F','D'], - OtherPars => 'double x0; double r;IV ns => n; gsl_rng *rng', - Code =>' -double xx=$COMP(x0); - -loop (n) %{ $x() = xx; xx = $COMP(r)*(1-xx)*xx; %}', - PMCode =>' - sub ran_ver { - my ($obj,$x0,$r,$n) = @_; - barf("Not enough parameters for ran_ver!") if $n<=0; - my $p = zeroes($n); - _ran_ver_int($p,$x0,$r,$n,$obj); - return $p; - } - ', - Doc => <<'EOF', -=for ref - -Returns an ndarray with C<$n> values generated by the Verhulst map from C<$x0> and -parameter C<$r>. - -Usage: - -=for usage - - $rng->ran_ver($x0, $r, $n); -EOF -); - -pp_def('ran_caos', - Pars => '[o]x(n)', - GenericTypes => ['F','D'], - OtherPars => 'double m; IV ns => n; gsl_rng *rng', - Code =>' -double xx=gsl_ran_gaussian($COMP(rng),0.1)+0.5; - -loop (n) %{ $x() = (xx-0.5)*$COMP(m); xx = 4.0*(1-xx)*xx; %}', - PMCode =>' - sub ran_caos { - my ($obj,$m,$n) = @_; - barf("Not enough parameters for ran_caos!") if $n<=0; - my $p = zeroes($n); - _ran_caos_int($p,$m,$n,$obj); - return $p; - } - ', - Doc => <<'EOF', -=for ref - -Returns values from Verhuls map with C<$r=4.0> and randomly chosen -C<$x0>. The values are scaled by C<$m>. - -Usage: - -=for usage - - $rng->ran_caos($m,$n); -EOF -); - -# XS function for the RNG object - -pp_addxs('',' -MODULE = PDL::GSL::RNG PACKAGE = PDL::GSL::RNG - -#define DEF_RNG(X) if (!strcmp(TYPE,#X)) rng=gsl_rng_alloc( gsl_rng_ ## X ); strcat(rngs,#X ", "); - -gsl_rng * -new (CLASS,TYPE) - char *CLASS - char *TYPE - CODE: - gsl_rng * rng = NULL; - char rngs[5000]; -strcpy(rngs,""); -DEF_RNG(borosh13) -DEF_RNG(coveyou) -DEF_RNG(cmrg) -DEF_RNG(fishman18) -DEF_RNG(fishman20) -DEF_RNG(fishman2x) -DEF_RNG(gfsr4) -DEF_RNG(knuthran) -DEF_RNG(knuthran2) -DEF_RNG(knuthran2002) -DEF_RNG(lecuyer21) -DEF_RNG(minstd) -DEF_RNG(mrg) -DEF_RNG(mt19937) -DEF_RNG(mt19937_1999) -DEF_RNG(mt19937_1998) -DEF_RNG(r250) -DEF_RNG(ran0) -DEF_RNG(ran1) -DEF_RNG(ran2) -DEF_RNG(ran3) -DEF_RNG(rand) -DEF_RNG(rand48) -DEF_RNG(random128_bsd) -DEF_RNG(random128_glibc2) -DEF_RNG(random128_libc5) -DEF_RNG(random256_bsd) -DEF_RNG(random256_glibc2) -DEF_RNG(random256_libc5) -DEF_RNG(random32_bsd) -DEF_RNG(random32_glibc2) -DEF_RNG(random32_libc5) -DEF_RNG(random64_bsd) -DEF_RNG(random64_glibc2) -DEF_RNG(random64_libc5) -DEF_RNG(random8_bsd) -DEF_RNG(random8_glibc2) -DEF_RNG(random8_libc5) -DEF_RNG(random_bsd) -DEF_RNG(random_glibc2) -DEF_RNG(random_libc5) -DEF_RNG(randu) -DEF_RNG(ranf) -DEF_RNG(ranlux) -DEF_RNG(ranlux389) -DEF_RNG(ranlxd1) -DEF_RNG(ranlxd2) -DEF_RNG(ranlxs0) -DEF_RNG(ranlxs1) -DEF_RNG(ranlxs2) -DEF_RNG(ranmar) -DEF_RNG(slatec) -DEF_RNG(taus) -DEF_RNG(taus2) -DEF_RNG(taus113) -DEF_RNG(transputer) -DEF_RNG(tt800) -DEF_RNG(uni) -DEF_RNG(uni32) -DEF_RNG(vax) -DEF_RNG(waterman14) -DEF_RNG(zuf) -DEF_RNG(default) - if (rng==NULL) { - barf("Unknown RNG, please use one of the following: %s", rngs); - } - else - RETVAL = rng; - OUTPUT: - RETVAL - -void -set_seed(rng, seed) - gsl_rng * rng - int seed - PPCODE: - gsl_rng_set(rng,seed); - XPUSHs(ST(0)); /* return self */ - -unsigned int -min(rng) - gsl_rng * rng - CODE: - RETVAL = gsl_rng_min(rng); - OUTPUT: - RETVAL - -unsigned int -max(rng) - gsl_rng * rng - CODE: - RETVAL = gsl_rng_max(rng); - OUTPUT: - RETVAL - -char* -name(rng) - gsl_rng * rng - CODE: - RETVAL = (char *) gsl_rng_name(rng); - OUTPUT: - RETVAL - -void -DESTROY(sv) - SV * sv - CODE: - gsl_rng *rng = INT2PTR(gsl_rng *, SvIV((SV*)SvRV(sv))); - gsl_rng_free(rng); - -gsl_ran_discrete_t * -ran_discrete_preproc(rng, p) - gsl_rng * rng - pdl * p - CODE: - IV n; - if (p->ndims!=1 || p->datatype!=PDL_D) { - barf("Bad input to ran_discrete_preproc!"); - } - PDL->barf_if_error(PDL->make_physical(p)); - n = p->dims[0]; - RETVAL = gsl_ran_discrete_preproc(n,(double *) p->data); - OUTPUT: - RETVAL - -void -ran_shuffle(rng, in) - gsl_rng * rng - pdl * in - CODE: - IV size, n; - PDL->barf_if_error(PDL->make_physical(in)); - n = in->nvals; - size = PDL->howbig(in->datatype); - gsl_ran_shuffle(rng,in->data,n,size); - PDL->changed(in, PDL_PARENTDATACHANGED, 0); - -void -ran_choose(rng, in, out) - gsl_rng * rng - pdl * in - pdl * out - CODE: - IV size, n,m; - - if (in->datatype != out->datatype) barf("Data Types must match for ran_chooser"); - PDL->barf_if_error(PDL->make_physical(in)); - PDL->barf_if_error(PDL->make_physical(out)); - n = in->nvals; - m = out->nvals; - size = PDL->howbig(in->datatype); - gsl_ran_choose(rng,out->data, m, in->data,n,size); - PDL->changed(out, PDL_PARENTDATACHANGED, 0); -'); - -pp_core_importList(' qw/ zeroes long barf /'); # import just a named list to our namespace, so we don't get warning - # messages like 'warning 'min' redefined at line ...' - -pp_export_nothing; # set to not export anything. (This is a OO package, it doesn't need to export any methods.) -pp_add_exported(@export_names); # ... except functions suitable for that - -pp_add_boot('gsl_set_error_handler_off(); -'); - -pp_done(); diff --git a/Libtmp/GSL/lib/PDL/GSL/SF.pd b/Libtmp/GSL/lib/PDL/GSL/SF.pd deleted file mode 100644 index fbcdf7847..000000000 --- a/Libtmp/GSL/lib/PDL/GSL/SF.pd +++ /dev/null @@ -1,2012 +0,0 @@ -use strict; -use warnings; - -pp_addpm({At=>'Top'},<<'EOD'); -use strict; -use warnings; - -=head1 NAME - -PDL::GSL::SF - PDL interface to GSL Special Functions - -=head1 DESCRIPTION - -This is an interface to the Special Function package present in the GNU Scientific Library. - -=cut - -EOD - -pp_add_boot("gsl_set_error_handler_off();\n"); - -pp_addhdr(' -#include -#include -#include "gslerr.h" -'); - -sub airy_func { - my ($which, $doc) = @_; - pp_def($which, - Doc => $doc, - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code => < 0.' -); -airy_func( - 'gsl_sf_airy_Bi_scaled', - 'Scaled Airy Function Bi(x). Bi(x) for x < 0 and exp(+2/3 x^{3/2}) Bi(x) for x > 0.' -); -airy_func('gsl_sf_airy_Ai_deriv', 'Derivative Airy Function Ai`(x).'); -airy_func('gsl_sf_airy_Bi_deriv', 'Derivative Airy Function Bi`(x).'); -airy_func( - 'gsl_sf_airy_Ai_deriv_scaled', - 'Derivative Scaled Airy Function Ai(x). Ai`(x) for x < 0 and exp(+2/3 x^{3/2}) Ai`(x) for x > 0.' -); -airy_func( - 'gsl_sf_airy_Bi_deriv_scaled', - 'Derivative Scaled Airy Function Bi(x). Bi`(x) for x < 0 and exp(+2/3 x^{3/2}) Bi`(x) for x > 0.' -); - -pp_addpm({At=>'Bot'},<<'EOD'); -=head1 AUTHOR - -This file copyright (C) 1999 Christian Pellegrin -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. - -The GSL SF modules were written by G. Jungman. - -=cut - -EOD - -use strict; -use warnings; - -pp_def('gsl_sf_bessel_Jn', - GenericTypes => ['D'], - OtherPars =>'int n', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_bessel_Jn_e,($COMP(n),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Regular Bessel Function J_n(x).' - ); - -pp_def('gsl_sf_bessel_Jn_array', - GenericTypes => ['D'], - OtherPars =>'int s; IV n=>num', - Pars=>'double x(); double [o]y(num)', - Code =>' -GSLERR(gsl_sf_bessel_Jn_array,($COMP(s),$COMP(s)+$COMP(n)-1,$x(),$P(y))) -', - Doc =>'Array of Regular Bessel Functions J_{s}(x) to J_{s+n-1}(x).' - ); - -pp_def('gsl_sf_bessel_Yn', - GenericTypes => ['D'], - OtherPars =>'int n', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_bessel_Yn_e,($COMP(n),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'IrRegular Bessel Function Y_n(x).' - ); - -pp_def('gsl_sf_bessel_Yn_array', - GenericTypes => ['D'], - OtherPars =>'int s; IV n=>num', - Pars=>'double x(); double [o]y(num)', - Code =>' -GSLERR(gsl_sf_bessel_Yn_array,($COMP(s),$COMP(s)+$COMP(n)-1,$x(),$P(y))) -', - Doc =>'Array of Regular Bessel Functions Y_{s}(x) to Y_{s+n-1}(x).' - ); - -pp_def('gsl_sf_bessel_In', - GenericTypes => ['D'], - OtherPars =>'int n', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_bessel_In_e,($COMP(n),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Regular Modified Bessel Function I_n(x).' - ); - -pp_def('gsl_sf_bessel_I_array', - GenericTypes => ['D'], - OtherPars =>'int s; IV n=>num', - Pars=>'double x(); double [o]y(num)', - Code =>' -GSLERR(gsl_sf_bessel_In_array,($COMP(s),$COMP(s)+$COMP(n)-1,$x(),$P(y))) -', - Doc =>'Array of Regular Modified Bessel Functions I_{s}(x) to I_{s+n-1}(x).' - ); - -pp_def('gsl_sf_bessel_In_scaled', - GenericTypes => ['D'], - OtherPars =>'int n', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_bessel_In_scaled_e,($COMP(n),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Scaled Regular Modified Bessel Function exp(-|x|) I_n(x).' - ); - -pp_def('gsl_sf_bessel_In_scaled_array', - GenericTypes => ['D'], - OtherPars =>'int s; IV n=>num', - Pars=>'double x(); double [o]y(num)', - Code =>' -GSLERR(gsl_sf_bessel_In_scaled_array,($COMP(s),$COMP(s)+$COMP(n)-1,$x(),$P(y))) -', - Doc =>'Array of Scaled Regular Modified Bessel Functions exp(-|x|) I_{s}(x) to exp(-|x|) I_{s+n-1}(x).' - ); - -pp_def('gsl_sf_bessel_Kn', - GenericTypes => ['D'], - OtherPars =>'int n', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_bessel_Kn_e,($COMP(n),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'IrRegular Modified Bessel Function K_n(x).' - ); - -pp_def('gsl_sf_bessel_K_array', - GenericTypes => ['D'], - OtherPars =>'int s; IV n=>num', - Pars=>'double x(); double [o]y(num)', - Code =>' -GSLERR(gsl_sf_bessel_Kn_array,($COMP(s),$COMP(s)+$COMP(n)-1,$x(),$P(y))) -', - Doc =>'Array of IrRegular Modified Bessel Functions K_{s}(x) to K_{s+n-1}(x).' - ); - -pp_def('gsl_sf_bessel_Kn_scaled', - GenericTypes => ['D'], - OtherPars =>'int n', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_bessel_Kn_scaled_e,($COMP(n),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Scaled IrRegular Modified Bessel Function exp(-|x|) K_n(x).' - ); - -pp_def('gsl_sf_bessel_Kn_scaled_array', - GenericTypes => ['D'], - OtherPars =>'int s; IV n=>num', - Pars=>'double x(); double [o]y(num)', - Code =>' -GSLERR(gsl_sf_bessel_Kn_scaled_array,($COMP(s),$COMP(s)+$COMP(n)-1,$x(),$P(y))) -', - Doc =>'Array of Scaled IrRegular Modified Bessel Functions exp(-|x|) K_{s}(x) to exp(-|x|) K_{s+n-1}(x).' - ); - -pp_def('gsl_sf_bessel_jl', - GenericTypes => ['D'], - OtherPars =>'int n', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_bessel_jl_e,($COMP(n),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Regular Sphericl Bessel Function J_n(x).' - ); - -pp_def('gsl_sf_bessel_jl_array', - GenericTypes => ['D'], - OtherPars =>'int n=>num', - Pars=>'double x(); double [o]y(num)', - Code =>' -GSLERR(gsl_sf_bessel_jl_array,($COMP(n)-1,$x(),$P(y))) -', - Doc =>'Array of Spherical Regular Bessel Functions J_{0}(x) to J_{n-1}(x).' - ); - -pp_def('gsl_sf_bessel_yl', - GenericTypes => ['D'], - OtherPars =>'int n', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_bessel_yl_e,($COMP(n),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'IrRegular Spherical Bessel Function y_n(x).' - ); - -pp_def('gsl_sf_bessel_yl_array', - GenericTypes => ['D'], - OtherPars =>'int n=>num', - Pars=>'double x(); double [o]y(num)', - Code =>' -GSLERR(gsl_sf_bessel_yl_array,($COMP(n)-1,$x(),$P(y))) -', - Doc =>'Array of Regular Spherical Bessel Functions y_{0}(x) to y_{n-1}(x).' - ); - -pp_def('gsl_sf_bessel_il_scaled', - GenericTypes => ['D'], - OtherPars =>'int n', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_bessel_il_scaled_e,($COMP(n),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Scaled Regular Modified Spherical Bessel Function exp(-|x|) i_n(x).' - ); - -pp_def('gsl_sf_bessel_il_scaled_array', - GenericTypes => ['D'], - OtherPars =>'int n=>num', - Pars=>'double x(); double [o]y(num)', - Code =>' -GSLERR(gsl_sf_bessel_il_scaled_array,($COMP(n)-1,$x(),$P(y))) -', - Doc =>'Array of Scaled Regular Modified Spherical Bessel Functions exp(-|x|) i_{0}(x) to exp(-|x|) i_{n-1}(x).' - ); - -pp_def('gsl_sf_bessel_kl_scaled', - GenericTypes => ['D'], - OtherPars =>'int n', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_bessel_kl_scaled_e,($COMP(n),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Scaled IrRegular Modified Spherical Bessel Function exp(-|x|) k_n(x).' - ); - -pp_def('gsl_sf_bessel_kl_scaled_array', - GenericTypes => ['D'], - OtherPars =>'int n=>num', - Pars=>'double x(); double [o]y(num)', - Code =>' -GSLERR(gsl_sf_bessel_kl_scaled_array,($COMP(n)-1,$x(),$P(y))) -', - Doc =>'Array of Scaled IrRegular Modified Spherical Bessel Functions exp(-|x|) k_{s}(x) to exp(-|x|) k_{s+n-1}(x).' - ); - -pp_def('gsl_sf_bessel_Jnu', - GenericTypes => ['D'], - OtherPars =>'double n', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_bessel_Jnu_e,($COMP(n),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Regular Cylindrical Bessel Function J_nu(x).' - ); - -pp_def('gsl_sf_bessel_Ynu', - GenericTypes => ['D'], - OtherPars =>'double n', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_bessel_Ynu_e,($COMP(n),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'IrRegular Cylindrical Bessel Function J_nu(x).' - ); - -pp_def('gsl_sf_bessel_Inu_scaled', - GenericTypes => ['D'], - OtherPars =>'double n', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_bessel_Inu_scaled_e,($COMP(n),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Scaled Modified Cylindrical Bessel Function exp(-|x|) I_nu(x).' - ); - -pp_def('gsl_sf_bessel_Inu', - GenericTypes => ['D'], - OtherPars =>'double n', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_bessel_Inu_e,($COMP(n),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Modified Cylindrical Bessel Function I_nu(x).' - ); - -pp_def('gsl_sf_bessel_Knu_scaled', - GenericTypes => ['D'], - OtherPars =>'double n', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_bessel_Knu_scaled_e,($COMP(n),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Scaled Modified Cylindrical Bessel Function exp(-|x|) K_nu(x).' - ); - -pp_def('gsl_sf_bessel_Knu', - GenericTypes => ['D'], - OtherPars =>'double n', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_bessel_Knu_e,($COMP(n),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Modified Cylindrical Bessel Function K_nu(x).' - ); - -pp_def('gsl_sf_bessel_lnKnu', - GenericTypes => ['D'], - OtherPars =>'double n', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_bessel_lnKnu_e,($COMP(n),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Logarithm of Modified Cylindrical Bessel Function K_nu(x).' - ); - -pp_def('gsl_sf_clausen', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_clausen_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Clausen Integral. Cl_2(x) := Integrate[-Log[2 Sin[t/2]], {t,0,x}]' - ); - -pp_def('gsl_sf_hydrogenicR', - GenericTypes => ['D'], - OtherPars =>'int n; int l; double z', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_hydrogenicR_e,($COMP(n),$COMP(l),$COMP(z),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Normalized Hydrogenic bound states. Radial dipendence.' - ); - -pp_def('gsl_sf_coulomb_wave_FGp_array', - GenericTypes => ['D'], - OtherPars =>'double lam_min; IV kmax=>n; double eta', - Pars=>'double x(); double [o]fc(n); double [o]fcp(n); double [o]gc(n); double [o]gcp(n); int [o]ovfw(); double [o]fe(n); double [o]ge(n);', - Code =>' -int s; -s = gsl_sf_coulomb_wave_FGp_array($COMP(lam_min),$COMP(kmax),$COMP(eta),$x(),$P(fc),$P(fcp),$P(gc),$P(gcp),$P(fe),$P(ge)); -if (s==GSL_EOVRFLW) { -$ovfw()=1; -} -else {if (s) $CROAK("Error in gsl_sf_coulomb_wave_FGp_array: %s",gsl_strerror(s)); -else {$ovfw()=0;}} -', - Doc =>' Coulomb wave functions F_{lam_F}(eta,x), G_{lam_G}(eta,x) and their derivatives; lam_G := lam_F - k_lam_G. if ovfw is signaled then F_L(eta,x) = fc[k_L] * exp(fe) and similar. ' - ); - - -pp_def('gsl_sf_coulomb_wave_sphF_array', - GenericTypes => ['D'], - OtherPars =>'double lam_min; IV kmax=>n; double eta', - Pars=>'double x(); double [o]fc(n); int [o]ovfw(); double [o]fe(n);', - Code =>' -int s; -s = gsl_sf_coulomb_wave_sphF_array($COMP(lam_min),$COMP(kmax),$COMP(eta),$x(),$P(fc),$P(fe)); -if (s==GSL_EOVRFLW) { -$ovfw()=1; -} -else {if (s) $CROAK("Error in gsl_sf_coulomb_wave_sphF_array: %s",gsl_strerror(s)); -else {$ovfw()=0;}} -', - Doc =>' Coulomb wave function divided by the argument, F(xi, eta)/xi. This is the function which reduces to spherical Bessel functions in the limit eta->0. ' - ); - -pp_def('gsl_sf_coulomb_CL_e', - GenericTypes => ['D'], - Pars=>'double L(); double eta(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_coulomb_CL_e,($L(),$eta(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Coulomb wave function normalization constant. [Abramowitz+Stegun 14.1.8, 14.1.9].' - ); - -pp_def('gsl_sf_coupling_3j', - GenericTypes => ['L'], - Pars=>'ja(); jb(); jc(); ma(); mb(); mc(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_coupling_3j_e,($ja(),$jb(),$jc(),$ma(),$mb(),$mc(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'3j Symbols: (ja jb jc) over (ma mb mc).' - ); - -pp_def('gsl_sf_coupling_6j', - GenericTypes => ['L'], - Pars=>'ja(); jb(); jc(); jd(); je(); jf(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_coupling_6j_e,($ja(),$jb(),$jc(),$jd(),$je(),$jf(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'6j Symbols: (ja jb jc) over (jd je jf).' - ); - -pp_def('gsl_sf_coupling_9j', - GenericTypes => ['L'], - Pars=>'ja(); jb(); jc(); jd(); je(); jf(); jg(); jh(); ji(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_coupling_9j_e,($ja(),$jb(),$jc(),$jd(),$je(),$jf(),$jg(),$jh(),$ji(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'9j Symbols: (ja jb jc) over (jd je jf) over (jg jh ji).' - ); - -pp_def('gsl_sf_dawson', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_dawson_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Dawsons integral: Exp[-x^2] Integral[ Exp[t^2], {t,0,x}]' - ); - -pp_def('gsl_sf_debye_1', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_debye_1_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'D_n(x) := n/x^n Integrate[t^n/(e^t - 1), {t,0,x}]' - ); - -pp_def('gsl_sf_debye_2', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_debye_2_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'D_n(x) := n/x^n Integrate[t^n/(e^t - 1), {t,0,x}]' - ); - -pp_def('gsl_sf_debye_3', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_debye_3_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'D_n(x) := n/x^n Integrate[t^n/(e^t - 1), {t,0,x}]' - ); - -pp_def('gsl_sf_debye_4', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_debye_4_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'D_n(x) := n/x^n Integrate[t^n/(e^t - 1), {t,0,x}]' - ); - -pp_def('gsl_sf_dilog', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_dilog_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'/* Real part of DiLogarithm(x), for real argument. In Lewins notation, this is Li_2(x). Li_2(x) = - Re[ Integrate[ Log[1-s] / s, {s, 0, x}] ]' - ); - -pp_def('gsl_sf_complex_dilog', - GenericTypes => ['D'], - Pars=>'double r(); double t(); double [o]re(); double [o]im(); double [o]ere(); double [o]eim()', - Code =>' -gsl_sf_result re; -gsl_sf_result im; -GSLERR(gsl_sf_complex_dilog_e,($r(),$t(),&re,&im)) -$re() = re.val; -$ere() = re.err; -$im() = im.val; -$eim() = im.err; -', - Doc =>'DiLogarithm(z), for complex argument z = r Exp[i theta].' - ); - -pp_def('gsl_sf_multiply', - GenericTypes => ['D'], - Pars=>'double x(); double xx(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_multiply_e,($x(),$xx(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Multiplication.' - ); - -pp_def('gsl_sf_multiply_err', - GenericTypes => ['D'], - Pars=>'double x(); double xe(); double xx(); double xxe(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_multiply_err_e,($x(),$xe(),$xx(),$xxe(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Multiplication with associated errors.' - ); - -pp_def('gsl_sf_ellint_Kcomp', - GenericTypes => ['D'], - Pars=>'double k(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_ellint_Kcomp_e,($k(),GSL_PREC_DOUBLE,&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Legendre form of complete elliptic integrals K(k) = Integral[1/Sqrt[1 - k^2 Sin[t]^2], {t, 0, Pi/2}].' - ); - -pp_def('gsl_sf_ellint_Ecomp', - GenericTypes => ['D'], - Pars=>'double k(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_ellint_Ecomp_e,($k(),GSL_PREC_DOUBLE,&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Legendre form of complete elliptic integrals E(k) = Integral[ Sqrt[1 - k^2 Sin[t]^2], {t, 0, Pi/2}]' - ); - -pp_def('gsl_sf_ellint_F', - GenericTypes => ['D'], - Pars=>'double phi(); double k(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_ellint_F_e,($phi(),$k(),GSL_PREC_DOUBLE,&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Legendre form of incomplete elliptic integrals F(phi,k) = Integral[1/Sqrt[1 - k^2 Sin[t]^2], {t, 0, phi}]' - ); - -pp_def('gsl_sf_ellint_E', - GenericTypes => ['D'], - Pars=>'double phi(); double k(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_ellint_E_e,($phi(),$k(),GSL_PREC_DOUBLE,&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Legendre form of incomplete elliptic integrals E(phi,k) = Integral[ Sqrt[1 - k^2 Sin[t]^2], {t, 0, phi}]' - ); - -pp_def('gsl_sf_ellint_P', - GenericTypes => ['D'], - Pars=>'double phi(); double k(); double n(); - double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_ellint_P_e,($phi(),$k(),$n(),GSL_PREC_DOUBLE,&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Legendre form of incomplete elliptic integrals P(phi,k,n) = Integral[(1 + n Sin[t]^2)^(-1)/Sqrt[1 - k^2 Sin[t]^2], {t, 0, phi}]' - ); - -pp_def('gsl_sf_ellint_D', - GenericTypes => ['D'], - Pars=>'double phi(); double k(); - double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_ellint_D_e,($phi(),$k(),GSL_PREC_DOUBLE,&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Legendre form of incomplete elliptic integrals D(phi,k)' - ); - -pp_def('gsl_sf_ellint_RC', - GenericTypes => ['D'], - Pars=>'double x(); double yy(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_ellint_RC_e,($x(),$yy(),GSL_PREC_DOUBLE,&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Carlsons symmetric basis of functions RC(x,y) = 1/2 Integral[(t+x)^(-1/2) (t+y)^(-1)], {t,0,Inf}' - ); - -pp_def('gsl_sf_ellint_RD', - GenericTypes => ['D'], - Pars=>'double x(); double yy(); double z(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_ellint_RD_e,($x(),$yy(),$z(),GSL_PREC_DOUBLE,&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Carlsons symmetric basis of functions RD(x,y,z) = 3/2 Integral[(t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-3/2), {t,0,Inf}]' - ); - -pp_def('gsl_sf_ellint_RF', - GenericTypes => ['D'], - Pars=>'double x(); double yy(); double z(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_ellint_RF_e,($x(),$yy(),$z(),GSL_PREC_DOUBLE,&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Carlsons symmetric basis of functions RF(x,y,z) = 1/2 Integral[(t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-1/2), {t,0,Inf}]' - ); - -pp_def('gsl_sf_ellint_RJ', - GenericTypes => ['D'], - Pars=>'double x(); double yy(); double z(); double p(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_ellint_RJ_e,($x(),$yy(),$z(),$p(),GSL_PREC_DOUBLE,&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Carlsons symmetric basis of functions RJ(x,y,z,p) = 3/2 Integral[(t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-1/2) (t+p)^(-1), {t,0,Inf}]' - ); - -pp_def('gsl_sf_elljac', - GenericTypes => ['D'], - Pars=>'double u(); double m(); double [o]sn(); double [o]cn(); double [o]dn()', - Code =>' -if (gsl_sf_elljac_e($u(),$m(),$P(sn),$P(cn),$P(dn))) {$CROAK("Error in gsl_sf_elljac");}; -', - Doc =>'Jacobian elliptic functions sn, dn, cn by descending Landen transformations' - ); - -pp_def('gsl_sf_erfc', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_erfc_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Complementary Error Function erfc(x) := 2/Sqrt[Pi] Integrate[Exp[-t^2], {t,x,Infinity}]' - ); - -pp_def('gsl_sf_log_erfc', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_log_erfc_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Log Complementary Error Function' - ); - -pp_def('gsl_sf_erf', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_erf_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Error Function erf(x) := 2/Sqrt[Pi] Integrate[Exp[-t^2], {t,0,x}]' - ); - -pp_def('gsl_sf_erf_Z', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_erf_Z_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Z(x) : Abramowitz+Stegun 26.2.1' - ); - -pp_def('gsl_sf_erf_Q', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_erf_Q_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Q(x) : Abramowitz+Stegun 26.2.1' - ); - -pp_def('gsl_sf_exp', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_exp_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Exponential' - ); - -pp_def('gsl_sf_exprel_n', - GenericTypes => ['D'], - OtherPars => 'int n', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_exprel_n_e,($COMP(n),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'N-relative Exponential. exprel_N(x) = N!/x^N (exp(x) - Sum[x^k/k!, {k,0,N-1}]) = 1 + x/(N+1) + x^2/((N+1)(N+2)) + ... = 1F1(1,1+N,x)' - ); - -pp_def('gsl_sf_exp_err', - GenericTypes => ['D'], - Pars=>'double x(); double dx(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_exp_err_e,($x(),$dx(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Exponential of a quantity with given error.' - ); - -pp_def('gsl_sf_expint_E1', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_expint_E1_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'E_1(x) := Re[ Integrate[ Exp[-xt]/t, {t,1,Infinity}] ]' - ); - -pp_def('gsl_sf_expint_E2', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_expint_E2_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'E_2(x) := Re[ Integrate[ Exp[-xt]/t^2, {t,1,Infity}] ]' - ); - -pp_def('gsl_sf_expint_Ei', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_expint_Ei_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Ei(x) := PV Integrate[ Exp[-t]/t, {t,-x,Infinity}]' - ); - -pp_def('gsl_sf_Shi', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_Shi_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Shi(x) := Integrate[ Sinh[t]/t, {t,0,x}]' - ); - -pp_def('gsl_sf_Chi', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_Chi_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Chi(x) := Re[ M_EULER + log(x) + Integrate[(Cosh[t]-1)/t, {t,0,x}] ]' - ); - -pp_def('gsl_sf_expint_3', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_expint_3_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Ei_3(x) := Integral[ Exp[-t^3], {t,0,x}]' - ); - -pp_def('gsl_sf_Si', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_Si_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Si(x) := Integrate[ Sin[t]/t, {t,0,x}]' - ); - -pp_def('gsl_sf_Ci', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_Ci_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Ci(x) := -Integrate[ Cos[t]/t, {t,x,Infinity}]' - ); - -pp_def('gsl_sf_atanint', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_atanint_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'AtanInt(x) := Integral[ Arctan[t]/t, {t,0,x}]' - ); - -pp_def('gsl_sf_fermi_dirac_int', - GenericTypes => ['D'], - OtherPars => 'int j', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_fermi_dirac_int_e,($COMP(j),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Complete integral F_j(x) for integer j' - ); - -pp_def('gsl_sf_fermi_dirac_mhalf', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_fermi_dirac_mhalf_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Complete integral F_{-1/2}(x)' - ); - -pp_def('gsl_sf_fermi_dirac_half', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_fermi_dirac_half_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Complete integral F_{1/2}(x)' - ); - -pp_def('gsl_sf_fermi_dirac_3half', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_fermi_dirac_3half_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Complete integral F_{3/2}(x)' - ); - -pp_def('gsl_sf_fermi_dirac_inc_0', - GenericTypes => ['D'], - OtherPars => 'double b', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_fermi_dirac_inc_0_e,($x(),$COMP(b),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Incomplete integral F_0(x,b) = ln(1 + e^(b-x)) - (b-x)' - ); - -pp_def('gsl_sf_lngamma', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]s(); double [o]e()', - Code =>' -gsl_sf_result r; -double sgn; -GSLERR(gsl_sf_lngamma_sgn_e,($x(),&r,&sgn)) -$y() = r.val; -$e() = r.err; -$s() = sgn; -', - Doc =>'Log[Gamma(x)], x not a negative integer Uses real Lanczos method. Determines the sign of Gamma[x] as well as Log[|Gamma[x]|] for x < 0. So Gamma[x] = sgn * Exp[result_lg].' - ); - -pp_def('gsl_sf_gamma', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_gamma_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Gamma(x), x not a negative integer' - ); - -pp_def('gsl_sf_gammastar', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_gammastar_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Regulated Gamma Function, x > 0 Gamma^*(x) = Gamma(x)/(Sqrt[2Pi] x^(x-1/2) exp(-x)) = (1 + 1/(12x) + ...), x->Inf' - ); - -pp_def('gsl_sf_gammainv', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_gammainv_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'1/Gamma(x)' - ); - -pp_def('gsl_sf_lngamma_complex', - GenericTypes => ['D'], - Pars=>'double zr(); double zi(); double [o]x(); double [o]y(); double [o]xe(); double [o]ye()', - Code =>' -gsl_sf_result r; -gsl_sf_result ri; -GSLERR(gsl_sf_lngamma_complex_e,($zr(),$zi(),&r,&ri)) -$x() = r.val; -$xe() = r.err; -$y() = ri.val; -$ye() = ri.err; -', - Doc =>'Log[Gamma(z)] for z complex, z not a negative integer. Calculates: lnr = log|Gamma(z)|, arg = arg(Gamma(z)) in (-Pi, Pi]' - ); - -pp_def('gsl_sf_taylorcoeff', - GenericTypes => ['D'], - OtherPars => 'int n', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_taylorcoeff_e,($COMP(n),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'x^n / n!' - ); - -pp_def('gsl_sf_fact', - GenericTypes => ['L'], - Pars=>'x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_fact_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'n!' - ); - -pp_def('gsl_sf_doublefact', - GenericTypes => ['L'], - Pars=>'x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_doublefact_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'n!! = n(n-2)(n-4)' - ); - -pp_def('gsl_sf_lnfact', - GenericTypes => ['L'], - Pars=>'x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_lnfact_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'ln n!' - ); - -pp_def('gsl_sf_lndoublefact', - GenericTypes => ['L'], - Pars=>'x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_lndoublefact_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'ln n!!' - ); - -pp_def('gsl_sf_lnchoose', - GenericTypes => ['L'], - Pars=>'n(); m(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_lnchoose_e,($n(), $m(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'log(n choose m)' - ); - -pp_def('gsl_sf_choose', - GenericTypes => ['L'], - Pars=>'n(); m(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_choose_e,($n(), $m(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'n choose m' - ); - -pp_def('gsl_sf_lnpoch', - GenericTypes => ['D'], - OtherPars => 'double a', - Pars=>'double x(); double [o]y(); double [o]s(); double [o]e()', - Code =>' -gsl_sf_result r; -double sgn; -GSLERR(gsl_sf_lnpoch_sgn_e,($COMP(a),$x(),&r,&sgn)) -$y() = r.val; -$e() = r.err; -$s() = sgn; -', - Doc =>'Logarithm of Pochammer (Apell) symbol, with sign information. result = log( |(a)_x| ), sgn = sgn( (a)_x ) where (a)_x := Gamma[a + x]/Gamma[a]' - ); - -pp_def('gsl_sf_poch', - GenericTypes => ['D'], - OtherPars => 'double a', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_poch_e,($COMP(a),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Pochammer (Apell) symbol (a)_x := Gamma[a + x]/Gamma[x]' - ); - -pp_def('gsl_sf_pochrel', - GenericTypes => ['D'], - OtherPars => 'double a', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_pochrel_e,($COMP(a),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Relative Pochammer (Apell) symbol ((a,x) - 1)/x where (a,x) = (a)_x := Gamma[a + x]/Gamma[a]' - ); - -pp_def('gsl_sf_gamma_inc_Q', - GenericTypes => ['D'], - OtherPars => 'double a', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_gamma_inc_Q_e,($COMP(a),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Normalized Incomplete Gamma Function Q(a,x) = 1/Gamma(a) Integral[ t^(a-1) e^(-t), {t,x,Infinity} ]' - ); - -pp_def('gsl_sf_gamma_inc_P', - GenericTypes => ['D'], - OtherPars => 'double a', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_gamma_inc_P_e,($COMP(a),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Complementary Normalized Incomplete Gamma Function P(a,x) = 1/Gamma(a) Integral[ t^(a-1) e^(-t), {t,0,x} ]' - ); - -pp_def('gsl_sf_lnbeta', - GenericTypes => ['D'], - Pars=>'double a(); double b(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_lnbeta_e,($a(),$b(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Logarithm of Beta Function Log[B(a,b)]' - ); - -pp_def('gsl_sf_beta', - GenericTypes => ['D'], - OtherPars => '', - Pars=>'double a(); double b();double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_beta_e,($a(),$b(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Beta Function B(a,b)' - ); - -pp_def('gsl_sf_gegenpoly_n', - GenericTypes => ['D'], - OtherPars =>'int n; double lambda', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_gegenpoly_n_e,($COMP(n),$COMP(lambda), $x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Evaluate Gegenbauer polynomials.' - ); - -pp_def('gsl_sf_gegenpoly_array', - GenericTypes => ['D'], - OtherPars =>'int n=>num; double lambda', - Pars=>'double x(); double [o]y(num)', - Code =>' -GSLERR(gsl_sf_gegenpoly_array,($COMP(n)-1,$COMP(lambda),$x(),$P(y))) -', - Doc =>'Calculate array of Gegenbauer polynomials from 0 to n-1.' - ); - -pp_def('gsl_sf_hyperg_0F1', - GenericTypes => ['D'], - OtherPars =>'double c', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_hyperg_0F1_e,($COMP(c), $x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'/* Hypergeometric function related to Bessel functions 0F1[c,x] = Gamma[c] x^(1/2(1-c)) I_{c-1}(2 Sqrt[x]) Gamma[c] (-x)^(1/2(1-c)) J_{c-1}(2 Sqrt[-x])' - ); - -pp_def('gsl_sf_hyperg_1F1', - GenericTypes => ['D'], - OtherPars =>'double a; double b', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_hyperg_1F1_e,($COMP(a),$COMP(b), $x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Confluent hypergeometric function for integer parameters. 1F1[a,b,x] = M(a,b,x)' - ); - -pp_def('gsl_sf_hyperg_U', - GenericTypes => ['D'], - OtherPars =>'double a; double b', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_hyperg_U_e,($COMP(a),$COMP(b), $x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Confluent hypergeometric function for integer parameters. U(a,b,x)' - ); - -pp_def('gsl_sf_hyperg_2F1', - GenericTypes => ['D'], - OtherPars =>'double a; double b; double c', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_hyperg_2F1_e,($COMP(a),$COMP(b), $COMP(c),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Confluent hypergeometric function for integer parameters. 2F1[a,b,c,x]' - ); - -pp_def('gsl_sf_hyperg_2F1_conj', - GenericTypes => ['D'], - OtherPars =>'double a; double b; double c', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_hyperg_2F1_conj_e,($COMP(a),$COMP(b), $COMP(c),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Gauss hypergeometric function 2F1[aR + I aI, aR - I aI, c, x]' - ); - -pp_def('gsl_sf_hyperg_2F1_renorm', - GenericTypes => ['D'], - OtherPars =>'double a; double b; double c', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_hyperg_2F1_renorm_e,($COMP(a),$COMP(b), $COMP(c),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Renormalized Gauss hypergeometric function 2F1[a,b,c,x] / Gamma[c]' - ); - -pp_def('gsl_sf_hyperg_2F1_conj_renorm', - GenericTypes => ['D'], - OtherPars =>'double a; double b; double c', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_hyperg_2F1_conj_renorm_e,($COMP(a),$COMP(b), $COMP(c),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Renormalized Gauss hypergeometric function 2F1[aR + I aI, aR - I aI, c, x] / Gamma[c]' - ); - -pp_def('gsl_sf_hyperg_2F0', - GenericTypes => ['D'], - OtherPars =>'double a; double b', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_hyperg_2F0_e,($COMP(a),$COMP(b), $x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Mysterious hypergeometric function. The series representation is a divergent hypergeometric series. However, for x < 0 we have 2F0(a,b,x) = (-1/x)^a U(a,1+a-b,-1/x)' - ); - -pp_def('gsl_sf_laguerre_n', - GenericTypes => ['D'], - OtherPars =>'int n; double a', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_laguerre_n_e,($COMP(n),$COMP(a), $x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Evaluate generalized Laguerre polynomials.' - ); - -pp_def('gsl_sf_legendre_Pl', - GenericTypes => ['D'], - OtherPars =>'int l', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_legendre_Pl_e,($COMP(l),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'P_l(x)' - ); - -pp_def('gsl_sf_legendre_Pl_array', - GenericTypes => ['D'], - OtherPars =>'int l=>num', - Pars=>'double x(); double [o]y(num)', - Code =>' -GSLERR(gsl_sf_legendre_Pl_array,($COMP(l)-1,$x(),$P(y))) -', - Doc =>'P_l(x) from 0 to n-1.' - ); - -pp_def('gsl_sf_legendre_Ql', - GenericTypes => ['D'], - OtherPars =>'int l', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_legendre_Ql_e,($COMP(l),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Q_l(x)' - ); - -pp_def('gsl_sf_legendre_Plm', - GenericTypes => ['D'], - OtherPars =>'int l; int m', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_legendre_Plm_e,($COMP(l),$COMP(m),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'P_lm(x)' - ); - -pp_def('gsl_sf_legendre_array', - GenericTypes => ['D'], - OtherPars => 'char norm; int lmax; int csphase', - Pars => 'double x(); double [o]y(n=CALC($COMP(lmax)*($COMP(lmax)+1)/2+$COMP(lmax)+1)); double [t]work(wn=CALC(gsl_sf_legendre_array_n($COMP(lmax))))', - HandleBad => 1, - Code => <<'EOBC', - PDL_IF_BAD(if ( $ISBAD( x() ) ) { - loop(n) %{ $SETBAD ( y() ); %} - } else,) { - if($x()<-1||$x()>1) $CROAK("The input to gsl_sf_legendre_array must be abs(x)<=1, and you input %f",$x()); - gsl_sf_legendre_t norm; - switch ($COMP(norm)) { - case 'S': norm = GSL_SF_LEGENDRE_SCHMIDT; break; - case 'Y': norm = GSL_SF_LEGENDRE_SPHARM; break; - case 'N': norm = GSL_SF_LEGENDRE_FULL; break; - default: norm = GSL_SF_LEGENDRE_NONE; break; - } - GSLERR(gsl_sf_legendre_array_e,(norm, $COMP(lmax), $x(), $COMP(csphase), $P(work))); - loop(n) %{ $y() = $work(wn=>n); %} - } -EOBC - Doc => <<'EOD', -=for ref - -Calculate all normalized associated Legendre polynomials. - -=for usage - -$Plm = gsl_sf_legendre_array($x,'P',4,-1); - -The calculation is done for degree 0 <= l <= lmax and order 0 <= m <= l on the range abs(x)<=1. - -The parameter norm should be: - -=over 3 - -=item 'S' for Schmidt semi-normalized associated Legendre polynomials S_l^m(x), - -=item 'Y' for spherical harmonic associated Legendre polynomials Y_l^m(x), or - -=item 'N' for fully normalized associated Legendre polynomials N_l^m(x). - -=item 'P' (or any other) for unnormalized associated Legendre polynomials P_l^m(x), - -=back - -lmax is the maximum degree l. -csphase should be (-1) to INCLUDE the Condon-Shortley phase factor (-1)^m, or (+1) to EXCLUDE it. - -See L to get the value of C and C in the returned vector. -EOD - ); - -pp_def('gsl_sf_legendre_array_index', - OtherPars => 'int lmax', - Pars => 'int [o]l(n=CALC($COMP(lmax)*($COMP(lmax)+1)/2+$COMP(lmax)+1)); int [o]m(n)', - Code => q/ - int ell, em, index; - for (ell=0; ell<=$COMP(lmax); ell++){ - for (em=0; em<=ell; em++){ - index = gsl_sf_legendre_array_index(ell,em); - $l(n=>index)=ell; - $m(n=>index)=em; - } - }/, - Doc =>'=for ref - -Calculate the relation between gsl_sf_legendre_arrays index and l and m values. - -=for usage -($l,$m) = gsl_sf_legendre_array_index($lmax); - -Note that this function is called differently than the corresponding GSL function, to make it more useful for PDL: here you just input the maximum l (lmax) that was used in C and it calculates all l and m values.' - ); - -pp_def('gsl_sf_legendre_sphPlm', - GenericTypes => ['D'], - OtherPars =>'int l; int m', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_legendre_sphPlm_e,($COMP(l),$COMP(m),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'P_lm(x), normalized properly for use in spherical harmonics' - ); - -pp_def('gsl_sf_conicalP_half', - GenericTypes => ['D'], - OtherPars =>'double lambda', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_conicalP_half_e,($COMP(lambda),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Irregular Spherical Conical Function P^{1/2}_{-1/2 + I lambda}(x)' - ); - -pp_def('gsl_sf_conicalP_mhalf', - GenericTypes => ['D'], - OtherPars =>'double lambda', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_conicalP_mhalf_e,($COMP(lambda),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Regular Spherical Conical Function P^{-1/2}_{-1/2 + I lambda}(x)' - ); - -pp_def('gsl_sf_conicalP_0', - GenericTypes => ['D'], - OtherPars =>'double lambda', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_conicalP_0_e,($COMP(lambda),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Conical Function P^{0}_{-1/2 + I lambda}(x)' - ); - -pp_def('gsl_sf_conicalP_1', - GenericTypes => ['D'], - OtherPars =>'double lambda', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_conicalP_1_e,($COMP(lambda),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Conical Function P^{1}_{-1/2 + I lambda}(x)' - ); - -pp_def('gsl_sf_conicalP_sph_reg', - GenericTypes => ['D'], - OtherPars =>'int l; double lambda', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_conicalP_sph_reg_e,($COMP(l),$COMP(lambda),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Regular Spherical Conical Function P^{-1/2-l}_{-1/2 + I lambda}(x)' - ); - -pp_def('gsl_sf_conicalP_cyl_reg_e', - GenericTypes => ['D'], - OtherPars =>'int m; double lambda', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_conicalP_cyl_reg_e,($COMP(m),$COMP(lambda),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Regular Cylindrical Conical Function P^{-m}_{-1/2 + I lambda}(x)' - ); - -pp_def('gsl_sf_legendre_H3d', - GenericTypes => ['D'], - OtherPars =>'int l; double lambda; double eta', - Pars=>'double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_legendre_H3d_e,($COMP(l),$COMP(lambda),$COMP(eta),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'lth radial eigenfunction of the Laplacian on the 3-dimensional hyperbolic space.' - ); - -pp_def('gsl_sf_legendre_H3d_array', - GenericTypes => ['D'], - OtherPars =>'int l=>num; double lambda; double eta', - Pars=>'double [o]y(num)', - Code =>' -GSLERR(gsl_sf_legendre_H3d_array,($COMP(l)-1,$COMP(lambda),$COMP(eta),$P(y))) -', - Doc =>'Array of H3d(ell), for l from 0 to n-1.' - ); - -pp_def('gsl_sf_log', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_log_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Provide a logarithm function with GSL semantics.' - ); - - -pp_def('gsl_sf_complex_log', - GenericTypes => ['D'], - Pars=>'double zr(); double zi(); double [o]x(); double [o]y(); double [o]xe(); double [o]ye()', - Code =>' -gsl_sf_result r; -gsl_sf_result ri; -GSLERR(gsl_sf_complex_log_e,($zr(),$zi(),&r,&ri)) -$x() = r.val; -$xe() = r.err; -$y() = ri.val; -$ye() = ri.err; -', - Doc =>'Complex Logarithm exp(lnr + I theta) = zr + I zi Returns argument in [-pi,pi].' - ); - -pp_def('gsl_poly_eval', - GenericTypes => ['D'], - Pars=>'double x(); double c(m); double [o]y()', - Code =>' -$y() = gsl_poly_eval($P(c),$SIZE(m),$x()); -', - Doc =>'c[0] + c[1] x + c[2] x^2 + ... + c[m-1] x^(m-1)' - ); - -pp_def('gsl_sf_pow_int', - GenericTypes => ['D'], - OtherPars => 'int n', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_pow_int_e,($x(),$COMP(n),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Calculate x^n.' - ); - -pp_def('gsl_sf_psi', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_psi_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Di-Gamma Function psi(x).' - ); - -pp_def('gsl_sf_psi_1piy', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_psi_1piy_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Di-Gamma Function Re[psi(1 + I y)]' - ); - -pp_def('gsl_sf_psi_n', - GenericTypes => ['D'], - OtherPars => 'int n', - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_psi_n_e,($COMP(n),$x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Poly-Gamma Function psi^(n)(x)' - ); - -pp_def('gsl_sf_synchrotron_1', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_synchrotron_1_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'First synchrotron function: synchrotron_1(x) = x Integral[ K_{5/3}(t), {t, x, Infinity}]' - ); - -pp_def('gsl_sf_synchrotron_2', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_synchrotron_2_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Second synchroton function: synchrotron_2(x) = x * K_{2/3}(x)' - ); - -pp_def('gsl_sf_transport_2', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_transport_2_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'J(2,x)' - ); - -pp_def('gsl_sf_transport_3', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_transport_3_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'J(3,x)' - ); - -pp_def('gsl_sf_transport_4', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_transport_4_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'J(4,x)' - ); - -pp_def('gsl_sf_transport_5', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_transport_5_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'J(5,x)' - ); - -pp_def('gsl_sf_sin', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_sin_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Sin(x) with GSL semantics.' - ); - -pp_def('gsl_sf_cos', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_cos_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Cos(x) with GSL semantics.' - ); - -pp_def('gsl_sf_hypot', - GenericTypes => ['D'], - Pars=>'double x(); double xx(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_hypot_e,($x(),$xx(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Hypot(x,xx) with GSL semantics.' - ); - -pp_def('gsl_sf_complex_sin', - GenericTypes => ['D'], - Pars=>'double zr(); double zi(); double [o]x(); double [o]y(); double [o]xe(); double [o]ye()', - Code =>' -gsl_sf_result r; -gsl_sf_result ri; -GSLERR(gsl_sf_complex_sin_e,($zr(),$zi(),&r,&ri)) -$x() = r.val; -$xe() = r.err; -$y() = ri.val; -$ye() = ri.err; -', - Doc =>'Sin(z) for complex z' - ); - -pp_def('gsl_sf_complex_cos', - GenericTypes => ['D'], - Pars=>'double zr(); double zi(); double [o]x(); double [o]y(); double [o]xe(); double [o]ye()', - Code =>' -gsl_sf_result r; -gsl_sf_result ri; -GSLERR(gsl_sf_complex_cos_e,($zr(),$zi(),&r,&ri)) -$x() = r.val; -$xe() = r.err; -$y() = ri.val; -$ye() = ri.err; -', - Doc =>'Cos(z) for complex z' - ); - -pp_def('gsl_sf_complex_logsin', - GenericTypes => ['D'], - Pars=>'double zr(); double zi(); double [o]x(); double [o]y(); double [o]xe(); double [o]ye()', - Code =>' -gsl_sf_result r; -gsl_sf_result ri; -GSLERR(gsl_sf_complex_logsin_e,($zr(),$zi(),&r,&ri)) -$x() = r.val; -$xe() = r.err; -$y() = ri.val; -$ye() = ri.err; -', - Doc =>'Log(Sin(z)) for complex z' - ); - -pp_def('gsl_sf_lnsinh', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_lnsinh_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Log(Sinh(x)) with GSL semantics.' - ); - -pp_def('gsl_sf_lncosh', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_lncosh_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Log(Cos(x)) with GSL semantics.' - ); - -pp_def('gsl_sf_polar_to_rect', - GenericTypes => ['D'], - Pars=>'double r(); double t(); double [o]x(); double [o]y(); double [o]xe(); double [o]ye()', - Code =>' -gsl_sf_result r; -gsl_sf_result ri; -GSLERR(gsl_sf_polar_to_rect,($r(),$t(),&r,&ri)) -$x() = r.val; -$xe() = r.err; -$y() = ri.val; -$ye() = ri.err; -', - Doc =>'Convert polar to rectlinear coordinates.' - ); - -pp_def('gsl_sf_rect_to_polar', - GenericTypes => ['D'], - Pars=>'double x(); double y(); double [o]r(); double [o]t(); double [o]re(); double [o]te()', - Code =>' -gsl_sf_result r; -gsl_sf_result ri; -GSLERR(gsl_sf_rect_to_polar,($x(),$y(),&r,&ri)) -$r() = r.val; -$re() = r.err; -$t() = ri.val; -$te() = ri.err; -', - Doc =>'Convert rectlinear to polar coordinates. return argument in range [-pi, pi].' - ); - -pp_def('gsl_sf_angle_restrict_symm', - GenericTypes => ['D'], - Pars=>'double [o]y();', - Code =>' -GSLERR(gsl_sf_angle_restrict_symm_e,($P(y))) -', - Doc =>'Force an angle to lie in the range (-pi,pi].' - ); - -pp_def('gsl_sf_angle_restrict_pos', - GenericTypes => ['D'], - Pars=>'double [o]y();', - Code =>' -GSLERR(gsl_sf_angle_restrict_pos_e,($P(y))) -', - Doc =>'Force an angle to lie in the range [0,2 pi).' - ); - -pp_def('gsl_sf_sin_err', - GenericTypes => ['D'], - Pars=>'double x(); double dx(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_sin_err_e,($x(),$dx(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Sin(x) for quantity with an associated error.' - ); - -pp_def('gsl_sf_cos_err', - GenericTypes => ['D'], - Pars=>'double x(); double dx(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_cos_err_e,($x(),$dx(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Cos(x) for quantity with an associated error.' - ); - -pp_def('gsl_sf_zeta', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_zeta_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Riemann Zeta Function zeta(x) = Sum[ k^(-s), {k,1,Infinity} ], s != 1.0' - ); - -pp_def('gsl_sf_hzeta', - GenericTypes => ['D'], - OtherPars =>'double q', - Pars=>'double s(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_hzeta_e,($s(), $COMP(q), &r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Hurwicz Zeta Function zeta(s,q) = Sum[ (k+q)^(-s), {k,0,Infinity} ]' - ); - -pp_def('gsl_sf_eta', - GenericTypes => ['D'], - Pars=>'double x(); double [o]y(); double [o]e()', - Code =>' -gsl_sf_result r; -GSLERR(gsl_sf_eta_e,($x(),&r)) -$y() = r.val; -$e() = r.err; -', - Doc =>'Eta Function eta(s) = (1-2^(1-s)) zeta(s)' - ); - -pp_done(); diff --git a/Libtmp/GSL/lib/PDL/Stats/Distr.pd b/Libtmp/GSL/lib/PDL/Stats/Distr.pd deleted file mode 100644 index 31fe195ea..000000000 --- a/Libtmp/GSL/lib/PDL/Stats/Distr.pd +++ /dev/null @@ -1,815 +0,0 @@ -use strict; -use warnings; -use PDL::Types 'types'; - -my $F = [map $_->ppsym, grep !$_->integer && $_->real, types()]; - -pp_addpm({At=>'Top'}, <<'EOD'); - -use strict; -use warnings; - -use Carp; -use PDL::LiteF; - -my $DEV = ($^O =~ /win/i)? '/png' : '/xs'; - -=head1 NAME - -PDL::Stats::Distr -- parameter estimations and probability density functions for distributions. - -=head1 DESCRIPTION - -Parameter estimate is maximum likelihood estimate when there is closed form estimate, otherwise it is method of moments estimate. - -=head1 SYNOPSIS - - use PDL::LiteF; - use PDL::Stats::Distr; - - # do a frequency (probability) plot with fitted normal curve - my $data = grandom(100)->abs; - - my ($xvals, $hist) = $data->hist; - - # turn frequency into probability - $hist /= $data->nelem; - - # get maximum likelihood estimates of normal curve parameters - my ($m, $v) = $data->mle_gaussian(); - - # fitted normal curve probabilities - my $p = $xvals->pdf_gaussian($m, $v); - - use PDL::Graphics::PGPLOT::Window; - my $win = pgwin( Dev=>"/xs" ); - - $win->bin( $hist ); - $win->hold; - $win->line( $p, {COLOR=>2} ); - $win->close; - -Or, play with different distributions with B :) - - $data->plot_distr( 'gaussian', 'lognormal' ); - -=cut - -EOD - -pp_addhdr(' -#include -#include -#include -'); - -pp_def('mme_beta', - Pars => 'a(n); float+ [o]alpha(); float+ [o]beta()', - GenericTypes => $F, - HandleBad => 1, - Code => ' - $GENERIC(alpha) sa = 0, a2 = 0; - PDL_Indx N = PDL_IF_BAD(0,$SIZE(n)); - loop (n) %{ - PDL_IF_BAD(if ($ISBAD($a())) continue;,) - sa += $a(); - a2 += pow($a(), 2); - PDL_IF_BAD(N++;,) - %} - if (N < 1) { $SETBAD(alpha()); $SETBAD(beta()); continue; } - $GENERIC(alpha) m = sa / N; - $GENERIC(alpha) v = a2 / N - pow(m, 2); - $alpha() = m * ( m * (1 - m) / v - 1 ); - $beta() = (1 - m) * ( m * (1 - m) / v - 1 ); - ', - Doc => ' -=for usage - - my ($a, $b) = $data->mme_beta(); - -=for ref - -beta distribution. pdf: f(x; a,b) = 1/B(a,b) x^(a-1) (1-x)^(b-1) - -=cut -', - -); - -pp_def('pdf_beta', - Pars => 'x(); a(); b(); float+ [o]p()', - GenericTypes => $F, - HandleBad => 1, - Code => ' -PDL_IF_BAD(if ($ISBAD($x()) || $ISBAD($a()) || $ISBAD($b())) { $SETBAD($p()); continue; },) -if (!($x()>=0 && $x()<=1)) $CROAK("x out of range [0,1]"); -double B_1 = 1 / gsl_sf_beta( $a(), $b() ); -$p() = B_1 * pow($x(), $a()-1) * pow(1-$x(), $b()-1); -', - Doc => ' -=for ref - -probability density function for beta distribution. x defined on [0,1]. - -=cut -', -); - -pp_def('mme_binomial', - Pars => 'a(n); int [o]n_(); float+ [o]p()', - GenericTypes => $F, - HandleBad => 1, - Code => ' - $GENERIC(p) sa, a2, m, v; - sa = 0; a2 = 0; - PDL_Indx N = PDL_IF_BAD(0,$SIZE(n)); - loop (n) %{ - PDL_IF_BAD(if ($ISBAD($a())) continue;,) - sa += $a(); - a2 += pow($a(), 2); - PDL_IF_BAD(N++;,) - %} - if (N < 1) { $SETBAD(n_()); $SETBAD(p()); continue; } - m = sa / N; - v = a2 / N - pow(m, 2); - $p() = 1 - v/m; - $n_() = m / $p() >= 0? (int) (m / $p() + .5) : (int) (m / $p() - .5); - $p() = m / $n_(); - ', - Doc => ' -=for usage - - my ($n, $p) = $data->mme_binomial; - -=for ref - -binomial distribution. pmf: f(k; n,p) = (n k) p^k (1-p)^(n-k) for k = 0,1,2..n -', -); - -pp_def('pmf_binomial', - Pars => 'ushort x(); ushort n(); p(); float+ [o]out()', - GenericTypes => $F, - HandleBad => 1, - Code => ' -PDL_IF_BAD(if ($ISBAD($x()) || $ISBAD($n()) || $ISBAD($p())) { - $SETBAD( $out() ); - continue; -},) -$GENERIC(out) bc = gsl_sf_choose($n(), $x()); -$out() = bc * pow($p(), $x()) * pow(1-$p(), $n() - $x()); -', - Doc => ' -=for ref - -probability mass function for binomial distribution. -', - -); - -pp_def('mle_exp', - Pars => 'a(n); float+ [o]l()', - GenericTypes => $F, - HandleBad => 1, - Code => ' - $GENERIC(l) sa = 0; - PDL_Indx N = PDL_IF_BAD(0,$SIZE(n)); - loop (n) %{ - PDL_IF_BAD(if ($ISBAD($a())) continue;,) - sa += $a(); - PDL_IF_BAD(N++;,) - %} - PDL_IF_BAD(if (sa <= 0) { $SETBAD(l()); continue; },) - $l() = N / sa; -', - Doc => ' -=for usage - - my $lamda = $data->mle_exp; - -=for ref - -exponential distribution. mle same as method of moments estimate. -', -); - -pp_def('pdf_exp', - Pars => 'x(); l(); float+ [o]p()', - GenericTypes => $F, - HandleBad => 1, - Code => ' -PDL_IF_BAD(if ( $ISBAD($x()) || $ISBAD($l()) ) { $SETBAD( $p() ); continue; },) -$p() = $l() * exp( -1 * $l() * $x() ); -', - Doc => ' -=for ref - -probability density function for exponential distribution. -', -); - -pp_def('mme_gamma', - Pars => 'a(n); float+ [o]shape(); float+ [o]scale()', - GenericTypes => $F, - HandleBad => 1, - Code => ' - $GENERIC(shape) sa = 0, a2 = 0; - PDL_Indx N = PDL_IF_BAD(0,$SIZE(n)); - loop (n) %{ - PDL_IF_BAD(if ($ISBAD($a())) continue;,) - sa += $a(); - a2 += pow($a(), 2); - PDL_IF_BAD(N++;,) - %} - PDL_IF_BAD(if (N < 1) { $SETBAD(shape()); $SETBAD(scale()); continue; },) - $GENERIC(shape) m = sa / N; - $GENERIC(shape) v = a2 / N - pow(m, 2); - $shape() = pow(m, 2) / v; - $scale() = v / m; -', - Doc => ' -=for usage - - my ($shape, $scale) = $data->mme_gamma(); - -=for ref - -two-parameter gamma distribution -', -); - -pp_def('pdf_gamma', - Pars => 'x(); a(); t(); float+ [o]p()', - GenericTypes => $F, - HandleBad => 1, - Code => ' -PDL_IF_BAD(if ($ISBAD($x()) || $ISBAD($a()) || $ISBAD($t())) { $SETBAD( $p() ); continue; },) -double g = gsl_sf_gamma( $a() ); -$p() = pow($x(), $a()-1) * exp(-1*$x() / $t()) / (pow($t(), $a()) * g); -', - Doc => ' -=for ref - -probability density function for two-parameter gamma distribution. -', -); - -pp_def('mle_gaussian', - Pars => 'a(n); float+ [o]m(); float+ [o]v()', - GenericTypes => $F, - HandleBad => 1, - Code => ' - $GENERIC(m) sa = 0, a2 = 0; - PDL_Indx N = PDL_IF_BAD(0,$SIZE(n)); - loop (n) %{ - PDL_IF_BAD(if ($ISBAD($a())) continue;,) - sa += $a(); - a2 += pow($a(), 2); - PDL_IF_BAD(N++;,) - %} - if (N < 1) { $SETBAD(m()); $SETBAD(v()); continue; } - $m() = sa / N; - $v() = a2 / N - pow($m(),2); -', - Doc => ' -=for usage - - my ($m, $v) = $data->mle_gaussian(); - -=for ref - -gaussian aka normal distribution. same results as $data->average and $data->var. mle same as method of moments estimate. -', -); - -pp_def('pdf_gaussian', - Pars => 'x(); m(); v(); float+ [o]p()', - GenericTypes => $F, - HandleBad => 1, - Code => ' -PDL_IF_BAD(if ($ISBAD($x()) || $ISBAD($m()) || $ISBAD($v())) { $SETBAD( $p() ); continue; },) -$p() = 1 / sqrt($v() * 2 * M_PI) - * exp( -1 * pow($x() - $m(), 2) / (2*$v()) ); -', - Doc => ' -=for ref - -probability density function for gaussian distribution. -', -); - -pp_def('mle_geo', - Pars => 'a(n); float+ [o]p();', - GenericTypes => $F, - HandleBad => 1, - Code => ' - $GENERIC(p) sa = 0; - PDL_Indx N = PDL_IF_BAD(0,$SIZE(n)); - loop (n) %{ - PDL_IF_BAD(if ($ISBAD($a())) continue;,) - sa += $a(); - PDL_IF_BAD(N++;,) - %} - if (N < 1) { $SETBAD(p()); continue; } - $p() = 1 / (1 + sa/N); -', - Doc => ' -=for ref - -geometric distribution. mle same as method of moments estimate. -', -); - -pp_def('pmf_geo', - Pars => 'ushort x(); p(); float+ [o]out()', - GenericTypes => $F, - HandleBad => 1, - Code => ' -PDL_IF_BAD(if ($ISBAD($x()) || $ISBAD($p())) { $SETBAD( $out() ); continue; },) -$out() = pow(1-$p(), $x()) * $p(); -', - Doc => ' -=for ref - -probability mass function for geometric distribution. x >= 0. -', -); - -pp_def('mle_geosh', - Pars => 'a(n); float+ [o]p();', - GenericTypes => $F, - HandleBad => 1, - Code => ' - $GENERIC(p) sa = 0; - PDL_Indx N = PDL_IF_BAD(0,$SIZE(n)); - loop (n) %{ - PDL_IF_BAD(if ($ISBAD($a())) continue;,) - sa += $a(); - PDL_IF_BAD(N++;,) - %} - if (sa <= 0) { $SETBAD(p()); continue; } - $p() = N / sa; - ', - Doc => ' -=for ref - -shifted geometric distribution. mle same as method of moments estimate. -', -); - -pp_def('pmf_geosh', - Pars => 'ushort x(); p(); float+ [o]out()', - GenericTypes => $F, - HandleBad => 1, - Code => ' -PDL_IF_BAD(if ($ISBAD($x()) || $ISBAD($p())) { $SETBAD( $out() ); continue; },) -if ( $x() < 1 ) { $CROAK( "x >= 1 please" ); } -$out() = pow(1-$p(), $x()-1) * $p(); -', - Doc => ' -=for ref - -probability mass function for shifted geometric distribution. x >= 1. -', -); - -pp_def('mle_lognormal', - Pars => 'a(n); float+ [o]m(); float+ [o]v()', - GenericTypes => $F, - HandleBad => 1, - Code => ' - $GENERIC(m) sa = 0, a2 = 0; - PDL_Indx N = PDL_IF_BAD(0,$SIZE(n)); - loop (n) %{ - PDL_IF_BAD(if ($ISBAD($a())) continue;,) - sa += log($a()); - PDL_IF_BAD(N++;,) - %} - if (N < 1) { $SETBAD(m()); $SETBAD(v()); continue; } - $m() = sa / N; - loop (n) %{ - PDL_IF_BAD(if ($ISBAD($a())) continue;,) - a2 += pow(log($a()) - $m(), 2); - %} - $v() = a2 / N; -', - Doc => ' -=for usage - - my ($m, $v) = $data->mle_lognormal(); - -=for ref - -lognormal distribution. maximum likelihood estimation. -', -); - -pp_def('mme_lognormal', - Pars => 'a(n); float+ [o]m(); float+ [o]v()', - GenericTypes => $F, - HandleBad => 1, - Code => ' - $GENERIC(m) sa = 0, a2 = 0; - PDL_Indx N = PDL_IF_BAD(0,$SIZE(n)); - loop (n) %{ - PDL_IF_BAD(if ($ISBAD($a())) continue;,) - sa += $a(); - a2 += pow($a(), 2); - PDL_IF_BAD(N++;,) - %} - if (N < 1) { $SETBAD(m()); $SETBAD(v()); continue; } - $m() = 2 * log(sa / N) - 1/2 * log( a2 / N ); - $v() = log( a2 / N ) - 2 * log( sa / N ); - ', - Doc => ' -=for usage - - my ($m, $v) = $data->mme_lognormal(); - -=for ref - -lognormal distribution. method of moments estimation. -', -); - -pp_def('pdf_lognormal', - Pars => 'x(); m(); v(); float+ [o]p()', - GenericTypes => $F, - HandleBad => 1, - Code => ' -PDL_IF_BAD(if ( $ISBAD($x()) || $ISBAD($m()) || $ISBAD($v()) ) { - $SETBAD( $p() ); - continue; -},) -if (!($x() > 0 && $v() > 0)) $CROAK( "x and v > 0 please" ); -$p() = 1 / ($x() * sqrt($v() * 2 * M_PI)) - * exp( -1 * pow(log($x()) - $m(), 2) / (2*$v()) ); -', - Doc => ' -=for ref - -probability density function for lognormal distribution. x > 0. v > 0. -', -); - -pp_def('mme_nbd', - Pars => 'a(n); float+ [o]r(); float+ [o]p()', - GenericTypes => $F, - HandleBad => 1, - Code => ' - $GENERIC(p) sa = 0, a2 = 0; - PDL_Indx N = PDL_IF_BAD(0,$SIZE(n)); - loop (n) %{ - PDL_IF_BAD(if ($ISBAD($a())) continue;,) - sa += $a(); - a2 += pow($a(), 2); - PDL_IF_BAD(N++;,) - %} - PDL_IF_BAD(if (N < 1) { $SETBAD(r()); $SETBAD(p()); continue; },) - $GENERIC(p) m = sa / N; - $GENERIC(p) v = a2 / N - pow(m, 2); - $r() = pow(m, 2) / (v - m); - $p() = m / v; -', - Doc => ' -=for usage - - my ($r, $p) = $data->mme_nbd(); - -=for ref - -negative binomial distribution. pmf: f(x; r,p) = (x+r-1 r-1) p^r (1-p)^x for x=0,1,2... -', -); - -pp_def('pmf_nbd', - Pars => 'ushort x(); r(); p(); float+ [o]out()', - GenericTypes => $F, - HandleBad => 1, - Code => ' -PDL_IF_BAD(if ( $ISBAD($x()) || $ISBAD($r()) || $ISBAD($p()) ) { - $SETBAD( $out() ); - continue; -},) - $GENERIC(out) nbc - = gsl_sf_gamma($x()+$r()) / (gsl_sf_fact($x()) * gsl_sf_gamma($r())); - $out() = nbc * pow($p(),$r()) * pow(1-$p(), $x()); -', - Doc => ' -=for ref - -probability mass function for negative binomial distribution. -', -); - -pp_def('mme_pareto', - Pars => 'a(n); float+ [o]k(); float+ [o]xm()', - GenericTypes => $F, - HandleBad => 1, - Code => ' - $GENERIC(xm) sa = 0, min = $a(n=>0); - PDL_Indx N = PDL_IF_BAD(0,$SIZE(n)); - loop (n) %{ - PDL_IF_BAD(if ($ISBAD($a())) continue;,) - sa += $a(); - if (min > $a()) - min = $a(); - PDL_IF_BAD(N++;,) - %} - if (min <= 0) $CROAK("min <= 0!"); - $k() = (sa - min) / ( N*( sa/N - min ) ); - $xm() = (N * $k() - 1) * min / ( N * $k() ); -', - Doc => ' -=for usage - - my ($k, $xm) = $data->mme_pareto(); - -=for ref - -pareto distribution. pdf: f(x; k,xm) = k xm^k / x^(k+1) for x >= xm > 0. -', -); - -pp_def('pdf_pareto', - Pars => 'x(); k(); xm(); float+ [o]p()', - GenericTypes => $F, - HandleBad => 1, - Code => ' -PDL_IF_BAD(if ( $ISBAD($x()) || $ISBAD($k()) || $ISBAD($xm()) ) { - $SETBAD( $p() ); - continue; -},) -if (!($xm() > 0 && $x() >= $xm() )) $CROAK("x >= xm > 0 please"); -$p() = $k() * pow($xm(),$k()) / pow($x(), $k()+1); -', - Doc => ' -=for ref - -probability density function for pareto distribution. x >= xm > 0. -', -); - -pp_def('mle_poisson', - Pars => 'a(n); float+ [o]l();', - GenericTypes => $F, - HandleBad => 1, - Code => ' - $GENERIC(l) sa = 0; - PDL_Indx N = PDL_IF_BAD(0,$SIZE(n)); - loop (n) %{ - PDL_IF_BAD(if ($ISBAD($a())) continue;,) - sa += $a(); - PDL_IF_BAD(N++;,) - %} - if (N < 1) { $SETBAD(l()); continue; } - $l() = sa / N; -', - Doc => ' -=for usage - - my $lamda = $data->mle_poisson(); - -=for ref - -poisson distribution. pmf: f(x;l) = e^(-l) * l^x / x! -', -); - -pp_def('pmf_poisson', - Pars => 'x(); l(); float+ [o]p()', - GenericTypes => $F, - HandleBad => 1, - Code => q{ -PDL_IF_BAD(if ($ISBAD($x()) || $ISBAD($l())) { $SETBAD($p()); continue; },) -if ($x() < 0) { - $p() = 0; -} -else if ($x() < GSL_SF_FACT_NMAX / 2) { - /* Exact formula */ - $p() = exp( -1 * $l()) * pow($l(),$x()) / gsl_sf_fact( (unsigned int) $x() ); -} -else { - /* Use Stirling's approximation. See - * http://en.wikipedia.org/wiki/Stirling%27s_approximation - */ - double log_p = $x() - $l() + $x() * log($l() / $x()) - - 0.5 * log(2*M_PI * $x()) - 1. / 12. / $x() - + 1 / 360. / $x()/$x()/$x() - 1. / 1260. / $x()/$x()/$x()/$x()/$x(); - $p() = exp(log_p); -} -}, - Doc => q{ -=for ref - -Probability mass function for poisson distribution. Uses Stirling's formula for x > 85. -}, -); - -pp_def('pmf_poisson_stirling', - Pars => 'x(); l(); [o]p()', - GenericTypes => $F, - HandleBad => 1, - Code => q{ -PDL_IF_BAD(if ($ISBAD($x()) || $ISBAD($l())) { $SETBAD($p()); continue; },) -if ($x() < 0) { - $p() = 0; -} -else if ($x() == 0) { - $p() = exp(-$l()); -} -else { - /* Use Stirling's approximation. See - * http://en.wikipedia.org/wiki/Stirling%27s_approximation - */ - double log_p = $x() - $l() + $x() * log($l() / $x()) - - 0.5 * log(2*M_PI * $x()) - 1. / 12. / $x() - + 1 / 360. / $x()/$x()/$x() - 1. / 1260. / $x()/$x()/$x()/$x()/$x(); - $p() = exp(log_p); -} -}, - Doc => q{ -=for ref - -Probability mass function for poisson distribution. Uses Stirling's formula for all values of the input. See http://en.wikipedia.org/wiki/Stirling's_approximation for more info. -}, -); - -pp_def('pmf_poisson_factorial', - Pars => 'ushort x(); l(); float+ [o]p()', - GenericTypes => $F, - HandleBad => 1, - Code => q{ -PDL_IF_BAD(if ($ISBAD($x()) || $ISBAD($l())) { $SETBAD($p()); continue; },) -if ($x() < GSL_SF_FACT_NMAX) { - $p() = exp( -1 * $l()) * pow($l(),$x()) / gsl_sf_fact( $x() ); -} else { - /* bail out */ - $p() = 0; -} -}, - PMCode => pp_line_numbers(__LINE__, <<'EOF'), -sub PDL::pmf_poisson_factorial { - my ($x, $l) = @_; - my $pdlx = PDL->topdl($x); - croak "Does not support input greater than 170. Please use pmf_poisson or pmf_poisson_stirling instead." - if any($pdlx >= 170); - PDL::_pmf_poisson_factorial_int($pdlx, $l, my $p = PDL->null); - $p; -} -EOF - Doc => <<'EOD', -=for ref - -Probability mass function for poisson distribution. Input is limited to x < 170 to avoid gsl_sf_fact() overflow. -EOD -); - -pp_addpm {At=>'Bot'}, pp_line_numbers(__LINE__, <<'EOD'); - -=head2 plot_distr - -=for ref - -Plots data distribution. When given specific distribution(s) to fit, returns % ref to sum log likelihood and parameter values under fitted distribution(s). See FUNCTIONS above for available distributions. - -=for options - -Default options (case insensitive): - - MAXBN => 20, - # see PDL::Graphics::PGPLOT::Window for next options - WIN => undef, # pgwin object. not closed here if passed - # allows comparing multiple distr in same plot - # set env before passing WIN - DEV => '/xs' , # open and close dev for plotting if no WIN - # defaults to '/png' in Windows - COLOR => 1, # color for data distr - -=for usage - -Usage: - - # yes it threads :) - my $data = grandom( 500, 3 )->abs; - # ll on plot is sum across 3 data curves - my ($ll, $pars) - = $data->plot_distr( 'gaussian', 'lognormal', {DEV=>'/png'} ); - - # pars are from normalized data (ie data / bin_size) - print "$_\t@{$pars->{$_}}\n" for (sort keys %$pars); - print "$_\t$ll->{$_}\n" for (sort keys %$ll); - -=cut - -*plot_distr = \&PDL::plot_distr; -sub PDL::plot_distr { - require PDL::Graphics::PGPLOT::Window; - my ($self, @distr) = @_; - - my %opt = ( - MAXBN => 20, - WIN => undef, # pgwin object. not closed here if passed - DEV => $DEV, # open and close default win if no WIN - COLOR => 1, # color for data distr - ); - my $opt = pop @distr - if ref $distr[-1] eq 'HASH'; - $opt and $opt{uc $_} = $opt->{$_} for (keys %$opt); - - $self = $self->squeeze; - - # use int range, step etc for int xvals--pmf compatible - my $INT = 1 - if grep { /(?:binomial)|(?:geo)|(?:nbd)|(?:poisson)/ } @distr; - - my ($range, $step, $step_int); - $range = $self->max->sclr - $self->min->sclr; - $step = $range / $opt{MAXBN}; - $step_int = ($range <= $opt{MAXBN})? 1 - : PDL::ceil( $range / $opt{MAXBN} ) - ; - $opt{MAXBN} = PDL::ceil( $range / $step )->min->sclr; - - my $hist = $self->double->histogram($step, $self->min->sclr, $opt{MAXBN}); - # turn fre into prob - $hist /= $self->dim(0); - - my $xvals = $self->min->sclr + sequence( $opt{MAXBN} ) * $step; - my $xvals_int - = PDL::ceil($self->min->sclr) + sequence( $opt{MAXBN} ) * $step_int; - $xvals_int = $xvals_int->where( $xvals_int <= $xvals->max )->sever; - - my $win = $opt{WIN}; - if (!$win) { - $win = PDL::Graphics::PGPLOT::Window::pgwin( Dev=>$opt{DEV} ); - $win->env($xvals->minmax,0,1, {XTitle=>'xvals', YTitle=>'probability'}); - } - - $win->line( $xvals, $hist, { COLOR=>$opt{COLOR} } ); - - if (!@distr) { - $win->close - unless defined $opt{WIN}; - return; - } - - my (%ll, %pars, @text, $c); - $c = $opt{COLOR}; # fitted lines start from ++$c - for my $distr ( @distr ) { - # find mle_ or mme_$distr; - my @funcs = grep { /_$distr$/ } (keys %PDL::Stats::Distr::); - if (!@funcs) { - carp "Do not recognize $distr distribution!"; - next; - } - # might have mle and mme for a distr. sort so mle comes first - @funcs = sort @funcs; - my ($f_para, $f_prob) = @funcs[0, -1]; - - my $nrmd = $self / $step; - eval { - my @paras = $nrmd->$f_para(); - $pars{$distr} = \@paras; - - @paras = map { $_->dummy(0) } @paras; - $ll{$distr} = $nrmd->$f_prob( @paras )->log->sumover; - push @text, sprintf "$distr LL = %.2f", $ll{$distr}->sum; - - if ($f_prob =~ /^pdf/) { - $win->line( $xvals, ($xvals/$step)->$f_prob(@paras), {COLOR=>++$c} ); - } - else { - $win->points( $xvals_int, ($xvals_int/$step_int)->$f_prob(@paras), {COLOR=>++$c} ); - } - }; - carp $@ if $@; - } - $win->legend(\@text, ($xvals->min->sclr + $xvals->max->sclr)/2, .95, - {COLOR=>[$opt{COLOR}+1 .. $c], TextFraction=>.75} ); - $win->close - unless defined $opt{WIN}; - return (\%ll, \%pars); -} - -=head1 DEPENDENCIES - -GSL - GNU Scientific Library - -=head1 SEE ALSO - -PDL::Graphics::PGPLOT - -PDL::GSL::CDF - -=head1 AUTHOR - -Copyright (C) 2009 Maggie J. Xiong , David Mertens - -All rights reserved. There is no warranty. You are allowed to redistribute this software / documentation as described in the file COPYING in the PDL distribution. - -=cut - -EOD - -pp_done(); diff --git a/Libtmp/GSL/t/cdf.t b/Libtmp/GSL/t/cdf.t deleted file mode 100644 index 60e10291b..000000000 --- a/Libtmp/GSL/t/cdf.t +++ /dev/null @@ -1,108 +0,0 @@ -use strict; -use warnings; -use Test::More; -use PDL::LiteF; -use PDL::GSL::CDF; -use Test::PDL; - -{ - my $a = sequence 5; - my $a_bad = $a->append(pdl 'BAD'); - is_pdl $a->gsl_cdf_tdist_P(1999), my $texp = pdl('0.5 0.841284230695725 0.977182329566921 0.998633419750658 0.999967176033837'); - is_pdl $a_bad->gsl_cdf_tdist_P(1999), $texp->append(pdl 'BAD'); - is_pdl $a->gsl_cdf_fdist_P(1,1999), my $fexp = pdl('0 0.682568461391451 0.842545068166059 0.916581225902023 0.95436465913384'); - is_pdl $a_bad->gsl_cdf_fdist_P(1,1999), $fexp->append(pdl 'BAD'); -} -{ - my $x = pdl('[0.185 0.242 0.385 0.528 0.671 0.814]'); - is_pdl gsl_cdf_gaussian_P ($x, 0.6) , pdl([ 0.621085647787498, 0.656648488459185, 0.739455181655133, 0.810570345223288, 0.868287672230057, 0.912556425780295 ]); - is_pdl gsl_cdf_gaussian_Q ($x, 0.6) , pdl([ 0.378914352212502, 0.343351511540815, 0.260544818344867, 0.189429654776712, 0.131712327769943, 0.0874435742197054 ]); - is_pdl gsl_cdf_gaussian_Pinv ($x, 0.6) , pdl([ -0.53788401840115, -0.419930160118405, -0.175424937736083, 0.04214598829315, 0.265605686522293, 0.535639994592513 ]); - is_pdl gsl_cdf_gaussian_Qinv ($x, 0.6) , pdl([ 0.53788401840115, 0.419930160118405, 0.175424937736083, -0.04214598829315, -0.265605686522293, -0.535639994592513 ]); - is_pdl gsl_cdf_ugaussian_P ($x) , pdl([ 0.573385482213374, 0.5956099183852, 0.649881292234724, 0.701250332027385, 0.748889735236971, 0.792177524517595 ]); - is_pdl gsl_cdf_ugaussian_Q ($x) , pdl([ 0.426614517786626, 0.4043900816148, 0.350118707765276, 0.298749667972615, 0.251110264763029, 0.207822475482405 ]); - is_pdl gsl_cdf_ugaussian_Pinv ($x) , pdl([ -0.896473364001916, -0.699883600197341, -0.292374896226804, 0.0702433138219167, 0.442676144203822, 0.892733324320856 ]); - is_pdl gsl_cdf_ugaussian_Qinv ($x) , pdl([ 0.896473364001916, 0.699883600197341, 0.292374896226804, -0.0702433138219167, -0.442676144203822, -0.892733324320856 ]); - is_pdl gsl_cdf_exponential_Pinv ($x, 0.25) , pdl([ 0.0511417914353186, 0.0692679733349413, 0.121533252793905, 0.187694073349145, 0.277924382054191, 0.420502151317234 ]); - is_pdl gsl_cdf_exponential_Q ($x, 0.25) , pdl([ 0.477113915521034, 0.379841962851379, 0.214381101426978, 0.120995732814878, 0.0682894493122845, 0.0385422591266925 ]); - is_pdl gsl_cdf_exponential_P ($x, 0.25) , pdl([ 0.522886084478966, 0.620158037148621, 0.785618898573022, 0.879004267185122, 0.931710550687715, 0.961457740873307 ]); - is_pdl gsl_cdf_exponential_Qinv ($x, 0.25) , pdl([ 0.421849863475953, 0.354704388206363, 0.238627986173588, 0.159664748818969, 0.0997465355026138, 0.0514487282448992 ]); - is_pdl gsl_cdf_laplace_P ($x, 100) , pdl([ 0.500924144902391, 0.501208537080327, 0.501921299125978, 0.502633042650321, 0.5033437691088, 0.504053479954779 ]); - is_pdl gsl_cdf_laplace_Q ($x, 100) , pdl([ 0.499075855097609, 0.498791462919674, 0.498078700874022, 0.497366957349679, 0.4966562308912, 0.495946520045221 ]); - is_pdl gsl_cdf_laplace_Pinv ($x, 100) , pdl([ -99.4252273343867, -72.5670372265505, -26.1364764134408, 5.76291128366364, 41.855034765682, 98.886142470899 ]); - is_pdl gsl_cdf_laplace_Qinv ($x, 100) , pdl([ 99.4252273343867, 72.5670372265505, 26.1364764134408, -5.76291128366364, -41.855034765682, -98.886142470899 ]); - is_pdl gsl_cdf_exppow_P ($x, 10, 20) , pdl([ 0.509501756003763, 0.512429324069787, 0.51977392465648, 0.527118525243172, 0.534463125829865, 0.541807726416557 ]); - is_pdl gsl_cdf_exppow_Q ($x, 10, 20) , pdl([ 0.490498243996237, 0.487570675930213, 0.48022607534352, 0.472881474756828, 0.465536874170135, 0.458192273583443 ]); - is_pdl gsl_cdf_cauchy_P ($x, 10) , pdl([ 0.505888061226044, 0.507701596026139, 0.512248881040388, 0.516791169875417, 0.521326624632281, 0.525853424029018 ]); - is_pdl gsl_cdf_cauchy_Q ($x, 10) , pdl([ 0.494111938773956, 0.492298403973861, 0.487751118959612, 0.483208830124583, 0.478673375367719, 0.474146575970982 ]); - is_pdl gsl_cdf_cauchy_Pinv ($x, 10) , pdl([ -15.2235450689613, -10.5157249987904, -3.77868511775821, 0.881921819968596, 5.95646629389323, 15.1198164986352 ]); - is_pdl gsl_cdf_cauchy_Qinv ($x, 10) , pdl([ 15.2235450689613, 10.5157249987904, 3.77868511775821, -0.881921819968596, -5.95646629389323, -15.1198164986352 ]); - is_pdl gsl_cdf_rayleigh_P ($x, 0.6) , pdl([ 0.0464226377672899, 0.078118776801827, 0.1860595358795, 0.321044709711432, 0.46491770539411, 0.601589429889673 ]); - is_pdl gsl_cdf_rayleigh_Q ($x, 0.6) , pdl([ 0.95357736223271, 0.921881223198173, 0.8139404641205, 0.678955290288568, 0.53508229460589, 0.398410570110327 ]); - is_pdl gsl_cdf_rayleigh_Pinv ($x, 0.6) , pdl([ 0.383781655806681, 0.44664500803729, 0.591621304591413, 0.735227129019012, 0.894663188197699, 1.10047544079531 ]); - is_pdl gsl_cdf_rayleigh_Qinv ($x, 0.6) , pdl([ 1.10223754554576, 1.01071689311811, 0.829004583931798, 0.678110961862902, 0.535975766474127, 0.384931600866063 ]); - is_pdl gsl_cdf_gamma_P ($x, 1, 2) , pdl([ 0.0883507890355383, 0.113966040407124, 0.175105681796397, 0.232026460343294, 0.285019493288257, 0.334355809698479 ]); - is_pdl gsl_cdf_gamma_Q ($x, 1, 2) , pdl([ 0.911649210964462, 0.886033959592876, 0.824894318203603, 0.767973539656706, 0.714980506711743, 0.665644190301521 ]); - is_pdl gsl_cdf_gamma_Pinv ($x, 1, 2) , pdl([ 0.409134331482549, 0.554143786679531, 0.972266022351238, 1.50155258679316, 2.22339505643353, 3.36401721053787 ]); - is_pdl gsl_cdf_gamma_Qinv ($x, 1, 2) , pdl([ 3.37479890780762, 2.8376351056509, 1.90902388938871, 1.27731799055175, 0.79797228402091, 0.411589825959194 ]); - is_pdl gsl_cdf_flat_P ($x, 0, 2) , pdl([0.0925, 0.121, 0.1925, 0.264, 0.3355, 0.407]); - is_pdl gsl_cdf_flat_Q ($x, 0, 2) , pdl([0.9075, 0.879, 0.8075, 0.736, 0.6645, 0.593]); - is_pdl gsl_cdf_flat_Pinv ($x, 1, 2) , pdl([1.185, 1.242, 1.385, 1.528, 1.671, 1.814]); - is_pdl gsl_cdf_flat_Qinv ($x, 1, 2) , pdl([1.815, 1.758, 1.615, 1.472, 1.329, 1.186]); - is_pdl gsl_cdf_lognormal_P ($x, 0.3, 0.6), pdl([ 0.000462607682986099, 0.00208704926075921, 0.0182706108740019, 0.0588581153586641, 0.122014171048952, 0.199616592770584 ]); - is_pdl gsl_cdf_lognormal_Q ($x, 0.3, 0.6), pdl([ 0.999537392317014, 0.997912950739241, 0.981729389125998, 0.941141884641336, 0.877985828951048, 0.800383407229416 ]); - is_pdl gsl_cdf_lognormal_Pinv ($x, 0.3, 0.6), pdl([ 0.788294113403652, 0.886982381298518, 1.13266703782153, 1.40796582903933, 1.76051377932779, 2.30628958893332 ]); - is_pdl gsl_cdf_lognormal_Qinv ($x, 0.3, 0.6), pdl([ 2.31147076885183, 2.05428973428196, 1.60869764860025, 1.29414987410153, 1.03499263782317, 0.790065050431614 ]); - is_pdl gsl_cdf_chisq_P ($x, 0.9) , pdl([ 0.376001891123672, 0.420680573652885, 0.507543979704043, 0.573017615826538, 0.625398070381443, 0.668713631235268 ]); - is_pdl gsl_cdf_chisq_Q ($x, 0.9) , pdl([ 0.623998108876328, 0.579319426347115, 0.492456020295957, 0.426982384173462, 0.374601929618557, 0.331286368764732 ]); - is_pdl gsl_cdf_chisq_Pinv ($x, 0.9) , pdl([ 0.0363724832898845, 0.0667572490273275, 0.195682928939501, 0.425980141237634, 0.822286592027187, 1.57020451220134 ]); - is_pdl gsl_cdf_chisq_Qinv ($x, 0.9) , pdl([ 1.57778431676383, 1.21061535673313, 0.640242963605861, 0.320808004604559, 0.135209087406343, 0.0368164405929002 ]); - is_pdl gsl_cdf_fdist_P ($x, 0.5, 0.8), pdl([ 0.395899685715852, 0.421622763059381, 0.468763623791231, 0.502470355270446, 0.528704043285766, 0.550125367078533 ]); - is_pdl gsl_cdf_fdist_Q ($x, 0.5, 0.8), pdl([ 0.604100314284148, 0.578377236940619, 0.531236376208769, 0.497529644729554, 0.471295956714234, 0.449874632921467 ]); - is_pdl gsl_cdf_fdist_Pinv ($x, 1.5, 1.8), pdl([ 0.159905620801234, 0.241769364468725, 0.535256991978199, 1.04012738093102, 2.03748438382058, 4.73267894511171 ]); - is_pdl gsl_cdf_fdist_Qinv ($x, 1.5, 1.8), pdl([ 4.76724629388056, 3.2748615203822, 1.55210103408465, 0.806318097041121, 0.402458766933549, 0.161204550680556 ]); - is_pdl gsl_cdf_tdist_P ($x, 0.7) , pdl([ 0.553772454538094, 0.569686093296725, 0.607329840803688, 0.641083174744614, 0.670718124690489, 0.696452396237433 ]); - is_pdl gsl_cdf_tdist_Q ($x, 0.7) , pdl([ 0.446227545461906, 0.430313906703275, 0.392670159196312, 0.358916825255386, 0.329281875309511, 0.303547603762567 ]); - is_pdl gsl_cdf_tdist_Pinv ($x, 2) , pdl([ -1.14725744068977, -0.851907453506246, -0.33422960488652, 0.079320431592384, 0.514697310821784, 1.1412373360793 ]); - is_pdl gsl_cdf_tdist_Qinv ($x, 2) , pdl([ 1.14725744068977, 0.851907453506246, 0.33422960488652, -0.079320431592384, -0.514697310821784, -1.1412373360793 ]); - is_pdl gsl_cdf_beta_P ($x, 1, 2) , pdl([ 0.335775, 0.425436, 0.621775, 0.777216, 0.891759, 0.965404 ]); - is_pdl gsl_cdf_beta_Q ($x, 1, 2) , pdl([ 0.664225, 0.574564, 0.378225, 0.222784, 0.108241, 0.034596 ]); - is_pdl gsl_cdf_beta_Pinv ($x, 1, 2) , pdl([ 0.0972264957366104, 0.129368045612843, 0.215780642932094, 0.312977438507293, 0.426414784012001, 0.568722826943044 ]); - is_pdl gsl_cdf_beta_Qinv ($x, 1, 2) , pdl([ 0.569883736647869, 0.508065044950046, 0.379516317700457, 0.273363915016602, 0.180854103349104, 0.097780514508803 ]); - is_pdl gsl_cdf_logistic_P ($x, 2) , pdl([ 0.523108525489089, 0.53021314643553, 0.547976937636796, 0.565619324932982, 0.583097006692731, 0.600368317448984 ]); - is_pdl gsl_cdf_logistic_Q ($x, 2) , pdl([ 0.476891474510911, 0.46978685356447, 0.452023062363204, 0.434380675067018, 0.416902993307269, 0.399631682551016 ]); - is_pdl gsl_cdf_logistic_Pinv ($x, 2) , pdl([ -2.96566457632508, -2.28349131897137, -0.936757867037467, 0.224234596241412, 1.42542277241262, 2.95242738457868 ]); - is_pdl gsl_cdf_logistic_Qinv ($x, 2) , pdl([ 2.96566457632508, 2.28349131897137, 0.936757867037467, -0.224234596241412, -1.42542277241262, -2.95242738457868 ]); - is_pdl gsl_cdf_pareto_P ($x, 0.1, .05), pdl([ 0.122635965393589, 0.145886740885255, 0.184637135967562, 0.209988145413383, 0.228697483330511, 0.243455357847612 ]); - is_pdl gsl_cdf_pareto_Q ($x, 0.1, .05), pdl([ 0.877364034606411, 0.854113259114745, 0.815362864032438, 0.790011854586617, 0.771302516669489, 0.756544642152388 ]); - is_pdl gsl_cdf_pareto_Pinv ($x, 11, 22) , pdl([ 22.4129623633136, 22.5611817388298, 22.994070157874, 23.553980832887, 24.3396294860749, 25.6348390467721 ]); - is_pdl gsl_cdf_pareto_Qinv ($x, 11, 22) , pdl([ 25.6474051745322, 25.0287675983915, 23.9942991451309, 23.3151266418161, 22.8126206663786, 22.4154640895281 ]); - is_pdl gsl_cdf_weibull_P ($x, 1, 2) , pdl([ 0.0336459494865304, 0.0568821210761311, 0.137762910230166, 0.243296666516025, 0.362525498247215, 0.484488671492325 ]); - is_pdl gsl_cdf_weibull_Q ($x, 1, 2) , pdl([ 0.96635405051347, 0.943117878923869, 0.862237089769834, 0.756703333483975, 0.637474501752785, 0.515511328507675 ]); - is_pdl gsl_cdf_weibull_Pinv ($x, 1, 2) , pdl([ 0.452291018859842, 0.526376189943813, 0.6972323939517, 0.866473481069433, 1.05437067875428, 1.29692274452603 ]); - is_pdl gsl_cdf_weibull_Qinv ($x, 1, 2) , pdl([ 1.29899940488971, 1.19114128163936, 0.976991271554845, 0.799161432550317, 0.631653498375854, 0.453646242108977 ]); - is_pdl gsl_cdf_gumbel1_P ($x, 1, 2) , pdl([ 0.189719508676777, 0.208021808732089, 0.256429559887203, 0.307411906647793, 0.359732618472373, 0.412233898397574 ]); - is_pdl gsl_cdf_gumbel1_Q ($x, 1, 2) , pdl([ 0.810280491323223, 0.791978191267911, 0.743570440112797, 0.692588093352207, 0.640267381527627, 0.587766101602426 ]); - is_pdl gsl_cdf_gumbel1_Pinv ($x, 1, 2) , pdl([ 0.169958621442342, 0.343323365119587, 0.73970230236935, 1.14153180132375, 1.61197577505875, 2.27402235471976 ]); - is_pdl gsl_cdf_gumbel1_Qinv ($x, 1, 2) , pdl([ 2.28000609915082, 1.9766254442873, 1.4144201875401, 0.979794730452658, 0.58725902873148, 0.173158502928973 ]); - is_pdl gsl_cdf_gumbel2_P ($x, 1, 2) , pdl([ 2.01801560369643e-05, 0.000257507217296498, 0.00554529646233709, 0.0226435827865344, 0.0507610509809621, 0.0856914301179353 ]); - is_pdl gsl_cdf_gumbel2_Q ($x, 1, 2) , pdl([ 0.999979819843963, 0.999742492782704, 0.994454703537663, 0.977356417213466, 0.949238949019038, 0.914308569882065 ]); - is_pdl gsl_cdf_gumbel2_Pinv ($x, 1, 2) , pdl([ 1.18525580612995, 1.40962451163448, 2.09531165232346, 3.13156162332933, 5.0127054286201, 9.71841320586135 ]); - is_pdl gsl_cdf_gumbel2_Qinv ($x, 1, 2) , pdl([ 9.77674003915904, 7.21834313792145, 4.11410036764091, 2.6639093663331, 1.79905050540872, 1.18905455877868 ]); - my $y = long('[7 8 9 10 11 12 14]'); - is_pdl gsl_cdf_poisson_P ($y, 10) , pdl([ 0.220220646601699, 0.332819678750718, 0.457929714471852, 0.583039750192985, 0.696776146303107, 0.791556476394874, 0.916541527065337 ]); - is_pdl gsl_cdf_poisson_Q ($y, 10) , pdl([ 0.779779353398301, 0.667180321249282, 0.542070285528148, 0.416960249807015, 0.303223853696893, 0.208443523605126, 0.0834584729346629 ]); - is_pdl gsl_cdf_binomial_P ($y, 0.3, 33) , pdl([ 0.182152437649663, 0.304315672104686, 0.449748094074962, 0.599335728101516, 0.733381789761681, 0.838703695351808, 0.956260798765594 ]); - is_pdl gsl_cdf_binomial_Q ($y, 0.3, 33) , pdl([ 0.817847562350337, 0.695684327895314, 0.550251905925038, 0.400664271898484, 0.266618210238319, 0.161296304648192, 0.0437392012344063 ]); - is_pdl gsl_cdf_negative_binomial_P ($y, 0.3, 3.3), pdl([ 0.559244604008049, 0.633219417269081, 0.69823505870183, 0.754213525975423, 0.801591665095168, 0.841112929477554, 0.900208316084222 ]); - is_pdl gsl_cdf_negative_binomial_Q ($y, 0.3, 3.3), pdl([ 0.440755395991951, 0.366780582730919, 0.30176494129817, 0.245786474024577, 0.198408334904832, 0.158887070522446, 0.0997916839157782 ]); - is_pdl gsl_cdf_pascal_P ($y, 0.3, 3) , pdl([ 0.617217213599996, 0.687259545750001, 0.747184652145002, 0.7975217415168, 0.839164242724381, 0.873172285377237, 0.922614747387929 ]); - is_pdl gsl_cdf_pascal_Q ($y, 0.3, 3) , pdl([ 0.382782786400004, 0.312740454249998, 0.252815347854998, 0.2024782584832, 0.160835757275619, 0.126827714622763, 0.077385252612071 ]); - is_pdl gsl_cdf_geometric_P ($y, 0.3) , pdl([ 0.9176457, 0.94235199, 0.959646393, 0.9717524751, 0.98022673257, 0.986158712799, 0.99321776927151 ]); - is_pdl gsl_cdf_geometric_Q ($y, 0.3) , pdl([ 0.0823543, 0.05764801, 0.040353607, 0.0282475249, 0.01977326743, 0.013841287201, 0.00678223072849001 ]); - my $z = long('[6 7 8 9 10]'); - is_pdl gsl_cdf_hypergeometric_P ($z, 25, 5, 11), pdl([ 0.00324196875921016, 0.0472401162056337, 0.24523177971454, 0.61921603300914, 0.918403435644814 ]); - is_pdl gsl_cdf_hypergeometric_Q ($z, 25, 5, 11), pdl([ 0.99675803124079, 0.952759883794366, 0.75476822028546, 0.38078396699086, 0.0815965643551856 ]); -} - -done_testing; diff --git a/Libtmp/GSL/t/gsl_diff.t b/Libtmp/GSL/t/gsl_diff.t deleted file mode 100644 index 44053a92f..000000000 --- a/Libtmp/GSL/t/gsl_diff.t +++ /dev/null @@ -1,34 +0,0 @@ -# Test Script for the PDL interface to the GSL library -# This tests mainly that the interface is working, i.e. that the -# functions can be called. -# The GSL library already has a extensive test suite, and we -# do not want to duplicate that effort here. - -use strict; -use warnings; -use Test::More; -use PDL::LiteF; -use PDL::GSL::DIFF; - -my @res = gsldiff(\&myf,1.5); - -ok(abs($res[0]- 28.4632075095177) < 1e-6 ); - -@res = gsldiff(\&myf,1.5,{Method => 'central'}); - -ok(abs($res[0]- 28.4632075095177) < 1e-6 ); - -@res = gsldiff(\&myf,1.5,{Method => 'forward'}); - -ok(abs($res[0]- 28.4632852673531) < 1e-6 ); - -@res = gsldiff(\&myf,1.5,{Method => 'backward'}); - -ok(abs($res[0]-28.4631297516823 ) < 1e-6 ); - -done_testing; - -sub myf{ - my ($x) = @_; - return exp($x**2); -} diff --git a/Libtmp/GSL/t/gsl_integ.t b/Libtmp/GSL/t/gsl_integ.t deleted file mode 100644 index f4a00c430..000000000 --- a/Libtmp/GSL/t/gsl_integ.t +++ /dev/null @@ -1,132 +0,0 @@ -# Test Script for the PDL interface to the GSL library -# This tests mainly that the interface is working, i.e. that the -# functions can be called. -# The GSL library already has a extensive test suite, and we -# do not want to duplicate that effort here. - -use strict; -use warnings; -use Test::More; -use PDL::LiteF; -use PDL::GSL::INTEG; - -my $alfa = 2.6; - -my ($res,$abserr,$neval,$ierr) = gslinteg_qng(\&f1,0,1,0,1e-9); -ok(abs($res - 0.0771604938270651) < 1e-6); -($res,$abserr,$neval,$ierr) = gslinteg_qng(\&f1,0,1,0,1e-9,{Warn => 'y'}); -ok(abs($res - 0.0771604938270651) < 1e-6); - -($res,$abserr,$ierr) = gslinteg_qag(\&f1,0,1,0,1e-10,1000,1); -ok(abs($res - 0.0771604938271586) < 1e-6); -($res,$abserr,$ierr) = gslinteg_qag(\&f1,0,1,0,1e-10,1000,1,{Warn => 'y'}); -ok(abs($res - 0.0771604938271586) < 1e-6); - -($res,$abserr,$ierr) = gslinteg_qags(\&f1,0,1,0,1e-10,1000); -ok(abs($res - 0.0771604938271579) < 1e-6); -($res,$abserr,$ierr) = gslinteg_qags(\&f1,0,1,0,1e-10,1000,{Warn => 'y'}); -ok(abs($res - 0.0771604938271579) < 1e-6); - -my $points = pdl(0,1,sqrt(2),3); -($res,$abserr,$ierr) = gslinteg_qagp(\&f454,$points,0,1e-3,1000); -ok(abs($res - 52.7408061167272) < 1e-6); -($res,$abserr,$ierr) = gslinteg_qagp(\&f454,$points,0,1e-3,1000,{Warn => 'y'}); -ok(abs($res - 52.7408061167272) < 1e-6); - -($res,$abserr,$ierr) = gslinteg_qagi(\&myfn1,1e-7,0,1000); -ok(abs($res -2.27587579446875 ) < 1e-6); -($res,$abserr,$ierr) = gslinteg_qagi(\&myfn1,1e-7,0,1000,{Warn => 'y'}); -ok(abs($res -2.27587579446875 ) < 1e-6); - -$alfa = 1; -($res,$abserr,$ierr) = gslinteg_qagiu(\&f16,99.9,1e-7,0,1000); -ok(abs($res -0.000100000000000671) < 1e-6); -($res,$abserr,$ierr) = gslinteg_qagiu(\&f16,99.9,1e-7,0,1000,{Warn => 'y'}); -ok(abs($res -0.000100000000000671) < 1e-6); - -($res,$abserr,$ierr) = gslinteg_qagil(\&myfn2,1.0,1e-7,0,1000); -ok(abs($res -2.71828182845905) < 1e-6); -($res,$abserr,$ierr) = gslinteg_qagil(\&myfn2,1.0,1e-7,0,1000,{Warn => 'y'}); -ok(abs($res -2.71828182845905) < 1e-6); - - -($res,$abserr,$ierr) = gslinteg_qawc(\&f459,-1,5,0,0,1e-3,1000); -ok(abs($res + 0.08994400695837) < 1e-6); -($res,$abserr,$ierr) = gslinteg_qawc(\&f459,-1,5,0,0,1e-3,1000,{Warn => 'y'}); -ok(abs($res + 0.08994400695837) < 1e-6); - - -($res,$abserr,$ierr) = gslinteg_qaws(\&f458,0,0,1,0,0,1,0,1e-7,1000); -ok(abs($res + 0.18927518534894) < 1e-6); -($res,$abserr,$ierr) = gslinteg_qaws(\&f458,0,0,1,0,0,1,0,1e-7,1000,{Warn => 'y'}); -ok(abs($res + 0.18927518534894) < 1e-6); - - -my $PI = 3.14159265358979323846264338328; -($res,$abserr,$ierr) = gslinteg_qawo(\&f456,10*$PI,'sin',0,1,0,1e-7,1000); -ok(abs($res + 0.128136848399167) < 1e-6); -($res,$abserr,$ierr) = gslinteg_qawo(\&f456,10*$PI,'sin',0,1,0,1e-7,1000,{Warn => 'y'}); -ok(abs($res + 0.128136848399167) < 1e-6); - -($res,$abserr,$ierr) = gslinteg_qawf(\&f457,$PI/2.0,'cos',0,1e-7,1000); -ok(abs($res -0.999999999927978) < 1e-6); -($res,$abserr,$ierr) = gslinteg_qawf(\&f457,$PI/2.0,'cos',0,1e-7,1000,{Warn => 'y'}); -ok(abs($res -0.999999999927978) < 1e-6); - -done_testing; - -sub f1{ - my ($x) = @_; - return ($x**$alfa)*log(1.0/$x); -} - -sub f454{ - my ($x) = @_; - my $x2 = $x**2; - my $x3 = $x**3; - return $x3 * log(abs(($x2-1.0)*($x2-2.0))); -} - -sub myfn1{ - my ($x) = @_; - return exp(-$x - $x*$x) ; -} - -sub f16 { - my ($x) = @_; - if (($x==0) && ($alfa == 1)) {return 1;} - if (($x==0) && ($alfa > 1)) {return 0;} - return ($x**($alfa-1))/((1+10*$x)**2); -} - -sub myfn2{ - my ($x) = @_; - return exp($alfa*$x); -} - -sub f459{ - my ($x) = @_; - return 1.0 / (5.0 * $x * $x * $x + 6.0) ; -} - -sub f458{ - my ($x) = @_; - if($x==0){return 0;} - else{ - my $u = log($x); - my $v = 1 + $u*$u; - return 1.0/($v*$v); - } -} - -sub f456{ - my ($x) = @_; - if($x==0){return 0;} - else{ return log($x);} -} - -sub f457{ - my ($x) = @_; - if ($x == 0){return 0;} - return 1.0/sqrt($x) -} diff --git a/Libtmp/GSL/t/gsl_interp.t b/Libtmp/GSL/t/gsl_interp.t deleted file mode 100644 index fc44b2efa..000000000 --- a/Libtmp/GSL/t/gsl_interp.t +++ /dev/null @@ -1,41 +0,0 @@ - -# Test Script for the PDL interface to the GSL library -# This tests mainly that the interface is working, i.e. that the -# functions can be called. -# The GSL library already has a extensive test suite, and we -# do not want to duplicate that effort here. - -use strict; -use warnings; -use PDL::LiteF; -use PDL::GSL::INTERP; -use Test::More; - -my $x = sequence(10); -my $y = exp($x); -my $spl = PDL::GSL::INTERP->init('cspline',$x,$y); - -ok(defined $spl, 'create cspline'); - -my $spl2 = PDL::GSL::INTERP->init('cspline',$x,$y,{Sort => 0}); - -ok(defined $spl2, 'create cspline no sort'); - -ok(abs($spl->eval(5.5)-237.641810667697) < 1e-6, 'eval 5.5' ); -ok(abs($spl->deriv(5.5)-237.669424604497) < 1e-6, 'deriv 5.5' ); -ok(abs($spl->deriv2(5.5)-306.23332503967) < 1e-6, 'deriv2 5.5' ); -ok(abs($spl->integ(3.2,8.5)-4925.23555581654) < 1e-6, 'integ 3.2 to 8.5' ); -eval { $spl->eval(10) }; isnt($@,'', 'eval 10 (extrapolation not allowed)') ; -eval { $spl->deriv(10) }; isnt($@,'', 'deriv 10 (extrapolation not allowed)' ); -eval { $spl->deriv2(10) }; isnt($@,'', 'deriv2 10 (extrapolation not allowed)'); -eval { $spl->integ(8.5,10)}; isnt($@,'', 'integ 8.5 to 10 (extrapolation not allowed)'); - -# Bad value test added 5/31/2005 D. Hunt - -ok ($spl->eval(pdl(0)->setbadat(0))->isbad, 'cspline eval w badvalue'); - -my $nx = ($x)*($x<=3) + ($x-1)*($x>3); # x value not monotonically increasing -my $i; eval { $i = PDL::GSL::INTERP->init('cspline',$nx, $y) }; -isnt($@,'',"module exception handling"); - -done_testing; diff --git a/Libtmp/GSL/t/gsl_linalg.t b/Libtmp/GSL/t/gsl_linalg.t deleted file mode 100644 index 77571d336..000000000 --- a/Libtmp/GSL/t/gsl_linalg.t +++ /dev/null @@ -1,40 +0,0 @@ -use strict; -use warnings; -use Test::More; -use Test::PDL; -use PDL::LiteF; -use PDL::GSL::LINALG; - -my $A = pdl [ - [0.18, 0.60, 0.57, 0.96], - [0.41, 0.24, 0.99, 0.58], - [0.14, 0.30, 0.97, 0.66], - [0.51, 0.13, 0.19, 0.85], -]; -my $B = sequence(2,4); # column vectors, but must transpose for GSL - -LU_decomp(my $lu=$A->copy, my $p=null, my $signum=null); -LU_solve($lu, $p, $B->transpose, my $x=null); -$x = $x->inplace->transpose; -is_pdl $A x $x, $B; -is_pdl LU_det($lu, $signum), $lu->diagonal(0,1)->prodover * $signum; - -my $D = 3; -#solve (1 + i)(b_{n + 1} - b_n)=1 - i with homogeneous BCs -my $c=zeroes($D); # subdiag -my $d=-ones($D); # diagonal -my $e=ones($D); $e->slice('(-1)').=0; # superdiag -my $b=ones($D); $b->slice('(-1)').=(1-$D); -($b, my $info) = solve_tridiag($d, $e, $c, $b, $x=null); -is_pdl $x, sequence($D), "tridiag"; - -# now complex -$A = czip($A, 1e-9); -$B = czip($B, 2); -LU_decomp($lu=$A->copy, $p=null, $signum=null); -LU_solve($lu, $p, $B->transpose, $x=null); -$x = $x->inplace->transpose; -is_pdl $A x $x, $B; -is_pdl LU_det($lu, $signum), $lu->diagonal(0,1)->prodover * $signum; - -done_testing; diff --git a/Libtmp/GSL/t/gsl_mroot.t b/Libtmp/GSL/t/gsl_mroot.t deleted file mode 100644 index 29546743a..000000000 --- a/Libtmp/GSL/t/gsl_mroot.t +++ /dev/null @@ -1,33 +0,0 @@ -# Test Script for the PDL interface to the GSL library -# This tests mainly that the interface is working, i.e. that the -# functions can be called. -# The GSL library already has a extensive test suite, and we -# do not want to duplicate that effort here. - -use strict; -use warnings; -use PDL::LiteF; -use PDL::GSL::MROOT; -use Test::More; - -my $init = pdl (-10.00, -5.0); -my $epsabs = 1e-7; - -my $res = gslmroot_fsolver($init, \&rosenbrock,{Method => 0, EpsAbs => $epsabs}); - -my @res = list ($res); - -ok(abs($res[0]- 1) < 1e-6 ); -ok(abs($res[1]- 1) < 1e-6 ); - -done_testing; - -sub rosenbrock{ - my ($x) = @_; - my $a1 = 1; - my $b1 = 10; - my $y = zeroes($x); - $y->slice(0) .= $a1 * (1 - $x->slice(0)); - $y->slice(1) .= $b1 * ($x->slice(1) - $x->slice(0)**2); - return $y; -} diff --git a/Libtmp/GSL/t/gsl_rng.t b/Libtmp/GSL/t/gsl_rng.t deleted file mode 100644 index 4b7cd6b42..000000000 --- a/Libtmp/GSL/t/gsl_rng.t +++ /dev/null @@ -1,138 +0,0 @@ -# Test Script for the PDL interface to the GSL library -# This tests only that the interface is working, i.e. that the -# functions can be called. The actual return values are not -# checked. -# The GSL library already has a extensive test suite, and we -# do not want to duplicate that effort here. - -use strict; -use warnings; -use Test::More; -use PDL::LiteF; -use PDL::GSL::RNG; - -my $image = zeroes(10,10); -my $ndim = 2; -my $name = ''; -my $sigma = 1; - -# new() function Test: -my $rng = PDL::GSL::RNG->new('taus'); - -pass('new() function'); - -# set_seed(); function Test: -$rng->set_seed(666); - -pass('set_seed(); function'); - -my $rng2 = PDL::GSL::RNG->new('taus')->set_seed(666); -is(ref $rng2, 'PDL::GSL::RNG', 'PDL::GSL::RNG->new(..)->set_seed(..)'); - -# min() function Test: -my $min = $rng->min(); my $max = $rng->max(); - -pass('min() function'); - -# rmax() function Test: -$min = $rng->min(); $max = $rng->max(); - -pass('rmax() function'); - -# name() function Test: -$name = $rng->name(); - -pass('name() function'); - -# get_uniform() function Test: -my $x = zeroes 5,6; $max=100; -my $o = $rng->get_uniform(10,10); $rng->get_uniform($x); - -pass('get_uniform() function'); - -# get_uniform_pos() function Test: -$x = zeroes 5,6; - -$o = $rng->get_uniform_pos(10,10); $rng->get_uniform_pos($x); - -pass('get_uniform_pos() function'); - -# get() function Test: -$x = zeroes 5,6; - -$o = $rng->get(10,10); $rng->get($x); - -pass('get() function'); - -# get_int() function Test: -$x = zeroes 5,6; $max=100; - -$o = $rng->get(10,10); $rng->get($x); - -pass('get_int() function'); - -# ran_gaussian() function Test: -$o = $rng->ran_gaussian($sigma,10,10); - -$rng->ran_gaussian($sigma,$x); - - -pass('ran_gaussian() function'); - -# $rng->ran_gaussian_var() function Test: -my $sigma_pdl = rvals zeroes 11,11; $o = $rng->ran_gaussian_var($sigma_pdl); - -pass('ran_gaussian_var() method'); - -# ran_additive_gaussian() function Test: -$rng->ran_additive_gaussian(1,$image); - -pass('ran_additive_gaussian() method'); - -# ran_additive_poisson() function Test: -$rng->ran_additive_poisson(1,$image); - -pass('ran_additive_poisson() method'); - -# ran_feed_poisson() function Test: -$rng->ran_feed_poisson($image); - -pass('ran_feed_poisson() method'); - -# ran_bivariate_gaussian() function Test: -$o = $rng->ran_bivariate_gaussian(1,2,0.5,1000); - -pass('ran_bivariate_gaussian() method'); - -# ran_dir() function Test: -$o = $rng->ran_dir($ndim,12); - -pass('ran_dir() method'); - -# ran_discrete_preproc() function Test: -my $prob = pdl [0.1,0.3,0.6]; - -my $discrete_dist_handle = $rng->ran_discrete_preproc($prob); - -$o = $rng->ran_discrete($discrete_dist_handle,100); - -pass('ran_discrete_preproc() method'); - -# ran_discrete() function Test: -$prob = pdl [0.1,0.3,0.6]; - -$discrete_dist_handle = $rng->ran_discrete_preproc($prob); - -$o = $rng->ran_discrete($discrete_dist_handle,100); - -pass('ran_discrete() method'); - -my $vec2d = sequence(10,10); -$rng->ran_shuffle($_) for $vec2d->dog; -ok any($vec2d != sequence(10,10)), 'ran_shuffle() method'; - -$vec2d = sequence(10,10); -$rng->ran_shuffle_1d($vec2d); -ok any($vec2d != sequence(10,10)), 'ran_shuffle_1d() method'; - -done_testing; diff --git a/Libtmp/GSL/t/gsl_sf.t b/Libtmp/GSL/t/gsl_sf.t deleted file mode 100644 index 49e8e6ff6..000000000 --- a/Libtmp/GSL/t/gsl_sf.t +++ /dev/null @@ -1,39 +0,0 @@ -# Test Script for the PDL interface to the GSL library -# This tests only that the interface is working, i.e. that the -# functions can be called. The actual return values are not -# checked. -# The GSL library already has a extensive test suite, and we -# do not want to duplicate that effort here. - -use strict; -use warnings; -use Test::More; -use Test::PDL; -use PDL::LiteF; -use PDL::GSL::SF; - -my ($y,$err) = gsl_sf_bessel_Jn(5.0, 0); -is_pdl $y, pdl(-0.17759677131433830434739701), "GSL SF Bessel function"; - -chomp(my $version = `gsl-config --version`); -if ((split /\./, $version)[0] >= 2) { - my $Ylm = gsl_sf_legendre_array(xvals(21)/10-1,'Y',4,-1); - ok($Ylm->slice("(0)")->uniq->nelem == 1, "Legendre Y00 is constant"); - is_pdl $Ylm->slice("(0),(0)"), pdl(0.5/sqrt(3.141592654)), "Y00 value is correct"; -} - -{ -# Check that the PDL error handler gets called instead of aborting -# Failure is an abort. -my @warning; -local $SIG{__WARN__} = sub { push @warning, @_ }; -my $err_test = eval {gsl_sf_lngamma(pdl(0))}; -isnt $@, '', "Got an error for invalid input"; -ok !@warning, 'no warnings' or diag explain \@warning; -} - -($y, my $e) = gsl_sf_airy_Ai(sequence(4)); -is_pdl $y, pdl([0.35502805, 0.13529242, 0.03492413, 0.0065911394]); -is_pdl $e, pdl([8.3366727e-17, 5.6151774e-17, 3.9261626e-17, 1.0852712e-17]); - -done_testing; diff --git a/Libtmp/GSL/t/stats_distr.t b/Libtmp/GSL/t/stats_distr.t deleted file mode 100644 index 28f3f2829..000000000 --- a/Libtmp/GSL/t/stats_distr.t +++ /dev/null @@ -1,138 +0,0 @@ -use strict; -use warnings; -use Test::More; -use PDL::Stats::Distr; -use PDL::LiteF; -use Test::PDL; - -{ - my $a = sequence(5) / 10; - is_pdl pdl($a->mme_beta), pdl(1.4, 5.6); - is_pdl $a->pdf_beta(1, 3), pdl '3 2.43 1.92 1.47 1.08'; -} -{ - my $a = sequence 5; - $a %= 2; - is_pdl pdl($a->mme_binomial), pdl(1, .4); - is_pdl $a->pmf_binomial(2,.4), pdl '0.36 0.48 0.36 0.48 0.36'; -} -{ - my $a = sequence 5; - is_pdl pdl($a->mle_exp), pdl .5; - is_pdl $a->pdf_exp(2.5), pdl '2.5 0.205212 0.016844 0.001382 0.000113'; -} -{ - my $a = sequence 5; - is_pdl pdl($a->mle_gaussian), pdl(2,2); - is_pdl $a->pdf_gaussian(1,2), pdl '0.219695644733861 0.282094791773878 0.219695644733861 0.103776874355149 0.0297325723059073'; -} -{ - my $a = sequence 5; - is_pdl pdl($a->mle_geo), pdl 0.333333333333333; - is_pdl $a->pmf_geo(.5), pdl '0.5 0.25 0.125 0.0625 0.03125'; -} -{ - my $a = sequence 5; - $a += 1; - is_pdl pdl($a->mle_geosh), pdl 0.333333333333333; - is_pdl $a->pmf_geosh(.5), pdl '0.5 0.25 0.125 0.0625 0.03125'; -} -{ - my $a = sequence(5) + 1; - is_pdl pdl($a->mle_lognormal), pdl(0.957498348556409, 0.323097797388514); - is_pdl pdl($a->mme_lognormal), pdl(2.19722457733622, 0.200670695462151); - is_pdl $a->pdf_lognormal(1,2), pdl '0.219695644733861 0.137765961149997 0.0938032750793072 0.0679412228146167 0.0514161127408299'; -} -{ - my $a = sequence 5; - $a *= $a; - is_pdl pdl($a->mme_nbd), pdl(1.25, 0.172413793103448); - is_pdl $a->pmf_nbd(2, .4), pdl '0.16 0.192 0.10368 0.0161243136 0.000767341894828032'; -} -{ - my $a = sequence 5; - $a += 1; - is_pdl pdl($a->mme_pareto), pdl(1.4, 0.857142857142857); - is_pdl $a->pdf_pareto(2, .4), pdl '0.32 0.04 0.0118518518518519 0.005 0.00256'; -} -{ - my $a = sequence 5; - $a %= 2; - is_pdl pdl($a->mle_poisson), pdl .4; - is_pdl $a->pmf_poisson(.4), pdl '0.670320046035639 0.268128018414256 0.670320046035639 0.268128018414256 0.670320046035639'; -} -{ - my $a = sequence 6; - $a->setbadat(-1); - $a /= 10; - is_pdl pdl($a->mme_beta), pdl(1.4, 5.6); - is_pdl $a->pdf_beta(1, 3), pdl '3 2.43 1.92 1.47 1.08 BAD'; -} -{ - my $a = sequence 6; - $a->setbadat(-1); - $a %= 2; - is_pdl pdl($a->mme_binomial), pdl(1, .4); - is_pdl $a->pmf_binomial(2,.4), pdl '0.36 0.48 0.36 0.48 0.36 BAD'; -} -{ - my $a = sequence 6; - $a->setbadat(-1); - is_pdl pdl($a->mle_exp), pdl .5; - is_pdl $a->pdf_exp(2.5), pdl '2.5 0.205212496559747 0.0168448674977137 0.00138271092536958 0.000113499824406212 BAD'; -} -{ - my $a = sequence 6; - $a->setbadat(-1); - is_pdl pdl($a->mle_gaussian), pdl(2,2); - is_pdl $a->pdf_gaussian(1,2), pdl '0.219695644733861 0.282094791773878 0.219695644733861 0.103776874355149 0.0297325723059073 BAD'; -} -{ - my $a = sequence 6; - $a->setbadat(-1); - is_pdl pdl($a->mle_geo), pdl 0.333333333333333; - is_pdl $a->pmf_geo(.5), pdl '0.5 0.25 0.125 0.0625 0.03125 BAD'; -} -{ - my $a = sequence 6; - $a->setbadat(-1); - $a += 1; - is_pdl pdl($a->mle_geosh), pdl 0.333333333333333; - is_pdl $a->pmf_geosh(.5), pdl '0.5 0.25 0.125 0.0625 0.03125 BAD'; -} -{ - my $a = sequence 6; - $a->setbadat(-1); - $a += 1; - is_pdl pdl($a->mle_lognormal), pdl(0.957498348556409, 0.323097797388514); - is_pdl pdl($a->mme_lognormal), pdl(2.19722457733622, 0.200670695462151); - is_pdl $a->pdf_lognormal(1,2), pdl '0.219695644733861 0.137765961149997 0.0938032750793072 0.0679412228146167 0.0514161127408299 BAD'; -} -{ - my $a = sequence 6; - $a->setbadat(-1); - $a *= $a; - is_pdl pdl($a->mme_nbd), pdl(1.25, 0.172413793103448); - is_pdl $a->pmf_nbd(2, .4), pdl '0.16 0.192 0.10368 0.0161243136 0.000767341894828032 BAD'; -} -{ - my $a = sequence 6; - $a->setbadat(-1); - $a += 1; - is_pdl pdl($a->mme_pareto), pdl(1.4, 0.857142857142857); - is_pdl $a->pdf_pareto(2, .4), pdl '0.32 0.04 0.0118518518518519 0.005 0.00256 BAD'; -} -{ - my $a = sequence 6; - $a->setbadat(-1); - $a %= 2; - is_pdl pdl($a->mle_poisson), pdl .4; - is_pdl $a->pmf_poisson(.4), pdl '0.670320046035639 0.268128018414256 0.670320046035639 0.268128018414256 0.670320046035639 BAD'; - is_pdl $a->pmf_poisson_factorial(.4), pdl '0.670320046035639 0.268128018414256 0.670320046035639 0.268128018414256 0.670320046035639 BAD'; - is_pdl $a->pmf_poisson_stirling(.4), pdl '0.670320046035639 0.268050878476493 0.670320046035639 0.268050878476493 0.670320046035639 BAD'; - - $a += 171; - is_pdl $a->pmf_poisson_stirling(10), pdl '0 0 0 0 0 BAD'; -} - -done_testing; diff --git a/Libtmp/GSL/typemap b/Libtmp/GSL/typemap deleted file mode 100644 index 93ee3929a..000000000 --- a/Libtmp/GSL/typemap +++ /dev/null @@ -1,15 +0,0 @@ -TYPEMAP -GslSpline * T_PTROBJ -GslAccel * T_PTROBJ -gsl_spline * T_PTR -gsl_interp_accel * T_PTR -gsl_rng * ANY_OBJ -gsl_ran_discrete_t * T_PTROBJ - -OUTPUT -ANY_OBJ - sv_setref_pv($arg, CLASS, (void *) $var); - -INPUT -ANY_OBJ - $var = INT2PTR($type, SvIV((SV*)SvRV($arg))) diff --git a/MANIFEST b/MANIFEST index 57d4901bb..4500af935 100644 --- a/MANIFEST +++ b/MANIFEST @@ -274,33 +274,6 @@ Libtmp/Fit/Polynomial.pm Libtmp/Fit/t/poly.t Libtmp/Func.pm Libtmp/Func_demo.pm -Libtmp/GSL/gslerr.h -Libtmp/GSL/lib/PDL/Demos/GSL_CDF.pm -Libtmp/GSL/lib/PDL/Demos/GSL_RNG.pm -Libtmp/GSL/lib/PDL/GSL/CDF.pd -Libtmp/GSL/lib/PDL/GSL/DIFF-FUNC.c -Libtmp/GSL/lib/PDL/GSL/DIFF.pd -Libtmp/GSL/lib/PDL/GSL/INTEG-FUNC.c -Libtmp/GSL/lib/PDL/GSL/INTEG.pd -Libtmp/GSL/lib/PDL/GSL/INTERP.pd -Libtmp/GSL/lib/PDL/GSL/LINALG.pd -Libtmp/GSL/lib/PDL/GSL/MROOT-FUNC.c -Libtmp/GSL/lib/PDL/GSL/MROOT.pd -Libtmp/GSL/lib/PDL/GSL/RNG.pd -Libtmp/GSL/lib/PDL/GSL/SF.pd -Libtmp/GSL/lib/PDL/Stats/Distr.pd -Libtmp/GSL/Makefile.PL -Libtmp/GSL/README -Libtmp/GSL/t/cdf.t -Libtmp/GSL/t/gsl_diff.t -Libtmp/GSL/t/gsl_integ.t -Libtmp/GSL/t/gsl_interp.t -Libtmp/GSL/t/gsl_linalg.t -Libtmp/GSL/t/gsl_mroot.t -Libtmp/GSL/t/gsl_rng.t -Libtmp/GSL/t/gsl_sf.t -Libtmp/GSL/t/stats_distr.t -Libtmp/GSL/typemap Libtmp/Makefile.PL Libtmp/Simplex/Demo.pm Libtmp/Simplex/Makefile.PL diff --git a/MANIFEST.SKIP b/MANIFEST.SKIP index 91bcdbce8..1ecda2ec7 100644 --- a/MANIFEST.SKIP +++ b/MANIFEST.SKIP @@ -57,15 +57,6 @@ Makefile\.old ^Basic/lib/PDL/Core/pdl\.h$ ^Basic/script/pdl$ ^Libtmp/Fit/Gaussian/Gaussian\..* -^Libtmp/GSL/lib/PDL/GSL/CDF(\.(pm|xs|c)$|-pp-) -^Libtmp/GSL/lib/PDL/GSL/DIFF(\.(pm|xs|c)$|-pp-) -^Libtmp/GSL/lib/PDL/GSL/INTEG(\.(pm|xs|c)$|-pp-) -^Libtmp/GSL/lib/PDL/GSL/INTERP(\.(pm|xs|c)$|-pp-) -^Libtmp/GSL/lib/PDL/GSL/LINALG(\.(pm|xs|c)$|-pp-) -^Libtmp/GSL/lib/PDL/GSL/MROOT(\.(pm|xs|c)$|-pp-) -^Libtmp/GSL/lib/PDL/GSL/RNG(\.(pm|xs|c)$|-pp-) -^Libtmp/GSL/lib/PDL/GSL/SF(\.(pm|xs|c)$|-pp-) -^Libtmp/GSL/lib/PDL/Stats/Distr(\.(pm|xs|c)$|-pp-) ^Libtmp/Slatec/Slatec..* ^Libtmp/Slatec/slatec/.*\.c$ \b[\._]Inline