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 #include <Standard_OStream.hxx>
  22
  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      throw Standard_RangeError("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 // ACos : Returns the value of the arc cosine of a real
  43 //-------------------------------------------------------------------
  44 Standard_Real ACos (const Standard_Real Value)
  45 {
  46   if ( (Value < -1.) || (Value > 1.) ){
  47     throw Standard_RangeError();
  48   }
  49   return acos(Value);
  50 }
  51
  52 //-------------------------------------------------------------------
  53 // ACosApprox : Returns the approximate value of the arc cosine of a real.
  54 //              The max error is about 1 degree near Value=0.
  55 //-------------------------------------------------------------------
  56
  57 inline Standard_Real apx_for_ACosApprox (const Standard_Real x)
  58 {
  59   return  (-0.000007239283986332 +
  60     x * (2.000291665285952400 +
  61     x * (0.163910606547823220 +
  62     x * (0.047654245891495528 -
  63     x * (0.005516443930088506 +
  64     0.015098965761299077 * x))))) / sqrt(2*x);
  65 }
  66
  67 Standard_Real ACosApprox (const Standard_Real Value)
  68 {
  69   double XX;
  70   if (Value < 0.) {
  71     XX = 1.+Value;
  72     if (XX < RealSmall())
  73       return 0.;
  74     return M_PI - apx_for_ACosApprox(XX);
  75   }
  76   XX = 1.-Value;
  77   if (XX < RealSmall())
  78     return 0.;
  79   return apx_for_ACosApprox(XX);
  80
  81 // The code above is the same but includes 2 comparisons instead of 3
  82 //   Standard_Real xn = 1.+Value;
  83 //   Standard_Real xp = 1.-Value;
  84 //   if (xp < RealSmall() || xn < RealSmall())
  85 //     return 0.;
  86 //   if (Value < 0.)
  87 //     return M_PI - apx_for_ACosApprox (xn);
  88 //   return apx_for_ACosApprox (xp);
  89 }
  90
  91 //-------------------------------------------------------------------
  92 // ASin : Returns the value of the arc sine of a real
  93 //-------------------------------------------------------------------
  94 Standard_Real ASin (const Standard_Real Value)
  95 {
  96   if ( Value < -1 || Value > 1 ){
  97     throw Standard_RangeError();
  98   }
  99   return asin(Value);
 100 }
 101
 102 //-------------------------------------------------------------------
 103 // ATan2 : Returns the arc tangent of a real divide by an another real
 104 //-------------------------------------------------------------------
 105 Standard_Real ATan2 (const Standard_Real Value, const Standard_Real Other)
 106 {
 107   if ( Value == 0. && Other == 0. ){
 108     throw Standard_NullValue();
 109   }
 110   return atan2(Value,Other);
 111 }
 112
 113 //-------------------------------------------------------------------
 114 // Sign : Returns |a| if B >= 0; -|a| if b < 0.
 115 //-------------------------------------------------------------------
 116 Standard_Real Sign(const Standard_Real a, const Standard_Real b)
 117 {
 118   if (b >= 0.0) {
 119     return Abs(a);
 120   } else {
 121     return (-1.0 * Abs(a));
 122   }
 123 }
 124
 125 //==========================================================================
 126 //===== The special routines for "IEEE" and differents hardwares ===========
 127 //==========================================================================
 128 union RealMap {
 129   double real;
 130   unsigned int map[2];
 131 };
 132
 133 //--------------------------------------------------------------------
 134 // HardwareHighBitsOfDouble :
 135 //    Returns 1 if the low bits are at end.   (exemple: decmips and ALPHA )
 136 //    Returns 0 if the low bits are at begin. (exemple: sun, sgi, ...)
 137 //--------------------------------------------------------------------
 138 static int HardwareHighBitsOfDouble()
 139 {
 140   RealMap MaxDouble;
 141   MaxDouble.real = DBL_MAX;
 142   //=========================================================
 143   // reperesentation of the max double in IEEE is
 144   //      "7fef ffff ffff ffff"   for the big indiens.
 145   //      "ffff ffff 7fef ffff"   for the littel indiens.
 146   //=========================================================
 147
 148   if(MaxDouble.map[1] != 0xffffffff){
 149     return 1;
 150   } else {
 151     return 0;
 152   }
 153 }
 154
 155 //--------------------------------------------------------------------
 156 // HardwareLowBitsOfDouble :
 157 //    Returns 0 if the low bits are at end.   (exemple: decmips )
 158 //    Returns 1 if the low bits are at begin. (exemple: sun, sgi, ...)
 159 //--------------------------------------------------------------------
 160 static int HardwareLowBitsOfDouble()
 161 {
 162   RealMap MaxDouble;
 163   MaxDouble.real = DBL_MAX;
 164   //=========================================================
 165   // reperesentation of the max double in IEEE is
 166   //      "7fef ffff ffff ffff"   for the big indiens.
 167   //      "ffff ffff 7fef ffff"   for the littel indiens.
 168   //=========================================================
 169
 170   if(MaxDouble.map[1] != 0xffffffff){
 171     return 0;
 172   } else {
 173     return 1;
 174   }
 175 }
 176
 177 static int HighBitsOfDouble = HardwareHighBitsOfDouble();
 178 static int LowBitsOfDouble = HardwareLowBitsOfDouble();
 179
 180 double NextAfter(const double x, const double y)
 181 {
 182   RealMap res;
 183
 184   res.real=x;
 185
 186   if (x == 0.0) {
 187         return DBL_MIN;
 188   }
 189   if(x==y) {
 190     //=========================================
 191     //   -oo__________0___________+oo
 192     //               x=y
 193     //  The direction is "Null", so there is nothing after
 194     //=========================================
 195
 196   } else if (((x<y) && (x>=0.0)) || ((x>y) && (x<0.0))) {
 197     //=========================================
 198     //   -oo__________0___________+oo
 199     //        y <- x     x -> y
 200     //
 201     //=========================================
 202     if (res.map[LowBitsOfDouble]==0xffffffff) {
 203       res.map[LowBitsOfDouble]=0;
 204       res.map[HighBitsOfDouble]++;
 205     } else {
 206       res.map[LowBitsOfDouble]++;
 207     }
 208   } else {
 209     //=========================================
 210     //   -oo__________0___________+oo
 211     //        x -> y     y <- x
 212     //
 213     //=========================================
 214     if (res.map[LowBitsOfDouble]==0) {
 215       if (res.map[HighBitsOfDouble]==0) {
 216         res.map[HighBitsOfDouble]=0x80000000;
 217         res.map[LowBitsOfDouble]=0x00000001;
 218       } else {
 219         res.map[LowBitsOfDouble]=0xffffffff;
 220         res.map[HighBitsOfDouble]--;
 221       }
 222     } else {
 223       res.map[LowBitsOfDouble]--;
 224     }
 225   }
 226   return res.real;
 227 }
 228
 229 //-------------------------------------------------------------------
 230 // ATanh : Returns the value of the hyperbolic arc tangent of a real
 231 //-------------------------------------------------------------------
 232 Standard_Real     ATanh(const Standard_Real Value)
 233 {
 234   if ( (Value <= -1.) || (Value >= 1.) ){
 235 #ifdef OCCT_DEBUG
 236     cout << "Illegal agument in ATanh" << endl ;
 237 #endif
 238     throw Standard_NumericError("Illegal agument in ATanh");
 239   }
 240 #if __QNX__
 241   return std::atanh(Value);
 242 #else
 243   return atanh(Value);
 244 #endif
 245 }
 246
 247 //-------------------------------------------------------------------
 248 // ACosh : Returns the hyperbolic Arc cosine of a real
 249 //-------------------------------------------------------------------
 250 Standard_Real     ACosh (const Standard_Real Value)
 251 {
 252   if ( Value < 1. ){
 253 #ifdef OCCT_DEBUG
 254     cout << "Illegal agument in ACosh" << endl ;
 255 #endif
 256     throw Standard_NumericError("Illegal agument in ACosh");
 257   }
 258 #if __QNX__
 259   return std::acosh(Value);
 260 #else
 261   return acosh(Value);
 262 #endif
 263 }
 264
 265 //-------------------------------------------------------------------
 266 // Cosh : Returns the hyperbolic cosine of a real
 267 //-------------------------------------------------------------------
 268 Standard_Real     Cosh (const Standard_Real Value)
 269 {
 270   if ( Abs(Value) > 0.71047586007394394e+03 ){
 271 #ifdef OCCT_DEBUG
 272     cout << "Result of Cosh exceeds the maximum value Standard_Real" << endl ;
 273 #endif
 274     throw Standard_NumericError("Result of Cosh exceeds the maximum value Standard_Real");
 275   }
 276   return cosh(Value);
 277 }
 278
 279 //-------------------------------------------------------------------
 280 // Sinh : Returns the hyperbolicsine of a real
 281 //-------------------------------------------------------------------
 282 Standard_Real     Sinh (const Standard_Real Value)
 283 {
 284   if ( Abs(Value) > 0.71047586007394394e+03 ){
 285 #ifdef OCCT_DEBUG
 286     cout << "Result of Sinh exceeds the maximum value Standard_Real" << endl ;
 287 #endif
 288     throw Standard_NumericError("Result of Sinh exceeds the maximum value Standard_Real");
 289   }
 290   return sinh(Value);
 291 }
 292
 293 //-------------------------------------------------------------------
 294 // Log : Returns the naturaOPl logarithm of a real
 295 //-------------------------------------------------------------------
 296 Standard_Real     Log (const Standard_Real Value)
 297 {   if ( Value <= 0. ){
 298 #ifdef OCCT_DEBUG
 299     cout << "Illegal agument in Log" << endl ;
 300 #endif
 301     throw Standard_NumericError("Illegal agument in Log");
 302   }
 303  return log(Value);
 304 }
 305 //-------------------------------------------------------------------
 306 // Sqrt : Returns the square root of a real
 307 //-------------------------------------------------------------------
 308 Standard_Real     Sqrt (const Standard_Real Value)
 309 {
 310   if (  Value < 0. ){
 311 #ifdef OCCT_DEBUG
 312     cout << "Illegal agument in Sqrt" << endl ;
 313 #endif
 314     throw Standard_NumericError("Illegal agument in Sqrt");
 315   }
 316  return sqrt(Value);
 317 }
 318