File : s-lidosq-ada.adb
1 ------------------------------------------------------------------------------
2 -- --
3 -- GNAT COMPILER COMPONENTS --
4 -- --
5 -- S Y S T E M . L I B M _ D O U B L E . S Q R T --
6 -- --
7 -- B o d y --
8 -- --
9 -- Copyright (C) 2014-2015, Free Software Foundation, Inc. --
10 -- --
11 -- GNAT is free software; you can redistribute it and/or modify it under --
12 -- terms of the GNU General Public License as published by the Free Soft- --
13 -- ware Foundation; either version 3, or (at your option) any later ver- --
14 -- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
15 -- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
16 -- or FITNESS FOR A PARTICULAR PURPOSE. --
17 -- --
18 -- --
19 -- --
20 -- --
21 -- --
22 -- You should have received a copy of the GNU General Public License and --
23 -- a copy of the GCC Runtime Library Exception along with this program; --
24 -- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see --
25 -- <http://www.gnu.org/licenses/>. --
26 -- --
27 -- GNAT was originally developed by the GNAT team at New York University. --
28 -- Extensive contributions were provided by Ada Core Technologies Inc. --
29 -- --
30 ------------------------------------------------------------------------------
31
32 -- This is the Ada Cert Math specific implementation of sqrt (powerpc)
33
34 with Ada.Unchecked_Conversion;
35
36 with System.Machine_Code;
37
38 package body System.Libm_Double.Squareroot is
39
40 function Rsqrt (X : Long_Float) return Long_Float;
41 -- Compute the reciprocal square root. There are two reasons for computing
42 -- the reciprocal square root instead of computing directly the square
43 -- root: PowerPc provides an instruction (fsqrte) to compute an estimate of
44 -- the reciprocal (with 5 bits of precision), and the Newton-Raphson method
45 -- is more efficient on the reciprocal than on the direct root (because the
46 -- direct root needs divisions, while the reciprocal does not). Note that
47 -- PowerPc core e300 doesn't support the direct square root operation.
48
49 -----------
50 -- Rsqrt --
51 -----------
52
53 function Rsqrt (X : Long_Float) return Long_Float is
54 X_Half : constant Long_Float := X * 0.5;
55 Y, Y1 : Long_Float;
56
57 begin
58 if Standard'Target_Name = "powerpc-elf" then
59
60 -- On powerpc, the precision of fsqrte is at least 5 binary digits
61
62 System.Machine_Code.Asm ("frsqrte %0,%1",
63 Outputs => Long_Float'Asm_Output ("=f", Y),
64 Inputs => Long_Float'Asm_Input ("f", X));
65 else
66 -- Provide the exact result for 1.0
67
68 if X = 1.0 then
69 return X;
70 end if;
71
72 -- Use the method described in Fast Inverse Square Root article by
73 -- Chris Lomont (http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf),
74 -- although the code was known before that article.
75
76 declare
77 type Unsigned_Long is mod 2**64;
78
79 function To_Unsigned_Long is new Ada.Unchecked_Conversion
80 (Long_Float, Unsigned_Long);
81 function From_Unsigned_Long is new Ada.Unchecked_Conversion
82 (Unsigned_Long, Long_Float);
83 U : Unsigned_Long;
84
85 begin
86 U := To_Unsigned_Long (X);
87 U := 16#5fe6ec85_e7de30da# - (U / 2);
88 Y := From_Unsigned_Long (U);
89
90 -- Precision is about 4 digits
91 end;
92 end if;
93
94 -- Newton iterations: X <- X - F(X)/F'(X)
95 -- Here F(X) = 1/X^2 - A, so F'(X) = -2/X^3
96 -- So: X <- X - (1/X^2 - A) / (-2/X^3)
97 -- <- X + .5(X - A*X^3)
98 -- <- X + .5*X*(1 - A*X^2)
99 -- <- X (1 + .5 - .5*A*X^2)
100 -- <- X(1.5 - .5*A*X^2)
101 -- Precision is doubled at each iteration.
102
103 -- Refine: 10 digits (PowerPc) or 8 digits (fast method)
104
105 Y := Y * (1.5 - X_Half * Y * Y);
106
107 -- Refine: 20 digits (PowerPc) or 16 digits (fast method)
108
109 Y := Y * (1.5 - X_Half * Y * Y);
110
111 -- Refine: 40 digits (PowerPc) or 32 digits (fast method)
112
113 Y := Y * (1.5 - X_Half * Y * Y);
114
115 -- Refine (beyond the precision of Long_Float)
116
117 Y1 := Y * (1.5 - X_Half * Y * Y);
118
119 if Y = Y1 then
120 return Y1;
121 else
122 Y := Y1;
123 end if;
124
125 -- Empirical tests show the above iterations are inadequate in some
126 -- cases and that two more iterations are needed to converge. Other
127 -- algorithms may need to be explored. ???
128
129 Y1 := Y * (1.5 - X_Half * Y * Y);
130
131 if Y = Y1 then
132 return Y1;
133 else
134 Y := Y1;
135 end if;
136
137 Y := Y * (1.5 - X_Half * Y * Y);
138
139 -- This algorithm doesn't always provide exact results. For example,
140 -- Sqrt (25.0) /= 5.0 exactly (it's wrong in the last bit).
141
142 return Y;
143 end Rsqrt;
144
145 ----------
146 -- Sqrt --
147 ----------
148
149 function Sqrt (X : Long_Float) return Long_Float is
150 begin
151 if X <= 0.0 then
152 if X = 0.0 then
153 return X;
154 else
155 return NaN;
156 end if;
157
158 elsif not Long_Float'Machine_Overflows and then X = Infinity then
159 -- Note that if Machine_Overflow is True Infinity won't return.
160 -- But in that case, we can assume that X is not infinity.
161 return X;
162
163 else
164 return X * Rsqrt (X);
165 end if;
166 end Sqrt;
167
168 end System.Libm_Double.Squareroot;