File: erfc_check

package info (click to toggle)
libmath-gsl-perl 0.45-2
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid, trixie
  • size: 192,156 kB
  • sloc: ansic: 895,524; perl: 24,682; makefile: 12
file content (35 lines) | stat: -rwxr-xr-x 1,033 bytes parent folder | download | duplicates (5)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
#!/usr/bin/perl
# interfc
#   v0.2   Keith Lofstrom   2008 Sept 18
#         with tweaks by Jonathan Leto
#
# Compare the results of numerical integration of
# (2.0/sqrt(PI)*integral[x,inf](e^(-x^2) ) with the
# erfc() function.
#
# The results are pretty good, within 1e-14 or better
# for the examples tried

use Math::GSL::SF qw/:all/;
use Math::GSL::Integration qw/:all/;
use Math::GSL::Const qw/:all/;
use strict;

my $workspace  = gsl_integration_workspace_alloc(101);
my $maxdiff = 0;

printf " arg         erfc             int()             difference\n";

for( my $i=-5.2 ; $i<=5.201 ; $i += 0.1 ) {
   my $e  = gsl_sf_erfc($i);
   my ($status, $result, $abserr ) = gsl_integration_qagiu(
                     sub { exp( -($_[0]**2) ) },
                     $i, 0, 1e-9, 100, $workspace );

   my $yy = $M_2_SQRTPI * $result ;
   my $diff=$yy-$e ;
   $maxdiff = abs($diff) > $maxdiff ? abs($diff) : $maxdiff;

   printf "%4.1f%18.12f%18.12f %.18g\n", $i, $e, $yy, $diff ;
}
printf "Maximum local error: %.12g\n", $maxdiff;