0032402: Coding Rules - eliminate msvc warning C4668 (symbol is not defined as a...
[occt.git] / src / Standard / Standard_Real.cxx
1 // Copyright (c) 1998-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
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
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.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
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>
20 #include <Standard_Stream.hxx>
21
22 static const Standard_Real ACosLimit = 1. + Epsilon(1.);
23
24 //============================================================================
25 // function : HashCode
26 // purpose  :
27 //============================================================================
28 Standard_Integer HashCode (const Standard_Real theReal, const Standard_Integer theUpperBound)
29 {
30   if (theUpperBound < 1)
31   {
32     throw Standard_RangeError ("Try to apply HashCode method with negative or null argument.");
33   }
34   union
35   {
36     Standard_Real    R;
37     Standard_Integer I[2];
38   } U;
39
40   //  U.R = Abs(me); // Treat me = -0.0 ADN 27/11/97
41   U.R = theReal;
42
43   return HashCode (U.I[0] ^ U.I[1], theUpperBound);
44 }
45
46 //-------------------------------------------------------------------
47 // ACos : Returns the value of the arc cosine of a real
48 //-------------------------------------------------------------------
49 Standard_Real ACos (const Standard_Real Value) 
50
51   if ((Value < -ACosLimit) || (Value > ACosLimit)){
52     throw Standard_RangeError();
53   }
54   else if (Value > 1.)
55   {
56     return 0.; //acos(1.)
57   }
58   else if (Value < -1.)
59   {
60     return M_PI; //acos(-1.)
61   }
62   return acos(Value);
63 }
64
65 //-------------------------------------------------------------------
66 // ACosApprox : Returns the approximate value of the arc cosine of a real.
67 //              The max error is about 1 degree near Value=0.
68 //-------------------------------------------------------------------
69
70 inline Standard_Real apx_for_ACosApprox (const Standard_Real x)
71 {
72   return  (-0.000007239283986332 +
73     x * (2.000291665285952400 +
74     x * (0.163910606547823220 +
75     x * (0.047654245891495528 -
76     x * (0.005516443930088506 +
77     0.015098965761299077 * x))))) / sqrt(2*x);
78 }
79
80 Standard_Real ACosApprox (const Standard_Real Value)
81 {
82   double XX;
83   if (Value < 0.) {
84     XX = 1.+Value;
85     if (XX < RealSmall())
86       return 0.;
87     return M_PI - apx_for_ACosApprox(XX);
88   }
89   XX = 1.-Value;
90   if (XX < RealSmall())
91     return 0.;
92   return apx_for_ACosApprox(XX);
93
94 // The code above is the same but includes 2 comparisons instead of 3
95 //   Standard_Real xn = 1.+Value;
96 //   Standard_Real xp = 1.-Value;
97 //   if (xp < RealSmall() || xn < RealSmall())
98 //     return 0.;
99 //   if (Value < 0.)
100 //     return M_PI - apx_for_ACosApprox (xn);
101 //   return apx_for_ACosApprox (xp);
102 }
103
104 //-------------------------------------------------------------------
105 // ASin : Returns the value of the arc sine of a real
106 //-------------------------------------------------------------------
107 Standard_Real ASin (const Standard_Real Value) 
108
109   if ((Value < -ACosLimit) || (Value > ACosLimit)){
110     throw Standard_RangeError();
111   }
112   else if (Value > 1.)
113   {
114     return M_PI_2; //asin(1.)
115   }
116   else if (Value < -1.)
117   {
118     return -M_PI_2; //asin(-1.)
119   }
120   return asin(Value);
121 }
122
123 //-------------------------------------------------------------------
124 // ATan2 : Returns the arc tangent of a real divide by an another real
125 //-------------------------------------------------------------------
126 Standard_Real ATan2 (const Standard_Real Value, const Standard_Real Other) 
127
128   if ( Value == 0. && Other == 0. ){
129     throw Standard_NullValue();
130   }
131   return atan2(Value,Other); 
132 }
133
134 //-------------------------------------------------------------------
135 // Sign : Returns |a| if B >= 0; -|a| if b < 0.
136 //-------------------------------------------------------------------
137 Standard_Real Sign(const Standard_Real a, const Standard_Real b)
138 {
139   if (b >= 0.0) {
140     return Abs(a);
141   } else {
142     return (-1.0 * Abs(a));
143   }
144 }
145
146 //==========================================================================
147 //===== The special routines for "IEEE" and different hardware =============
148 //==========================================================================
149 union RealMap {
150   double real;
151   unsigned int map[2];
152 };
153
154 //--------------------------------------------------------------------
155 // HardwareHighBitsOfDouble :  
156 //    Returns 1 if the low bits are at end.   (example: decmips and ALPHA )
157 //    Returns 0 if the low bits are at begin. (example: sun, sgi, ...)
158 //--------------------------------------------------------------------
159 static int HardwareHighBitsOfDouble()
160 {
161   RealMap MaxDouble;
162   MaxDouble.real = DBL_MAX;
163   //=========================================================
164   // representation of the max double in IEEE is
165   //      "7fef ffff ffff ffff"   for the big endians.
166   //      "ffff ffff 7fef ffff"   for the little endians.
167   //=========================================================
168
169   if(MaxDouble.map[1] != 0xffffffff){
170     return 1;
171   } else {
172     return 0;
173   }
174 }
175
176 //--------------------------------------------------------------------
177 // HardwareLowBitsOfDouble :  
178 //    Returns 0 if the low bits are at end.   (example: decmips )
179 //    Returns 1 if the low bits are at begin. (example: sun, sgi, ...)
180 //--------------------------------------------------------------------
181 static int HardwareLowBitsOfDouble()
182 {
183   RealMap MaxDouble;
184   MaxDouble.real = DBL_MAX;
185   //=========================================================
186   // representation of the max double in IEEE is
187   //      "7fef ffff ffff ffff"   for the big endians.
188   //      "ffff ffff 7fef ffff"   for the little endians.
189   //=========================================================
190
191   if(MaxDouble.map[1] != 0xffffffff){
192     return 0;
193   } else {
194     return 1;
195   }
196 }
197
198 static int HighBitsOfDouble = HardwareHighBitsOfDouble();
199 static int LowBitsOfDouble = HardwareLowBitsOfDouble();
200
201 double NextAfter(const double x, const double y)
202 {
203   RealMap res;
204
205   res.real=x;
206   
207   if (x == 0.0) {
208         return DBL_MIN;
209   }
210   if(x==y) {
211     //=========================================
212     //   -oo__________0___________+oo
213     //               x=y
214     //  The direction is "Null", so there is nothing after
215     //=========================================
216
217   } else if (((x<y) && (x>=0.0)) || ((x>y) && (x<0.0))) {
218     //=========================================
219     //   -oo__________0___________+oo
220     //        y <- x     x -> y
221     //
222     //=========================================
223     if (res.map[LowBitsOfDouble]==0xffffffff) {
224       res.map[LowBitsOfDouble]=0;
225       res.map[HighBitsOfDouble]++;
226     } else {
227       res.map[LowBitsOfDouble]++;
228     }
229   } else {
230     //=========================================
231     //   -oo__________0___________+oo
232     //        x -> y     y <- x
233     //
234     //=========================================
235     if (res.map[LowBitsOfDouble]==0) {
236       if (res.map[HighBitsOfDouble]==0) {
237         res.map[HighBitsOfDouble]=0x80000000;
238         res.map[LowBitsOfDouble]=0x00000001;
239       } else {
240         res.map[LowBitsOfDouble]=0xffffffff;
241         res.map[HighBitsOfDouble]--;
242       }
243     } else {
244       res.map[LowBitsOfDouble]--;
245     }
246   }
247   return res.real;
248 }
249
250 //-------------------------------------------------------------------
251 // ATanh : Returns the value of the hyperbolic arc tangent of a real
252 //-------------------------------------------------------------------
253 Standard_Real     ATanh(const Standard_Real Value) 
254
255   if ( (Value <= -1.) || (Value >= 1.) ){
256 #ifdef OCCT_DEBUG
257     std::cout << "Illegal argument in ATanh" << std::endl ;
258 #endif
259     throw Standard_NumericError("Illegal argument in ATanh");
260   }
261 #if defined(__QNX__)
262   return std::atanh(Value);
263 #else
264   return atanh(Value);
265 #endif
266 }
267
268 //-------------------------------------------------------------------
269 // ACosh : Returns the hyperbolic Arc cosine of a real
270 //-------------------------------------------------------------------
271 Standard_Real     ACosh (const Standard_Real Value) 
272
273   if ( Value < 1. ){
274 #ifdef OCCT_DEBUG
275     std::cout << "Illegal argument in ACosh" << std::endl ;
276 #endif
277     throw Standard_NumericError("Illegal argument in ACosh");
278   }
279 #if defined(__QNX__)
280   return std::acosh(Value);
281 #else
282   return acosh(Value);
283 #endif
284 }
285
286 //-------------------------------------------------------------------
287 // Cosh : Returns the hyperbolic cosine of a real
288 //-------------------------------------------------------------------
289 Standard_Real     Cosh (const Standard_Real Value) 
290
291   if ( Abs(Value) > 0.71047586007394394e+03 ){
292 #ifdef OCCT_DEBUG
293     std::cout << "Result of Cosh exceeds the maximum value Standard_Real" << std::endl ;
294 #endif
295     throw Standard_NumericError("Result of Cosh exceeds the maximum value Standard_Real");
296   } 
297   return cosh(Value); 
298 }
299
300 //-------------------------------------------------------------------
301 // Sinh : Returns the hyperbolicsine of a real
302 //-------------------------------------------------------------------
303 Standard_Real     Sinh (const Standard_Real Value) 
304
305   if ( Abs(Value) > 0.71047586007394394e+03 ){
306 #ifdef OCCT_DEBUG
307     std::cout << "Result of Sinh exceeds the maximum value Standard_Real" << std::endl ;
308 #endif
309     throw Standard_NumericError("Result of Sinh exceeds the maximum value Standard_Real");
310   } 
311   return sinh(Value); 
312 }
313
314 //-------------------------------------------------------------------
315 // Log : Returns the naturaOPl logarithm of a real
316 //-------------------------------------------------------------------
317 Standard_Real     Log (const Standard_Real Value) 
318 {   if ( Value <= 0. ){
319 #ifdef OCCT_DEBUG
320     std::cout << "Illegal argument in Log" << std::endl ;
321 #endif
322     throw Standard_NumericError("Illegal argument in Log");
323   } 
324  return log(Value); 
325 }
326 //-------------------------------------------------------------------
327 // Sqrt : Returns the square root of a real
328 //-------------------------------------------------------------------
329 Standard_Real     Sqrt (const Standard_Real Value) 
330
331   if (  Value < 0. ){
332 #ifdef OCCT_DEBUG
333     std::cout << "Illegal argument in Sqrt" << std::endl ;
334 #endif
335     throw Standard_NumericError("Illegal argument in Sqrt");
336   } 
337  return sqrt(Value); 
338 }
339