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