tools/process-amath.pl

   1 #!/usr/bin/perl -w
   2
   3 # This script processes the test cases from amath, see
   4 # http://www.wolfgang-ehrhardt.de/amath_functions.html
   5
   6 use strict;
   7
   8 my $debug_underflow = 0;
   9 my $debug_overflow = 0;
  10 my $debug_arguments = 1;
  11 my $dir = $ARGV[0];
  12 die "$0: missing amath directory\n" unless (defined $dir) && -d $dir;
  13
  14 my @test_files =
  15     ('t_sfd1a.pas',
  16      't_sfd1.pas',
  17      't_sfd3a.pas',
  18      't_sfd3b.pas',
  19      't_sfd3c.pas',
  20      't_sfd3.pas',
  21      't_sfd4.pas',
  22      't_sfd6.pas',
  23      't_amath1.pas',
  24      't_amathm.pas',
  25     );
  26
  27 my %name_map =
  28     ('lnbeta' => 'betaln',
  29      'beta' => 'beta',
  30      'lngamma' => 'gammaln',
  31      'gamma' => 'gamma',
  32      'fac' => 'fact',           # no actual tests
  33      'dfac' => 'factdouble',
  34      'pochhammer' => 'pochhammer',
  35      'binomial' => 'combin',
  36      'cauchy_cdf' => 'r.pcauchy',
  37      'cauchy_inv' => 'r.qcauchy',
  38      'cauchy_pdf' => 'r.dcauchy',
  39      'chi2_cdf' => 'r.pchisq',
  40      'chi2_inv' => 'r.qchisq',
  41      'chi2_pdf' => 'r.dchisq',
  42      'exp_cdf' => 'r.pexp',
  43      'exp_inv' => 'r.qexp',
  44      'exp_pdf' => 'r.dexp',
  45      'gamma_cdf' => 'r.pgamma',
  46      'gamma_inv' => 'r.qgamma',
  47      'gamma_pdf' => 'r.dgamma',
  48      'laplace_pdf' => 'laplace',
  49      'logistic_pdf' => 'logistic',
  50      'lognormal_cdf' => 'r.plnorm',
  51      'lognormal_inv' => 'r.qlnorm',
  52      'lognormal_pdf' => 'r.dlnorm',
  53      'pareto_pdf' => 'pareto',
  54      'weibull_cdf' => 'r.pweibull',
  55      'weibull_inv' => 'r.qweibull',
  56      'weibull_pdf' => 'r.dweibull',
  57      'binomial_pmf' => 'r.dbinom',
  58      'binomial_cdf' => 'r.pbinom',
  59      'poisson_pmf' => 'r.dpois',
  60      'poisson_cdf' => 'r.ppois',
  61      'negbinom_pmf' => 'r.dnbinom',
  62      'negbinom_cdf' => 'r.pnbinom',
  63      'hypergeo_pmf' => 'r.dhyper',
  64      'hypergeo_cdf' => 'r.phyper',
  65      'rayleigh_pdf' => 'rayleigh',
  66      'normal_cdf' => 'r.pnorm',
  67      'normal_inv' => 'r.qnorm',
  68      'normal_pdf' => 'r.dnorm',
  69      'beta_cdf' => 'r.pbeta',
  70      'beta_inv' => 'r.qbeta',
  71      'beta_pdf' => 'r.dbeta',
  72      't_cdf' => 'r.pt',
  73      't_inv' => 'r.qt',
  74      't_pdf' => 'r.dt',
  75      'f_cdf' => 'r.pf',
  76      'f_inv' => 'r.qf',
  77      'f_pdf' => 'r.df',
  78      'erf' => 'erf',
  79      'erfc' => 'erfc',
  80      'bessel_j0' => 'besselj0', # Really named besselj
  81      'bessel_j1' => 'besselj1', # Really named besselj
  82      'bessel_jv' => 'besselj',
  83      'bessel_y0' => 'bessely0', # Really named bessely
  84      'bessel_y1' => 'bessely1', # Really named bessely
  85      'bessel_yv' => 'bessely',
  86      'bessel_i0' => 'besseli0', # Really named besseli
  87      'bessel_i1' => 'besseli1', # Really named besseli
  88      'bessel_iv' => 'besseli',
  89      'bessel_k0' => 'besselk0', # Really named besselk
  90      'bessel_k1' => 'besselk1', # Really named besselk
  91      'bessel_kv' => 'besselk',
  92      'exp' => 'exp',
  93      'exp2' => 'exp2',
  94      'exp10' => 'exp10',
  95      'expm1' => 'expm1',
  96      'ln' => 'ln',
  97      'ln1p' => 'ln1p',
  98      'log10' => 'log10',
  99      'log2' => 'log2',
 100      'arccos' => 'acos',
 101      'arccosh' => 'acosh',
 102      'arcsin' => 'asin',
 103      'arcsinh' => 'asinh',
 104      'arccot' => 'acot',
 105      'arccoth' => 'acoth',
 106      'arctan' => 'atan',
 107      'arctanh' => 'atanh',
 108      'cos' => 'cos',
 109      'cosh' => 'cosh',
 110      'cot' => 'cot',
 111      'coth' => 'coth',
 112      'csc' => 'csc',
 113      'csch' => 'csch',
 114      'sec' => 'sec',
 115      'sech' => 'sech',
 116      'sin' => 'sin',
 117      'sinh' => 'sinh',
 118      'tan' => 'tan',
 119      'tanh' => 'tanh',
 120      'gd' => 'gd',
 121     );
 122
 123 my %invalid_tests =
 124     (# Magically changed to something else
 125      'cos(1.0)' => 1,
 126      'cos(0.0)' => 1,
 127      'cos(1e26)' => 1, # 1e26 not representable
 128      'cot(1e26)' => 1, # 1e26 not representable
 129      'sin(1e26)' => 1, # 1e26 not representable
 130      'tan(1e26)' => 1, # 1e26 not representable
 131
 132      # Just plain wrong (and would depend on representation anyway)
 133      'besselj(11.791534439014281614,0)' => 1,
 134      'besselj(13.323691936314223032,1)' => 1,
 135      'bessely(13.36109747387276348,0)' => 1,
 136      'bessely(14.89744212833672538,1)' => 1,
 137
 138      # Overflow, not zero
 139      'bessely(1.5,-1700.5)' => 1,
 140     );
 141
 142 sub def_expr_handler {
 143     my ($f,$pa) = @_;
 144     my $expr = "$f(" . join (",", @$pa) . ")";
 145     return undef if exists $invalid_tests{$expr};
 146     return $expr;
 147 }
 148
 149 my %expr_handlers =
 150     ('beta' => \&non_negative_handler,
 151      'gammaln' => \&non_negative_handler,
 152      'factdouble' => \&non_negative_handler,
 153      'combin' => \&non_negative_handler,
 154      'r.dcauchy' => sub { &reorder_handler ("3,1,2", @_); },
 155      'r.pcauchy' => sub { &reorder_handler ("3,1,2", @_); },
 156      'r.qcauchy' => sub { &reorder_handler ("3,1,2", @_); },
 157      'r.dchisq' => sub { &reorder_handler ("2,1", @_); },
 158      'r.pchisq' => sub { &reorder_handler ("2,1", @_); },
 159      'r.qchisq' => sub { &reorder_handler ("2,1", @_); },
 160      'r.dexp' => sub { my ($f,$pa) = @_; &def_expr_handler ($f,["$pa->[2]-$pa->[0]","1/$pa->[1]"]); },
 161      'r.pexp' => sub { my ($f,$pa) = @_; &def_expr_handler ($f,["$pa->[2]-$pa->[0]","1/$pa->[1]"]); },
 162      'r.qexp' => sub { my ($f,$pa) = @_; &def_expr_handler ($f,[$pa->[2],"1/$pa->[1]"]) . "+$pa->[0]"; },
 163      'r.dgamma' => sub { &reorder_handler ("3,1,2", @_); },
 164      'r.pgamma' => sub { &reorder_handler ("3,1,2", @_); },
 165      'r.qgamma' => sub { &reorder_handler ("3,1,2", @_); },
 166      'laplace' => sub { my ($f,$pa) = @_; &def_expr_handler ($f,["$pa->[2]-$pa->[0]",$pa->[1]]); },
 167      'logistic' => sub { my ($f,$pa) = @_; &def_expr_handler ($f,["$pa->[2]-$pa->[0]",$pa->[1]]); },
 168      'r.dlnorm' => sub { &reorder_handler ("3,1,2", @_); },
 169      'r.plnorm' => sub { &reorder_handler ("3,1,2", @_); },
 170      'r.qlnorm' => sub { &reorder_handler ("3,1,2", @_); },
 171      'pareto' => sub { &reorder_handler ("3,2,1", @_); },
 172      'r.dweibull' => sub { &reorder_handler ("3,1,2", @_); },
 173      'r.pweibull' => sub { &reorder_handler ("3,1,2", @_); },
 174      'r.qweibull' => sub { &reorder_handler ("3,1,2", @_); },
 175      'r.dbinom' => sub { &reorder_handler ("3,2,1", @_); },
 176      'r.pbinom' => sub { &reorder_handler ("3,2,1", @_); },
 177      'r.dpois' => sub { &reorder_handler ("2,1", @_); },
 178      'r.ppois' => sub { &reorder_handler ("2,1", @_); },
 179      'r.dnbinom' => sub { &reorder_handler ("3,2,1", @_); },
 180      'r.pnbinom' => sub { &reorder_handler ("3,2,1", @_); },
 181      'r.dhyper' => sub { &reorder_handler ("4,1,2,3", @_); },
 182      'r.phyper' => sub { &reorder_handler ("4,1,2,3", @_); },
 183      'rayleigh' => sub { &reorder_handler ("2,1", @_); },
 184      'r.dnorm' => sub { &reorder_handler ("3,1,2", @_); },
 185      'r.pnorm' => sub { &reorder_handler ("3,1,2", @_); },
 186      'r.qnorm' => sub { &reorder_handler ("3,1,2", @_); },
 187      'r.dbeta' => sub { &reorder_handler ("3,1,2", @_); },
 188      'r.pbeta' => sub { &reorder_handler ("3,1,2", @_); },
 189      'r.qbeta' => sub { &reorder_handler ("3,1,2", @_); },
 190      'r.dt' => sub { &reorder_handler ("2,1", @_); },
 191      'r.pt' => sub { &reorder_handler ("2,1", @_); },
 192      'r.qt' => sub { &reorder_handler ("2,1", @_); },
 193      'r.df' => sub { &reorder_handler ("3,1,2", @_); },
 194      'r.pf' => sub { &reorder_handler ("3,1,2", @_); },
 195      'r.qf' => sub { &reorder_handler ("3,1,2", @_); },
 196      'besselj0' => sub { my ($f,$pa) = @_; &def_expr_handler ('besselj',[@$pa,0]); },
 197      'besselj1' => sub { my ($f,$pa) = @_; &def_expr_handler ('besselj',[@$pa,1]); },
 198      'besselj' => sub { &reorder_handler ("2,1", @_); },
 199      'bessely0' => sub { my ($f,$pa) = @_; &def_expr_handler ('bessely',[@$pa,0]); },
 200      'bessely1' => sub { my ($f,$pa) = @_; &def_expr_handler ('bessely',[@$pa,1]); },
 201      'bessely' => sub { &reorder_handler ("2,1", @_); },
 202      'besseli0' => sub { my ($f,$pa) = @_; &def_expr_handler ('besseli',[@$pa,0]); },
 203      'besseli1' => sub { my ($f,$pa) = @_; &def_expr_handler ('besseli',[@$pa,1]); },
 204      'besseli' => sub { &reorder_handler ("2,1", @_); },
 205      'besselk0' => sub { my ($f,$pa) = @_; &def_expr_handler ('besselk',[@$pa,0]); },
 206      'besselk1' => sub { my ($f,$pa) = @_; &def_expr_handler ('besselk',[@$pa,1]); },
 207      'besselk' => sub { &reorder_handler ("2,1", @_); },
 208      'exp2' => sub { my ($f,$pa) = @_; &def_expr_handler ('power',[2,@$pa]); },
 209      'exp10' => sub { my ($f,$pa) = @_; &def_expr_handler ('power',[10,@$pa]); },
 210      'ln' => \&positive_handler,
 211      'log10' => \&positive_handler,
 212      'log2' => \&positive_handler,
 213     );
 214
 215 my %constants =
 216     # Use lower case.
 217     ('pi_1' => 3.1415926535897932385,
 218      'pi_2' => 1.5707963267948966192,
 219      'pi_3' => 1.0471975511965977462,
 220      'pi_4' => 0.78539816339744830962,
 221      'pi_6' => 0.52359877559829887308,
 222      'sqrt2' => 1.4142135623730950488,
 223      '-sqrt2' => -1.4142135623730950488,
 224      'sqrt3' => 1.7320508075688772935,
 225      '-sqrt3' => 1.7320508075688772935,
 226      'sqrt_5' => 0.7071067811865475244,
 227     );
 228
 229 # -----------------------------------------------------------------------------
 230
 231 my $last_func = '';
 232 my @test_lines = ();
 233
 234 sub output_test {
 235     my ($gfunc,$expr,$res) = @_;
 236
 237     my $gfunc0 = ($gfunc eq $last_func) ? '' : $gfunc;
 238     $res = "=$res" if $res =~ m{[*/]};
 239
 240     my $N = 1 + @test_lines;
 241     push @test_lines, "\"$gfunc0\",\"=$expr\",\"$res\",\"=IF(B$N=C$N,\"\"\"\",IF(C$N=0,-LOG10(ABS(B$N)),-LOG10(ABS((B$N-C$N)/C$N))))\"";
 242
 243     $last_func = $gfunc;
 244 }
 245
 246 # -----------------------------------------------------------------------------
 247
 248 sub interpret_number {
 249     my ($s) = @_;
 250
 251     if ($s =~ /^[-+]?(\d+\.?|\d*\.\d+)([eE][-+]?\d+)?$/) {
 252         return $s;
 253     } else {
 254         return undef;
 255     }
 256 }
 257
 258 # -----------------------------------------------------------------------------
 259
 260 sub reorder_handler {
 261     my ($order,$f,$pargs) = @_;
 262
 263     my @res;
 264     foreach (split (',',$order)) {
 265         push @res, $pargs->[$_ - 1];
 266     }
 267
 268     return &def_expr_handler ($f,\@res);
 269 }
 270
 271 sub non_negative_handler {
 272     my ($f,$pargs) = @_;
 273
 274     foreach (@$pargs) {
 275         my $x = &interpret_number ($_);
 276         return undef unless defined ($x) && $x >= 0;
 277     }
 278
 279     return &def_expr_handler ($f,$pargs);
 280 }
 281
 282 sub positive_handler {
 283     my ($f,$pargs) = @_;
 284
 285     foreach (@$pargs) {
 286         my $x = &interpret_number ($_);
 287         return undef unless defined ($x) && $x > 0;
 288     }
 289
 290     return &def_expr_handler ($f,$pargs);
 291 }
 292
 293 # -----------------------------------------------------------------------------
 294
 295 sub simplify_val {
 296     my ($val,$pvars) = @_;
 297
 298     $val =~ s/^\s+//;
 299     $val =~ s/\s+$//;
 300
 301     # Avoid a perl bug that underflows 0.153e-305
 302     while ($val =~ /^(.*)\b0\.(\d)(\d*)[eE]-(\d+)\b(.*)$/) {
 303         $val = "$1$2.$3e-" . ($4 + 1) . $5;
 304     }
 305
 306     $val =~ s/\bldexp\s*\(\s*([-+.eE0-9_]+)\s*[,;]\s*([-+]?\d+)\s*\)/($1*2^$2)/g;
 307
 308     if ($val =~ m{^[-+*/^() .eE0-9]+$}) {
 309         if ($val =~ /^[-+]?[0-9.]+[eE][-+]?\d+$/) {
 310             if ($val == 0) {
 311                 print STDERR "DEBUG: $val --> 0\n" if $debug_underflow;
 312                 return 0;
 313             }
 314             if (($val + 0) =~ /inf/ ) {
 315                 print STDERR "DEBUG: $val --> inf\n" if $debug_overflow;
 316                 return undef;
 317             }
 318         }
 319
 320         return $val;
 321     } elsif (exists $pvars->{lc $val}) {
 322         return $pvars->{lc $val};
 323     } else {
 324         print STDERR "DEBUG: Argument $val unresolved.\n" if $debug_arguments;
 325         return undef;
 326     }
 327 }
 328
 329 # -----------------------------------------------------------------------------
 330
 331 push @test_lines, ("") x 100;
 332
 333 my $func_no = 0;
 334 foreach my $f (@test_files) {
 335     my $fn = "$dir/tests/$f";
 336
 337     my ($afunc,$gfunc);
 338
 339     my %vars;
 340     my $expr;
 341
 342     my $first_row = 1 + @test_lines;
 343
 344     open (my $src, "<", $fn) or die "$0: Cannot read $fn: $!\n";
 345     while (<$src>) {
 346         last if /^implementation\b/;
 347     }
 348
 349
 350     while (<$src>) {
 351         if (/^procedure\s+test_([a-zA-Z0-9_]+)\s*;/) {
 352             $afunc = $1;
 353             $gfunc = $name_map{$afunc};
 354             printf STDERR "Reading tests for $gfunc\n" if $gfunc;
 355             %vars = %constants;
 356             next;
 357         }
 358
 359         next unless defined $gfunc;
 360
 361         if (/^end;/i) {
 362             my $last_row = @test_lines;
 363             if ($last_row >= $first_row) {
 364                 my $count = $last_row - $first_row + 1;
 365                 $test_lines[$func_no + 2] =
 366                     "$gfunc,$count,\"=min(D${first_row}:D${last_row},99)\"";
 367                 $func_no++;
 368                 $first_row = $last_row + 1;
 369             }
 370         }
 371
 372         if (s/^\s*y\s*:=\s*([a-zA-Z0-9_]+)\s*\(([^;{}]+)\)\s*;// &&
 373             $1 eq $afunc) {
 374             my $argtxt = $2;
 375
 376             $argtxt =~ s/\bldexp\s*\(\s*([-+.eE0-9_]+)\s*,\s*([-+]?\d+)\s*\)/ldexp($1;$2)/;
 377             my @args = split (',',$argtxt);
 378             my $ok = 1;
 379
 380             foreach (@args) {
 381                 $_ = &simplify_val ($_, \%vars);
 382                 if (!defined $_) {
 383                     $ok = 0;
 384                     last;
 385                 }
 386             }
 387             next unless $ok;
 388
 389             my $h = $expr_handlers{$gfunc} || \&def_expr_handler;
 390             $expr = &$h ($gfunc,\@args);
 391         }
 392
 393         while (s/^\s*([a-zA-Z0-9]+)\s*:=\s*([^;]+)\s*;//) {
 394             my $var = lc $1;
 395             my $val = $2;
 396             $val = &simplify_val ($val, \%vars);
 397             if (defined $val) {
 398                 $vars{$var} = $val;
 399             } else {
 400                 delete $vars{$var};
 401             }
 402         }
 403
 404         if (/^\s*test(rel|rele|abs)\s*/ && exists $vars{'f'} && defined ($expr)) {
 405             &output_test ($gfunc, $expr, $vars{'f'});
 406             $expr = undef;
 407         }
 408
 409         if (/^\s*TData:\s*array/ ... /;\s*(\{.*\}\s*)*$/) {
 410             if (/^\s*\(\s*tx\s*:([^;]+);\s*ty\s*:([^\)]+)\)\s*,?\s*$/) {
 411                 my $tx = $1;
 412                 my $ty = $2;
 413                 my $x = &simplify_val ($tx, \%constants);
 414                 my $y = &simplify_val ($ty, \%constants);
 415                 my $h = $expr_handlers{$gfunc} || \&def_expr_handler;
 416                 my $expr = (defined $x) ? &$h ($gfunc,[$x]) : undef;
 417                 if (defined ($expr) && defined ($y)) {
 418                     &output_test ($gfunc, $expr, $y);
 419                 }
 420             }
 421         }
 422     }
 423 }
 424 {
 425     my $r0 = 3;
 426     my $r1 = $func_no + 2;
 427     $test_lines[0] = "\"Function\",\"Number of Tests\",\"Accuracy\",\"=min(C${r0}:C${r1})\"";
 428 }
 429
 430 foreach (@test_lines) {
 431     print "$_\n";
 432 }
 433
 434 # -----------------------------------------------------------------------------