root/lib/lua/lfmathlib.c

DEFINITIONS

1. arg
2. newval
3. fmath_new
4. fmath_toStr
5. cordic
6. fmod
7. cathetus
8. sincosCordic
9. atanhypCordic
10. twoargfn
11. oneargfn
12. boolfn
13. intfn
14. trigfn
15. atrigfn
16. fmath_rec
17. fmath_pol
18. fmath_mul
19. fmath_div
20. fmath_add
21. fmath_sub
22. fmath_pow
23. fmath_mod
24. fmath_eq
25. fmath_lt
26. fmath_le
27. fmath_neg
28. fmath_log
29. fmath_log2
30. fmath_log10
31. fmath_sqrt
32. fmath_deg
33. fmath_rad
34. fmath_int
35. fmath_ceil
36. fmath_floor
37. fmath_round
38. fmath_sin
39. fmath_cos
40. fmath_tan
41. fmath_asin
42. fmath_acos
43. fmath_atan
44. luaopen_fmath

   1 /*
   2 lfmathlib.c,v 1.0
   3 double precision floating point based mathematics for CHDK Lua scripts
   4 
   5 (c)2021 philmoz
   6 License GPL 2.0
   7 
   8 fmath library usage in Lua:
   9 
  10 Create a new fmath double variable:
  11     x = fmath.new(n,d)
  12 Set x to n / d.
  13 E.G.
  14     x = fmath.new(45,10) - x = 4.5
  15     y = fmath.new(25,100) = y = 0.25
  16 Parameter d is optional and defaults to 1. x = fmath.new(5) is equivalent to x = fmath.new(5,1)
  17 
  18 Arithmetic - use standard operators.
  19 E.G.
  20     z = x + y
  21     z = x - y
  22     z = x * y
  23     z = x / y
  24     z = -x      (unary minus)
  25     z = x ^ y   (power function)
  26     z = x % y   (modulus)
  27 So long as either parameter is an fmath double the result will be an fmath double. The following are both valid:
  28     z = x + 1
  29     z = 1 + x
  30 
  31 Comparison - as above use standard operators ==, ~=, <, <=, >, >=
  32 Note: for comparison operators both sides must be fmath doubles. You cannot mix fmath double and integer values in comparisons.
  33 
  34 Math functions:
  35     log, log2, log10, sqrt
  36 Can be called via fmath library or object oriented syntax. The following are equivalent:
  37     y = fmath.log(x)    ==      y = x:log()
  38     y = fmath.sqrt(x)   ==      y = x:sqrt();
  39 
  40 Conversion:
  41     int, ceil, floor, round, deg, rad
  42 These can also be called via the fmath library or object oriented syntax as above.
  43 The int, ceil, floor and round functions return integer values not fmath doubles.
  44 The deg and rad functions convert between degrees and radians and return an fmath double value.
  45 
  46 The frac function from imathlib is not implemented. The imathlib function is mathematically wrong for negative values.
  47 You can get the equivalent of the imathlib frac function using:
  48     f = (x - x:int()) * 1000
  49 The correct value would be:
  50     f = (x - x:floor()) * 1000
  51 
  52 Trigonometry:
  53     sin, cos, tan, asin, acos, atan, pol, rec
  54 The same algorithms as used by the imathlib library are used here.
  55 Note: all functions operate on or return values in radians.
  56 */
  58 #include "stdlib.h"
  59 #include "math.h"
  61 #define lfmathlib_c
  62 #define LUA_LIB
  64 #include "lua.h"
  66 #include "lauxlib.h"
  67 #include "lualib.h"
  69 #define fabs(x)     ((((x) < 0.0) ? -(x) : (x)))
  71 // get function argument, if not valid fmath 'double' check for integer
  72 static double arg(lua_State *L, int n) {
  73     double* v = (double*)lua_touserdata(L, n);
  74     if (v == 0) {
  75         int i = luaL_checknumber(L, n);
  76         return (double)i;
  77     }
  78     return *v;
  79 }
  81 // create a new fmath 'double' with given value
  82 static double* newval(lua_State *L, double v) {
  83     double* r = (double*)lua_newuserdata(L, sizeof(double));
luaL_getmetatable(L, "fmathmeta");
lua_setmetatable(L, -2);
  86      *r = v;
  87     return r;
  88 }
  90 static int fmath_new (lua_State *L) {
  91     int n = luaL_checknumber(L, 1);
  92     int d = luaL_optnumber(L, 2, 1);
newval(L, (double)n / (double)d);
  94     return 1;
  95 }
  97 // Due to lack of support for %f format string on some cameras the toStr function uses a simple integer conversion for formatting
  98 // The optional 2nd parameter specifies the number of decimal digits to output - default is 6.
  99 static int fmath_toStr (lua_State *L) {
 100     char buf[100];
 101     // value
 102     double a = arg(L, 1);
 103     // precision
 104     int p = luaL_optnumber(L, 2, 6);
 105     if (p < 0) p = 0;
 106     if (p > 9) p = 9;
 107     // decimal scaling / rounding factor
 108     double m = pow(10.0, p);
 109     // sign
 110     const char* s = (a < 0) ? "-" : "";
 111     // round absolute value
a = fabs(a) + (0.5 / m);
 113     // whole part
 114     unsigned int w = (int)a;
 115     if (p > 0) {
 116         // decimal part
 117         unsigned int d = (int)((a - w) * m);
 118         // format string
 119         char fmt[20];
sprintf(fmt, "%s%%u.%%0%du", s, p);
 121         // value string
sprintf(buf, fmt, w, d);
 123     } else {
 124         // value string
sprintf(buf, "%s%d", s, w);
 126     }
lua_pushstring(L, buf);
 128     return 1;
 129 }
 131 enum {
FN_MUL = 0,
FN_DIV,
FN_ADD,
FN_SUB,
FN_POW,
FN_MOD,
FN_LOG,
FN_LOG2,
FN_LOG10,
FN_SQRT,
FN_NEG,
FN_DEG,
FN_RAD,
FN_EQ,
FN_LT,
FN_LE,
FN_INT,
FN_CEIL,
FN_FLOOR,
FN_ROUND,
FN_SIN,
FN_COS,
FN_TAN,
FN_ASIN,
FN_ACOS,
FN_ATAN,
 158 };
 160 enum {
ROTATE = 0,
VECTOR
 163 };
 165 /*
 166     Cordic trig functions.
 167     Copied from cordic_math.c and converted to use double data type.
 168  */
 170 #define N   17
 171 #define M   9
 173 static double atan_tab[M] = {
 174     0.78539816339745, 0.46364760900081, 0.24497866312686, 0.12435499454676,
 175     0.06241880999596, 0.03123983343027, 0.01562372862048, 0.00781234106010,
 176     0.00390623013197
 177 };
 178 #define INV_GAIN_CIRCLE     0.60725293501477
 180 // CORDIC main routine
 181 static void cordic(int f, double *x, double *y, double *z) {
 182     double *patan_tab = atan_tab;
 183     double xstep, ystep, zstep = 1;
 184     double div = 1;
 185     int i;
 186     for (i = 0; i < N; i += 1) {
xstep = *x / div;
ystep = *y / div;
 189         if (i < M) {
zstep = *patan_tab;
 191             ++patan_tab;
 192         } else zstep = zstep / 2.0;
 193         int f1 = (f) ? *y >= 0 : *z < 0;
 194         *x = (f1) ? *x + ystep : *x - ystep;
 195         *y = (f1) ? *y - xstep : *y + xstep;
 196         *z = (f1) ? *z + zstep : *z - zstep;
div = div * 2.0;
 198     }
 199 }
 201 #define trunc(x)    ((double)((int)(x)))
 203 static double fmod(double x, double y) {
 204     return x - trunc(x / y) * y;
 205 }
 207 static double cathetus(double x) {
 208     return sqrt((1 + x) * (1 - x));
 209 }
 211 // base CIRCULAR mode, ROTATE
 212 static void sincosCordic(double phi, double *sinphi, double *cosphi) {
 213     // convert to quadrant 1
phi = fmod(fmod(phi, M_PI * 2.0) + (M_PI * 2.0), M_PI * 2.0);
 215     int q = (int)(phi / (M_PI / 2.0)) + 1;
 216     switch(q) {
 217         case 2: phi = M_PI - phi; break;
 218         case 3: phi = phi - M_PI; break;
 219         case 4: phi = M_PI * 2.0 - phi; break;
 220     }
 222     double x = INV_GAIN_CIRCLE, y = 0, z = phi;
cordic(ROTATE, &x, &y, &z);
 225     // convert to original quadrant
 226     if ((q == 2) || (q == 3)) { x = -x; }
 227     if ((q == 3) || (q == 4)) { y = -y; }
 229     *sinphi = y;
 230     *cosphi = x;
 231 }
 233 // base CIRCULAR mode, VECTOR
 234 static void atanhypCordic(double x, double y, double *phi, double *hyp) {
 235     if (x == 0 && y == 0) {
 236         // error input vales
 237         *phi = 0;
 238         *hyp = 0;
 239         return;
 240     }
 242     // convert to quadrant 1
 243     int fy = (y >= 0);
 244     int q = 1;
 245     if (x < 0) {
 246         if (fy)
q = 2;
 248         else
q = 3;
 250     } else if (!fy) {
q = 4;
 252     }
x = fabs(x);
y = fabs(y);
 256     double z = 0;
cordic(VECTOR, &x, &y, &z);
 259     // convert to original quadrant
 260     switch(q) {
 261         case 2: z = M_PI - z; break;
 262         case 3: z = z - M_PI; break;
 263         case 4: z = -z; break;
 264     }
 266     *phi = z;
 267     *hyp = x;
 268 }
 270 /*
 271     A note on the implementation below.
 272     Floating point math operations require long calls to core CHDK functions.
 273     GCC currently does not merge the long call offsets to the same core function across different functions here.
 274     So two local functions that perform floating point division will end up with two copies of the core division function offset.
 275     The implementation here with switch statements is a bit messy; but helps reduce the overhead and compiled size.
 276 */
 278 static int twoargfn (lua_State *L, int fn) {
 279     double a = arg(L, 1);
 280     double b = arg(L, 2);
 281     double r = 0;
 282     switch(fn) {
 283     case FN_MUL:
r = a * b;
 285         break;
 286     case FN_DIV:
r = a / b;
 288         break;
 289     case FN_ADD:
r = a + b;
 291         break;
 292     case FN_SUB:
r = a - b;
 294         break;
 295     case FN_POW:
r = pow(a, b);
 297         break;
 298     case FN_MOD:
r = fmod(a, b);
 300         break;
 301     default:
 302         break;
 303     }
newval(L, r);
 305     return 1;
 306 }
 308 static int oneargfn (lua_State *L, int fn) {
 309     double a = arg(L, 1);
 310     double r = 0;
 311     switch(fn) {
 312     case FN_LOG:
r = log(a);
 314         break;
 315     case FN_LOG2:
r = log2(a);
 317         break;
 318     case FN_LOG10:
r = log10(a);
 320         break;
 321     case FN_SQRT:
r = sqrt(a);
 323         break;
 324     case FN_NEG:
r = -a;
 326         break;
 327     case FN_DEG:
r =  a * 180.0 / M_PI;
 329         break;
 330     case FN_RAD:
r =  a * M_PI / 180.0;
 332         break;
 333     default:
 334         break;
 335     }
newval(L, r);
 337     return 1;
 338 }
 340 static int boolfn (lua_State *L, int fn) {
 341     double a = arg(L, 1);
 342     double b = arg(L, 2);
 343     int r = 0;
 344     switch(fn) {
 345     case FN_EQ:
r = (a == b);
 347         break;
 348     case FN_LT:
r = (a < b);
 350         break;
 351     case FN_LE:
r = (a <= b);
 353         break;
 354     default:
 355         break;
 356     }
lua_pushboolean(L, r);
 358     return 1;
 359 }
 361 static int intfn (lua_State *L, int fn) {
 362     double a = arg(L, 1);
 363     double frac = a - (int)a;
 364     int r = 0;
 365     switch(fn) {
 366     case FN_INT:
r = (int)a;
 368         break;
 369     case FN_CEIL:
 370         if ((a > 0.0) && (frac != 0.0))
r = (int)a + 1;
 372         else
r = (int)a;
 374         break;
 375     case FN_FLOOR:
 376         if ((a < 0.0) && (frac != 0.0))
r = (int)a - 1;
 378         else
r = (int)a;
 380         break;
 381     case FN_ROUND:
 382         if (a < 0.0)
r = (int)(a - 0.5);
 384         else
r = (int)(a + 0.5);
 386         break;
 387     default:
 388         break;
 389     }
lua_pushnumber(L, r);
 391     return 1;
 392 }
 394 static int trigfn(lua_State *L, int fn) {
 395     double a = arg(L, 1);
 396     double _sin, _cos;
 397     double r = 0;
sincosCordic(a, &_sin, &_cos);
 399     switch(fn) {
 400     case FN_SIN:
r = _sin;
 402         break;
 403     case FN_COS:
r = _cos;
 405         break;
 406     case FN_TAN:
r = _sin / _cos;
 408         break;
 409     default:
 410         break;
 411     }
newval(L, r);
 413     return 1;
 414 }
 416 static int atrigfn(lua_State *L, int fn) {
 417     double x = arg(L, 1);
 418     double y = 0;
 419     switch(fn) {
 420     case FN_ASIN:
y = x;
x = cathetus(x);
 423         break;
 424     case FN_ACOS:
y = cathetus(x);
 426         break;
 427     case FN_ATAN:
y = x;
x = 1.0;
 430         break;
 431     default:
 432         break;
 433     }
 434     double phi, hyp;
atanhypCordic(x, y, &phi, &hyp);
newval(L, phi);
 437     return 1;
 438 }
 440 static int fmath_rec(lua_State *L) {
 441     double r = arg(L, 1);
 442     double theta = arg(L, 2);
 443     double _sin, _cos;
sincosCordic(theta, &_sin, &_cos);
newval(L, r * _cos);
newval(L, r * _sin);
 447     return 2;
 448 }
 450 static int fmath_pol (lua_State *L) {
 451     double px = arg(L, 1);
 452     double py = arg(L, 2);
 453     double phi, hyp;
atanhypCordic(px, py, &phi, &hyp);
newval(L, hyp * INV_GAIN_CIRCLE);
newval(L, phi);
 457     return 2;
 458 }
 460 static int fmath_mul (lua_State *L) { return twoargfn(L, FN_MUL); }
 461 static int fmath_div (lua_State *L) { return twoargfn(L, FN_DIV); }
 462 static int fmath_add (lua_State *L) { return twoargfn(L, FN_ADD); }
 463 static int fmath_sub (lua_State *L) { return twoargfn(L, FN_SUB); }
 464 static int fmath_pow (lua_State *L) { return twoargfn(L, FN_POW); }
 465 static int fmath_mod (lua_State *L) { return twoargfn(L, FN_MOD); }
 467 static int fmath_eq (lua_State *L) { return boolfn(L, FN_EQ); }
 468 static int fmath_lt (lua_State *L) { return boolfn(L, FN_LT); }
 469 static int fmath_le (lua_State *L) { return boolfn(L, FN_LE); }
 471 static int fmath_neg (lua_State *L) { return oneargfn(L, FN_NEG); }
 472 static int fmath_log (lua_State *L) { return oneargfn(L, FN_LOG); }
 473 static int fmath_log2 (lua_State *L) { return oneargfn(L, FN_LOG2); }
 474 static int fmath_log10 (lua_State *L) { return oneargfn(L, FN_LOG10); }
 475 static int fmath_sqrt (lua_State *L) { return oneargfn(L, FN_SQRT); }
 476 static int fmath_deg (lua_State *L) { return oneargfn(L, FN_DEG); }
 477 static int fmath_rad (lua_State *L) { return oneargfn(L, FN_RAD); }
 479 static int fmath_int (lua_State *L) { return intfn(L, FN_INT); }
 480 static int fmath_ceil (lua_State *L) { return intfn(L, FN_CEIL); }
 481 static int fmath_floor (lua_State *L) { return intfn(L, FN_FLOOR); }
 482 static int fmath_round (lua_State *L) { return intfn(L, FN_ROUND); }
 484 static int fmath_sin (lua_State *L) { return trigfn(L, FN_SIN); }
 485 static int fmath_cos (lua_State *L) { return trigfn(L, FN_COS); }
 486 static int fmath_tan (lua_State *L) { return trigfn(L, FN_TAN); }
 488 static int fmath_asin (lua_State *L) { return atrigfn(L, FN_ASIN); }
 489 static int fmath_acos (lua_State *L) { return atrigfn(L, FN_ACOS); }
 490 static int fmath_atan (lua_State *L) { return atrigfn(L, FN_ATAN); }
 492 static const luaL_Reg fmathlib_m[] = {
 493     {"__add",       fmath_add},
 494     {"__sub",       fmath_sub},
 495     {"__mul",       fmath_mul},
 496     {"__div",       fmath_div},
 497     {"__pow",       fmath_pow},
 498     {"__mod",       fmath_mod},
 499     {"__unm",       fmath_neg},
 500     {"__eq",        fmath_eq},
 501     {"__le",        fmath_le},
 502     {"__lt",        fmath_lt},
 503 // If any new entries are added above then adjust the offset in the luaL_register call below
 504     {"new",         fmath_new},
 505     {"tostr",       fmath_toStr},
 506     {"log",         fmath_log},
 507     {"log2",        fmath_log2},
 508     {"log10",       fmath_log10},
 509     {"sqrt",        fmath_sqrt},
 510     {"int",         fmath_int},
 511     {"ceil",        fmath_ceil},
 512     {"floor",       fmath_floor},
 513     {"round",       fmath_round},
 514     {"deg",         fmath_deg},
 515     {"rad",         fmath_rad},
 516     {"sin",         fmath_sin},
 517     {"cos",         fmath_cos},
 518     {"tan",         fmath_tan},
 519     {"asin",        fmath_asin},
 520     {"acos",        fmath_acos},
 521     {"atan",        fmath_atan},
 522     {"pol",         fmath_pol},
 523     {"rec",         fmath_rec},
 524     {NULL, NULL}
 525 };
 527 /*
 528 ** Open fmath library
 529 */
LUALIB_API int luaopen_fmath (lua_State *L) {
luaL_newmetatable(L, "fmathmeta");
lua_pushstring(L, "__index");
lua_pushvalue(L, -2);   /* pushes the metatable */
lua_settable(L, -3);    /* metatable.__index = metatable */
luaL_register(L, NULL, fmathlib_m);
luaL_register(L, LUA_FMATHLIBNAME, &fmathlib_m[7]);     // adjust offset if table is changed
newval(L, M_PI * 2.0);
lua_setfield(L, -2, "pi2");
newval(L, M_PI);
lua_setfield(L, -2, "pi");
newval(L, M_PI / 2.0);
lua_setfield(L, -2, "pi_2");
 546     return 1;
 547 }