platform/win32/MathUtilsWin.cpp

   1 /* -*- Mode: C++; c-basic-offset: 4; indent-tabs-mode: nil; tab-width: 4 -*- */
   2 /* vi: set ts=4 sw=4 expandtab: (add to ~/.vimrc: set modeline modelines=5) */
   3 /* ***** BEGIN LICENSE BLOCK *****
   4  * Version: MPL 1.1/GPL 2.0/LGPL 2.1
   5  *
   6  * The contents of this file are subject to the Mozilla Public License Version
   7  * 1.1 (the "License"); you may not use this file except in compliance with
   8  * the License. You may obtain a copy of the License at
   9  * http://www.mozilla.org/MPL/
  10  *
  11  * Software distributed under the License is distributed on an "AS IS" basis,
  12  * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
  13  * for the specific language governing rights and limitations under the
  14  * License.
  15  *
  16  * The Original Code is [Open Source Virtual Machine.].
  17  *
  18  * The Initial Developer of the Original Code is
  19  * Adobe System Incorporated.
  20  * Portions created by the Initial Developer are Copyright (C) 2004-2006
  21  * the Initial Developer. All Rights Reserved.
  22  *
  23  * Contributor(s):
  24  *   Adobe AS3 Team
  25  *
  26  * Alternatively, the contents of this file may be used under the terms of
  27  * either the GNU General Public License Version 2 or later (the "GPL"), or
  28  * the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
  29  * in which case the provisions of the GPL or the LGPL are applicable instead
  30  * of those above. If you wish to allow use of your version of this file only
  31  * under the terms of either the GPL or the LGPL, and not to allow others to
  32  * use your version of this file under the terms of the MPL, indicate your
  33  * decision by deleting the provisions above and replace them with the notice
  34  * and other provisions required by the GPL or the LGPL. If you do not delete
  35  * the provisions above, a recipient may use your version of this file under
  36  * the terms of any one of the MPL, the GPL or the LGPL.
  37  *
  38  * ***** END LICENSE BLOCK ***** */
  39
  40 #include "avmplus.h"
  41 #include <math.h>
  42
  43 #ifdef AVMPLUS_IA32
  44 #define X86_MATH
  45 #endif
  46
  47 // warning this code is used by amd64 and arm builds
  48
  49 namespace avmplus
  50 {
  51 #ifdef AVMPLUS_ARM
  52
  53     const static double PI = 3.141592653589793;
  54     const static double PI3_BY_4 = 3*PI/4;
  55     const static double PI_BY_4 = PI/4;
  56     const static double PI2 = 2*PI;
  57
  58     // 0=no, 1=+0, -1=-0
  59     static int32_t isZero(double v)
  60     {
  61         int32_t r = (MathUtils::isNegZero(v)) ? -1 : (v==0.0)? 1 : 0;
  62         return r;
  63     }
  64
  65     // sin, cos, tan all function incorrectly when called with really large values on windows mobile
  66     // they all start failing at different values, but all start failing somewhere with values
  67     // greater than 210 million.
  68     #define AVMPLUS_TRIG_FUNC_MAX 210000000
  69
  70     const static bool broken_trig_funcs = MathUtils::isNaN(MathUtils::cos(250000000));
  71
  72     // Helper function to adjust a value for sin, cos, or tan into an equivalent value
  73     // in the range that works correctly.  This works because these functions all have a period
  74     // of 2*PI or PI, so there are many equivalent values.
  75     static double adjustValueForTrigFuncs(double v)
  76     {
  77         bool negate = false;
  78         if( v < 0 ) {
  79             v = -v;
  80             negate = true;
  81         }
  82
  83         int temp = (int)((v - AVMPLUS_TRIG_FUNC_MAX)/PI2);
  84         double offset = PI2*(temp+1);
  85         v -= offset;
  86
  87         if( negate )
  88             v = -v;
  89
  90         return v;
  91     }
  92
  93 #endif /* AVMPLUS_ARM */
  94
  95 #ifdef AVMPLUS_ARM
  96     double MathUtils::atan2(double y, double x)
  97     {
  98         int32_t zx = isZero(x);
  99         int32_t zy = isZero(y);
 100         if (zx==-1 && zy!=0)
 101             return zy*PI;  // +-0,-0 case
 102         else if (zy==-1 && (x==1.0 || x==-1.0))
 103             return -(::atan2(y,x));  // negate result
 104
 105         double r = ::atan2(y, x);
 106         if (MathUtils::isNaN(r)) {
 107             int32_t s = MathUtils::isInfinite(x);
 108             if (s==1)
 109                 r = MathUtils::isInfinite(y) * PI_BY_4;
 110             else if (s==-1)
 111                 r = MathUtils::isInfinite(y) * PI3_BY_4;
 112         }
 113         return r;
 114     }
 115 #endif /* AVMPLUS_ARM */
 116
 117 #ifdef X86_MATH
 118     double MathUtils::ceil(double value)
 119     {
 120         // todo avoid control word modification
 121         short oldcw, newcw;
 122         _asm fnstcw [oldcw];
 123         _asm mov ax, [oldcw];
 124         _asm and ax, 0xf3ff; // Set to round down.
 125         _asm or  ax, 0x800;
 126         _asm mov [newcw], ax;
 127         _asm fldcw [newcw];
 128         _asm fld [value];
 129         _asm frndint;
 130         _asm fldcw [oldcw];
 131     }
 132 #endif /* X86_MATH */
 133
 134 #ifdef AVMPLUS_ARM
 135     double MathUtils::cos(double value)
 136     {
 137         if( broken_trig_funcs && (value > AVMPLUS_TRIG_FUNC_MAX || value < -AVMPLUS_TRIG_FUNC_MAX) )
 138         {
 139             return ::cos(adjustValueForTrigFuncs(value));
 140         }
 141         else
 142         {
 143             return ::cos(value);
 144         }
 145     }
 146 #endif AVMPLUS_ARM
 147
 148 #ifdef X86_MATH
 149     // Utility function, this module only.
 150     REALLY_INLINE static double expInternal(double x)
 151     {
 152         double value, exponent;
 153         _asm fld [x];
 154         _asm fldl2e;
 155         _asm _emit 0xD8; // fmul st(1);
 156         _asm _emit 0xC9;
 157         _asm _emit 0xDD; // fst st(1);
 158         _asm _emit 0xD1;
 159         _asm frndint;
 160         _asm fxch;
 161         _asm _emit 0xD8; // fsub st1;
 162         _asm _emit 0xE1;
 163         _asm f2xm1;
 164         _asm fstp [value];
 165         _asm fstp [exponent];
 166
 167         value += 1.0;
 168
 169         _asm fld [exponent];
 170         _asm fld [value];
 171         _asm fscale;
 172
 173         _asm fxch;
 174         _asm _emit 0xDD; // fstp st(0);
 175         _asm _emit 0xD8;
 176     }
 177
 178     // Inlined on other architectures
 179     double MathUtils::exp(double value)
 180     {
 181         switch (isInfinite(value)) {
 182         case 1:
 183             return kInfinity;
 184         case -1:
 185             return +0;
 186         default:
 187             return expInternal(value);
 188         }
 189     }
 190 #endif /* X86_MATH */
 191
 192 #ifdef X86_MATH
 193     double MathUtils::floor(double value)
 194     {
 195         // todo avoid control word modification
 196         short oldcw, newcw;
 197         _asm fnstcw [oldcw];
 198         _asm mov ax, [oldcw];
 199         _asm and ax, 0xf3ff; // Set to round down.
 200         _asm or  ax, 0x400;
 201         _asm mov [newcw], ax;
 202         _asm fldcw [newcw];
 203         _asm fld [value];
 204         _asm frndint;
 205         _asm fldcw [oldcw];
 206     }
 207 #endif /* X86_MATH */
 208
 209     /* @(#)s_frexp.c 5.1 93/09/24 */
 210     /*
 211      * ====================================================
 212      * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 213      *
 214      * Developed at SunPro, a Sun Microsystems, Inc. business.
 215      * Permission to use, copy, modify, and distribute this
 216      * software is freely granted, provided that this notice
 217      * is preserved.
 218      * ====================================================
 219      */
 220
 221     /*
 222      * for non-zero x
 223      *  x = frexp(arg,&exp);
 224      * return a double fp quantity x such that 0.5 <= |x| <1.0
 225      * and the corresponding binary exponent "exp". That is
 226      *  arg = x*2^exp.
 227      * If arg is inf, 0.0, or NaN, then frexp(arg,&exp) returns arg
 228      * with *exp=0.
 229      */
 230
 231     /*
 232      * NOTE: This is little-endian, must be adjusted to work for
 233      * big-endian systems.
 234      */
 235 #define EXTRACT_WORDS(hx, lx, x) {DWORD *ptr = (DWORD*)&x; hx=ptr[1]; lx=ptr[0];}
 236 #define SET_HIGH_WORD(x, hx) {DWORD *ptr = (DWORD*)&x; ptr[1]=hx;}
 237 #define GET_HIGH_WORD(hx, x) {DWORD *ptr = (DWORD*)&x; hx=ptr[1];}
 238     static const double two54 =  1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */
 239
 240     REALLY_INLINE static double ExtractFraction(double x, int *eptr)
 241     {
 242         DWORD hx, ix, lx;
 243         EXTRACT_WORDS(hx,lx,x);
 244         ix = 0x7fffffff&hx;
 245         *eptr = 0;
 246         if(ix>=0x7ff00000||((ix|lx)==0)) return x;  /* 0,inf,nan */
 247         if (ix<0x00100000) {        /* subnormal */
 248             x *= two54;
 249             GET_HIGH_WORD(hx,x);
 250             ix = hx&0x7fffffff;
 251             *eptr = -54;
 252         }
 253         *eptr += (ix>>20)-1022;
 254         hx = (hx&0x800fffff)|0x3fe00000;
 255         SET_HIGH_WORD(x,hx);
 256         return x;
 257     }
 258
 259     uint64_t MathUtils::frexp(double x, int *eptr)
 260     {
 261         double fracMantissa = ExtractFraction(x, eptr);
 262         // correct mantissa and eptr to get integer values
 263         //  for both
 264         *eptr -= 53; // 52 mantissa bits + the hidden bit
 265         return (uint64_t)((fracMantissa) * (double)(1LL << 53));
 266     }
 267
 268 #ifdef X86_MATH
 269     // VC++ 2008 refuses to inline this, issues warning on _forceinline
 270     #pragma warning ( disable : 4740 ) // flow in or out of inline asm code suppresses global optimization
 271     double MathUtils::mod(double x, double y)
 272     {
 273         if (!y) {
 274             return kNaN;
 275         }
 276         _asm    fld [y];
 277         _asm    fld [x];
 278       ModLoop:
 279         _asm    fprem;
 280         _asm    fnstsw ax;
 281         _asm    sahf;
 282         _asm    jp ModLoop;
 283         _asm _emit 0xDD; // fstp st(1);
 284         _asm _emit 0xD9;
 285     }
 286     #pragma warning ( default : 4740 )
 287 #endif /* X86_MATH */
 288
 289 #ifdef X86_MATH
 290     // Std. library pow()
 291     double MathUtils::powInternal(double x, double y)
 292     {
 293         double value, exponent;
 294
 295         _asm fld1;
 296         _asm fld [x];
 297         _asm fyl2x;
 298         _asm fstp [value];
 299
 300         _asm fld [value];
 301         _asm fld [y];
 302         _asm _emit 0xD8; // fmul st(1);
 303         _asm _emit 0xC9;
 304         _asm _emit 0xDD; // fst st(1);
 305         _asm _emit 0xD1;
 306         _asm frndint;
 307         _asm fxch;
 308         _asm _emit 0xD8; // fsub st1;
 309         _asm _emit 0xE1;
 310         _asm f2xm1;
 311         _asm fstp [value];
 312         _asm fstp [exponent];
 313
 314         value += 1.0;
 315
 316         _asm fld [exponent];
 317         _asm fld [value];
 318         _asm fscale;
 319
 320         _asm fxch;
 321         _asm _emit 0xDD; // fstp st(0);
 322         _asm _emit 0xD8;
 323     }
 324 #endif /* X86_MATH */
 325
 326 #ifdef AVMPLUS_ARM
 327     double MathUtils::sin(double value)
 328     {
 329         if( broken_trig_funcs && (value > AVMPLUS_TRIG_FUNC_MAX || value < -AVMPLUS_TRIG_FUNC_MAX) )
 330         {
 331             return ::sin(adjustValueForTrigFuncs(value));
 332         }
 333         else
 334         {
 335             return ::sin(value);
 336         }
 337     }
 338 #endif /* AVMPLUS_ARM */
 339
 340 #ifdef X86_MATH
 341     double MathUtils::tan(double value)
 342     {
 343         // This is a good candidate for inlining, but VC++ 2008 chokes on it.
 344         _asm fld [value];
 345         _asm fptan;
 346         _asm _emit 0xDD; // fstp st(0);
 347         _asm _emit 0xD8;
 348     }
 349 #elif defined(AVMPLUS_ARM)
 350     double MathUtils::tan(double value)
 351     {
 352         if( broken_trig_funcs && (value > AVMPLUS_TRIG_FUNC_MAX || value < -AVMPLUS_TRIG_FUNC_MAX) )
 353         {
 354             return ::tan(adjustValueForTrigFuncs(value));
 355         }
 356         else
 357         {
 358             return ::tan(value);
 359         }
 360     }
 361 #endif /* X86_MATH */
 362 }