## root/lib/lua/lfmathlib.c

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

## DEFINITIONS

This source file includes following definitions.
```   1 /*
2 lfmathlib.c,v 1.0
3 double precision floating point based mathematics for CHDK Lua scripts
4
5 (c)2021 philmoz
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,
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,
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;
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;
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); }
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); }
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[] = {
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},
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] */