src/lib_math.c

   1 /*
   2 ** Math library.
   3 ** Copyright (C) 2005-2016 Mike Pall. See Copyright Notice in luajit.h
   4 */
   5
   6 #include <math.h>
   7
   8 #define lib_math_c
   9 #define LUA_LIB
  10
  11 #include "lua.h"
  12 #include "lauxlib.h"
  13 #include "lualib.h"
  14
  15 #include "lj_obj.h"
  16 #include "lj_lib.h"
  17 #include "lj_vm.h"
  18
  19 /* ------------------------------------------------------------------------ */
  20
  21 #define LJLIB_MODULE_math
  22
  23 LJLIB_ASM(math_abs)             LJLIB_REC(.)
  24 {
  25   lj_lib_checknumber(L, 1);
  26   return FFH_RETRY;
  27 }
  28 LJLIB_ASM_(math_floor)          LJLIB_REC(math_round IRFPM_FLOOR)
  29 LJLIB_ASM_(math_ceil)           LJLIB_REC(math_round IRFPM_CEIL)
  30
  31 LJLIB_ASM(math_sqrt)            LJLIB_REC(math_unary IRFPM_SQRT)
  32 {
  33   lj_lib_checknum(L, 1);
  34   return FFH_RETRY;
  35 }
  36 LJLIB_ASM_(math_log10)          LJLIB_REC(math_unary IRFPM_LOG10)
  37 LJLIB_ASM_(math_exp)            LJLIB_REC(math_unary IRFPM_EXP)
  38 LJLIB_ASM_(math_sin)            LJLIB_REC(math_unary IRFPM_SIN)
  39 LJLIB_ASM_(math_cos)            LJLIB_REC(math_unary IRFPM_COS)
  40 LJLIB_ASM_(math_tan)            LJLIB_REC(math_unary IRFPM_TAN)
  41 LJLIB_ASM_(math_asin)           LJLIB_REC(math_atrig FF_math_asin)
  42 LJLIB_ASM_(math_acos)           LJLIB_REC(math_atrig FF_math_acos)
  43 LJLIB_ASM_(math_atan)           LJLIB_REC(math_atrig FF_math_atan)
  44 LJLIB_ASM_(math_sinh)           LJLIB_REC(math_htrig IRCALL_sinh)
  45 LJLIB_ASM_(math_cosh)           LJLIB_REC(math_htrig IRCALL_cosh)
  46 LJLIB_ASM_(math_tanh)           LJLIB_REC(math_htrig IRCALL_tanh)
  47 LJLIB_ASM_(math_frexp)
  48 LJLIB_ASM_(math_modf)           LJLIB_REC(.)
  49
  50 LJLIB_PUSH(57.29577951308232)
  51 LJLIB_ASM_(math_deg)            LJLIB_REC(math_degrad)
  52
  53 LJLIB_PUSH(0.017453292519943295)
  54 LJLIB_ASM_(math_rad)            LJLIB_REC(math_degrad)
  55
  56 LJLIB_ASM(math_log)             LJLIB_REC(math_log)
  57 {
  58   double x = lj_lib_checknum(L, 1);
  59   if (L->base+1 < L->top) {
  60     double y = lj_lib_checknum(L, 2);
  61 #ifdef LUAJIT_NO_LOG2
  62     x = log(x); y = 1.0 / log(y);
  63 #else
  64     x = lj_vm_log2(x); y = 1.0 / lj_vm_log2(y);
  65 #endif
  66     setnumV(L->base-1, x*y);  /* Do NOT join the expression to x / y. */
  67     return FFH_RES(1);
  68   }
  69   return FFH_RETRY;
  70 }
  71
  72 LJLIB_ASM(math_atan2)           LJLIB_REC(.)
  73 {
  74   lj_lib_checknum(L, 1);
  75   lj_lib_checknum(L, 2);
  76   return FFH_RETRY;
  77 }
  78 LJLIB_ASM_(math_pow)            LJLIB_REC(.)
  79 LJLIB_ASM_(math_fmod)
  80
  81 LJLIB_ASM(math_ldexp)           LJLIB_REC(.)
  82 {
  83   lj_lib_checknum(L, 1);
  84 #if LJ_DUALNUM && !LJ_TARGET_X86ORX64
  85   lj_lib_checkint(L, 2);
  86 #else
  87   lj_lib_checknum(L, 2);
  88 #endif
  89   return FFH_RETRY;
  90 }
  91
  92 LJLIB_ASM(math_min)             LJLIB_REC(math_minmax IR_MIN)
  93 {
  94   int i = 0;
  95   do { lj_lib_checknumber(L, ++i); } while (L->base+i < L->top);
  96   return FFH_RETRY;
  97 }
  98 LJLIB_ASM_(math_max)            LJLIB_REC(math_minmax IR_MAX)
  99
 100 LJLIB_PUSH(3.14159265358979323846) LJLIB_SET(pi)
 101 LJLIB_PUSH(1e310) LJLIB_SET(huge)
 102
 103 /* ------------------------------------------------------------------------ */
 104
 105 /* This implements a Tausworthe PRNG with period 2^223. Based on:
 106 **   Tables of maximally-equidistributed combined LFSR generators,
 107 **   Pierre L'Ecuyer, 1991, table 3, 1st entry.
 108 ** Full-period ME-CF generator with L=64, J=4, k=223, N1=49.
 109 */
 110
 111 /* PRNG state. */
 112 struct RandomState {
 113   uint64_t gen[4];      /* State of the 4 LFSR generators. */
 114   int valid;            /* State is valid. */
 115 };
 116
 117 /* Union needed for bit-pattern conversion between uint64_t and double. */
 118 typedef union { uint64_t u64; double d; } U64double;
 119
 120 /* Update generator i and compute a running xor of all states. */
 121 #define TW223_GEN(i, k, q, s) \
 122   z = rs->gen[i]; \
 123   z = (((z<<q)^z) >> (k-s)) ^ ((z&((uint64_t)(int64_t)-1 << (64-k)))<<s); \
 124   r ^= z; rs->gen[i] = z;
 125
 126 /* PRNG step function. Returns a double in the range 1.0 <= d < 2.0. */
 127 LJ_NOINLINE uint64_t LJ_FASTCALL lj_math_random_step(RandomState *rs)
 128 {
 129   uint64_t z, r = 0;
 130   TW223_GEN(0, 63, 31, 18)
 131   TW223_GEN(1, 58, 19, 28)
 132   TW223_GEN(2, 55, 24,  7)
 133   TW223_GEN(3, 47, 21,  8)
 134   return (r & U64x(000fffff,ffffffff)) | U64x(3ff00000,00000000);
 135 }
 136
 137 /* PRNG initialization function. */
 138 static void random_init(RandomState *rs, double d)
 139 {
 140   uint32_t r = 0x11090601;  /* 64-k[i] as four 8 bit constants. */
 141   int i;
 142   for (i = 0; i < 4; i++) {
 143     U64double u;
 144     uint32_t m = 1u << (r&255);
 145     r >>= 8;
 146     u.d = d = d * 3.14159265358979323846 + 2.7182818284590452354;
 147     if (u.u64 < m) u.u64 += m;  /* Ensure k[i] MSB of gen[i] are non-zero. */
 148     rs->gen[i] = u.u64;
 149   }
 150   rs->valid = 1;
 151   for (i = 0; i < 10; i++)
 152     lj_math_random_step(rs);
 153 }
 154
 155 /* PRNG extract function. */
 156 LJLIB_PUSH(top-2)  /* Upvalue holds userdata with RandomState. */
 157 LJLIB_CF(math_random)           LJLIB_REC(.)
 158 {
 159   int n = (int)(L->top - L->base);
 160   RandomState *rs = (RandomState *)(uddata(udataV(lj_lib_upvalue(L, 1))));
 161   U64double u;
 162   double d;
 163   if (LJ_UNLIKELY(!rs->valid)) random_init(rs, 0.0);
 164   u.u64 = lj_math_random_step(rs);
 165   d = u.d - 1.0;
 166   if (n > 0) {
 167 #if LJ_DUALNUM
 168     int isint = 1;
 169     double r1;
 170     lj_lib_checknumber(L, 1);
 171     if (tvisint(L->base)) {
 172       r1 = (lua_Number)intV(L->base);
 173     } else {
 174       isint = 0;
 175       r1 = numV(L->base);
 176     }
 177 #else
 178     double r1 = lj_lib_checknum(L, 1);
 179 #endif
 180     if (n == 1) {
 181       d = lj_vm_floor(d*r1) + 1.0;  /* d is an int in range [1, r1] */
 182     } else {
 183 #if LJ_DUALNUM
 184       double r2;
 185       lj_lib_checknumber(L, 2);
 186       if (tvisint(L->base+1)) {
 187         r2 = (lua_Number)intV(L->base+1);
 188       } else {
 189         isint = 0;
 190         r2 = numV(L->base+1);
 191       }
 192 #else
 193       double r2 = lj_lib_checknum(L, 2);
 194 #endif
 195       d = lj_vm_floor(d*(r2-r1+1.0)) + r1;  /* d is an int in range [r1, r2] */
 196     }
 197 #if LJ_DUALNUM
 198     if (isint) {
 199       setintV(L->top-1, lj_num2int(d));
 200       return 1;
 201     }
 202 #endif
 203   }  /* else: d is a double in range [0, 1] */
 204   setnumV(L->top++, d);
 205   return 1;
 206 }
 207
 208 /* PRNG seed function. */
 209 LJLIB_PUSH(top-2)  /* Upvalue holds userdata with RandomState. */
 210 LJLIB_CF(math_randomseed)
 211 {
 212   RandomState *rs = (RandomState *)(uddata(udataV(lj_lib_upvalue(L, 1))));
 213   random_init(rs, lj_lib_checknum(L, 1));
 214   return 0;
 215 }
 216
 217 /* ------------------------------------------------------------------------ */
 218
 219 #include "lj_libdef.h"
 220
 221 LUALIB_API int luaopen_math(lua_State *L)
 222 {
 223   RandomState *rs;
 224   rs = (RandomState *)lua_newuserdata(L, sizeof(RandomState));
 225   rs->valid = 0;  /* Use lazy initialization to save some time on startup. */
 226   LJ_LIB_REG(L, LUA_MATHLIBNAME, math);
 227 #if defined(LUA_COMPAT_MOD) && !LJ_52
 228   lua_getfield(L, -1, "fmod");
 229   lua_setfield(L, -2, "mod");
 230 #endif
 231   return 1;
 232 }
 233