File : s-lisisq-ada.adb
1 ------------------------------------------------------------------------------
2 -- --
3 -- GNAT COMPILER COMPONENTS --
4 -- --
5 -- S Y S T E M . L I B M _ S I N G 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_Single.Squareroot is
39
40 function Rsqrt (X : Float) return 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 : Float) return Float is
54 X_Half : constant Float := X * 0.5;
55 Y, Y1 : 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 => Float'Asm_Output ("=f", Y),
64 Inputs => 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 declare
73 type Unsigned is mod 2**32;
74
75 function To_Unsigned is new Ada.Unchecked_Conversion
76 (Float, Unsigned);
77 function From_Unsigned is new Ada.Unchecked_Conversion
78 (Unsigned, Float);
79 U : Unsigned;
80
81 begin
82 U := To_Unsigned (X);
83 U := 16#5f3759df# - (U / 2);
84 Y := From_Unsigned (U);
85
86 -- Precision is 4 binary digits (but the next iteration is
87 -- much better)
88 end;
89 end if;
90
91 -- Refine: 10 digits (PowerPc) or 8 digits (fast method)
92
93 Y := Y * (1.5 - X_Half * Y * Y);
94
95 -- Refine: 20 digits (PowerPc) or 16 digits (fast method)
96
97 Y := Y * (1.5 - X_Half * Y * Y);
98
99 -- Refine (beyond the precision of Float)
100
101 Y1 := Y * (1.5 - X_Half * Y * Y);
102
103 if Y = Y1 then
104 return Y1;
105 else
106 Y := Y1;
107 end if;
108
109 -- Empirical tests show the above iterations are inadequate in some
110 -- cases and that two more iterations are needed to converge. Other
111 -- algorithms may need to be explored. ???
112
113 Y1 := Y * (1.5 - X_Half * Y * Y);
114
115 if Y = Y1 then
116 return Y1;
117 else
118 Y := Y1;
119 end if;
120
121 Y := Y * (1.5 - X_Half * Y * Y);
122
123 -- This algorithm doesn't always provide exact results. For example,
124 -- Sqrt (25.0) /= 5.0 exactly (it's wrong in the last bit).
125
126 return Y;
127 end Rsqrt;
128
129 ----------
130 -- Sqrt --
131 ----------
132
133 function Sqrt (X : Float) return Float is
134 begin
135 if X <= 0.0 then
136 if X = 0.0 then
137 return X;
138 else
139 return NaN;
140 end if;
141
142 elsif not Float'Machine_Overflows and then X = Infinity then
143 -- Note that if Machine_Overflow is True Infinity won't return.
144 -- But in that case, we can assume that X is not infinity.
145 return X;
146
147 else
148 return X * Rsqrt (X);
149 end if;
150 end Sqrt;
151 end System.Libm_Single.Squareroot;