3 # This script processes the test cases from amath, see
4 # http://www.wolfgang-ehrhardt.de/amath_functions.html
8 my $debug_underflow = 0;
9 my $debug_overflow = 0;
10 my $debug_arguments = 1;
12 die "$0: missing amath directory\n" unless (defined $dir) && -d
$dir;
28 ('lnbeta' => 'betaln',
30 'lngamma' => 'gammaln',
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',
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',
101 'arccosh' => 'acosh',
103 'arcsinh' => 'asinh',
105 'arccoth' => 'acoth',
107 'arctanh' => 'atanh',
124 (# Magically changed to something else
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
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,
139 'bessely(1.5,-1700.5)' => 1,
142 sub def_expr_handler
{
144 my $expr = "$f(" . join (",", @
$pa) . ")";
145 return undef if exists $invalid_tests{$expr};
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
,
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,
229 # -----------------------------------------------------------------------------
235 my ($gfunc,$expr,$res) = @_;
237 my $gfunc0 = ($gfunc eq $last_func) ?
'' : $gfunc;
238 $res = "=$res" if $res =~ m{[*/]};
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))))\"";
246 # -----------------------------------------------------------------------------
248 sub interpret_number
{
251 if ($s =~ /^[-+]?(\d+\.?|\d*\.\d+)([eE][-+]?\d+)?$/) {
258 # -----------------------------------------------------------------------------
260 sub reorder_handler
{
261 my ($order,$f,$pargs) = @_;
264 foreach (split (',',$order)) {
265 push @res, $pargs->[$_ - 1];
268 return &def_expr_handler
($f,\
@res);
271 sub non_negative_handler
{
275 my $x = &interpret_number
($_);
276 return undef unless defined ($x) && $x >= 0;
279 return &def_expr_handler
($f,$pargs);
282 sub positive_handler
{
286 my $x = &interpret_number
($_);
287 return undef unless defined ($x) && $x > 0;
290 return &def_expr_handler
($f,$pargs);
293 # -----------------------------------------------------------------------------
296 my ($val,$pvars) = @_;
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;
306 $val =~ s/\bldexp\s*\(\s*([-+.eE0-9_]+)\s*[,;]\s*([-+]?\d+)\s*\)/($1*2^$2)/g;
308 if ($val =~ m{^[-+*/^() .eE0-9]+$}) {
309 if ($val =~ /^[-+]?[0-9.]+[eE][-+]?\d+$/) {
311 print STDERR
"DEBUG: $val --> 0\n" if $debug_underflow;
314 if (($val + 0) =~ /inf/ ) {
315 print STDERR
"DEBUG: $val --> inf\n" if $debug_overflow;
321 } elsif (exists $pvars->{lc $val}) {
322 return $pvars->{lc $val};
324 print STDERR
"DEBUG: Argument $val unresolved.\n" if $debug_arguments;
329 # -----------------------------------------------------------------------------
331 push @test_lines, ("") x
100;
334 foreach my $f (@test_files) {
335 my $fn = "$dir/tests/$f";
342 my $first_row = 1 + @test_lines;
344 open (my $src, "<", $fn) or die "$0: Cannot read $fn: $!\n";
346 last if /^implementation\b/;
351 if (/^procedure\s+test_([a-zA-Z0-9_]+)\s*;/) {
353 $gfunc = $name_map{$afunc};
354 printf STDERR
"Reading tests for $gfunc\n" if $gfunc;
359 next unless defined $gfunc;
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)\"";
368 $first_row = $last_row + 1;
372 if (s/^\s*y\s*:=\s*([a-zA-Z0-9_]+)\s*\(([^;{}]+)\)\s*;// &&
376 $argtxt =~ s/\bldexp\s*\(\s*([-+.eE0-9_]+)\s*,\s*([-+]?\d+)\s*\)/ldexp($1;$2)/;
377 my @args = split (',',$argtxt);
381 $_ = &simplify_val
($_, \
%vars);
389 my $h = $expr_handlers{$gfunc} || \
&def_expr_handler
;
390 $expr = &$h ($gfunc,\
@args);
393 while (s/^\s*([a-zA-Z0-9]+)\s*:=\s*([^;]+)\s*;//) {
396 $val = &simplify_val
($val, \
%vars);
404 if (/^\s*test(rel|rele|abs)\s*/ && exists $vars{'f'} && defined ($expr)) {
405 &output_test
($gfunc, $expr, $vars{'f'});
409 if (/^\s*TData:\s*array/ ... /;\s*(\{.*\}\s*)*$/) {
410 if (/^\s*\(\s*tx\s*:([^;]+);\s*ty\s*:([^\)]+)\)\s*,?\s*$/) {
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);
426 my $r1 = $func_no + 2;
427 $test_lines[0] = "\"Function\",\"Number of Tests\",\"Accuracy\",\"=min(C${r0}:C${r1})\"";
430 foreach (@test_lines) {
434 # -----------------------------------------------------------------------------