| 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__ |