1# -*- mode: perl; -*-
2
3use strict;
4use warnings;
5
6use Test::More tests => 3528;
7
8use Math::Complex ();
9
10use Math::BigFloat;
11
12my $inf = $Math::Complex::Inf;
13my $nan = $inf - $inf;
14
15my $class = 'Math::BigFloat';
16
17# isnan X
18
19sub isnan { $_[0] != $_[0] }
20
21# linspace MIN, MAX, N
22#
23# Returns N linearly spaced elements from MIN to MAX.
24
25sub linspace {
26    my ($xmin, $xmax, $n) = @_;
27
28    if ($n == 0) {
29        return ();
30    } elsif ($n == 1) {
31        return ($xmin);
32    } else {
33        my $c = ($xmax - $xmin) / ($n - 1);
34        return map { $xmin + $c * $_  } 0 .. ($n - 1);
35    }
36}
37
38# logspace MIN, MAX, N
39#
40# Returns N logarithmically spaced elements from MIN to MAX.
41
42sub logspace {
43    my ($xmin, $xmax, $n) = @_;
44    if ($n == 0) {
45        return ();
46    } elsif ($n == 1) {
47        return ($xmin);
48    } else {
49        my @lin = linspace(log($xmin), log($xmax), $n);
50        my @log = map { exp } @lin;
51        $log[   0   ] = $xmin;
52        $log[ $#log ] = $xmax;
53        return @log;
54    }
55}
56
57my @x;
58@x = logspace(0.01, 12, 20);
59@x = map { sprintf "%.3g", $_ } @x;
60@x = (reverse(map( { -$_ } @x)), 0, @x, $nan);
61
62my $accu       = 16;
63my $tol        = 1e-13;
64my $max_relerr = 0;
65
66for my $ply (@x) {
67    for my $plx (@x) {
68        my $plz = CORE::atan2($ply, $plx);
69
70        # $y -> batan2($x) where $x is a scalar
71
72        {
73            my $y = $class -> new($ply);
74            $y -> batan2($plx, $accu);
75
76            my $desc = qq|\$y = $class->new("$ply");|
77                     . qq| \$y->batan2("$plx", $accu)|
78                     . qq| vs. CORE::atan2($ply, $plx)|;
79
80            if (isnan($plz)) {
81                is($y, "NaN", $desc);
82            } elsif ($plz == 0) {
83                cmp_ok($y, '==', $plz, $desc);
84            } else {
85                my $relerr = abs(($y - $plz) / $plz);
86                if (!cmp_ok($relerr, '<', $tol, "relative error of $desc")) {
87                    diag(sprintf("             CORE::atan2(...): %.15g\n" .
88                                 "  Math::BigFloat->batan2(...): %.15g\n",
89                                 $plz, $y));
90                }
91                $max_relerr = $relerr if $relerr > $max_relerr;
92            }
93        }
94
95        # $y -> batan2($x) where $x is an object
96
97        {
98            my $x = $class -> new($plx);
99            my $y = $class -> new($ply);
100            $y -> batan2($plx, $accu);
101
102            my $desc = qq|\$y = $class->new("$ply");|
103                     . qq| \$x = $class->new("$plx");|
104                     . qq| \$y->batan2(\$x, $accu)|
105                     . qq| vs. CORE::atan2($ply, $plx)|;
106
107            if (isnan($plz)) {
108                is($y, "NaN", $desc);
109            } elsif ($plz == 0) {
110                cmp_ok($y, '==', $plz, $desc);
111            } else {
112                my $relerr = abs(($y - $plz) / $plz);
113                if (!cmp_ok($relerr, '<', $tol, "relative error of $desc")) {
114                    diag(sprintf("             CORE::atan2(...): %.15g\n" .
115                                 "  Math::BigFloat->batan2(...): %.15g\n",
116                                  $plz, $y));
117                }
118                $max_relerr = $relerr if $relerr > $max_relerr;
119            }
120        }
121
122    }
123}
124
125diag("Maximum relative error = ", $max_relerr -> numify(), "\n");
126