root/lib/lua/lfmathlib.c

/* [<][>][^][v][top][bottom][index][help] */

DEFINITIONS

This source file includes following 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 */
  57 
  58 #include "stdlib.h"
  59 #include "math.h"
  60 
  61 #define lfmathlib_c
  62 #define LUA_LIB
  63 
  64 #include "lua.h"
  65 
  66 #include "lauxlib.h"
  67 #include "lualib.h"
  68 
  69 #define fabs(x)     ((((x) < 0.0) ? -(x) : (x)))
  70 
  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 }
  80 
  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));
  84     luaL_getmetatable(L, "fmathmeta");
  85     lua_setmetatable(L, -2);
  86      *r = v;
  87     return r;
  88 }
  89 
  90 static int fmath_new (lua_State *L) {
  91     int n = luaL_checknumber(L, 1);
  92     int d = luaL_optnumber(L, 2, 1);
  93     newval(L, (double)n / (double)d);
  94     return 1;
  95 }
  96 
  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
 112     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];
 120         sprintf(fmt, "%s%%u.%%0%du", s, p);
 121         // value string
 122         sprintf(buf, fmt, w, d);
 123     } else {
 124         // value string
 125         sprintf(buf, "%s%d", s, w);
 126     }
 127     lua_pushstring(L, buf);
 128     return 1;
 129 }
 130 
 131 enum {
 132     FN_MUL = 0,
 133     FN_DIV,
 134     FN_ADD,
 135     FN_SUB,
 136     FN_POW,
 137     FN_MOD,
 138     FN_LOG,
 139     FN_LOG2,
 140     FN_LOG10,
 141     FN_SQRT,
 142     FN_NEG,
 143     FN_DEG,
 144     FN_RAD,
 145     FN_EQ,
 146     FN_LT,
 147     FN_LE,
 148     FN_INT,
 149     FN_CEIL,
 150     FN_FLOOR,
 151     FN_ROUND,
 152     FN_SIN,
 153     FN_COS,
 154     FN_TAN,
 155     FN_ASIN,
 156     FN_ACOS,
 157     FN_ATAN,
 158 };
 159 
 160 enum {
 161     ROTATE = 0,
 162     VECTOR
 163 };
 164 
 165 /*
 166     Cordic trig functions.
 167     Copied from cordic_math.c and converted to use double data type.
 168  */
 169  
 170 #define N   17
 171 #define M   9
 172 
 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
 179 
 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) {
 187         xstep = *x / div;
 188         ystep = *y / div;
 189         if (i < M) {
 190             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;
 197         div = div * 2.0;
 198     }
 199 }
 200 
 201 #define trunc(x)    ((double)((int)(x)))
 202 
 203 static double fmod(double x, double y) {
 204     return x - trunc(x / y) * y;
 205 }
 206 
 207 static double cathetus(double x) {
 208     return sqrt((1 + x) * (1 - x));
 209 }
 210 
 211 // base CIRCULAR mode, ROTATE
 212 static void sincosCordic(double phi, double *sinphi, double *cosphi) {
 213     // convert to quadrant 1
 214     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     }
 221 
 222     double x = INV_GAIN_CIRCLE, y = 0, z = phi;
 223     cordic(ROTATE, &x, &y, &z);
 224 
 225     // convert to original quadrant
 226     if ((q == 2) || (q == 3)) { x = -x; }
 227     if ((q == 3) || (q == 4)) { y = -y; }
 228 
 229     *sinphi = y;
 230     *cosphi = x;
 231 }
 232 
 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     }
 241 
 242     // convert to quadrant 1
 243     int fy = (y >= 0);
 244     int q = 1;
 245     if (x < 0) {
 246         if (fy)
 247             q = 2;
 248         else
 249             q = 3;
 250     } else if (!fy) {
 251         q = 4;
 252     }
 253     x = fabs(x);
 254     y = fabs(y);
 255 
 256     double z = 0;
 257     cordic(VECTOR, &x, &y, &z);
 258 
 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     }
 265 
 266     *phi = z;
 267     *hyp = x;
 268 }
 269 
 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 */
 277 
 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:
 284         r = a * b;
 285         break;
 286     case FN_DIV:
 287         r = a / b;
 288         break;
 289     case FN_ADD:
 290         r = a + b;
 291         break;
 292     case FN_SUB:
 293         r = a - b;
 294         break;
 295     case FN_POW:
 296         r = pow(a, b);
 297         break;
 298     case FN_MOD:
 299         r = fmod(a, b);
 300         break;
 301     default:
 302         break;
 303     }
 304     newval(L, r);
 305     return 1;
 306 }
 307 
 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:
 313         r = log(a);
 314         break;
 315     case FN_LOG2:
 316         r = log2(a);
 317         break;
 318     case FN_LOG10:
 319         r = log10(a);
 320         break;
 321     case FN_SQRT:
 322         r = sqrt(a);
 323         break;
 324     case FN_NEG:
 325         r = -a;
 326         break;
 327     case FN_DEG:
 328         r =  a * 180.0 / M_PI;
 329         break;
 330     case FN_RAD:
 331         r =  a * M_PI / 180.0;
 332         break;
 333     default:
 334         break;
 335     }
 336     newval(L, r);
 337     return 1;
 338 }
 339 
 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:
 346         r = (a == b);
 347         break;
 348     case FN_LT:
 349         r = (a < b);
 350         break;
 351     case FN_LE:
 352         r = (a <= b);
 353         break;
 354     default:
 355         break;
 356     }
 357     lua_pushboolean(L, r);
 358     return 1;
 359 }
 360 
 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:
 367         r = (int)a;
 368         break;
 369     case FN_CEIL:
 370         if ((a > 0.0) && (frac != 0.0))
 371             r = (int)a + 1;
 372         else
 373             r = (int)a;
 374         break;
 375     case FN_FLOOR:
 376         if ((a < 0.0) && (frac != 0.0))
 377             r = (int)a - 1;
 378         else
 379             r = (int)a;
 380         break;
 381     case FN_ROUND:
 382         if (a < 0.0)
 383             r = (int)(a - 0.5);
 384         else
 385             r = (int)(a + 0.5);
 386         break;
 387     default:
 388         break;
 389     }
 390     lua_pushnumber(L, r);
 391     return 1;
 392 }
 393 
 394 static int trigfn(lua_State *L, int fn) {
 395     double a = arg(L, 1);
 396     double _sin, _cos;
 397     double r = 0;
 398     sincosCordic(a, &_sin, &_cos);
 399     switch(fn) {
 400     case FN_SIN:
 401         r = _sin;
 402         break;
 403     case FN_COS:
 404         r = _cos;
 405         break;
 406     case FN_TAN:
 407         r = _sin / _cos;
 408         break;
 409     default:
 410         break;
 411     }
 412     newval(L, r);
 413     return 1;
 414 }
 415 
 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:
 421         y = x;
 422         x = cathetus(x);
 423         break;
 424     case FN_ACOS:
 425         y = cathetus(x);
 426         break;
 427     case FN_ATAN:
 428         y = x;
 429         x = 1.0;
 430         break;
 431     default:
 432         break;
 433     }
 434     double phi, hyp;
 435     atanhypCordic(x, y, &phi, &hyp);
 436     newval(L, phi);
 437     return 1;
 438 }
 439 
 440 static int fmath_rec(lua_State *L) {
 441     double r = arg(L, 1);
 442     double theta = arg(L, 2);
 443     double _sin, _cos;
 444     sincosCordic(theta, &_sin, &_cos);
 445     newval(L, r * _cos);
 446     newval(L, r * _sin);
 447     return 2;
 448 }
 449 
 450 static int fmath_pol (lua_State *L) {
 451     double px = arg(L, 1);
 452     double py = arg(L, 2);
 453     double phi, hyp;
 454     atanhypCordic(px, py, &phi, &hyp);
 455     newval(L, hyp * INV_GAIN_CIRCLE);
 456     newval(L, phi);
 457     return 2;
 458 }
 459 
 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); }
 466 
 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); }
 470 
 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); }
 478 
 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); }
 483 
 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); }
 487 
 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); }
 491 
 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 };
 526 
 527 /*
 528 ** Open fmath library
 529 */
 530 LUALIB_API int luaopen_fmath (lua_State *L) {
 531     luaL_newmetatable(L, "fmathmeta");
 532     lua_pushstring(L, "__index");
 533     lua_pushvalue(L, -2);   /* pushes the metatable */
 534     lua_settable(L, -3);    /* metatable.__index = metatable */
 535     luaL_register(L, NULL, fmathlib_m);
 536 
 537     luaL_register(L, LUA_FMATHLIBNAME, &fmathlib_m[7]);     // adjust offset if table is changed
 538 
 539     newval(L, M_PI * 2.0);
 540     lua_setfield(L, -2, "pi2");
 541     newval(L, M_PI);
 542     lua_setfield(L, -2, "pi");
 543     newval(L, M_PI / 2.0);
 544     lua_setfield(L, -2, "pi_2");
 545 
 546     return 1;
 547 }

/* [<][>][^][v][top][bottom][index][help] */