File : s-exnllf.adb
1 ------------------------------------------------------------------------------
2 -- --
3 -- GNAT RUN-TIME COMPONENTS --
4 -- --
5 -- S Y S T E M . E X N _ L L F --
6 -- --
7 -- B o d y --
8 -- --
9 -- Copyright (C) 1992-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 -- Note: the reason for treating exponents in the range 0 .. 4 specially is
33 -- to ensure identical results to the static inline expansion in the case of
34 -- a compile time known exponent in this range. The use of Float'Machine and
35 -- Long_Float'Machine is to avoid unwanted extra precision in the results.
36
37 package body System.Exn_LLF is
38
39 function Exp
40 (Left : Long_Long_Float;
41 Right : Integer) return Long_Long_Float;
42 -- Common routine used if Right not in 0 .. 4
43
44 ---------------
45 -- Exn_Float --
46 ---------------
47
48 function Exn_Float
49 (Left : Float;
50 Right : Integer) return Float
51 is
52 Temp : Float;
53 begin
54 case Right is
55 when 0 =>
56 return 1.0;
57 when 1 =>
58 return Left;
59 when 2 =>
60 return Float'Machine (Left * Left);
61 when 3 =>
62 return Float'Machine (Left * Left * Left);
63 when 4 =>
64 Temp := Float'Machine (Left * Left);
65 return Float'Machine (Temp * Temp);
66 when others =>
67 return
68 Float'Machine
69 (Float (Exp (Long_Long_Float (Left), Right)));
70 end case;
71 end Exn_Float;
72
73 --------------------
74 -- Exn_Long_Float --
75 --------------------
76
77 function Exn_Long_Float
78 (Left : Long_Float;
79 Right : Integer) return Long_Float
80 is
81 Temp : Long_Float;
82 begin
83 case Right is
84 when 0 =>
85 return 1.0;
86 when 1 =>
87 return Left;
88 when 2 =>
89 return Long_Float'Machine (Left * Left);
90 when 3 =>
91 return Long_Float'Machine (Left * Left * Left);
92 when 4 =>
93 Temp := Long_Float'Machine (Left * Left);
94 return Long_Float'Machine (Temp * Temp);
95 when others =>
96 return
97 Long_Float'Machine
98 (Long_Float (Exp (Long_Long_Float (Left), Right)));
99 end case;
100 end Exn_Long_Float;
101
102 -------------------------
103 -- Exn_Long_Long_Float --
104 -------------------------
105
106 function Exn_Long_Long_Float
107 (Left : Long_Long_Float;
108 Right : Integer) return Long_Long_Float
109 is
110 Temp : Long_Long_Float;
111 begin
112 case Right is
113 when 0 =>
114 return 1.0;
115 when 1 =>
116 return Left;
117 when 2 =>
118 return Left * Left;
119 when 3 =>
120 return Left * Left * Left;
121 when 4 =>
122 Temp := Left * Left;
123 return Temp * Temp;
124 when others =>
125 return Exp (Left, Right);
126 end case;
127 end Exn_Long_Long_Float;
128
129 ---------
130 -- Exp --
131 ---------
132
133 function Exp
134 (Left : Long_Long_Float;
135 Right : Integer) return Long_Long_Float
136 is
137 Result : Long_Long_Float := 1.0;
138 Factor : Long_Long_Float := Left;
139 Exp : Integer := Right;
140
141 begin
142 -- We use the standard logarithmic approach, Exp gets shifted right
143 -- testing successive low order bits and Factor is the value of the
144 -- base raised to the next power of 2. If the low order bit or Exp is
145 -- set, multiply the result by this factor. For negative exponents,
146 -- invert result upon return.
147
148 if Exp >= 0 then
149 loop
150 if Exp rem 2 /= 0 then
151 Result := Result * Factor;
152 end if;
153
154 Exp := Exp / 2;
155 exit when Exp = 0;
156 Factor := Factor * Factor;
157 end loop;
158
159 return Result;
160
161 -- Here we have a negative exponent, and we compute the result as:
162
163 -- 1.0 / (Left ** (-Right))
164
165 -- Note that the case of Left being zero is not special, it will
166 -- simply result in a division by zero at the end, yielding a
167 -- correctly signed infinity, or possibly generating an overflow.
168
169 -- Note on overflow: The coding of this routine assumes that the
170 -- target generates infinities with standard IEEE semantics. If this
171 -- is not the case, then the code below may raise Constraint_Error.
172 -- This follows the implementation permission given in RM 4.5.6(12).
173
174 else
175 begin
176 loop
177 if Exp rem 2 /= 0 then
178 Result := Result * Factor;
179 end if;
180
181 Exp := Exp / 2;
182 exit when Exp = 0;
183 Factor := Factor * Factor;
184 end loop;
185
186 return 1.0 / Result;
187 end;
188 end if;
189 end Exp;
190
191 end System.Exn_LLF;