line |
stmt |
bran |
cond |
sub |
pod |
time |
code |
1
|
|
|
|
|
|
|
## Math/MatrixDecomposition/Util.pm --- utility functions. |
2
|
|
|
|
|
|
|
|
3
|
|
|
|
|
|
|
# Copyright (C) 2010 Ralph Schleicher. All rights reserved. |
4
|
|
|
|
|
|
|
|
5
|
|
|
|
|
|
|
# This program is free software; you can redistribute it and/or |
6
|
|
|
|
|
|
|
# modify it under the same terms as Perl itself. |
7
|
|
|
|
|
|
|
|
8
|
|
|
|
|
|
|
## Code: |
9
|
|
|
|
|
|
|
|
10
|
|
|
|
|
|
|
package Math::MatrixDecomposition::Util; |
11
|
|
|
|
|
|
|
|
12
|
3
|
|
|
3
|
|
63432
|
use strict; |
|
3
|
|
|
|
|
7
|
|
|
3
|
|
|
|
|
94
|
|
13
|
3
|
|
|
3
|
|
16
|
use warnings; |
|
3
|
|
|
|
|
5
|
|
|
3
|
|
|
|
|
93
|
|
14
|
3
|
|
|
3
|
|
16
|
use Exporter qw(import); |
|
3
|
|
|
|
|
12
|
|
|
3
|
|
|
|
|
112
|
|
15
|
3
|
|
|
3
|
|
594
|
use POSIX qw(fmod); |
|
3
|
|
|
|
|
6502
|
|
|
3
|
|
|
|
|
17
|
|
16
|
|
|
|
|
|
|
|
17
|
|
|
|
|
|
|
BEGIN |
18
|
|
|
|
|
|
|
{ |
19
|
3
|
|
|
3
|
|
1694
|
our $VERSION = '1.04'; |
20
|
3
|
|
|
|
|
19
|
our @EXPORT_OK = qw(eps mod min max sign hypot cdiv); |
21
|
3
|
|
|
|
|
1291
|
our %EXPORT_TAGS = (all => [@EXPORT_OK]); |
22
|
|
|
|
|
|
|
} |
23
|
|
|
|
|
|
|
|
24
|
|
|
|
|
|
|
# Machine precision. |
25
|
|
|
|
|
|
|
my $epsilon = 1.0; |
26
|
|
|
|
|
|
|
|
27
|
38
|
|
|
38
|
|
216
|
*eps = sub () { $epsilon; }; |
28
|
|
|
|
|
|
|
|
29
|
|
|
|
|
|
|
INIT |
30
|
|
|
|
|
|
|
{ |
31
|
3
|
|
|
3
|
|
235
|
my $tem; |
32
|
|
|
|
|
|
|
|
33
|
3
|
|
|
|
|
6
|
while (1) |
34
|
|
|
|
|
|
|
{ |
35
|
159
|
|
|
|
|
172
|
$tem = 1.0 + $epsilon / 2.0; |
36
|
159
|
100
|
|
|
|
233
|
last if $tem == 1.0; |
37
|
156
|
|
|
|
|
159
|
$epsilon /= 2.0; |
38
|
|
|
|
|
|
|
} |
39
|
|
|
|
|
|
|
} |
40
|
|
|
|
|
|
|
|
41
|
|
|
|
|
|
|
# Remainder of floating-point division. |
42
|
|
|
|
|
|
|
*mod = \&fmod; |
43
|
|
|
|
|
|
|
|
44
|
|
|
|
|
|
|
# Minimum value. |
45
|
|
|
|
|
|
|
sub min ($$) |
46
|
|
|
|
|
|
|
{ |
47
|
22
|
|
|
22
|
1
|
39
|
my ($a, $b) = @_; |
48
|
|
|
|
|
|
|
|
49
|
22
|
100
|
|
|
|
60
|
$a < $b ? $a : $b; |
50
|
|
|
|
|
|
|
} |
51
|
|
|
|
|
|
|
|
52
|
|
|
|
|
|
|
# Maximum value. |
53
|
|
|
|
|
|
|
sub max ($$) |
54
|
|
|
|
|
|
|
{ |
55
|
3
|
|
|
3
|
1
|
6
|
my ($a, $b) = @_; |
56
|
|
|
|
|
|
|
|
57
|
3
|
100
|
|
|
|
28
|
$a > $b ? $a : $b; |
58
|
|
|
|
|
|
|
} |
59
|
|
|
|
|
|
|
|
60
|
|
|
|
|
|
|
# Transfer sign. |
61
|
|
|
|
|
|
|
sub sign ($$) |
62
|
|
|
|
|
|
|
{ |
63
|
34
|
|
|
34
|
1
|
60
|
my ($a, $b) = @_; |
64
|
|
|
|
|
|
|
|
65
|
34
|
100
|
|
|
|
112
|
($a < 0) == ($b < 0) ? $a : -$a; |
66
|
|
|
|
|
|
|
} |
67
|
|
|
|
|
|
|
|
68
|
|
|
|
|
|
|
# Length of the hypotenuse of a right triangle. |
69
|
|
|
|
|
|
|
sub hypot ($$) |
70
|
|
|
|
|
|
|
{ |
71
|
15
|
|
|
15
|
1
|
26
|
my $a = abs (shift); |
72
|
15
|
|
|
|
|
21
|
my $b = abs (shift); |
73
|
|
|
|
|
|
|
|
74
|
|
|
|
|
|
|
# Work variables. |
75
|
15
|
|
|
|
|
16
|
my ($s, $t); |
76
|
|
|
|
|
|
|
|
77
|
15
|
100
|
|
|
|
26
|
if ($a >= $b) |
78
|
|
|
|
|
|
|
{ |
79
|
12
|
|
|
|
|
13
|
$s = $a; |
80
|
|
|
|
|
|
|
# Avoid division by zero. |
81
|
12
|
50
|
|
|
|
21
|
$t = ($a == $b ? 1.0 : $b / $a); |
82
|
|
|
|
|
|
|
} |
83
|
|
|
|
|
|
|
else |
84
|
|
|
|
|
|
|
{ |
85
|
3
|
|
|
|
|
4
|
$s = $b; |
86
|
3
|
|
|
|
|
5
|
$t = $a / $b; |
87
|
|
|
|
|
|
|
} |
88
|
|
|
|
|
|
|
|
89
|
15
|
|
|
|
|
40
|
$s * sqrt (1.0 + $t * $t); |
90
|
|
|
|
|
|
|
} |
91
|
|
|
|
|
|
|
|
92
|
|
|
|
|
|
|
# Complex division. |
93
|
|
|
|
|
|
|
sub cdiv ($$$$) |
94
|
|
|
|
|
|
|
{ |
95
|
1
|
|
|
1
|
1
|
4
|
my ($a_re, $a_im, $b_re, $b_im) = @_; |
96
|
|
|
|
|
|
|
|
97
|
|
|
|
|
|
|
# Work variables. |
98
|
1
|
|
|
|
|
1
|
my ($r, $d, @z); |
99
|
|
|
|
|
|
|
|
100
|
1
|
50
|
|
|
|
5
|
if (abs ($b_re) > abs ($b_im)) |
101
|
|
|
|
|
|
|
{ |
102
|
0
|
|
|
|
|
0
|
$r = $b_im / $b_re; |
103
|
0
|
|
|
|
|
0
|
$d = $b_re + $r * $b_im; |
104
|
0
|
|
|
|
|
0
|
@z = (($a_re + $r * $a_im) / $d, |
105
|
|
|
|
|
|
|
($a_im - $r * $a_re) / $d); |
106
|
|
|
|
|
|
|
} |
107
|
|
|
|
|
|
|
else |
108
|
|
|
|
|
|
|
{ |
109
|
1
|
|
|
|
|
2
|
$r = $b_re / $b_im; |
110
|
1
|
|
|
|
|
3
|
$d = $b_im + $r * $b_re; |
111
|
1
|
|
|
|
|
6
|
@z = (($r * $a_re + $a_im) / $d, |
112
|
|
|
|
|
|
|
($r * $a_im - $a_re) / $d); |
113
|
|
|
|
|
|
|
} |
114
|
|
|
|
|
|
|
|
115
|
1
|
|
|
|
|
4
|
@z; |
116
|
|
|
|
|
|
|
} |
117
|
|
|
|
|
|
|
|
118
|
|
|
|
|
|
|
1; |
119
|
|
|
|
|
|
|
|
120
|
|
|
|
|
|
|
__END__ |