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