scid.ports.minpack.enorm source code

1 // This code has been mechanically translated from the original FORTRAN
2 // code at http://netlib.org/minpack.
3 
4 /** Authors:    Lars Tandle Kyllingstad
5     Copyright:  Copyright (c) 2009, Lars T. Kyllingstad. All rights reserved.
6     License:    Boost Licence 1.0
7 */
8 module scid.ports.minpack.enorm;
9 
10 
11 import std.math;
12 
13 
14 
15 /** Given an n-vector x, this function calculates the
16     Euclidean norm of x.
17 
18     The Euclidean norm is computed by accumulating the sum of
19     squares in three different sums. The sums of squares for the
20     small and large components are scaled so that no overflows
21     occur. Non-destructive underflows are permitted. Underflows
22     and overflows do not occur in the computation of the unscaled
23     sum of squarees for the intermediate components.
24     The definitions of small, intermediate and large components
25     depend on two constants, rdwarf and rgiant. The main
26     restrictions on these constants are that rdwarf^2 not
27     underflow and rgiant^2 not overflow. The constants
28     given here are suitable for every known computer.
29 
30     Parameters:
31         x = an input array.
32 */
33 Real enorm(Real)(size_t n, Real* x)
34 {
35     enum : Real
36     {
37         one = 1.0,
38         zero = 0.0,
39         rdwarf = 3.834e-20,
40         rgiant = 1.304e19
41     }
42 
43     Real s1 = zero;
44     Real s2 = zero;
45     Real s3 = zero;
46     Real x1max = zero;
47     Real x3max = zero;
48     immutable Real floatn = n;
49     immutable Real agiant = rgiant/floatn;
50 
51     Real xabs;
52     for (size_t i=0; i<n; i++)
53     {
54         xabs = abs(x[i]);
55         if (xabs <= rdwarf  ||  xabs >= agiant)
56         {
57             if (xabs > rdwarf)
58             {
59                 // Sum for large components.
60                 if (xabs > x1max)
61                 {
62                     s1 = one + s1*((x1max/xabs)^^2);
63                     x1max = xabs;
64                 }
65                 else
66                 {
67                     s1 += (xabs/x1max)^^2;
68                 }
69             }
70             else
71             {
72                 // Sum for small components.
73                 if (xabs > x3max)
74                 {
75                     s3 = one + s3*((x3max/xabs)^^2);
76                     x3max = xabs;
77                 }
78                 else
79                 {
80                     if (xabs != zero) s3 += (xabs/x3max)^^2;
81                 }
82             }
83         }
84         else
85         {
86             // Sum for intermediate components.
87             s2 += xabs*xabs;
88         }
89     }
90 
91     
92     // Calculation of norm.
93     Real result;
94     if (s1 != zero)
95     {
96         result = x1max*sqrt(s1+(s2/x1max)/x1max);
97     }
98     else
99     {
100         if (s2 != zero)
101         {
102             if (s2 >= x3max)
103                 result = sqrt(s2*(one+(x3max/s2)*(x3max*s3)));
104             else
105                 result = sqrt(x3max*((s2/x3max)+(x3max*s3)));
106         }
107         else
108         {
109             result = x3max*sqrt(s3);
110         }
111     }
112     
113     return result;
114 }
115 
116 
117 unittest
118 {
119     double[] v = [ 1.0, 2.0, 3.0 ];
120     double norm = enorm(3, v.ptr);
121     assert (abs(norm-sqrt(14.0)) < 1e-10);
122 }
123