mono/metadata/sysmath.c

   1 /*
   2  * sysmath.c: these are based on bob smith's csharp routines
   3  *
   4  * Author:
   5  *      Mono Project (http://www.mono-project.com)
   6  *      Ludovic Henry (ludovic@xamarin.com)
   7  *
   8  * Copyright 2001-2003 Ximian, Inc (http://www.ximian.com)
   9  * Copyright 2004-2009 Novell, Inc (http://www.novell.com)
  10  * Copyright 2015 Xamarin, Inc (https://www.xamarin.com)
  11  */
  12
  13 //
  14 // Copyright (c) Microsoft. All rights reserved.
  15 // Licensed under the MIT license. See LICENSE file in the project root for full license information.
  16 //
  17 // Files:
  18 //  - src/classlibnative/float/floatnative.cpp
  19 //  - src/pal/src/cruntime/floatnative.cpp
  20 //
  21 // Ported from C++ to C and adjusted to Mono runtime
  22
  23 #define __USE_ISOC99
  24
  25 #include <math.h>
  26 #include <mono/metadata/sysmath.h>
  27
  28 #include "number-ms.h"
  29 #include "utils/mono-compiler.h"
  30
  31 static const MonoDouble_double NaN = { .s = { .sign = 0x0, .exp = 0x7FF, .mantHi = 0x80000, .mantLo = 0x0 } };
  32
  33 /* +Infinity */
  34 static const MonoDouble_double PInfinity = { .s = { .sign = 0x0, .exp = 0x7FF, .mantHi = 0x0, .mantLo = 0x0 } };
  35
  36 /* -Infinity */
  37 static const MonoDouble_double MInfinity = { .s = { .sign = 0x1, .exp = 0x7FF, .mantHi = 0x0, .mantLo = 0x0 } };
  38
  39 /* +1 */
  40 static const MonoDouble_double POne = { .s = { .sign = 0x0, .exp = 0x3FF, .mantHi = 0x0, .mantLo = 0x0 } };
  41
  42 /* -1 */
  43 static const MonoDouble_double MOne = { .s = { .sign = 0x1, .exp = 0x3FF, .mantHi = 0x0, .mantLo = 0x0 } };
  44
  45 static MONO_ALWAYS_INLINE gboolean
  46 isplusinfinity (gdouble d)
  47 {
  48         return d == PInfinity.d;
  49 }
  50
  51 static MONO_ALWAYS_INLINE gboolean
  52 isminusinfinity (gdouble d)
  53 {
  54         return d == MInfinity.d;
  55 }
  56
  57 static MONO_ALWAYS_INLINE gboolean
  58 isinfinity (gdouble d)
  59 {
  60         return isplusinfinity (d) || isminusinfinity (d);
  61 }
  62
  63 static MONO_ALWAYS_INLINE gboolean
  64 isplusone (gdouble d)
  65 {
  66         return d == POne.d;
  67 }
  68
  69 static MONO_ALWAYS_INLINE gboolean
  70 isminusone (gdouble d)
  71 {
  72         return d == MOne.d;
  73 }
  74
  75 gdouble
  76 ves_icall_System_Math_Floor (gdouble x)
  77 {
  78         return floor(x);
  79 }
  80
  81 gdouble
  82 ves_icall_System_Math_Round (gdouble x)
  83 {
  84         gdouble tmp, floor_tmp;
  85
  86         /* If the number has no fractional part do nothing This shortcut is necessary
  87          * to workaround precision loss in borderline cases on some platforms */
  88         if (x == (gdouble)(gint64) x)
  89                 return x;
  90
  91         tmp = x + 0.5;
  92         floor_tmp = floor (tmp);
  93
  94         if (floor_tmp == tmp) {
  95                 if (fmod (tmp, 2.0) != 0)
  96                         floor_tmp -= 1.0;
  97         }
  98
  99         return copysign (floor_tmp, x);
 100 }
 101
 102 gdouble
 103 ves_icall_System_Math_Sin (gdouble x)
 104 {
 105         return sin (x);
 106 }
 107
 108 gdouble
 109 ves_icall_System_Math_Cos (gdouble x)
 110 {
 111         return cos (x);
 112 }
 113
 114 gdouble
 115 ves_icall_System_Math_Tan (gdouble x)
 116 {
 117         return tan (x);
 118 }
 119
 120 gdouble
 121 ves_icall_System_Math_Sinh (gdouble x)
 122 {
 123         return sinh (x);
 124 }
 125
 126 gdouble
 127 ves_icall_System_Math_Cosh (gdouble x)
 128 {
 129         return cosh (x);
 130 }
 131
 132 gdouble
 133 ves_icall_System_Math_Tanh (gdouble x)
 134 {
 135         return tanh (x);
 136 }
 137
 138 gdouble
 139 ves_icall_System_Math_Acos (gdouble x)
 140 {
 141         if (x < -1 || x > 1)
 142                 return NaN.d;
 143
 144         return acos (x);
 145 }
 146
 147 gdouble
 148 ves_icall_System_Math_Asin (gdouble x)
 149 {
 150         if (x < -1 || x > 1)
 151                 return NaN.d;
 152
 153         return asin (x);
 154 }
 155
 156 gdouble
 157 ves_icall_System_Math_Atan (gdouble x)
 158 {
 159         return atan (x);
 160 }
 161
 162 gdouble
 163 ves_icall_System_Math_Atan2 (gdouble y, gdouble x)
 164 {
 165         gdouble result;
 166
 167         if (isinfinity (x) && isinfinity (y))
 168                 return NaN.d;
 169
 170         result = atan2 (y, x);
 171         return result == -0.0 ? 0.0: result;
 172 }
 173
 174 gdouble
 175 ves_icall_System_Math_Exp (gdouble x)
 176 {
 177         if (isinfinity (x))
 178                 return x < 0 ? 0.0 : x;
 179
 180         return exp (x);
 181 }
 182
 183 gdouble
 184 ves_icall_System_Math_Log (gdouble x)
 185 {
 186         if (x == 0)
 187                 return MInfinity.d;
 188         else if (x < 0)
 189                 return NaN.d;
 190
 191         return log (x);
 192 }
 193
 194 gdouble
 195 ves_icall_System_Math_Log10 (gdouble x)
 196 {
 197         if (x == 0)
 198                 return MInfinity.d;
 199         else if (x < 0)
 200                 return NaN.d;
 201
 202         return log10 (x);
 203 }
 204
 205 gdouble
 206 ves_icall_System_Math_Pow (gdouble x, gdouble y)
 207 {
 208         gdouble result;
 209
 210         if (isnan (y))
 211                 return y;
 212         if (isnan (x))
 213                 return x;
 214
 215         if (isinfinity (y)) {
 216                 if (isplusone (x))
 217                         return x;
 218                 if (isminusone (x))
 219                         return NaN.d;
 220         }
 221
 222         /* following are cases from PAL_pow which abstract the implementation of pow for posix and win32 platforms
 223          * (https://github.com/dotnet/coreclr/blob/master/src/pal/src/cruntime/finite.cpp#L331) */
 224
 225         if (isplusinfinity (y) && !isnan (x)) {
 226                 if (isplusone (x) || isminusone (x))
 227                         result = NaN.d;
 228                 else if (x > MOne.d && x < POne.d)
 229                         result = 0.0;
 230                 else
 231                         result = PInfinity.d;
 232         } else if (isminusinfinity (y) && !isnan (x)) {
 233                 if (isplusone (x) || isminusone (x))
 234                         result = NaN.d;
 235                 if (x > MOne.d && x < POne.d)
 236                         result = PInfinity.d;
 237                 else
 238                         result = 0.0;
 239         } else if (x == 0.0 && y < 0.0) {
 240                 result = PInfinity.d;
 241         } else if (y == 0.0 && isnan (x)) {
 242                 /* Windows returns NaN for pow(NaN, 0), but POSIX specifies
 243                  * a return value of 1 for that case.  We need to return
 244                  * the same result as Windows. */
 245                 result = NaN.d;
 246         } else {
 247                 result = pow (x, y);
 248         }
 249
 250         if (result == PInfinity.d && x < 0.0 && isfinite (x) && ceil (y / 2) != floor (y / 2))
 251                 result = MInfinity.d;
 252
 253         /*
 254          * The even/odd test in the if (this one and the one above) used to be ((long long) y % 2 == 0)
 255          * on SPARC (long long) y for large y (>2**63) is always 0x7fffffff7fffffff, which
 256          * is an odd number, so the test ((long long) y % 2 == 0) will always fail for
 257          * large y. Since large double numbers are always even (e.g., the representation of
 258          * 1E20+1 is the same as that of 1E20, the last .+1. is too insignificant to be part
 259          * of the representation), this test will always return the wrong result for large y.
 260          *
 261          * The (ceil(y/2) == floor(y/2)) test is slower, but more robust.
 262          */
 263         if (result == MInfinity.d && x < 0.0 && isfinite (x) && ceil (y / 2) == floor (y / 2))
 264                 result = PInfinity.d;
 265
 266 #if defined (__linux__) && SIZEOF_VOID_P == 4
 267         /* On Linux 32bits, some tests erroneously return NaN */
 268         if (isnan (result)) {
 269                 if (isminusone (x) && (y > 9007199254740991.0 || y < -9007199254740991.0)) {
 270                         /* Math.Pow (-1, Double.MaxValue) and Math.Pow (-1, Double.MinValue) should return 1 */
 271                         result = POne.d;
 272                 } else if (x < -9007199254740991.0 && y < -9007199254740991.0) {
 273                         /* Math.Pow (Double.MinValue, Double.MinValue) should return 0 */
 274                         result = 0.0;
 275                 } else if (x < -9007199254740991.0 && y > 9007199254740991.0) {
 276                         /* Math.Pow (Double.MinValue, Double.MaxValue) should return Double.PositiveInfinity */
 277                         result = PInfinity.d;
 278                 }
 279         }
 280 #endif
 281
 282         return result == -0.0 ? 0 : result;
 283 }
 284
 285 gdouble
 286 ves_icall_System_Math_Sqrt (gdouble x)
 287 {
 288         if (x < 0)
 289                 return NaN.d;
 290
 291         return sqrt (x);
 292 }
 293
 294 gdouble
 295 ves_icall_System_Math_Abs_double (gdouble v)
 296 {
 297         return fabs (v);
 298 }
 299
 300 gfloat
 301 ves_icall_System_Math_Abs_single (gfloat v)
 302 {
 303         return fabsf (v);
 304 }
 305
 306 gdouble
 307 ves_icall_System_Math_Ceiling (gdouble v)
 308 {
 309         return ceil (v);
 310 }
 311
 312 gdouble
 313 ves_icall_System_Math_SplitFractionDouble (gdouble *v)
 314 {
 315         return modf (*v, v);
 316 }