0023024: Update headers of OCCT files
[occt.git] / src / Standard / Standard_Real.cxx
1 // Copyright (c) 1998-1999 Matra Datavision
2 // Copyright (c) 1999-2012 OPEN CASCADE SAS
3 //
4 // The content of this file is subject to the Open CASCADE Technology Public
5 // License Version 6.5 (the "License"). You may not use the content of this file
6 // except in compliance with the License. Please obtain a copy of the License
7 // at http://www.opencascade.org and read it completely before using this file.
8 //
9 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
10 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
11 //
12 // The Original Code and all software distributed under the License is
13 // distributed on an "AS IS" basis, without warranty of any kind, and the
14 // Initial Developer hereby disclaims all such warranties, including without
15 // limitation, any warranties of merchantability, fitness for a particular
16 // purpose or non-infringement. Please see the License for the specific terms
17 // and conditions governing the rights and limitations under the License.
18
19
20 #include <float.h>
21 #include <Standard_Real.hxx>
22 #include <Standard_RangeError.hxx>
23 #include <Standard_NumericError.hxx>
24 #include <Standard_NullValue.hxx>
25 #ifndef _Standard_Stream_HeaderFile
26 #include <Standard_Stream.hxx>
27 #endif
28 #ifndef _Standard_OStream_HeaderFile
29 #include <Standard_OStream.hxx>
30 #endif
31
32 const Handle_Standard_Type& Standard_Real_Type_() 
33 {
34   static Handle_Standard_Type _aType = 
35     new Standard_Type("Standard_Real",sizeof(Standard_Real),0,NULL);
36   
37   return _aType;
38 }
39
40 // ------------------------------------------------------------------
41 // Hascode : Computes a hascoding value for a given real
42 // ------------------------------------------------------------------
43 Standard_Integer HashCode(const Standard_Real me, const Standard_Integer Upper)
44 {
45   if (Upper < 1){
46      Standard_RangeError::
47       Raise("Try to apply HashCode method with negative or null argument.");
48   }
49   union 
50     {
51     Standard_Real R;
52     Standard_Integer I[2];
53     } U;
54 //  U.R = Abs(me); // Treat me = -0.0 ADN 27/11/97
55   U.R = me ;
56   return HashCode( ( U.I[0] ^ U.I[1] ) , Upper ) ;
57   }
58
59 // ------------------------------------------------------------------
60 // ShallowCopy : Makes a copy of a real value
61 // ------------------------------------------------------------------
62 Standard_Real ShallowCopy (const Standard_Real me) 
63 {
64   return me;
65 }
66
67 //-------------------------------------------------------------------
68 // ACos : Returns the value of the arc cosine of a real
69 //-------------------------------------------------------------------
70 Standard_Real ACos (const Standard_Real Value) 
71
72   if ( (Value < -1.) || (Value > 1.) ){
73     Standard_RangeError::Raise();
74   } 
75   return acos(Value); 
76 }
77
78 //-------------------------------------------------------------------
79 // ACosApprox : Returns the approximate value of the arc cosine of a real.
80 //              The max error is about 1 degree near Value=0.
81 //-------------------------------------------------------------------
82
83 inline Standard_Real apx_for_ACosApprox (const Standard_Real x)
84 {
85   return  (-0.000007239283986332 +
86     x * (2.000291665285952400 +
87     x * (0.163910606547823220 +
88     x * (0.047654245891495528 -
89     x * (0.005516443930088506 +
90     0.015098965761299077 * x))))) / sqrt(2*x);
91 }
92
93 Standard_Real ACosApprox (const Standard_Real Value)
94 {
95   double XX;
96   if (Value < 0.) {
97     XX = 1.+Value;
98     if (XX < RealSmall())
99       return 0.;
100     return M_PI - apx_for_ACosApprox(XX);
101   }
102   XX = 1.-Value;
103   if (XX < RealSmall())
104     return 0.;
105   return apx_for_ACosApprox(XX);
106
107 // The code above is the same but includes 2 comparisons instead of 3
108 //   Standard_Real xn = 1.+Value;
109 //   Standard_Real xp = 1.-Value;
110 //   if (xp < RealSmall() || xn < RealSmall())
111 //     return 0.;
112 //   if (Value < 0.)
113 //     return M_PI - apx_for_ACosApprox (xn);
114 //   return apx_for_ACosApprox (xp);
115 }
116
117 //-------------------------------------------------------------------
118 // ASin : Returns the value of the arc sine of a real
119 //-------------------------------------------------------------------
120 Standard_Real ASin (const Standard_Real Value) 
121
122   if ( Value < -1 || Value > 1 ){
123     Standard_RangeError::Raise();
124   }
125   return asin(Value); 
126 }
127
128 //-------------------------------------------------------------------
129 // ATan2 : Returns the arc tangent of a real divide by an another real
130 //-------------------------------------------------------------------
131 Standard_Real ATan2 (const Standard_Real Value, const Standard_Real Other) 
132
133   if ( Value == 0. && Other == 0. ){
134     Standard_NullValue::Raise();
135   }
136   return atan2(Value,Other); 
137 }
138
139 //-------------------------------------------------------------------
140 // Sign : Returns |a| if B >= 0; -|a| if b < 0.
141 //             from x in the direction y
142 //-------------------------------------------------------------------
143 Standard_Real Sign(const Standard_Real a, const Standard_Real b)
144 {
145   //==== We use the function "nextafter()" fom library "math.h" ==============
146   if (b >= 0.0) {
147     return Abs(a);
148   } else {
149     return (-1.0 * Abs(a));
150   }
151 }
152
153 //==========================================================================
154 //===== The special routines for "IEEE" and differents hardwares ===========
155 //==========================================================================
156 union RealMap {
157   double real;
158   unsigned int map[2];
159 };
160
161 //--------------------------------------------------------------------
162 // HardwareHighBitsOfDouble :  
163 //    Returns 1 if the low bits are at end.   (exemple: decmips and ALPHA )
164 //    Returns 0 if the low bits are at begin. (exemple: sun, sgi, ...)
165 //--------------------------------------------------------------------
166 static int HardwareHighBitsOfDouble()
167 {
168   RealMap MaxDouble;
169   MaxDouble.real = DBL_MAX;
170   //=========================================================
171   // reperesentation of the max double in IEEE is 
172   //      "7fef ffff ffff ffff"   for the big indiens.
173   //      "ffff ffff 7fef ffff"   for the littel indiens.
174   //=========================================================
175
176   if(MaxDouble.map[1] != 0xffffffff){
177     return 1;
178   } else {
179     return 0;
180   }
181 }
182
183 //--------------------------------------------------------------------
184 // HardwareLowBitsOfDouble :  
185 //    Returns 0 if the low bits are at end.   (exemple: decmips )
186 //    Returns 1 if the low bits are at begin. (exemple: sun, sgi, ...)
187 //--------------------------------------------------------------------
188 static int HardwareLowBitsOfDouble()
189 {
190   RealMap MaxDouble;
191   MaxDouble.real = DBL_MAX;
192   //=========================================================
193   // reperesentation of the max double in IEEE is 
194   //      "7fef ffff ffff ffff"   for the big indiens.
195   //      "ffff ffff 7fef ffff"   for the littel indiens.
196   //=========================================================
197
198   if(MaxDouble.map[1] != 0xffffffff){
199     return 0;
200   } else {
201     return 1;
202   }
203 }
204
205 static int HighBitsOfDouble = HardwareHighBitsOfDouble();
206 static int LowBitsOfDouble = HardwareLowBitsOfDouble();
207
208 double NextAfter(const double x, const double y)
209 {
210   RealMap res;
211
212   res.real=x;
213   
214   if (x == 0.0) {
215         return DBL_MIN;
216   }
217   if(x==y) {
218     //=========================================
219     //   -oo__________0___________+oo
220     //               x=y
221     //  The direction is "Null", so there is nothing after 
222     //=========================================
223
224   } else if (((x<y) && (x>=0.0)) || ((x>y) && (x<0.0))) {
225     //=========================================
226     //   -oo__________0___________+oo
227     //        y <- x     x -> y
228     //
229     //=========================================
230     if (res.map[LowBitsOfDouble]==0xffffffff) {
231       res.map[LowBitsOfDouble]=0;
232       res.map[HighBitsOfDouble]++;
233     } else {
234       res.map[LowBitsOfDouble]++;
235     }
236   } else {
237     //=========================================
238     //   -oo__________0___________+oo
239     //        x -> y     y <- x
240     //
241     //=========================================
242     if (res.map[LowBitsOfDouble]==0) {
243       if (res.map[HighBitsOfDouble]==0) {
244         res.map[HighBitsOfDouble]=0x80000000;
245         res.map[LowBitsOfDouble]=0x00000001;
246       } else {
247         res.map[LowBitsOfDouble]=0xffffffff;
248         res.map[HighBitsOfDouble]--;
249       }
250     } else {
251       res.map[LowBitsOfDouble]--;
252     }
253   }
254   return res.real;
255 }
256
257 // ------------------------------------------------------------------
258 // ShallowDump : Writes a real  value
259 // ------------------------------------------------------------------
260 Standard_EXPORT void              ShallowDump(const Standard_Real Value, 
261                                               Standard_OStream& s)
262 { s << Value << " Standard_Real" << "\n"; }
263
264
265 //-------------------------------------------------------------------
266 // ATanh : Returns the value of the hyperbolic arc tangent of a real
267 //-------------------------------------------------------------------
268 Standard_Real     ATanh(const Standard_Real Value) 
269
270   if ( (Value <= -1.) || (Value >= 1.) ){
271     Standard_NumericError::Raise("Illegal agument in ATanh");
272     cout << "Illegal agument in ATanh" << endl ;
273   } 
274   return atanh(Value); 
275 }
276
277 //-------------------------------------------------------------------
278 // ACosh : Returns the hyperbolic Arc cosine of a real
279 //-------------------------------------------------------------------
280 Standard_Real     ACosh (const Standard_Real Value) 
281
282   if ( Value < 1. ){
283     Standard_NumericError::Raise("Illegal agument in ACosh");
284     cout << "Illegal agument in ACosh" << endl ;
285   } 
286   return acosh(Value); 
287 }
288
289 //-------------------------------------------------------------------
290 // Log : Returns the naturaOPl logarithm of a real
291 //-------------------------------------------------------------------
292 Standard_Real     Log (const Standard_Real Value) 
293 {   if ( Value <= 0. ){
294     Standard_NumericError::Raise("Illegal agument in Log");
295     cout << "Illegal agument in Log" << endl ;
296   } 
297  return log(Value); 
298 }
299 //-------------------------------------------------------------------
300 // Sqrt : Returns the square root of a real
301 //-------------------------------------------------------------------
302 Standard_Real     Sqrt (const Standard_Real Value) 
303
304   if (  Value < 0. ){
305     Standard_NumericError::Raise("Illegal agument in Sqrt");
306     cout << "Illegal agument in Sqrt" << endl ;
307   } 
308  return sqrt(Value); 
309 }
310