b311480e |
1 | // Copyright (c) 1998-1999 Matra Datavision |
973c2be1 |
2 | // Copyright (c) 1999-2014 OPEN CASCADE SAS |
b311480e |
3 | // |
973c2be1 |
4 | // This file is part of Open CASCADE Technology software library. |
b311480e |
5 | // |
d5f74e42 |
6 | // This library is free software; you can redistribute it and/or modify it under |
7 | // the terms of the GNU Lesser General Public License version 2.1 as published |
973c2be1 |
8 | // by the Free Software Foundation, with special exception defined in the file |
9 | // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT |
10 | // distribution for complete text of the license and disclaimer of any warranty. |
b311480e |
11 | // |
973c2be1 |
12 | // Alternatively, this file may be used under the terms of Open CASCADE |
13 | // commercial license or contractual agreement. |
7fd59977 |
14 | |
15 | #include <float.h> |
16 | #include <Standard_Real.hxx> |
17 | #include <Standard_RangeError.hxx> |
18 | #include <Standard_NumericError.hxx> |
19 | #include <Standard_NullValue.hxx> |
7fd59977 |
20 | #include <Standard_Stream.hxx> |
7fd59977 |
21 | #include <Standard_OStream.hxx> |
7fd59977 |
22 | |
7fd59977 |
23 | // ------------------------------------------------------------------ |
24 | // Hascode : Computes a hascoding value for a given real |
25 | // ------------------------------------------------------------------ |
26 | Standard_Integer HashCode(const Standard_Real me, const Standard_Integer Upper) |
27 | { |
28 | if (Upper < 1){ |
29 | Standard_RangeError:: |
30 | Raise("Try to apply HashCode method with negative or null argument."); |
31 | } |
32 | union |
33 | { |
34 | Standard_Real R; |
35 | Standard_Integer I[2]; |
36 | } U; |
37 | // U.R = Abs(me); // Treat me = -0.0 ADN 27/11/97 |
38 | U.R = me ; |
39 | return HashCode( ( U.I[0] ^ U.I[1] ) , Upper ) ; |
40 | } |
41 | |
7fd59977 |
42 | //------------------------------------------------------------------- |
43 | // ACos : Returns the value of the arc cosine of a real |
44 | //------------------------------------------------------------------- |
45 | Standard_Real ACos (const Standard_Real Value) |
46 | { |
47 | if ( (Value < -1.) || (Value > 1.) ){ |
48 | Standard_RangeError::Raise(); |
49 | } |
50 | return acos(Value); |
51 | } |
52 | |
53 | //------------------------------------------------------------------- |
54 | // ACosApprox : Returns the approximate value of the arc cosine of a real. |
55 | // The max error is about 1 degree near Value=0. |
56 | //------------------------------------------------------------------- |
57 | |
58 | inline Standard_Real apx_for_ACosApprox (const Standard_Real x) |
59 | { |
60 | return (-0.000007239283986332 + |
61 | x * (2.000291665285952400 + |
62 | x * (0.163910606547823220 + |
63 | x * (0.047654245891495528 - |
64 | x * (0.005516443930088506 + |
65 | 0.015098965761299077 * x))))) / sqrt(2*x); |
66 | } |
67 | |
68 | Standard_Real ACosApprox (const Standard_Real Value) |
69 | { |
70 | double XX; |
71 | if (Value < 0.) { |
72 | XX = 1.+Value; |
73 | if (XX < RealSmall()) |
74 | return 0.; |
c6541a0c |
75 | return M_PI - apx_for_ACosApprox(XX); |
7fd59977 |
76 | } |
77 | XX = 1.-Value; |
78 | if (XX < RealSmall()) |
79 | return 0.; |
80 | return apx_for_ACosApprox(XX); |
81 | |
82 | // The code above is the same but includes 2 comparisons instead of 3 |
83 | // Standard_Real xn = 1.+Value; |
84 | // Standard_Real xp = 1.-Value; |
85 | // if (xp < RealSmall() || xn < RealSmall()) |
86 | // return 0.; |
87 | // if (Value < 0.) |
c6541a0c |
88 | // return M_PI - apx_for_ACosApprox (xn); |
7fd59977 |
89 | // return apx_for_ACosApprox (xp); |
90 | } |
91 | |
92 | //------------------------------------------------------------------- |
93 | // ASin : Returns the value of the arc sine of a real |
94 | //------------------------------------------------------------------- |
95 | Standard_Real ASin (const Standard_Real Value) |
96 | { |
97 | if ( Value < -1 || Value > 1 ){ |
98 | Standard_RangeError::Raise(); |
99 | } |
100 | return asin(Value); |
101 | } |
102 | |
103 | //------------------------------------------------------------------- |
104 | // ATan2 : Returns the arc tangent of a real divide by an another real |
105 | //------------------------------------------------------------------- |
106 | Standard_Real ATan2 (const Standard_Real Value, const Standard_Real Other) |
107 | { |
108 | if ( Value == 0. && Other == 0. ){ |
109 | Standard_NullValue::Raise(); |
110 | } |
111 | return atan2(Value,Other); |
112 | } |
113 | |
114 | //------------------------------------------------------------------- |
115 | // Sign : Returns |a| if B >= 0; -|a| if b < 0. |
116 | // from x in the direction y |
117 | //------------------------------------------------------------------- |
118 | Standard_Real Sign(const Standard_Real a, const Standard_Real b) |
119 | { |
120 | //==== We use the function "nextafter()" fom library "math.h" ============== |
121 | if (b >= 0.0) { |
122 | return Abs(a); |
123 | } else { |
124 | return (-1.0 * Abs(a)); |
125 | } |
126 | } |
127 | |
128 | //========================================================================== |
129 | //===== The special routines for "IEEE" and differents hardwares =========== |
130 | //========================================================================== |
131 | union RealMap { |
132 | double real; |
133 | unsigned int map[2]; |
134 | }; |
135 | |
136 | //-------------------------------------------------------------------- |
137 | // HardwareHighBitsOfDouble : |
138 | // Returns 1 if the low bits are at end. (exemple: decmips and ALPHA ) |
139 | // Returns 0 if the low bits are at begin. (exemple: sun, sgi, ...) |
140 | //-------------------------------------------------------------------- |
141 | static int HardwareHighBitsOfDouble() |
142 | { |
143 | RealMap MaxDouble; |
144 | MaxDouble.real = DBL_MAX; |
145 | //========================================================= |
146 | // reperesentation of the max double in IEEE is |
147 | // "7fef ffff ffff ffff" for the big indiens. |
148 | // "ffff ffff 7fef ffff" for the littel indiens. |
149 | //========================================================= |
150 | |
151 | if(MaxDouble.map[1] != 0xffffffff){ |
152 | return 1; |
153 | } else { |
154 | return 0; |
155 | } |
156 | } |
157 | |
158 | //-------------------------------------------------------------------- |
159 | // HardwareLowBitsOfDouble : |
160 | // Returns 0 if the low bits are at end. (exemple: decmips ) |
161 | // Returns 1 if the low bits are at begin. (exemple: sun, sgi, ...) |
162 | //-------------------------------------------------------------------- |
163 | static int HardwareLowBitsOfDouble() |
164 | { |
165 | RealMap MaxDouble; |
166 | MaxDouble.real = DBL_MAX; |
167 | //========================================================= |
168 | // reperesentation of the max double in IEEE is |
169 | // "7fef ffff ffff ffff" for the big indiens. |
170 | // "ffff ffff 7fef ffff" for the littel indiens. |
171 | //========================================================= |
172 | |
173 | if(MaxDouble.map[1] != 0xffffffff){ |
174 | return 0; |
175 | } else { |
176 | return 1; |
177 | } |
178 | } |
179 | |
180 | static int HighBitsOfDouble = HardwareHighBitsOfDouble(); |
181 | static int LowBitsOfDouble = HardwareLowBitsOfDouble(); |
182 | |
183 | double NextAfter(const double x, const double y) |
184 | { |
185 | RealMap res; |
186 | |
187 | res.real=x; |
188 | |
189 | if (x == 0.0) { |
190 | return DBL_MIN; |
191 | } |
192 | if(x==y) { |
193 | //========================================= |
194 | // -oo__________0___________+oo |
195 | // x=y |
196 | // The direction is "Null", so there is nothing after |
197 | //========================================= |
198 | |
199 | } else if (((x<y) && (x>=0.0)) || ((x>y) && (x<0.0))) { |
200 | //========================================= |
201 | // -oo__________0___________+oo |
202 | // y <- x x -> y |
203 | // |
204 | //========================================= |
205 | if (res.map[LowBitsOfDouble]==0xffffffff) { |
206 | res.map[LowBitsOfDouble]=0; |
207 | res.map[HighBitsOfDouble]++; |
208 | } else { |
209 | res.map[LowBitsOfDouble]++; |
210 | } |
211 | } else { |
212 | //========================================= |
213 | // -oo__________0___________+oo |
214 | // x -> y y <- x |
215 | // |
216 | //========================================= |
217 | if (res.map[LowBitsOfDouble]==0) { |
218 | if (res.map[HighBitsOfDouble]==0) { |
219 | res.map[HighBitsOfDouble]=0x80000000; |
220 | res.map[LowBitsOfDouble]=0x00000001; |
221 | } else { |
222 | res.map[LowBitsOfDouble]=0xffffffff; |
223 | res.map[HighBitsOfDouble]--; |
224 | } |
225 | } else { |
226 | res.map[LowBitsOfDouble]--; |
227 | } |
228 | } |
229 | return res.real; |
230 | } |
231 | |
7fd59977 |
232 | //------------------------------------------------------------------- |
233 | // ATanh : Returns the value of the hyperbolic arc tangent of a real |
234 | //------------------------------------------------------------------- |
235 | Standard_Real ATanh(const Standard_Real Value) |
236 | { |
237 | if ( (Value <= -1.) || (Value >= 1.) ){ |
238 | Standard_NumericError::Raise("Illegal agument in ATanh"); |
0797d9d3 |
239 | #ifdef OCCT_DEBUG |
7fd59977 |
240 | cout << "Illegal agument in ATanh" << endl ; |
63c629aa |
241 | #endif |
7fd59977 |
242 | } |
243 | return atanh(Value); |
244 | } |
245 | |
246 | //------------------------------------------------------------------- |
247 | // ACosh : Returns the hyperbolic Arc cosine of a real |
248 | //------------------------------------------------------------------- |
249 | Standard_Real ACosh (const Standard_Real Value) |
250 | { |
251 | if ( Value < 1. ){ |
252 | Standard_NumericError::Raise("Illegal agument in ACosh"); |
0797d9d3 |
253 | #ifdef OCCT_DEBUG |
7fd59977 |
254 | cout << "Illegal agument in ACosh" << endl ; |
63c629aa |
255 | #endif |
7fd59977 |
256 | } |
257 | return acosh(Value); |
258 | } |
259 | |
260 | //------------------------------------------------------------------- |
261 | // Log : Returns the naturaOPl logarithm of a real |
262 | //------------------------------------------------------------------- |
263 | Standard_Real Log (const Standard_Real Value) |
264 | { if ( Value <= 0. ){ |
265 | Standard_NumericError::Raise("Illegal agument in Log"); |
0797d9d3 |
266 | #ifdef OCCT_DEBUG |
7fd59977 |
267 | cout << "Illegal agument in Log" << endl ; |
63c629aa |
268 | #endif |
7fd59977 |
269 | } |
270 | return log(Value); |
271 | } |
272 | //------------------------------------------------------------------- |
273 | // Sqrt : Returns the square root of a real |
274 | //------------------------------------------------------------------- |
275 | Standard_Real Sqrt (const Standard_Real Value) |
276 | { |
277 | if ( Value < 0. ){ |
278 | Standard_NumericError::Raise("Illegal agument in Sqrt"); |
0797d9d3 |
279 | #ifdef OCCT_DEBUG |
7fd59977 |
280 | cout << "Illegal agument in Sqrt" << endl ; |
63c629aa |
281 | #endif |
7fd59977 |
282 | } |
283 | return sqrt(Value); |
284 | } |
285 | |