Mercurial > ~darius > hgwebdir.cgi > modulator
changeset 14:388074ff9474
Add fixed point code
| author | Daniel O'Connor <darius@dons.net.au> |
|---|---|
| date | Tue, 25 Feb 2025 13:28:29 +1030 |
| parents | 032acb7fbc04 |
| children | bf483219cb12 |
| files | q/LICENSE q/expr.c q/makefile q/q.c q/q.h q/readme.md q/t.c |
| diffstat | 7 files changed, 3157 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/q/LICENSE Tue Feb 25 13:28:29 2025 +1030 @@ -0,0 +1,24 @@ +This is free and unencumbered software released into the public domain. + +Anyone is free to copy, modify, publish, use, compile, sell, or +distribute this software, either in source code form or as a compiled +binary, for any purpose, commercial or non-commercial, and by any +means. + +In jurisdictions that recognize copyright laws, the author or authors +of this software dedicate any and all copyright interest in the +software to the public domain. We make this dedication for the benefit +of the public at large and to the detriment of our heirs and +successors. We intend this dedication to be an overt act of +relinquishment in perpetuity of all present and future rights to this +software under copyright law. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. +IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR +OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, +ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +OTHER DEALINGS IN THE SOFTWARE. + +For more information, please refer to <http://unlicense.org/>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/q/expr.c Tue Feb 25 13:28:29 2025 +1030 @@ -0,0 +1,245 @@ +/* Project: Small Expression Evaluator for Q library + * Author: Richard James Howe + * License: The Unlicense + * See: <https://en.wikipedia.org/wiki/Shunting-yard_algorithm> */ + +#include <assert.h> +#include <ctype.h> +#include <limits.h> +#include <stdarg.h> +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include "q.h" + +#define BUILD_BUG_ON(condition) ((void)sizeof(char[1 - 2*!!(condition)])) +#define DEFAULT_STACK_SIZE (64) + +enum { ASSOCIATE_NONE, ASSOCIATE_LEFT, ASSOCIATE_RIGHT, }; +enum { LEX_NUMBER, LEX_OPERATOR, LEX_END, }; + +static void expr_delete(qexpr_t *e) { + if (!e) + return; + free(e->ops); + free(e->numbers); + for (size_t i = 0; i < e->vars_max; i++) { + free(e->vars[i]->name); + free(e->vars[i]); + } + free(e->vars); + free(e); +} + +static qexpr_t *expr_new(size_t max) { + max = max ? max : 64; + qexpr_t *e = calloc(sizeof(*e), 1); + if (!e) + goto fail; + e->ops = calloc(sizeof(*e->ops), max); + e->numbers = calloc(sizeof(*(e->numbers)), max); + if (!(e->ops) || !(e->numbers)) + goto fail; + e->ops_max = max; + e->numbers_max = max; + qexpr_init(e); + return e; +fail: + expr_delete(e); + return NULL; +} + +static qvariable_t *variable_lookup(qexpr_t *e, const char *name) { + assert(e); + for (size_t i = 0; i < e->vars_max; i++) { + qvariable_t *v = e->vars[i]; + if (!strcmp(v->name, name)) + return v; + } + return NULL; +} + +static char *estrdup(const char *s) { + assert(s); + const size_t l = strlen(s) + 1; + char *r = malloc(l); + return memcpy(r, s, l); +} + +static int variable_name_is_valid(const char *n) { + assert(n); + if (!isalpha(*n) && !(*n == '_')) + return 0; + for (n++; *n; n++) + if (!isalnum(*n) && !(*n == '_')) + return 0; + return 1; +} + +static qvariable_t *variable_add(qexpr_t *e, const char *name, q_t value) { + assert(e); + assert(name); + qvariable_t *v = variable_lookup(e, name), **vs = e->vars; + if (v) { + v->value = value; + return v; + } + if (!variable_name_is_valid(name)) + return NULL; + char *s = estrdup(name); + vs = realloc(e->vars, (e->vars_max + 1) * sizeof(*vs)); + v = calloc(1, sizeof(*v)); + if (!vs || !v || !s) + goto fail; + v->name = s; + v->value = value; + vs[e->vars_max++] = v; + e->vars = vs; + return v; +fail: + free(v); + free(s); + free(vs); + if (vs != e->vars) + free(e->vars); + return NULL; +} + +static inline int tests(FILE *out) { + assert(out); + int report = 0; + static const struct test { + int r; + q_t result; + const char *expr; + } tests[] = { // NB. Variables defined later. + { -1, QINT( 0), "" }, + { -1, QINT( 0), "(" }, + { -1, QINT( 0), ")" }, + { -1, QINT( 0), "2**3" }, + { 0, QINT( 0), "0" }, + { 0, QINT( 2), "1+1" }, + { 0, -QINT( 1), "-1" }, + { 0, QINT( 1), "--1" }, + { 0, QINT(14), "2+(3*4)" }, + { 0, QINT(23), "a+(b*5)" }, + { -1, QINT(14), "(2+(3* 4)" }, + { -1, QINT(14), "2+(3*4)(" }, + { 0, QINT(14), "2+3*4" }, + { 0, QINT( 0), " 2==3 " }, + { 0, QINT( 1), "2 ==2" }, + { 0, QINT( 1), "2== (1+1)" }, + //{ 0, QINT( 8), "2 pow 3" }, + //{ -1, QINT( 0), "2pow3" }, + { 0, QINT(20), "(2+3)*4" }, + { 0, -QINT( 4), "(2+(-3))*4" }, + { -1, QINT( 0), "1/0" }, + { -1, QINT( 0), "1%0" }, + { 0, QINT(50), "100/2" }, + { 0, QINT( 2), "1--1", }, + { 0, QINT( 0), "1---1", }, + }; + + if (fputs("Running Built In Self Tests:\n", out) < 0) + report = -1; + const size_t length = sizeof (tests) / sizeof (tests[0]); + for (size_t i = 0; i < length; i++) { + qexpr_t *e = expr_new(64); + const struct test *test = &tests[i]; + if (!e) { + (void)fprintf(out, "test failed (unable to allocate)\n"); + report = -1; + goto end; + } + + qvariable_t *v1 = variable_add(e, "a", QINT(3)); + qvariable_t *v2 = variable_add(e, "b", QINT(4)); + qvariable_t *v3 = variable_add(e, "c", -QINT(5)); + if (!v1 || !v2 || !v3) { + (void)fprintf(out, "test failed (unable to assign variable)\n"); + report = -1; + goto end; + } + + const int r = qexpr(e, test->expr); + const q_t tos = e->numbers[0]; + const int pass = (r == test->r) && (r != 0 || tos == test->result); + if (fprintf(out, "%s: r(%2d), eval(\"%s\") = %lg \n", + pass ? " ok" : " FAIL", r, test->expr, (double)tos) < 0) + report = -1; + if (!pass) { + (void)fprintf(out, "\tExpected: r(%2d), %lg\n", + test->r, (double)(test->result)); + report = -1; + } + expr_delete(e); + } +end: + if (fprintf(out, "Tests Complete: %s\n", report == 0 ? "pass" : "FAIL") < 0) + report = -1; + return report; +} + +static qexpr_t *expr_new_with_vars(size_t max) { + qexpr_t *e = expr_new(max); + if (!e) return NULL; + if (!variable_add(e, "whole", qinfo.whole)) goto fail; + if (!variable_add(e, "fractional", qinfo.fractional)) goto fail; + if (!variable_add(e, "bit", qinfo.bit)) goto fail; + if (!variable_add(e, "smallest", qinfo.min)) goto fail; + if (!variable_add(e, "biggest", qinfo.max)) goto fail; + if (!variable_add(e, "pi", qinfo.pi)) goto fail; + if (!variable_add(e, "e", qinfo.e)) goto fail; + if (!variable_add(e, "sqrt2", qinfo.sqrt2)) goto fail; + if (!variable_add(e, "sqrt3", qinfo.sqrt3)) goto fail; + if (!variable_add(e, "ln2", qinfo.ln2)) goto fail; + if (!variable_add(e, "ln10", qinfo.ln10)) goto fail; + return e; +fail: + expr_delete(e); + return NULL; +} + +static int usage(FILE *out, const char *arg0) { + assert(out); + assert(arg0); + return fprintf(out, "usage: %s expr\n", arg0); +} + +int main(int argc, char *argv[]) { + + if (argc == 1) { + if (usage(stderr, argv[0]) < 0) { return 1; } + return tests(stderr); + } + + if (argc < 2) { + (void)fprintf(stderr, "usage: %s expr\n", argv[0]); + return 1; + } + + int r = 0; + + for (int i = 1; i < argc; i++) { + qexpr_t *e = expr_new_with_vars(0); + if (!e) { + (void)fprintf(stderr, "allocate failed\n"); + r = 1; + break; + } + if (qexpr(e, argv[i]) == 0) { + char n[64 + 1] = { 0, }; + qsprint(e->numbers[0], n, sizeof n); + if (printf("%s\n", n) < 0) + r = 1; + r = 0; + } else { + (void)fprintf(stderr, "error: %s\n", e->error_string); + r = 1; + } + expr_delete(e); + } + return r; +} + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/q/makefile Tue Feb 25 13:28:29 2025 +1030 @@ -0,0 +1,34 @@ +QVERSION=0x000902 +CFLAGS=-std=c99 -Wall -Wextra -O2 -pedantic -fwrapv -DQVERSION=${QVERSION} -Wmissing-prototypes +CC=cc +TARGET=q +RM=rm -fv +AR=ar +RANLIB=ranlib + +.PHONY: all run clean + +all: ${TARGET} expr + +run: test ${TARGET} t.q + ./${TARGET} t.q + +test: ${TARGET} + ./${TARGET} -t + + +q.o: q.c q.h + +lib${TARGET}.a: q.o + ${AR} rcs $@ $< + ${RANLIB} $@ + +check: *.c *.h + cppcheck --enable=all *.c *.h + +${TARGET}: lib${TARGET}.a t.o + +expr: lib${TARGET}.a expr.o + +clean: + git clean -dffx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/q/q.c Tue Feb 25 13:28:29 2025 +1030 @@ -0,0 +1,1839 @@ +/* Project: Q-Number (Q16.16, signed) library + * Author: Richard James Howe + * License: The Unlicense + * Email: howe.r.j.89@gmail.com + * Repo: <https://github.com/q> + * + * + * A Q32.32 version would be useful. + * + * The following should be changed/done for this library: + * + * - Moving towards a header-only model. + * - Removal of dependencies such as 'isalpha', 'tolower' + * as they are locale dependent. + * - Make components optional (filters, expression parser, ...) + * - Make hyperbolic arc sin/cos/tan functions. + * - Fix bugs / inaccuracies in CORDIC code. + * - Improve accuracy of all the functions and quantify error and + * their limits. + * + * BUG: Enter: 2.71791, get 2.0625, 2.7179 works fine. (Need to + * limit decimal places). + */ + +#include "q.h" +#include <assert.h> +#include <ctype.h> +#include <inttypes.h> +#include <limits.h> +#include <stdarg.h> /* for expression evaluator error handling */ +#include <stdio.h> /* vsnprintf, for expression evaluator */ +#include <string.h> + +#define UNUSED(X) ((void)(X)) +#define BOOLIFY(X) (!!(X)) +#define BUILD_BUG_ON(condition) ((void)sizeof(char[1 - 2*!!(condition)])) +#define MULTIPLIER (INT16_MAX) +#define DMIN (INT32_MIN) +#define DMAX (INT32_MAX) +#define MIN(X, Y) ((X) < (Y) ? (X) : (Y)) +#define MAX(X, Y) ((X) < (Y) ? (Y) : (X)) + +#ifndef CONFIG_Q_HIDE_FUNCS /* 1 = hide hidden (testing) functions, 0 = enable them */ +#define CONFIG_Q_HIDE_FUNCS (0) +#endif + +typedef int16_t hd_t; /* half Q width, signed */ +typedef uint64_t lu_t; /* double Q width, unsigned */ + +const qinfo_t qinfo = { + .whole = QBITS, + .fractional = QBITS, + .zero = (u_t)0uL << QBITS, + .bit = 1uL, + .one = (u_t)1uL << QBITS, + .min = (u_t)(QHIGH << QBITS), + .max = (u_t)((QHIGH << QBITS) - 1uL), + + .pi = QPI, /* 3.243F6 A8885 A308D 31319 8A2E0... */ + .e = QMK(0x2, 0xB7E1, 16), /* 2.B7E1 5162 8A... */ + .sqrt2 = QMK(0x1, 0x6A09, 16), /* 1.6A09 E667 F3... */ + .sqrt3 = QMK(0x1, 0xBB67, 16), /* 1.BB67 AE85 84... */ + .ln2 = QMK(0x0, 0xB172, 16), /* 0.B172 17F7 D1... */ + .ln10 = QMK(0x2, 0x4D76, 16), /* 2.4D76 3776 AA... */ + + .version = QVERSION, +}; + +qconf_t qconf = { /* Global Configuration Options */ + .bound = qbound_saturate, + .dp = 4, + .base = 10, +}; + +/********* Basic Library Routines ********************************************/ + + +static inline void implies(const int x, const int y) { + assert(!x || y); +} + +static inline void mutual(const int x, const int y) { /* mutual implication */ + assert(BOOLIFY(x) == BOOLIFY(y)); +} + +static inline void exclusive(const int x, const int y) { + assert(BOOLIFY(x) != BOOLIFY(y)); +} + +static inline void static_assertions(void) { + BUILD_BUG_ON(CHAR_BIT != 8); + // BUILD_BUG_ON((sizeof(q_t)*CHAR_BIT) != (QBITS * 2)); + BUILD_BUG_ON( sizeof(q_t) != sizeof(u_t)); + BUILD_BUG_ON( sizeof(u_t) != sizeof(d_t)); + BUILD_BUG_ON(sizeof(lu_t) != sizeof(ld_t)); + BUILD_BUG_ON(sizeof(d_t) != (sizeof(hd_t) * 2)); + BUILD_BUG_ON(sizeof(lu_t) != (sizeof(u_t) * 2)); +} + +q_t qbound_saturate(const ld_t s) { /**< default saturation handler */ + assert(s > DMAX || s < DMIN); + if (s > DMAX) return DMAX; + return DMIN; +} + +q_t qbound_wrap(const ld_t s) { /**< wrap numbers on overflow */ + assert(s > DMAX || s < DMIN); + if (s > DMAX) return DMIN + (s % DMAX); + return DMAX - ((-s) % DMAX); +} + +static inline q_t qsat(const ld_t s) { + static_assertions(); + if (s > DMAX || s < DMIN) return qconf.bound(s); + return s; +} + +d_t arshift(const d_t v, const unsigned p) { + u_t vn = v; + if (v >= 0l) + return vn >> p; + const u_t leading = ((u_t)(-1l)) << ((sizeof(v) * CHAR_BIT) - p - 1); + return leading | (vn >> p); +} + +static inline d_t divn(const d_t v, const unsigned p) { + /* return v / (1l << p); */ + const u_t shifted = ((u_t)v) >> p; + if (qispositive(v)) + return shifted; + const u_t leading = ((u_t)(-1l)) << ((sizeof(v)*CHAR_BIT) - p - 1); + return leading | shifted; +} + +/* These really all should be moved the header for efficiency reasons */ +static inline u_t qhigh(const q_t q) { return ((u_t)q) >> QBITS; } +static inline u_t qlow(const q_t q) { return ((u_t)q) & QMASK; } +static inline q_t qcons(const u_t hi, const u_t lo) { return (hi << QBITS) | (lo & QMASK); } + +int qtoi(const q_t toi) { return ((lu_t)((ld_t)toi)) >> QBITS; } +q_t qint(const int toq) { return ((u_t)((d_t)toq)) << QBITS; } +signed char qtoc(const q_t q) { return qtoi(q); } +q_t qchar(signed char c) { return qint(c); } +short qtoh(const q_t q) { return qtoi(q); } +q_t qshort(short s) { return qint(s); } +long qtol(const q_t q) { return qtoi(q); } +q_t qlong(long l) { return qint(l); } +long long qtoll(const q_t q) { return qtoi(q); } +q_t qvlong(long long ll) { return qint(ll); } + +q_t qisnegative(const q_t a) { return QINT(BOOLIFY(qhigh(a) & QHIGH)); } +q_t qispositive(const q_t a) { return QINT(!(qhigh(a) & QHIGH)); } +q_t qisinteger(const q_t a) { return QINT(!qlow(a)); } +q_t qisodd(const q_t a) { return QINT(qisinteger(a) && (qhigh(a) & 1)); } +q_t qiseven(const q_t a) { return QINT(qisinteger(a) && !(qhigh(a) & 1)); } +q_t qless(const q_t a, const q_t b) { return QINT(a < b); } +q_t qeqless(const q_t a, const q_t b) { return QINT(a <= b); } +q_t qmore(const q_t a, const q_t b) { return QINT(a > b); } +q_t qeqmore(const q_t a, const q_t b) { return QINT(a >= b); } +q_t qequal(const q_t a, const q_t b) { return QINT(a == b); } +q_t qunequal(const q_t a, const q_t b) { return QINT(a != b); } + +q_t qnegate(const q_t a) { return (~(u_t)a) + 1ULL; } +q_t qmin(const q_t a, const q_t b) { return qless(a, b) ? a : b; } +q_t qmax(const q_t a, const q_t b) { return qmore(a, b) ? a : b; } +q_t qabs(const q_t a) { return qisnegative(a) ? qnegate(a) : a; } +q_t qadd(const q_t a, const q_t b) { return qsat((ld_t)a + (ld_t)b); } +q_t qsub(const q_t a, const q_t b) { return qsat((ld_t)a - (ld_t)b); } +q_t qcopysign(const q_t a, const q_t b) { return qisnegative(b) ? qnegate(qabs(a)) : qabs(a); } +q_t qand(const q_t a, const q_t b) { return a & b; } +q_t qxor(const q_t a, const q_t b) { return a ^ b; } +q_t qor(const q_t a, const q_t b) { return a | b; } +q_t qinvert(const q_t a) { return ~a; } +q_t qnot(const q_t a) { return QINT(!a); } +q_t qlogical(const q_t a) { return QINT(BOOLIFY(a)); } + +q_t qlrs(const q_t a, const q_t b) { /* assert low bits == 0? */ return (u_t)a >> (u_t)qtoi(b); } +q_t qlls(const q_t a, const q_t b) { return (u_t)a << b; } +q_t qars(const q_t a, const q_t b) { return arshift(a, qtoi(b)); } +q_t qals(const q_t a, const q_t b) { return qsat((lu_t)a << b); } +q_t qsign(const q_t a) { return qisnegative(a) ? -QINT(1) : QINT(1); } +q_t qsignum(const q_t a) { return a ? qsign(a) : QINT(0); } + +q_t qapproxequal(const q_t a, const q_t b, const q_t epsilon) { + assert(qeqmore(epsilon, qint(0))); + return QINT(qless(qabs(qsub(a, b)), epsilon)); +} + +q_t qapproxunequal(const q_t a, const q_t b, const q_t epsilon) { + return QINT(!qapproxequal(a, b, epsilon)); +} + +q_t qwithin(q_t v, q_t b1, q_t b2) { + const q_t hi = qmax(b1, b2); + const q_t lo = qmin(b1, b2); + if (qequal(v, b1) || qequal(v, b2)) + return 1; + return qless(v, hi) && qmore(v, lo) ? QINT(1) : QINT(0); +} + +q_t qwithin_interval(q_t v, q_t expected, q_t allowance) { + const q_t b1 = qadd(expected, allowance); + const q_t b2 = qsub(expected, allowance); + return qwithin(v, b1, b2); +} + +q_t qfloor(const q_t q) { + return q & ~QMASK; +} + +q_t qceil(q_t q) { + const q_t adj = qisinteger(q) ? QINT(0) : QINT(1); + q = qadd(q, adj); + return ((u_t)q) & (QMASK << QBITS); +} + +q_t qtrunc(q_t q) { + const q_t adj = qisnegative(q) && qlow(q) ? QINT(1) : QINT(0); + q = qadd(q, adj); + return ((u_t)q) & (QMASK << QBITS); +} + +q_t qround(q_t q) { + const int negative = qisnegative(q); + q = qabs(q); + const q_t adj = (qlow(q) & QHIGH) ? QINT(1) : QINT(0); + q = qadd(q, adj); + q = ((u_t)q) & (QMASK << QBITS); + return negative ? qnegate(q) : q; +} + +int qpack(const q_t *q, char *buffer, const size_t length) { + assert(buffer); + if (length < sizeof(*q)) + return -1; + q_t qn = *q; + uint8_t *b = (uint8_t*)buffer; + for (size_t i = 0; i < sizeof(qn); i++) { + b[i] = qn; + qn = (u_t)qn >> CHAR_BIT; + } + return sizeof(qn); +} + +int qunpack(q_t *q, const char *buffer, const size_t length) { + assert(q); + assert(buffer); + if (length < sizeof(*q)) + return -1; + uint8_t *b = (uint8_t*)buffer; + u_t nq = 0; + for (size_t i = 0; i < sizeof(*q); i++) { + nq <<= CHAR_BIT; + nq |= b[sizeof(*q)-i-1]; + } + *q = nq; + return sizeof(*q); +} + +static inline ld_t multiply(const q_t a, const q_t b) { + const ld_t dd = ((ld_t)a * (ld_t)b) + (lu_t)QHIGH; + /* N.B. portable version of "dd >> QBITS", for double width signed values */ + return dd < 0 ? (-1ull << (2 * QBITS)) | ((lu_t)dd >> QBITS) : ((lu_t)dd) >> QBITS; +} + +q_t qmul(const q_t a, const q_t b) { + return qsat(multiply(a, b)); +} + +q_t qfma(const q_t a, const q_t b, const q_t c) { + return qsat(multiply(a, b) + (ld_t)c); +} + +q_t qdiv(const q_t a, const q_t b) { + assert(b); + const ld_t dd = ((ld_t)a) << QBITS; + ld_t bd2 = divn(b, 1); + if (!((dd >= 0 && b > 0) || (dd < 0 && b < 0))) + bd2 = -bd2; + /* Overflow not checked! */ + /*return (dd/b) + (bd2/b);*/ + return (dd + bd2) / b; +} + +q_t qrem(const q_t a, const q_t b) { + return qsub(a, qmul(qtrunc(qdiv(a, b)), b)); +} + +q_t qmod(q_t a, q_t b) { + return qsub(a, qmul(qfloor(qdiv(a, b)), b)); +} + +static char itoch(const unsigned ch) { + assert(ch < 36); + if (ch <= 9) + return ch + '0'; + return ch + 'A' - 10; +} + +static inline void swap(char *a, char *b) { + assert(a); + assert(b); + const int c = *a; + *a = *b; + *b = c; +} + +static void reverse(char *s, const size_t length) { + assert(s); + for (size_t i = 0; i < length/2; i++) + swap(&s[i], &s[length - i - 1]); +} + +static int uprint(u_t p, char *s, const size_t length, const d_t base) { + assert(s); + assert(base >= 2 && base <= 36); + if (length < 2) + return -1; + size_t i = 0; + do { + unsigned ch = p % base; + p /= base; + s[i++] = itoch(ch); + } while (p && i < length); + if (p && i >= length) + return -1; + reverse(s, i); + return i; +} + +/* <https://codereview.stackexchange.com/questions/109212> */ +int qsprintbdp(q_t p, char *s, size_t length, const u_t base, const d_t idp) { + assert(s); + const int negative = BOOLIFY(qisnegative(p)); + if (negative) + p = qnegate(p); + const d_t hi = qhigh(p); + char frac[QBITS + 2] = { '.', }; + memset(s, 0, length); + assert(base >= 2 && base <= 36); + u_t lo = qlow(p); + size_t i = 1; + for (i = 1; lo; i++) { + if (idp >= 0 && (int)i > idp) + break; + lo *= base; + assert(i < (QBITS + 2)); + frac[i] = itoch(lo >> QBITS); + lo &= QMASK; + } + if (negative) + s[0] = '-'; + const int hisz = uprint(hi, s + negative, length - (1 + negative), base); + if (hisz < 0 || (hisz + i + negative + 1) > length) + return -1; + memcpy(s + hisz + negative, frac, i); + return i + hisz; +} + +int qsprintb(q_t p, char *s, size_t length, const u_t base) { + return qsprintbdp(p, s, length, base, qconf.dp); +} + +int qsprint(const q_t p, char *s, const size_t length) { + return qsprintb(p, s, length, qconf.base); +} + +static inline int extract(unsigned char c, const int radix) { + c = tolower(c); + if (c >= '0' && c <= '9') + c -= '0'; + else if (c >= 'a' && c <= 'z') + c -= ('a' - 10); + else + return -1; + if (c < radix) + return c; + return -1; +} + +static inline q_t qmk(d_t integer, u_t fractional) { + const int negative = integer < 0; + integer = negative ? -integer : integer; + const q_t r = qcons((d_t)integer, fractional); + return negative ? qnegate(r) : r; +} + +static inline u_t integer_logarithm(u_t num, const u_t base) { + assert(num > 0 && base >= 2 && base <= 36); + u_t r = -1; + do r++; while (num /= base); + return r; +} + +int qnconvbdp(q_t *q, const char *s, size_t length, const d_t base, const u_t idp) { + assert(q); + assert(s); + assert(base >= 2 && base <= 36); + *q = QINT(0); + if (length < 1) + return -1; + d_t hi = 0, lo = 0, places = 1, negative = 0, overflow = 0; + size_t sidx = 0; + + if (s[sidx] == '-') { + if (length < 2) + return -1; + negative = 1; + sidx++; + } + + for (; sidx < length && s[sidx]; sidx++) { + const d_t e = extract(s[sidx], base); + if (e < 0) + break; + if (hi > MULTIPLIER) { /* continue on with conversion, do not accumulate */ + overflow = 1; + } else { + hi = (hi * base) + e; + } + } + if (sidx >= length || !s[sidx]) + goto done; + if (s[sidx] != '.') + return -2; + sidx++; + + const u_t ilog = integer_logarithm(0x10000, base); + const u_t max = MIN(idp, ilog); /* Calculate maximum decimal places given base */ + + for (u_t dp = 0; sidx < length && s[sidx]; sidx++, dp++) { + const int ch = extract(s[sidx], base); + if (ch < 0) + return -3; + if (dp < max) { /* continue on with conversion , do not accumulate */ + /* We could get more accuracy by looking at one digit + * passed the maximum digits allowed and rounding if + * that digit exists in the input. */ + lo = (lo * base) + ch; + if (places >= (DMAX / base)) + return -4; + places *= base; + } + assert((dp + 1) > dp); + } + if (!places) + return -5; + lo = ((d_t)((u_t)lo << QBITS) / places); +done: + if (overflow) { + *q = negative ? qinfo.min : qinfo.max; + return -6; + } else { + const q_t nq = qmk(hi, lo); + *q = negative ? qnegate(nq) : nq; + + } + return 0; +} + +int qnconvb(q_t *q, const char *s, size_t length, const d_t base) { + return qnconvbdp(q, s, length, base, qconf.dp); +} + +int qnconv(q_t *q, const char *s, size_t length) { + return qnconvb(q, s, length, qconf.base); +} + +int qconv(q_t *q, const char * const s) { + assert(s); + return qnconv(q, s, strlen(s)); +} + +int qconvb(q_t *q, const char * const s, const d_t base) { + assert(s); + return qnconvb(q, s, strlen(s), base); +} + +typedef enum { + CORDIC_MODE_VECTOR_E/* = 'VECT'*/, + CORDIC_MODE_ROTATE_E/* = 'ROT'*/, +} cordic_mode_e; + +typedef enum { + CORDIC_COORD_HYPERBOLIC_E = -1, + CORDIC_COORD_LINEAR_E = 0, + CORDIC_COORD_CIRCULAR_E = 1, +} cordic_coordinates_e; + +static const d_t cordic_circular_inverse_scaling = 0x9B74; /* 1/scaling-factor */ +static const d_t cordic_hyperbolic_inverse_scaling = 0x13520; /* 1/scaling-factor */ + +static inline int mulsign(d_t a, d_t b) { /* sign(a*b) */ + const int aneg = a < 0; + const int bneg = b < 0; + return aneg ^ bneg ? -QINT(1) : QINT(1); +} + +/* Universal CORDIC <https://en.wikibooks.org/wiki/Digital_Circuits/CORDIC> + * + * x(i+1) = x(i) - u.d(i).y(i).pow(2, -i) + * y(i+1) = y(i) + d(i).x(i).pow(2, -i) + * z(i+1) = z(i) - d(i).a(i) + * + * d(i) = sgn(z(i)) (rotation) + * d(i) = -sgn(x(i).y(i)) (vectoring) + * + * hyperbolic linear circular + * u = -1 0 1 + * a = atanh(pow(2, -i)) pow(2, -i) atan(pow(2, -i)) + * + * linear shift sequence: i = 0, 1, 2, 3, ... + * circular shift sequence: i = 1, 2, 3, 4, ... + * hyperbolic shift sequence: i = 1, 2, 3, 4, 4, 5, ... */ +static int cordic(const cordic_coordinates_e coord, const cordic_mode_e mode, int iterations, d_t *x0, d_t *y0, d_t *z0) { + assert(x0); + assert(y0); + assert(z0); + if (mode != CORDIC_MODE_VECTOR_E && mode != CORDIC_MODE_ROTATE_E) + return -1; + + BUILD_BUG_ON(sizeof(d_t) != sizeof(uint32_t)); + BUILD_BUG_ON(sizeof(u_t) != sizeof(uint32_t)); + + static const u_t arctans[] = { /* atan(2^0), atan(2^-1), atan(2^-2), ... */ + 0xC90FuL, 0x76B1uL, 0x3EB6uL, 0x1FD5uL, + 0x0FFAuL, 0x07FFuL, 0x03FFuL, 0x01FFuL, + 0x00FFuL, 0x007FuL, 0x003FuL, 0x001FuL, + 0x000FuL, 0x0007uL, 0x0003uL, 0x0001uL, + 0x0000uL, // 0x0000uL, + }; + static const size_t arctans_length = sizeof arctans / sizeof arctans[0]; + + static const u_t arctanhs[] = { /* atanh(2^-1), atanh(2^-2), ... */ + 0x8c9fuL, 0x4162uL, 0x202buL, 0x1005uL, + 0x0800uL, 0x0400uL, 0x0200uL, 0x0100uL, + 0x0080uL, 0x0040uL, 0x0020uL, 0x0010uL, + 0x0008uL, 0x0004uL, 0x0002uL, 0x0001uL, + 0x0000uL, // 0x0000uL, + }; + static const size_t arctanhs_length = sizeof arctanhs / sizeof arctanhs[0]; + + static const u_t halfs[] = { /* 2^0, 2^-1, 2^-2, ..*/ + 0x10000uL, + 0x8000uL, 0x4000uL, 0x2000uL, 0x1000uL, + 0x0800uL, 0x0400uL, 0x0200uL, 0x0100uL, + 0x0080uL, 0x0040uL, 0x0020uL, 0x0010uL, + 0x0008uL, 0x0004uL, 0x0002uL, 0x0001uL, + //0x0000uL, // 0x0000uL, + }; + static const size_t halfs_length = sizeof halfs / sizeof halfs[0]; + + const u_t *lookup = NULL; + size_t i = 0, j = 0, k = 0, length = 0; + const size_t *shiftx = NULL, *shifty = NULL; + int hyperbolic = 0; + + switch (coord) { + case CORDIC_COORD_CIRCULAR_E: + lookup = arctans; + length = arctans_length; + i = 0; + shifty = &i; + shiftx = &i; + break; + case CORDIC_COORD_HYPERBOLIC_E: + lookup = arctanhs; + length = arctanhs_length; + hyperbolic = 1; + i = 1; + shifty = &i; + shiftx = &i; + break; + case CORDIC_COORD_LINEAR_E: + lookup = halfs; + length = halfs_length; + shifty = &j; + shiftx = NULL; + i = 1; + break; + default: /* not implemented */ + return -2; + } + + iterations = iterations > (int)length ? (int)length : iterations; + iterations = iterations < 0 ? (int)length : iterations; + + d_t x = *x0, y = *y0, z = *z0; + + /* rotation mode: z determines direction, + * vector mode: y determines direction */ + for (; j < (unsigned)iterations; i++, j++) { + again: + { + const d_t m = mode == CORDIC_MODE_ROTATE_E ? z : -y /*-mulsign(x, y)*/; + const d_t d = -!!(m < 0); + const d_t xs = ((((shiftx ? divn(y, *shiftx) : 0)) ^ d) - d); + const d_t ys = (divn(x, *shifty) ^ d) - d; + const d_t xn = x - (hyperbolic ? -xs : xs); + const d_t yn = y + ys; + const d_t zn = z - ((lookup[j] ^ d) - d); + x = xn; /* cosine, in circular, rotation mode */ + y = yn; /* sine, in circular, rotation mode */ + z = zn; + } + if (hyperbolic) { /* Experimental/Needs bug fixing */ + switch (1) { // TODO: Correct hyperbolic redo of iteration + case 0: break; + case 1: if (k++ >= 3) { k = 0; goto again; } break; + case 2: { + assert(j <= 120); + size_t cmp = j + 1; + if (cmp == 4 || cmp == 13 /*|| cmp == 40 || cmp == 121 || cmp == floor(pow(3,i-1)/2) */) { + if (k) { + k = 0; + } else { + k = 1; + goto again; + } + } + break; + } + } + } + } + *x0 = x; + *y0 = y; + *z0 = z; + + return iterations; +} + +/* See: - <https://dspguru.com/dsp/faqs/cordic/> + * - <https://en.wikipedia.org/wiki/CORDIC> */ +static int qcordic(q_t theta, const int iterations, q_t *sine, q_t *cosine) { + assert(sine); + assert(cosine); + + static const q_t pi = QPI, npi = -QPI; + static const q_t hpi = QPI/2, hnpi = -(QPI/2); + static const q_t qpi = QPI/4, qnpi = -(QPI/4); + static const q_t dpi = QPI*2, dnpi = -(QPI*2); + + /* Convert to range -pi to pi, we could use qmod, + * however that uses multiplication and division, and + * if we can use those operators freely then there are + * other, better algorithms we can use instead of CORDIC + * for sine/cosine calculation. */ + while (qless(theta, npi)) theta = qadd(theta, dpi); + while (qmore(theta, pi)) theta = qadd(theta, dnpi); + + int negate = 0, shift = 0; + + /* convert to range -pi/2 to pi/2 */ + if (qless(theta, hnpi)) { + theta = qadd(theta, pi); + negate = 1; + } else if (qmore(theta, hpi)) { + theta = qadd(theta, npi); + negate = 1; + } + + /* convert to range -pi/4 to pi/4 */ + if (qless(theta, qnpi)) { + theta = qadd(theta, hpi); + shift = -1; + } else if (qmore(theta, qpi)) { + theta = qadd(theta, hnpi); + shift = 1; + } + + d_t x = cordic_circular_inverse_scaling, y = 0, z = theta /* no theta scaling needed */; + + /* CORDIC in Q2.16 format */ + if (cordic(CORDIC_COORD_CIRCULAR_E, CORDIC_MODE_ROTATE_E, iterations, &x, &y, &z) < 0) + return -1; + + /* undo shifting and quadrant changes */ + if (shift > 0) { + const d_t yt = y; + y = x; + x = -yt; + } else if (shift < 0) { + const d_t yt = y; + y = -x; + x = yt; + } + + if (negate) { + x = -x; + y = -y; + } + /* set output; no scaling needed */ + *cosine = x; + *sine = y; + return 0; +} + +q_t qatan(const q_t t) { + q_t x = qint(1), y = t, z = QINT(0); + cordic(CORDIC_COORD_CIRCULAR_E, CORDIC_MODE_VECTOR_E, -1, &x, &y, &z); + return z; +} + +q_t qatan2(const q_t a, const q_t b) { + q_t x = b, y = a, z = QINT(0); + if (qequal(b, QINT(0))) { + assert(qunequal(a, QINT(0))); + if (qmore(a, QINT(0))) + return QPI/2; + return -(QPI/2); + } else if (qless(b, QINT(0))) { + if (qeqmore(a, QINT(0))) + return qadd(qatan(qdiv(a, b)), QPI); + return qsub(qatan(qdiv(a, b)), QPI); + } + cordic(CORDIC_COORD_CIRCULAR_E, CORDIC_MODE_VECTOR_E, -1, &x, &y, &z); + return z; +} + +void qsincos(q_t theta, q_t *sine, q_t *cosine) { + assert(sine); + assert(cosine); + const int r = qcordic(theta, -1, sine, cosine); + assert(r >= 0); +} + +q_t qsin(const q_t theta) { + q_t sine = QINT(0), cosine = QINT(0); + qsincos(theta, &sine, &cosine); + return sine; +} + +q_t qcos(const q_t theta) { + q_t sine = QINT(0), cosine = QINT(0); + qsincos(theta, &sine, &cosine); + return cosine; +} + +q_t qtan(const q_t theta) { + q_t sine = QINT(0), cosine = QINT(0); + qsincos(theta, &sine, &cosine); + return qdiv(sine, cosine); /* can use qcordic_div, with range limits it imposes */ +} + +q_t qcot(const q_t theta) { + q_t sine = QINT(0), cosine = QINT(0); + qsincos(theta, &sine, &cosine); + return qdiv(cosine, sine); /* can use qcordic_div, with range limits it imposes */ +} + +q_t qcordic_mul(const q_t a, const q_t b) { /* works for small values; result < 4 */ + q_t x = a, y = QINT(0), z = b; + const int r = cordic(CORDIC_COORD_LINEAR_E, CORDIC_MODE_ROTATE_E, -1, &x, &y, &z); + assert(r >= 0); + return y; +} + +q_t qcordic_div(const q_t a, const q_t b) { + q_t x = b, y = a, z = QINT(0); + const int r = cordic(CORDIC_COORD_LINEAR_E, CORDIC_MODE_VECTOR_E, -1, &x, &y, &z); + assert(r >= 0); + return z; +} + +void qsincosh(const q_t a, q_t *sinh, q_t *cosh) { + assert(sinh); + assert(cosh); + q_t x = cordic_hyperbolic_inverse_scaling, y = QINT(0), z = a; /* (e^2x - 1) / (e^2x + 1) */ + const int r = cordic(CORDIC_COORD_HYPERBOLIC_E, CORDIC_MODE_ROTATE_E, -1, &x, &y, &z); + assert(r >= 0); + *sinh = y; + *cosh = x; +} + +q_t qtanh(const q_t a) { + q_t sinh = QINT(0), cosh = QINT(0); + qsincosh(a, &sinh, &cosh); + return qdiv(sinh, cosh); +} + +q_t qcosh(const q_t a) { + q_t sinh = QINT(0), cosh = QINT(0); + qsincosh(a, &sinh, &cosh); + return cosh; +} + +q_t qsinh(const q_t a) { + q_t sinh = QINT(0), cosh = QINT(0); + qsincosh(a, &sinh, &cosh); + return sinh; +} + +q_t qcordic_exp(const q_t e) { + q_t s = QINT(0), h = QINT(0); + qsincosh(e, &s, &h); + return qadd(s, h); +} + +q_t qcordic_ln(const q_t d) { + q_t x = qadd(d, QINT(1)), y = qsub(d, QINT(1)), z = QINT(0); + const int r = cordic(CORDIC_COORD_HYPERBOLIC_E, CORDIC_MODE_VECTOR_E, -1, &x, &y, &z); + assert(r >= 0); + return qadd(z, z); +} + +q_t qcordic_sqrt(const q_t n) { /* testing only; works for 0 < x < 2 */ + const q_t quarter = 1uLL << (QBITS - 2); /* 0.25 */ + q_t x = qadd(n, quarter), + y = qsub(n, quarter), + z = 0; + const int r = cordic(CORDIC_COORD_HYPERBOLIC_E, CORDIC_MODE_VECTOR_E, -1, &x, &y, &z); + assert(r >= 0); + return qmul(x, cordic_hyperbolic_inverse_scaling); +} + +q_t qhypot(const q_t a, const q_t b) { + q_t x = qabs(a), y = qabs(b), z = QINT(0); /* abs() should not be needed? */ + const int r = cordic(CORDIC_COORD_CIRCULAR_E, CORDIC_MODE_VECTOR_E, -1, &x, &y, &z); + assert(r >= 0); + return qmul(x, cordic_circular_inverse_scaling); +} + +q_t qatanh(q_t x) { + assert(qabs(qless(x, QINT(1)))); + return qmul(qlog(qdiv(qadd(QINT(1), x), qsub(QINT(1), x))), QMK(0, 0x8000, 16)); +} + +q_t qasinh(q_t x) { + return qlog(qadd(x, qsqrt(qadd(qmul(x, x), QINT(1))))); +} + +q_t qacosh(q_t x) { + assert(qeqmore(x, QINT(1))); + return qlog(qadd(x, qsqrt(qsub(qmul(x, x), QINT(1))))); +} + +void qpol2rec(const q_t magnitude, const q_t theta, q_t *i, q_t *j) { + assert(i); + assert(j); + q_t sin = QINT(0), cos = QINT(0); + qsincos(theta, &sin, &cos); + *i = qmul(sin, magnitude); + *j = qmul(cos, magnitude); +} + +void qrec2pol(const q_t i, const q_t j, q_t *magnitude, q_t *theta) { + assert(magnitude); + assert(theta); + const int is = qisnegative(i), js = qisnegative(j); + q_t x = qabs(i), y = qabs(j), z = QINT(0); + const int r = cordic(CORDIC_COORD_CIRCULAR_E, CORDIC_MODE_VECTOR_E, -1, &x, &y, &z); + assert(r >= 0); + *magnitude = qmul(x, cordic_circular_inverse_scaling); + if (is && js) + z = qadd(z, QPI); + else if (js) + z = qadd(z, QPI/2l); + else if (is) + z = qadd(z, (3l*QPI)/2l); + *theta = z; +} + +q_t qcordic_hyperbolic_gain(const int n) { + q_t x = QINT(1), y = QINT(0), z = QINT(0); + const int r = cordic(CORDIC_COORD_HYPERBOLIC_E, CORDIC_MODE_ROTATE_E, n, &x, &y, &z); + assert(r >= 0); + return x; +} + +q_t qcordic_circular_gain(const int n) { + q_t x = QINT(1), y = QINT(0), z = QINT(0); + const int r = cordic(CORDIC_COORD_CIRCULAR_E, CORDIC_MODE_ROTATE_E, n, &x, &y, &z); + assert(r >= 0); + return x; +} + +static inline int isodd(const unsigned n) { + return n & 1; +} + +d_t dpower(d_t b, unsigned e) { /* https://stackoverflow.com/questions/101439 */ + d_t result = 1; + for (;;) { + if (isodd(e)) + result *= b; + e >>= 1; + if (!e) + break; + b *= b; + } + return result; +} + +d_t dlog(d_t x, const unsigned base) { /* rounds up, look at remainder to round down */ + d_t b = 0; + assert(x && base > 1); + while ((x /= (d_t)base)) /* can use >> for base that are powers of two */ + b++; + return b; +} + +q_t qlog(q_t x) { + q_t logs = 0; + assert(qmore(x, 0)); + static const q_t lmax = QMK(9, 0x8000, 16); /* 9.5, lower limit needs checking */ + for (; qmore(x, lmax); x = divn(x, 1)) + logs = qadd(logs, qinfo.ln2); + return qadd(logs, qcordic_ln(x)); +} + +q_t qsqr(const q_t x) { + return qmul(x, x); +} + +q_t qexp(const q_t e) { /* exp(e) = exp(e/2)*exp(e/2) */ + if (qless(e, QINT(1))) /* 1.1268 is approximately the limit for qcordic_exp */ + return qcordic_exp(e); + return qsqr(qexp(divn(e, 1))); +} + +q_t qpow(q_t n, q_t exp) { + implies(qisnegative(n), qisinteger(exp)); + implies(qequal(n, QINT(0)), qunequal(exp, QINT(0))); + if (qequal(QINT(0), n)) + return QINT(1); + if (qisnegative(n)) { + const q_t abspow = qpow(qabs(n), exp); + return qisodd(exp) ? qnegate(abspow) : abspow; + } + if (qisnegative(exp)) + return qdiv(QINT(1), qpow(n, qabs(exp))); + return qexp(multiply(qlog(n), exp)); +} + +q_t qsqrt(const q_t x) { /* Newton-Rhaphson method */ + assert(qeqmore(x, 0)); + const q_t difference = qmore(x, QINT(100)) ? 0x0100 : 0x0010; + if (qequal(QINT(0), x)) + return QINT(0); + q_t guess = qmore(x, qinfo.sqrt2) ? divn(x, 1) : QINT(1); + while (qmore(qabs(qsub(qmul(guess, guess), x)), difference)) + guess = divn(qadd(qdiv(x, guess), guess), 1); + return qabs(guess); /* correct for overflow int very large numbers */ +} + +q_t qasin(const q_t t) { + assert(qless(qabs(t), QINT(1))); + /* can also use: return qatan(qdiv(t, qsqrt(qsub(QINT(1), qmul(t, t))))); */ + return qatan2(t, qsqrt(qsub(QINT(1), qmul(t, t)))); +} + +q_t qacos(const q_t t) { + assert(qeqless(qabs(t), QINT(1))); + /* can also use: return qatan(qdiv(qsqrt(qsub(QINT(1), qmul(t, t))), t)); */ + return qatan2(qsqrt(qsub(QINT(1), qmul(t, t))), t); +} + +q_t qdeg2rad(const q_t deg) { + return qdiv(qmul(QPI, deg), QINT(180)); +} + +q_t qrad2deg(const q_t rad) { + return qdiv(qmul(QINT(180), rad), QPI); +} + +void qfilter_init(qfilter_t *f, const q_t time, const q_t rc, const q_t seed) { + assert(f); + memset(f, 0, sizeof(*f)); + f->time = time; + f->rc = rc; + f->filtered = seed; /* alpha * seed for LPF */ + f->raw = seed; +} + +q_t qfilter_low_pass(qfilter_t *f, const q_t time, const q_t data) { + assert(f); + /* If the calling rate is constant (for example the function is + * guaranteed to be always called at a rate of 5 milliseconds) we + * can avoid the costly alpha calculation! */ + const q_t dt = (u_t)time - (u_t)f->time; + const q_t alpha = qdiv(dt, qadd(f->rc, dt)); + f->filtered = qfma(alpha, qsub(data, f->filtered), f->filtered); + f->time = time; + f->raw = data; + return f->filtered; +} + +q_t qfilter_high_pass(qfilter_t *f, const q_t time, const q_t data) { + assert(f); + const q_t dt = (u_t)time - (u_t)f->time; + const q_t alpha = qdiv(f->rc, qadd(f->rc, dt)); + f->filtered = qmul(alpha, qadd(f->filtered, qsub(data, f->raw))); + f->time = time; + f->raw = data; + return f->filtered; +} + +q_t qfilter_value(const qfilter_t *f) { + assert(f); + return f->filtered; +} + +/* Must be called at a constant rate; perhaps a PID which takes call time + * into account could be made, but that would complicate things. Differentiator + * term needs filtering also. It would be nice to create a version that took + * into account the time delta, see + * <https://www.quora.com/Do-I-need-to-sample-at-a-constant-rate-for-PID-control-or-is-it-sufficient-to-know-the-time-at-which-my-sample-was-taken-even-if-the-increment-varies> + * */ +q_t qpid_update(qpid_t *pid, const q_t error, const q_t position) { + assert(pid); + const q_t p = qmul(pid->p_gain, error); + pid->i_state = qadd(pid->i_state, error); + pid->i_state = qmax(pid->i_state, pid->i_min); + pid->i_state = qmin(pid->i_state, pid->i_max); + const q_t i = qmul(pid->i_state, pid->i_gain); + const q_t d = qmul(pid->d_gain, qsub(position, pid->d_state)); + pid->d_state = position; + return qsub(qadd(p, i), d); +} + +/* Simpsons method for numerical integration, from "Math Toolkit for + * Real-Time Programming" by Jack Crenshaw */ +q_t qsimpson(q_t (*f)(q_t), const q_t x1, const q_t x2, const unsigned n) { + assert(f); + assert((n & 1) == 0); + const q_t h = qdiv(qsub(x2, x1), QINT(n)); + q_t sum = 0, x = x1; + for (unsigned i = 0; i < (n / 2u); i++){ + sum = qadd(sum, qadd(f(x), qmul(QINT(2), f(qadd(x,h))))); + x = qadd(x, qmul(QINT(2), h)); + } + sum = qsub(qmul(QINT(2), sum), qadd(f(x1), f(x2))); + return qdiv(qmul(h, sum), QINT(3)); +} + +/* The matrix meta-data field is not used at the moment, but could be + * used for things like versioning, determining whether the matrix is + * all zeros, or is the identify matrix, whether it contains valid data, + * and more. Some common matrix operations are missing, such as factorization + * + * A function for image kernels might be useful. */ + +enum { METADATA, LENGTH, ROW, COLUMN, DATA, }; + +int qmatrix_is_valid(const q_t *m) { + const size_t size = m[LENGTH], row = m[ROW], column = m[COLUMN]; + const size_t elements = row * column; + if (elements < row || elements < column) /* overflow */ + return 0; + if (elements > size) + return 0; + return 1; +} + +int qmatrix_resize(q_t *m, const size_t row, const size_t column) { + const size_t rc = row * column; + const size_t sz = m[LENGTH]; + if ((row && column) && (rc < row || rc < column)) /* overflow */ + return -1; + if (rc > sz) + return -1; + m[ROW] = row; + m[COLUMN] = column; + return 0; +} + +int qmatrix_apply_unary(q_t *r, const q_t *a, q_t (*func)(q_t)) { + assert(r); + assert(qmatrix_is_valid(r)); + assert(a); + assert(qmatrix_is_valid(a)); + assert(func); + const q_t *ma = &a[DATA]; + q_t *mr = &r[DATA]; + const size_t arows = a[ROW], acolumns = a[COLUMN]; + if (qmatrix_resize(r, arows, acolumns) < 0) + return -1; + for (size_t i = 0; i < arows; i++) + for (size_t j = 0; j < acolumns; j++) + mr[i*acolumns + j] = func(ma[i*acolumns + j]); + return 0; +} + +int qmatrix_apply_scalar(q_t *r, const q_t *a, q_t (*func)(q_t, q_t), const q_t c) { + assert(r); + assert(qmatrix_is_valid(r)); + assert(a); + assert(qmatrix_is_valid(a)); + assert(func); + const q_t *ma = &a[DATA]; + q_t *mr = &r[DATA]; + const size_t arows = a[ROW], acolumns = a[COLUMN]; + if (qmatrix_resize(r, arows, acolumns) < 0) + return -1; + for (size_t i = 0; i < arows; i++) + for (size_t j = 0; j < acolumns; j++) + mr[i*acolumns + j] = func(ma[i*acolumns + j], c); + return 0; +} + +int qmatrix_apply_binary(q_t *r, const q_t *a, const q_t *b, q_t (*func)(q_t, q_t)) { + assert(a); + assert(qmatrix_is_valid(a)); + assert(b); + assert(qmatrix_is_valid(b)); + assert(r); + assert(qmatrix_is_valid(r)); + assert(func); + const q_t *ma = &a[DATA], *mb = &b[DATA]; + q_t *mr = &r[DATA]; + const size_t arows = a[ROW], acolumns = a[COLUMN]; + const size_t brows = b[ROW], bcolumns = b[COLUMN]; + const size_t rrows = r[ROW], rcolumns = r[COLUMN]; + if (arows != brows || acolumns != bcolumns) + return -1; + if (arows != rrows || acolumns != rcolumns) + return -1; + for (size_t i = 0; i < arows; i++) + for (size_t j = 0; j < acolumns; j++) { + const size_t idx = (i*acolumns) + j; + mr[idx] = func(ma[idx], mb[idx]); + } + return 0; +} + +static q_t qfz(q_t a) { UNUSED(a); return QINT(0); } +static q_t qf1(q_t a) { UNUSED(a); return QINT(1); } + +int qmatrix_zero(q_t *r) { return qmatrix_apply_unary(r, r, qfz); } +int qmatrix_one(q_t *r) { return qmatrix_apply_unary(r, r, qf1); } +int qmatrix_logical(q_t *r, const q_t *a) { return qmatrix_apply_unary(r, a, qlogical); } +int qmatrix_not(q_t *r, const q_t *a) { return qmatrix_apply_unary(r, a, qnot); } +int qmatrix_signum(q_t *r, const q_t *a) { return qmatrix_apply_unary(r, a, qsignum); } +int qmatrix_invert(q_t *r, const q_t *a) { return qmatrix_apply_unary(r, a, qinvert); } +int qmatrix_add(q_t *r, const q_t *a, const q_t *b) { return qmatrix_apply_binary(r, a, b, qadd); } +int qmatrix_sub(q_t *r, const q_t *a, const q_t *b) { return qmatrix_apply_binary(r, a, b, qsub); } +int qmatrix_and(q_t *r, const q_t *a, const q_t *b) { return qmatrix_apply_binary(r, a, b, qand); } +int qmatrix_or (q_t *r, const q_t *a, const q_t *b) { return qmatrix_apply_binary(r, a, b, qor); } +int qmatrix_xor(q_t *r, const q_t *a, const q_t *b) { return qmatrix_apply_binary(r, a, b, qxor); } + +int qmatrix_scalar_add(q_t *r, const q_t *a, const q_t scalar) { return qmatrix_apply_scalar(r, a, qadd, scalar); } +int qmatrix_scalar_sub(q_t *r, const q_t *a, const q_t scalar) { return qmatrix_apply_scalar(r, a, qsub, scalar); } +int qmatrix_scalar_mul(q_t *r, const q_t *a, const q_t scalar) { return qmatrix_apply_scalar(r, a, qmul, scalar); } +int qmatrix_scalar_div(q_t *r, const q_t *a, const q_t scalar) { return qmatrix_apply_scalar(r, a, qdiv, scalar); } +int qmatrix_scalar_mod(q_t *r, const q_t *a, const q_t scalar) { return qmatrix_apply_scalar(r, a, qmod, scalar); } +int qmatrix_scalar_rem(q_t *r, const q_t *a, const q_t scalar) { return qmatrix_apply_scalar(r, a, qrem, scalar); } +int qmatrix_scalar_and(q_t *r, const q_t *a, const q_t scalar) { return qmatrix_apply_scalar(r, a, qand, scalar); } +int qmatrix_scalar_or (q_t *r, const q_t *a, const q_t scalar) { return qmatrix_apply_scalar(r, a, qor, scalar); } +int qmatrix_scalar_xor(q_t *r, const q_t *a, const q_t scalar) { return qmatrix_apply_scalar(r, a, qxor, scalar); } + +int qmatrix_is_square(const q_t *m) { + assert(m); + assert(qmatrix_is_valid(m)); + return m[COLUMN] == m[ROW]; +} + +int qmatrix_identity(q_t *r) { + assert(r); + assert(qmatrix_is_valid(r)); + if (!qmatrix_is_square(r)) + return -1; + q_t *mr = &r[DATA]; + const size_t length = r[ROW]; + for (size_t i = 0; i < length; i++) + for (size_t j = 0; j < length; j++) + mr[i*length + j] = i == j ? QINT(1) : QINT(0); + return 0; +} + +int qmatrix_copy(q_t *r, const q_t *a) { + assert(r); + assert(qmatrix_is_valid(r)); + assert(a); + assert(qmatrix_is_valid(a)); + const size_t arows = a[ROW], acolumns = a[COLUMN]; + const size_t copy = arows * acolumns * sizeof (q_t); + if ((arows && acolumns) && (copy < arows || copy < acolumns)) + return -1; + if (qmatrix_resize(r, arows, acolumns) < 0) + return -1; + memcpy(&r[DATA], &a[DATA], copy); + return 0; +} + +q_t qmatrix_trace(const q_t *m) { + assert(m); + assert(qmatrix_is_square(m)); + const size_t length = m[ROW]; + const q_t *mm = &m[DATA]; + q_t tr = QINT(0); + for (size_t i = 0; i < length; i++) + for (size_t j = 0; j < length; j++) + if (i == j) + tr = qadd(tr, mm[i*length + j]); + return tr; +} + +q_t qmatrix_equal(const q_t *a, const q_t *b) { + assert(a); + assert(qmatrix_is_valid(a)); + assert(b); + assert(qmatrix_is_valid(b)); + const size_t arow = a[ROW], acolumn = a[COLUMN]; + const size_t brow = b[ROW], bcolumn = b[COLUMN]; + const q_t *ma = &a[DATA]; + const q_t *mb = &a[DATA]; + if (a == b) + return QINT(1); + if (arow != brow && acolumn != bcolumn) + return QINT(0); + return !memcmp(ma, mb, sizeof(q_t) * arow * brow); +} + +static q_t determine(const q_t *m, const size_t length) { + assert(m); + if (length == 1) + return m[0]; + if (length == 2) + return qsub(qmul(m[0], m[3]), qmul(m[1], m[2])); + size_t co1 = 0, co2 = 0; + q_t det = QINT(0), sgn = QINT(1); + q_t co[length*length]; /* This should really be passed in */ + for (size_t i = 0; i < length; i++) { + for (size_t j = 0; j < length; j++) + for (size_t k = 0; k < length; k++) + if (j && k != i) { + co[co1*length + co2] = m[j*length + k]; + if (++co2 > (length - 2)) { + co1++; + co2 = 0; + } + } + det = qadd(det, qcopysign(qmul(m[(0*length) + i], determine(co, length - 1)), sgn)); + sgn = qnegate(sgn); + } + return det; +} + +q_t qmatrix_determinant(const q_t *m) { + assert(m); + assert(qmatrix_is_square(m)); + assert(m[ROW] < 16); + const size_t length = m[ROW]; + const q_t *mm = &m[DATA]; + return determine(mm, length); +} + +int qmatrix_transpose(q_t *r, const q_t *m) { + assert(r); + assert(qmatrix_is_valid(r)); + assert(m); + assert(qmatrix_is_valid(m)); + q_t *mr = &r[DATA]; + const q_t *mm = &m[DATA]; + const size_t mrows = m[ROW], mcolumns = m[COLUMN]; + const size_t msize = mrows * mcolumns; + const size_t rsize = r[LENGTH]; + if (msize > rsize) + return -1; + for (size_t i = 0; i < mrows; i++) + for (size_t j = 0; j < mcolumns; j++) + mr[i*mcolumns + j] = mm[j*mcolumns + i]; + r[ROW] = mcolumns; + r[COLUMN] = mrows; + return 0; +} + +int qmatrix_mul(q_t *r, const q_t *a, const q_t *b) { + assert(a); + assert(qmatrix_is_valid(a)); + assert(b); + assert(qmatrix_is_valid(b)); + assert(r); + assert(qmatrix_is_valid(r)); + q_t *mr = &r[DATA]; + const q_t *ma = &a[DATA], *mb = &b[DATA]; + const size_t arows = a[ROW], acolumns = a[COLUMN]; + const size_t brows = b[ROW], bcolumns = b[COLUMN]; + if (acolumns != brows) + return -1; + if (qmatrix_resize(r, arows, bcolumns) < 0) + return -1; + for (size_t i = 0; i < arows; i++) + for (size_t j = 0; j < bcolumns; j++) { + q_t s = QINT(0); + for (size_t k = 0; k < brows; k++) + s = qadd(s, qmul(ma[i*acolumns + k], mb[k*bcolumns + j])); + mr[i*arows + j] = s; + } + return 0; +} + +static int addchar(char **str, size_t *length, const int ch) { + assert(str && *str); + assert(length); + if (!length) + return -1; + char *s = *str; + *s++ = ch; + *str = s; + *length -= 1; + return 0; +} + +static int addstr(char **str, size_t *length, char *addme) { + assert(str && *str); + assert(length); + assert(addme); + const size_t sz = strlen(addme); + for (size_t i = 0; i < sz; i++) + if (addchar(str, length, addme[i]) < 0) + return -1; + return 0; +} + +int qmatrix_sprintb(const q_t *m, char *str, size_t length, unsigned base) { + assert(str); + assert(m); + const q_t *mm = &m[DATA]; + const size_t rows = m[ROW], columns = m[COLUMN]; + if (base < 2 || base > 36) + return -1; + if (!qmatrix_is_valid(m)) + return addstr(&str, &length, "[ INVALID ]"); + if (addstr(&str, &length, "[ ") < 0) + return -1; + for (size_t i = 0; i < rows; i++) { + for (size_t j = 0; j < columns; j++) { + const int r = qsprintb(mm[i*columns + j], str, length, base); + if (r < 0) + return -1; + if ((length - r) > length) + return -1; + length -= r; + str += r; + if (rows) + if (addchar(&str, &length, columns && j < (columns - 1) ? ',' : i < rows - 1 ? ';' : ' ') < 0) + return -1; + if ((columns && j < (columns - 1)) || (i < (rows - 1))) + if (addchar(&str, &length, ' ') < 0) + return -1; + } + } + if (addchar(&str, &length, ']') < 0) + return -1; + return 0; +} + +size_t qmatrix_string_length(const q_t *m) { + assert(m); + if (!qmatrix_is_valid(m)) + return 128; /* space for invalid matrix message */ + const size_t msize = m[LENGTH]; + const size_t r = (msize * + (32 /*max length if base 2 used)*/ + + 2 /* '-' and '.' */ + + 2 /* space and comma/semi colon separator */ + )) + 16 /* space for extra formatting */; + return r; +} + +/* See <https://github.com/jamesbowman/sincos> + * and "Math Toolkit for Real-Time Programming" by Jack Crenshaw + * + * The naming of these functions ('furman_') is incorrect, they do their + * computation on numbers represented in Furmans but they do not use a 'Furman + * algorithm'. As I do not have a better name, the name shall stick. */ +static int16_t _sine(const int16_t y) { + const int16_t s1 = 0x6487, s3 = -0x2953, s5 = 0x04f8; + const int16_t z = arshift((int32_t)y * y, 12); + int16_t prod = arshift((int32_t)z * s5, 16); + int16_t sum = s3 + prod; + prod = arshift((int32_t)z * sum, 16); + sum = s1 + prod; + return arshift((int32_t)y * sum, 13); +} + +static int16_t _cosine(int16_t y) { + const int16_t c0 = 0x7fff, c2 = -0x4ee9, c4 = 0x0fbd; + const int16_t z = arshift((int32_t)y * y, 12); + int16_t prod = arshift((int32_t)z * c4, 16); + const int16_t sum = c2 + prod; + prod = arshift((int32_t)z * sum, 15); + return c0 + prod; +} + +int16_t furman_sin(int16_t x) { + const int16_t n = 3 & arshift(x + 0x2000, 14); + x -= n << 14; + const int16_t r = (n & 1) ? _cosine(x) : _sine(x); + return (n & 2) ? -r : r; +} + +int16_t furman_cos(int16_t x) { + return furman_sin(x + 0x4000); +} + +/* expression evaluator */ + +enum { ASSOCIATE_NONE, ASSOCIATE_LEFT, ASSOCIATE_RIGHT, }; +enum { LEX_NUMBER, LEX_OPERATOR, LEX_END, }; + +int qexpr_init(qexpr_t *e) { + assert(e); + e->lpar = qop("("); + e->rpar = qop(")"); + e->negate = qop("negate"); + e->minus = qop("-"); + e->initialized = 1; + assert(e->lpar && e->rpar && e->negate && e->minus); + return 0; +} + +static int error(qexpr_t *e, const char *fmt, ...) { + assert(e); + assert(fmt); + if (e->error) + return 0; + va_list ap; + va_start(ap, fmt); + (void)vsnprintf(e->error_string, sizeof (e->error_string), fmt, ap); + va_end(ap); + e->error = -1; + return -QINT(1); +} + +static q_t numberify(const char *s) { + assert(s); + q_t q = 0; + (void) qconv(&q, s); + return q; +} + +static q_t qbase(q_t b) { + int nb = qtoi(b); + if (nb < 2 || nb > 36) + return -QINT(1); + qconf.base = nb; + return b; +} + +static q_t qplaces(q_t places) { + /* TODO: Bounds checks given base */ + qconf.dp = qtoi(places); + return places; +} + +static q_t check_div0(qexpr_t *e, q_t a, q_t b) { + assert(e); + UNUSED(a); + if (!b) + return error(e, "division by zero"); + return QINT(0); +} + +static q_t check_nlz(qexpr_t *e, q_t a) { // Not Less Zero + assert(e); + if (qless(a, QINT(0))) + return error(e, "negative argument"); + return QINT(0); +} + +static q_t check_nlez(qexpr_t *e, q_t a) { // Not Less Equal Zero + assert(e); + if (qeqless(a, QINT(0))) + return error(e, "negative or zero argument"); + return QINT(0); +} + +static q_t check_nlo(qexpr_t *e, q_t a) { // Not less than one + assert(e); + if (qless(a, QINT(1))) + return error(e, "out of range [1, INF]"); + return QINT(0); +} + +static q_t check_alo(qexpr_t *e, q_t a) { + assert(e); + if (qmore(qabs(a), QINT(1))) + return error(e, "out of range [-1, 1]"); + return QINT(0); +} + +const qoperations_t *qop(const char *op) { + assert(op); + static const qoperations_t ops[] = { + /* Binary Search Table: Use 'LC_ALL="C" sort -k 2 < table' to sort this */ + /* name function check function precedence arity left/right-assoc hidden */ + { "!", .eval.unary = qnot, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "!=", .eval.binary = qunequal, .check.binary = NULL, 2, 2, ASSOCIATE_LEFT, 0, }, + { "%", .eval.binary = qrem,/*!*/ .check.binary = check_div0, 3, 2, ASSOCIATE_LEFT, 0, }, + { "&", .eval.binary = qand, .check.binary = NULL, 2, 2, ASSOCIATE_LEFT, 0, }, + { "(", .eval.unary = NULL, .check.unary = NULL, 0, 0, ASSOCIATE_NONE, 0, }, + { ")", .eval.unary = NULL, .check.unary = NULL, 0, 0, ASSOCIATE_NONE, 0, }, + { "*", .eval.binary = qmul, .check.binary = NULL, 3, 2, ASSOCIATE_LEFT, 0, }, + { "+", .eval.binary = qadd, .check.binary = NULL, 2, 2, ASSOCIATE_LEFT, 0, }, + { "-", .eval.binary = qsub, .check.binary = NULL, 2, 2, ASSOCIATE_LEFT, 0, }, + { "/", .eval.binary = qdiv, .check.binary = check_div0, 3, 2, ASSOCIATE_LEFT, 0, }, + { "<", .eval.binary = qless, .check.binary = NULL, 2, 2, ASSOCIATE_LEFT, 0, }, + { "<<", .eval.binary = qlls, .check.binary = NULL, 4, 2, ASSOCIATE_RIGHT, 0, }, + { "<=", .eval.binary = qeqless, .check.binary = NULL, 2, 2, ASSOCIATE_LEFT, 0, }, + { "==", .eval.binary = qequal, .check.binary = NULL, 2, 2, ASSOCIATE_LEFT, 0, }, + { ">", .eval.binary = qmore, .check.binary = NULL, 2, 2, ASSOCIATE_LEFT, 0, }, + { ">=", .eval.binary = qeqmore, .check.binary = NULL, 2, 2, ASSOCIATE_LEFT, 0, }, + { ">>", .eval.binary = qlrs, .check.binary = NULL, 4, 2, ASSOCIATE_RIGHT, 0, }, + { "^", .eval.binary = qxor, .check.binary = NULL, 2, 2, ASSOCIATE_LEFT, 0, }, + { "_div", .eval.binary = qcordic_div, .check.binary = NULL, 5, 2, ASSOCIATE_RIGHT, 1, }, + { "_exp", .eval.unary = qcordic_exp, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 1, }, + { "_ln", .eval.unary = qcordic_ln, .check.unary = check_nlez, 5, 1, ASSOCIATE_RIGHT, 1, }, + { "_mul", .eval.binary = qcordic_mul, .check.binary = NULL, 5, 2, ASSOCIATE_RIGHT, 1, }, + { "_sqrt", .eval.unary = qcordic_sqrt, .check.unary = check_nlz, 5, 1, ASSOCIATE_RIGHT, 1, }, + { "abs", .eval.unary = qabs, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "acos", .eval.unary = qacos, .check.unary = check_alo, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "acosh", .eval.unary = qacosh, .check.unary = check_nlo, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "arshift", .eval.binary = qars, .check.binary = NULL, 4, 2, ASSOCIATE_RIGHT, 1, }, + { "asin", .eval.unary = qasin, .check.unary = check_alo, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "asinh", .eval.unary = qasinh, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "atan", .eval.unary = qatan, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "atan2", .eval.binary = qatan2, .check.binary = NULL, 5, 2, ASSOCIATE_RIGHT, 1, }, + { "atanh", .eval.unary = qatanh, .check.unary = check_alo, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "base", .eval.unary = qbase, .check.unary = NULL, 2, 1, ASSOCIATE_RIGHT, 0, }, + { "ceil", .eval.unary = qceil, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "copysign", .eval.binary = qcopysign, .check.binary = NULL, 4, 2, ASSOCIATE_RIGHT, 1, }, + { "cos", .eval.unary = qcos, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "cosh", .eval.unary = qcosh, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "cot", .eval.unary = qcot, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "deg2rad", .eval.unary = qdeg2rad, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "even?", .eval.unary = qiseven, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "exp", .eval.unary = qexp, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "floor", .eval.unary = qfloor, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "hypot", .eval.binary = qhypot, .check.binary = NULL, 5, 2, ASSOCIATE_RIGHT, 0, }, + { "int?", .eval.unary = qisinteger, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "log", .eval.unary = qlog, .check.unary = check_nlez, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "lshift", .eval.binary = qlls, .check.binary = NULL, 4, 2, ASSOCIATE_RIGHT, 1, }, + { "max", .eval.binary = qmax, .check.binary = NULL, 5, 2, ASSOCIATE_RIGHT, 1, }, + { "min", .eval.binary = qmin, .check.binary = NULL, 5, 2, ASSOCIATE_RIGHT, 1, }, + { "mod", .eval.binary = qmod, .check.binary = check_div0, 3, 2, ASSOCIATE_LEFT, 0, }, + { "neg?", .eval.unary = qisnegative, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "negate", .eval.unary = qnegate, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "odd?", .eval.unary = qisodd, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "places", .eval.unary = qplaces, .check.unary = NULL, 2, 1, ASSOCIATE_RIGHT, 0, }, + { "pos?", .eval.unary = qispositive, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "pow", .eval.binary = qpow, .check.binary = NULL, 5, 2, ASSOCIATE_RIGHT, 0, }, + { "rad2deg", .eval.unary = qrad2deg, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "rem", .eval.binary = qrem, .check.binary = check_div0, 3, 2, ASSOCIATE_LEFT, 0, }, + { "round", .eval.unary = qround, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "rshift", .eval.binary = qlrs, .check.binary = NULL, 4, 2, ASSOCIATE_RIGHT, 1, }, + { "sign", .eval.unary = qsign, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "signum", .eval.unary = qsignum, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "sin", .eval.unary = qsin, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "sinh", .eval.unary = qsinh, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "sqrt", .eval.unary = qsqrt, .check.unary = check_nlz, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "tan", .eval.unary = qtan, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "tanh", .eval.unary = qtanh, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "trunc", .eval.unary = qtrunc, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + { "|", .eval.binary = qor, .check.binary = NULL, 2, 2, ASSOCIATE_LEFT, 0, }, + { "~", .eval.unary = qinvert, .check.unary = NULL, 5, 1, ASSOCIATE_RIGHT, 0, }, + }; + const size_t length = (sizeof ops / sizeof ops[0]); + size_t l = 0, r = length - 1; + while (l <= r) { // Iterative Binary Search + size_t m = l + ((r - l)/2u); + assert (m < length); + const int comp = strcmp(ops[m].name, op); + if (comp == 0) + return &ops[m]; + if (comp < 0) + l = m + 1; + else + r = m - 1; + } + return NULL; +} + +static int number_push(qexpr_t *e, q_t num) { + assert(e); + if (e->error) + return -1; + if (e->numbers_count > (e->numbers_max - 1)) { + error(e, "number stack overflow"); + return -1; + } + e->numbers[e->numbers_count++] = num; + return 0; +} + +static q_t number_pop(qexpr_t *e) { + assert(e); + if (e->error) + return -1; + if (!(e->numbers_count)) { + error(e, "number stack empty"); + return -1; /* error handled elsewhere */ + } + return e->numbers[--(e->numbers_count)]; +} + +static int op_push(qexpr_t *e, const qoperations_t *op) { + assert(e); + assert(op); + if (e->error) + return -1; + if (e->ops_count > (e->ops_max - 1)) { + error(e, "operator stack overflow"); + return -1; + } + e->ops[e->ops_count++] = op; + return 0; +} + +int qexpr_error(qexpr_t *e) { + assert(e); + assert(e->initialized); + return e->error; +} + +q_t qexpr_result(qexpr_t *e) { + assert(e); + assert(e->initialized); + assert(e->error == 0); + assert(e->numbers_count == 1); + return e->numbers[0]; +} + +static const qoperations_t *op_pop(qexpr_t *e) { + assert(e); + if (e->error) + return NULL; + if (!(e->ops_count)) { + error(e, "operator stack empty"); + return NULL; + } + return e->ops[--(e->ops_count)]; +} + +static int op_eval(qexpr_t *e) { + assert(e); + const qoperations_t *pop = op_pop(e); + if (!pop) + return -1; + const q_t a = number_pop(e); + const int exists = pop->arity == 1 ? BOOLIFY(pop->eval.unary) : BOOLIFY(pop->eval.binary); + if (!exists) { + error(e, "syntax error"); + return -1; + } + if (pop->arity == 1) { + if (pop->check.unary && pop->check.unary(e, a) < 0) { + error(e, "unary check failed"); + return -1; + } + return number_push(e, pop->eval.unary(a)); + } + const q_t b = number_pop(e); + if (pop->check.binary && pop->check.binary(e, b, a)) { + error(e, "binary check failed"); + return -1; + } + + return number_push(e, pop->eval.binary(b, a)); +} + +static int shunt(qexpr_t *e, const qoperations_t *op) { + assert(e); + assert(op); + if (op == e->lpar) { + return op_push(e, op); + } else if (op == e->rpar) { + while (e->ops_count && e->ops[e->ops_count - 1] != e->lpar) + if (op_eval(e) < 0 || e->error) + break; + const qoperations_t *pop = op_pop(e); + if (!pop || (pop != e->lpar)) { + e->error = 0; /* clear error so following error is printed */ + error(e, "expected \"(\""); + return -1; + } + return 0; + } else if (op->assocativity == ASSOCIATE_RIGHT) { + while (e->ops_count && op->precedence < e->ops[e->ops_count - 1]->precedence) + if (op_eval(e) < 0 || e->error) + break; + } else { + while (e->ops_count && op->precedence <= e->ops[e->ops_count - 1]->precedence) + if (op_eval(e) < 0 || e->error) + break; + } + return op_push(e, op); +} + +static int variable_name_is_valid(const char *n) { + assert(n); + if (!isalpha(*n) && !(*n == '_')) + return 0; + for (n++; *n; n++) + if (!isalnum(*n) && !(*n == '_')) + return 0; + return 1; +} + +static qvariable_t *variable_lookup(qexpr_t *e, const char *name) { + assert(e); + assert(name); + for (size_t i = 0; i < e->vars_max; i++) { + qvariable_t *v = e->vars[i]; + assert(v->name); + assert(variable_name_is_valid(v->name)); + if (!strcmp(v->name, name)) + return v; + } + return NULL; +} + +static int lex(qexpr_t *e, const char **expr) { + assert(e); + assert(expr && *expr); + int r = 0; + const char *s = *expr; + qvariable_t *v = NULL; + e->id_count = 0; + e->number = 0; + e->op = NULL; + memset(e->id, 0, sizeof (e->id)); + for (; *s && isspace(*s); s++) + ; + if (!(*s)) + return LEX_END; + if (isalpha(*s) || *s == '_') { + for (; e->id_count < sizeof(e->id) && *s && (isalnum(*s) || *s == '_');) + e->id[e->id_count++] = *s++; + if ((v = variable_lookup(e, e->id))) { + e->number = v->value; + r = LEX_NUMBER; + } else if ((e->op = qop(e->id))) { + r = LEX_OPERATOR; + } else { + r = -1; + } + } else { + if (ispunct(*s)) { + const qoperations_t *op1 = NULL, *op2 = NULL; + int set = 0; + e->id[e->id_count++] = *s++; + op1 = qop(e->id); + if (*s && ispunct(*s)) { + set = 1; + e->id[e->id_count++] = *s++; + op2 = qop(e->id); + } + r = (op1 || op2) ? LEX_OPERATOR : -1; + e->op = op2 ? op2 : op1; + if (e->op == op1 && set) { + s--; + e->id_count--; + e->id[1] = 0; + } + } else if (isdigit(*s)) { + r = LEX_NUMBER; + int dot = 0; + for (; e->id_count < sizeof(e->id) && *s; s++) { + const int ch = *s; + if (!(isdigit(ch) || (ch == '.' && !dot))) + break; + e->id[e->id_count++] = ch; + if (ch == '.') + dot = 1; + } + e->number = numberify(e->id); + } else { + r = -1; + } + } + /*printf("id(%d) %d => %s\n", (int)(s - *expr), r, e->id);*/ + *expr = s; + return r; +} + +int qexpr(qexpr_t *e, const char *expr) { + assert(e); + assert(expr); + int firstop = 1; + const qoperations_t *previous = NULL; + if (e->initialized) { + memset(e->error_string, 0, sizeof (e->error_string)); + e->error = 0; + e->ops_count = 0; + e->numbers_count = 0; + e->initialized = 1; + } + for (int l = 0; l != LEX_END && !(e->error);) { + switch ((l = lex(e, &expr))) { + case LEX_NUMBER: + number_push(e, e->number); + previous = NULL; + firstop = 0; + break; + case LEX_OPERATOR: { + const qoperations_t *op = e->op; + if (CONFIG_Q_HIDE_FUNCS && op->hidden) { + error(e, "unknown operator \"%s\"", op->name); + goto end; + } + if (firstop || (previous && previous != e->rpar)) { + if (e->op == e->minus) { + op = e->negate; + } else if (e->op->arity == 1) { + /* do nothing */ + } else if (e->op != e->lpar) { + assert(e->op); + error(e, "invalid use of \"%s\"", e->op->name); + goto end; + } + } + shunt(e, op); + previous = op; + firstop = 0; + break; + } + case LEX_END: break; + default: + error(e, "invalid symbol: %s", e->id); + l = LEX_END; + } + } + while (e->ops_count) + if (op_eval(e) < 0 || e->error) + break; + if (e->numbers_count != 1) { + error(e, "invalid expression: %d", e->numbers_count); + return -1; + } + implies(e->error == 0, e->numbers_count == 1); +end: + return e->error == 0 ? 0 : -1; +} + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/q/q.h Tue Feb 25 13:28:29 2025 +1030 @@ -0,0 +1,329 @@ +/* Project: Q-Number (Q16.16, signed) library + * Author: Richard James Howe + * License: The Unlicense + * Email: howe.r.j.89@gmail.com + * Repo: <https://github.com/q> */ + +#ifndef Q_H +#define Q_H + +#ifdef __cplusplus +extern "C" { +#endif + +#include <stdint.h> +#include <stddef.h> + +#define QMAX_ID (32) +#define QMAX_ERROR (256) + +#define QBITS (12) +#define QMASK ((1ULL << QBITS) - 1ULL) +#define QHIGH (1ULL << (QBITS - 1ULL)) + +#define QMK(QHIGH, LOW, SF) ((ld_t)((((lu_t)QHIGH) << QBITS) | (QMASK & ((((lu_t)LOW) << QBITS) >> (SF))))) +#define QINT(INT) ((q_t)((u_t)(INT) << QBITS)) +#define QPI (QMK(0x3, 0x243F, 16)) + +#ifndef PREPACK +#define PREPACK +#endif + +#ifndef POSTPACK +#define POSTPACK +#endif + +#ifndef RESTRICT +#ifdef __cplusplus +#define RESTRICT +#else +#define RESTRICT restrict +#endif +#endif + +typedef int32_t q_t; /* Q Fixed Point Number, (Q16.16, Signed) */ +typedef int64_t ld_t; /* Double width of Q, signed, for internal calculations, not in Q format */ +typedef int32_t d_t; /* same width as Q, signed, but not in Q format */ +typedef uint32_t u_t; /* same width as Q, unsigned, but not in Q format */ + +typedef PREPACK struct { + size_t whole, /* number of bits for whole, or integer, part of Q number */ + fractional; /* number of bits for fractional part of Q number */ + const q_t zero, /* the constant '0' */ + bit, /* smallest 'q' number representable */ + one, /* the constant '1' */ + pi, /* the constant 'pi' */ + e, /* the constant 'e' */ + sqrt2, /* the square root of 2 */ + sqrt3, /* the square root of 3 */ + ln2, /* the natural logarithm of 2 */ + ln10, /* the natural logarithm of 10 */ + min, /* most negative 'q' number */ + max; /* most positive 'q' number */ + const uint32_t version; /* version in X.Y.Z format (Z = lowest 8 bits) */ +} POSTPACK qinfo_t; + +typedef PREPACK struct { + q_t rc, /* time constant */ + time, /* time of previous measurement */ + raw, /* previous raw value */ + filtered; /* filtered value */ +} POSTPACK qfilter_t; /* High/Low Pass Filter */ + +typedef PREPACK struct { + q_t d_gain, d_state; /* differentiator; gain, state */ + q_t i_gain, i_state, i_min, i_max; /* integrator; gain, state, minimum and maximum */ + q_t p_gain; /* proportional gain */ +} POSTPACK qpid_t; /* PID Controller <https://en.wikipedia.org/wiki/PID_controller> */ + +struct qexpr; +typedef struct qexpr qexpr_t; + +typedef PREPACK struct { + char *name; + union { + q_t (*unary) (q_t a); + q_t (*binary) (q_t a1, q_t a2); + } eval; + union { + q_t (*unary) (qexpr_t *e, q_t a); + q_t (*binary) (qexpr_t *e, q_t a1, q_t a2); + } check; + int precedence, arity, assocativity, hidden; +} POSTPACK qoperations_t; /* use in the expression evaluator */ + +typedef PREPACK struct { + char *name; + q_t value; +} POSTPACK qvariable_t; /* Variable which can be used with the expression evaluator */ + +struct PREPACK qexpr { + const qoperations_t **ops, *lpar, *rpar, *negate, *minus; + qvariable_t **vars; + char id[QMAX_ID]; + char error_string[QMAX_ERROR]; + q_t number; + const qoperations_t *op; + q_t *numbers; + size_t ops_count, ops_max; + size_t numbers_count, numbers_max; + size_t id_count; + size_t vars_max; + int error; + int initialized; +} POSTPACK; /* An expression evaluator for the Q library */ + +typedef q_t (*qbounds_t)(ld_t s); + +q_t qbound_saturate(ld_t s); /* default over/underflow behavior, saturation */ +q_t qbound_wrap(ld_t s); /* over/underflow behavior, wrap around */ + +typedef PREPACK struct { + qbounds_t bound; /* handles saturation when a number over or underflows */ + int dp; /* decimal points to print, negative specifies maximum precision */ + unsigned base; /* base to use for numeric number conversion */ +} POSTPACK qconf_t; /* Q format configuration options */ + +extern const qinfo_t qinfo; /* information about the format and constants */ +extern qconf_t qconf; /* @warning GLOBAL Q configuration options */ + +int qtoi(q_t toi); +q_t qint(int toq); +signed char qtoc(const q_t q); +q_t qchar(signed char c); +short qtoh(const q_t q); +q_t qshort(short s); +long qtol(const q_t q); +q_t qlong(long l); +long long qtoll(const q_t q); +q_t qvlong(long long ll); + +q_t qisnegative(q_t a); +q_t qispositive(q_t a); +q_t qisinteger(q_t a); +q_t qisodd(q_t a); +q_t qiseven(q_t a); + +q_t qless(q_t a, q_t b); +q_t qmore(q_t a, q_t b); +q_t qeqless(q_t a, q_t b); +q_t qeqmore(q_t a, q_t b); +q_t qequal(q_t a, q_t b); +q_t qunequal(q_t a, q_t b); +q_t qapproxequal(q_t a, q_t b, q_t epsilon); +q_t qapproxunequal(q_t a, q_t b, q_t epsilon); +q_t qwithin(q_t v, q_t b1, q_t b2); +q_t qwithin_interval(q_t v, q_t expected, q_t allowance); + +q_t qnegate(q_t a); +q_t qmin(q_t a, q_t b); +q_t qmax(q_t a, q_t b); +q_t qabs(q_t a); +q_t qcopysign(q_t a, q_t b); +q_t qsign(q_t a); +q_t qsignum(q_t a); + +q_t qadd(q_t a, q_t b); +q_t qsub(q_t a, q_t b); +q_t qmul(q_t a, q_t b); +q_t qdiv(q_t a, q_t b); +q_t qrem(q_t a, q_t b); +q_t qmod(q_t a, q_t b); +q_t qfma(q_t a, q_t b, q_t c); +q_t qsqr(q_t x); +q_t qexp(q_t e); +q_t qlog(q_t n); +q_t qsqrt(q_t x); + +q_t qround(q_t q); +q_t qceil(q_t q); +q_t qtrunc(q_t q); +q_t qfloor(q_t q); + +q_t qand(q_t a, q_t b); +q_t qxor(q_t a, q_t b); +q_t qor(q_t a, q_t b); +q_t qinvert(q_t a); +q_t qnot(q_t a); +q_t qlogical(q_t a); + +q_t qlls(q_t a, q_t b); +q_t qlrs(q_t a, q_t b); +q_t qals(q_t a, q_t b); +q_t qars(q_t a, q_t b); + +q_t qpow(q_t n, q_t exp); + +int qsprint(q_t p, char *s, size_t length); +int qsprintb(q_t p, char *s, size_t length, u_t base); +int qsprintbdp(q_t p, char *s, size_t length, u_t base, d_t idp); +int qnconv(q_t *q, const char *s, size_t length); +int qnconvb(q_t *q, const char *s, size_t length, d_t base); +int qconv(q_t *q, const char *s); +int qconvb(q_t *q, const char * const s, d_t base); +int qnconvbdp(q_t *q, const char *s, size_t length, d_t base, u_t idp); + +void qsincos(q_t theta, q_t *sine, q_t *cosine); +q_t qsin(q_t theta); +q_t qcos(q_t theta); +q_t qtan(q_t theta); +q_t qcot(q_t theta); +q_t qhypot(q_t a, q_t b); + +q_t qatan(q_t t); +q_t qatan2(q_t x, q_t y); +q_t qasin(q_t t); +q_t qacos(q_t t); +q_t qsinh(q_t a); +q_t qcosh(q_t a); +q_t qtanh(q_t a); + +q_t qatanh(q_t t); +q_t qasinh(q_t t); +q_t qacosh(q_t t); + +void qsincosh(q_t a, q_t *sinh, q_t *cosh); + +q_t qcordic_ln(q_t d); /* CORDIC testing only */ +q_t qcordic_exp(q_t e); /* CORDIC testing only; useless for large values */ +q_t qcordic_sqrt(q_t a); /* CORDIC testing only; do not use, a <= 2, a >= 0 */ +q_t qcordic_mul(q_t a, q_t b); /* CORDIC testing only; do not use */ +q_t qcordic_div(q_t a, q_t b); /* CORDIC testing only; do not use */ +q_t qcordic_circular_gain(int n); +q_t qcordic_hyperbolic_gain(int n); + +void qpol2rec(q_t magnitude, q_t theta, q_t *i, q_t *j); +void qrec2pol(q_t i, q_t j, q_t *magnitude, q_t *theta); + +q_t qrad2deg(q_t rad); +q_t qdeg2rad(q_t deg); + +d_t dpower(d_t b, unsigned e); +d_t dlog(d_t n, unsigned base); +d_t arshift(d_t v, unsigned p); +int qpack(const q_t *q, char *buffer, size_t length); +int qunpack(q_t *q, const char *buffer, size_t length); + +q_t qsimpson(q_t (*f)(q_t), q_t x1, q_t x2, unsigned n); /* numerical integrator of f, between x1, x2, for n steps */ + +void qfilter_init(qfilter_t *f, q_t time, q_t rc, q_t seed); +q_t qfilter_low_pass(qfilter_t *f, q_t time, q_t data); +q_t qfilter_high_pass(qfilter_t *f, q_t time, q_t data); +q_t qfilter_value(const qfilter_t *f); + +q_t qpid_update(qpid_t *pid, const q_t error, const q_t position); + +/* A matrix consists of at least four elements, a meta data field, + * the length of the array (which must be big enough to store + * row*column, but may be * larger) and a row and a column count + * in unsigned integer format, and the array elements in Q format. + * This simplifies storage and declaration of matrices. + * + * An example, the 2x3 matrix: + * + * [ 1, 2, 3; 4, 5, 6 ] + * + * Should be defined as: + * + * q_t m[] = { 0, 2*3, 2, 3, QINT(1), QINT(2), QINT(3), QINT(4), QINT(5), QINT(6) }; + * + */ +int qmatrix_apply_unary(q_t *r, const q_t *a, q_t (*func)(q_t)); +int qmatrix_apply_scalar(q_t *r, const q_t *a, q_t (*func)(q_t, q_t), const q_t c); +int qmatrix_apply_binary(q_t * RESTRICT r, const q_t *a, const q_t *b, q_t (*func)(q_t, q_t)); +int qmatrix_sprintb(const q_t *m, char *str, size_t length, unsigned base); +int qmatrix_resize(q_t *m, const size_t row, const size_t column); +int qmatrix_copy(q_t *r, const q_t *a); +size_t qmatrix_string_length(const q_t *m); + +q_t qmatrix_trace(const q_t *m); +q_t qmatrix_determinant(const q_t *m); +q_t qmatrix_equal(const q_t *a, const q_t *b); + +int qmatrix_zero(q_t *r); +int qmatrix_one(q_t *r); +int qmatrix_identity(q_t *r); /* turn into identity matrix, r must be square */ + +int qmatrix_logical(q_t *r, const q_t *a); +int qmatrix_not(q_t *r, const q_t *a); +int qmatrix_signum(q_t *r, const q_t *a); +int qmatrix_invert(q_t *r, const q_t *a); + +int qmatrix_is_valid(const q_t *m); +int qmatrix_is_square(const q_t *m); + +int qmatrix_transpose(q_t * RESTRICT r, const q_t * RESTRICT m); +int qmatrix_add(q_t * RESTRICT r, const q_t *a, const q_t *b); +int qmatrix_sub(q_t * RESTRICT r, const q_t *a, const q_t *b); +int qmatrix_mul(q_t * RESTRICT r, const q_t *a, const q_t *b); +int qmatrix_and(q_t *r, const q_t *a, const q_t *b); +int qmatrix_or (q_t *r, const q_t *a, const q_t *b); +int qmatrix_xor(q_t *r, const q_t *a, const q_t *b); + +int qmatrix_scalar_add(q_t *r, const q_t *a, const q_t scalar); +int qmatrix_scalar_sub(q_t *r, const q_t *a, const q_t scalar); +int qmatrix_scalar_mul(q_t *r, const q_t *a, const q_t scalar); +int qmatrix_scalar_div(q_t *r, const q_t *a, const q_t scalar); +int qmatrix_scalar_mod(q_t *r, const q_t *a, const q_t scalar); +int qmatrix_scalar_rem(q_t *r, const q_t *a, const q_t scalar); +int qmatrix_scalar_and(q_t *r, const q_t *a, const q_t scalar); +int qmatrix_scalar_or (q_t *r, const q_t *a, const q_t scalar); +int qmatrix_scalar_xor(q_t *r, const q_t *a, const q_t scalar); + +/* Expression evaluator */ + +int qexpr(qexpr_t *e, const char *expr); +int qexpr_init(qexpr_t *e); +int qexpr_error(qexpr_t *e); +q_t qexpr_result(qexpr_t *e); +const qoperations_t *qop(const char *op); + +/* A better cosine/sine, not in Q format */ + +int16_t furman_sin(int16_t x); /* SINE: 1 Furman = 1/65536 of a circle */ +int16_t furman_cos(int16_t x); /* COSINE: 1 Furman = 1/65536 of a circle */ + +#ifdef __cplusplus +} +#endif +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/q/readme.md Tue Feb 25 13:28:29 2025 +1030 @@ -0,0 +1,105 @@ +# Q: Fixed Point Number Library + +| Project | Q: Fixed Point Number Library | +| --------- | --------------------------------- | +| Author | Richard James Howe | +| Copyright | 2018 Richard James Howe | +| License | MIT | +| Email | howe.r.j.89@gmail.com | +| Website | <https://github.com/howerj/q> | + +This is a small fixed point number library designed for embedded use. It +implements all of your favorite transcendental functions, plus only the best +basic operators, selected for your calculating pleasure. The format is a signed +[Q16.16][], which is good enough for [Doom][] and good enough for you. + +The default [makefile][] target builds the library and a test program, which will +require a [C][] compiler (and Make). The library is small enough that is can be +easily modified to suite your purpose. The dependencies are minimal; they are few +string handling functions, and '[tolower][]' for numeric input. This should allow +the code to ported to the platform of your choice. The 'run' make target builds +the test program (called 'q') and runs it on some input. The '-h' option will +spit out a more detailed help. + +I would compile the library with the '-fwrapv' option enabled, you might +be some kind of Maverick who doesn't play by no rules however. + +The trigonometric functions, and some others, are implemented internally with +[CORDIC][]. + +Of note, all operators are bounded by minimum and maximum values which are not +shown in the following table and by default all arithmetic is saturating. The +effective input range of a number might lower than what is possible given a +mathematical functions definition - either because of the limited range of the +[Q16.16][] type, or because the implementation of a function introduces too +much error along some part of its' input range. Caveat Emptor (although you're +not exactly paying for this library now, are you? Caveat lector perhaps). + +( This table needs completing, specifically the input ranges... ) + +| C Function | Operator | Input Range | Asserts | Notes | +| ------------- | ----------- | ----------- | -------- | ----------------------------------------------- | +| qadd(a, b) | a + b | | | Addition | +| qsub(a, b) | a \- b | | | Subtraction | +| qdiv(a, b) | a / b | b != 0 | Yes | Division | +| qmul(a, b) | a \* b | | | Multiplication | +| qrem(a, b) | a rem b | b != 0 | Yes | Remainder: remainder after division | +| qmod(a, b) | a mod b | b != 0 | Yes | Modulo | +| qsin(theta) | sin(theta) | | | Sine | +| qcos(theta) | cos(theta) | | | Cosine | +| qtan(theta) | tan(theta) | | | Tangent | +| qcot(theta) | cot(theta) | | | Cotangent | +| qhypot(a, b) | hypot(a, b) | | | Hypotenuse; sqrt(a\*a + b\*b) | +| qasin(x) | asin(x) | abs(x) <= 1 | Yes | Arcsine | +| qacos(x) | acos(x) | abs(x) <= 1 | Yes | Arccosine | +| qatan(t) | atan(t) | | | Arctangent | +| qsinh(a) | sinh(a) | | | Hyperbolic Sine | +| qcosh(a) | cosh(a) | | | Hyperbolic Cosine | +| qtanh(a) | tanh(a) | | | Hyperbolic Tangent | +| qasinh(a) | asinh(a) | | | Inverse Hyperbolic Sine | +| qacosh(a) | acosh(a) | | | Inverse Hyperbolic Cosine | +| qatanh(a) | atanh(a) | | | Inverse Hyperbolic Tangent | +| qexp(e) | exp(e) | e < ln(MAX) | No | Exponential function, High error for 'e' > ~7. | +| qlog(n) | log(n) | n > 0 | Yes | Natural Logarithm | +| qsqrt(x) | sqrt(x) | n >= 0 | Yes | Square Root | +| qround(q) | round(q) | | | Round to nearest value | +| qceil(q) | ceil(q) | | | Round up value | +| qtrunc(q) | trunc(q) | | | Truncate value | +| qfloor(q) | floor(q) | | | Round down value | +| qnegate(a) | -a | | | Negate a number | +| qabs(a) | abs(a) | | | Absolute value of a number | +| qfma(a, b, c) | (a\*b)+c | | | Fused Multiply Add, uses Q32.32 internally | +| qequal(a, b) | a == b | | | | +| qsignum(a) | signum(a) | | | Signum function | +| qsign(a) | sgn(a) | | | Sign function | +| | | | | | + +For the round/ceil/trunc/floor functions the following table from the +[cplusplus.com][] helps: + +| value | round | floor | ceil | trunc | +| ----- | ----- | ----- | ---- | ----- | +| 2.3 | 2.0 | 2.0 | 3.0 | 2.0 | +| 3.8 | 4.0 | 3.0 | 4.0 | 3.0 | +| 5.5 | 6.0 | 5.0 | 6.0 | 5.0 | +| -2.3 | -2.0 | -3.0 | -2.0 | -2.0 | +| -3.8 | -4.0 | -4.0 | -3.0 | -3.0 | +| -5.5 | -6.0 | -6.0 | -5.0 | -5.0 | + +Have fun with the adding, and subtracting, and stuff, I hope it goes well. It +would be cool to make an [APL][] interpreter built around this library. Testing +would become much easier as you could use programming language constructs to +create new tests over larger ranges of numbers. + +[APL]: https://en.wikipedia.org/wiki/APL_(programming_language) +[Doom]: https://en.wikipedia.org/wiki/Doom_(1993_video_game) +[tolower]: http://www.cplusplus.com/reference/cctype/tolower/ +[makefile]: https://en.wikipedia.org/wiki/Make_(software) +[C]: https://en.wikipedia.org/wiki/C_%28programming_language%29 +[cplusplus.com]: http://www.cplusplus.com/reference/cmath/round/ +[Q16.16]: https://en.wikipedia.org/wiki/Fixed-point_arithmetic +[CORDIC]: https://en.wikipedia.org/wiki/CORDIC +[VHDL]: https://en.wikipedia.org/wiki/VHDL +[FPGA]: https://en.wikipedia.org/wiki/Field-programmable_gate_array + +<style type="text/css">body{margin:40px auto;max-width:850px;line-height:1.6;font-size:16px;color:#444;padding:0 10px}h1,h2,h3{line-height:1.2}</style>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/q/t.c Tue Feb 25 13:28:29 2025 +1030 @@ -0,0 +1,581 @@ +/**@file t.c + * @brief Q-Number (Q16.16, signed) library test bench + * @copyright Richard James Howe (2018) + * @license MIT + * @email howe.r.j.89@gmail.com + * @site <https://github.com/howerj/q> + * + * This file contains a wrapper for testing the Q Number library, + * it includes a command processor and some built in tests. View + * the help string later on for more information */ +#include "q.h" +#include <assert.h> +#include <ctype.h> +#include <errno.h> +#include <stdbool.h> +#include <stdio.h> +#include <stdlib.h> +#include <string.h> + +#define N (16) +#define BUILD_BUG_ON(condition) ((void)sizeof(char[1 - 2*!!(condition)])) + +static int qprint(FILE *out, const q_t p) { + assert(out); + char buf[64 + 1] = { 0, }; + const int r = qsprint(p, buf, sizeof buf); + return r < 0 ? r : fprintf(out, "%s", buf); +} + +static int printq(FILE *out, q_t q, const char *msg) { + assert(out); + assert(msg); + if (fprintf(out, "%s = ", msg) < 0) + return -1; + if (qprint(out, q) < 0) + return -1; + if (fputc('\n', out) != '\n') + return -1; + return 0; +} + +static int print_sincos(FILE *out, const q_t theta) { + assert(out); + q_t sine = qinfo.zero, cosine = qinfo.zero; + qsincos(theta, &sine, &cosine); + if (qprint(out, theta) < 0) + return -1; + if (fputc(',', out) != ',') + return -1; + if (qprint(out, sine) < 0) + return -1; + if (fputc(',', out) != ',') + return -1; + if (qprint(out, cosine) < 0) + return -1; + if (fputc('\n', out) != '\n') + return -1; + return 0; +} + +static int print_sincos_table(FILE *out) { + assert(out); + const q_t tpi = qdiv(qinfo.pi, qint(20)); + const q_t end = qmul(qinfo.pi, qint(2)); + const q_t start = qnegate(end); + if (fprintf(out, "theta,sine,cosine\n") < 0) + return -1; + for (q_t i = start; qless(i, end); i = qadd(i, tpi)) + if (print_sincos(out, i) < 0) + return -1; + return 0; +} + +static int qinfo_print(FILE *out, const qinfo_t *qi) { + assert(out); + assert(qi); + if (fprintf(out, "Q%u.%u Info\n", (unsigned)qi->whole, (unsigned)qi->fractional) < 0) return -1; + if (printq(out, qi->bit, "bit") < 0) return -1; + if (printq(out, qi->one, "one") < 0) return -1; + if (printq(out, qi->zero, "zero") < 0) return -1; + if (printq(out, qi->pi, "pi") < 0) return -1; + if (printq(out, qi->e, "e") < 0) return -1; + if (printq(out, qi->sqrt2, "sqrt2") < 0) return -1; + if (printq(out, qi->sqrt3, "sqrt3") < 0) return -1; + if (printq(out, qi->ln2, "ln2") < 0) return -1; + if (printq(out, qi->ln10, "ln10") < 0) return -1; + if (printq(out, qi->min, "min") < 0) return -1; + if (printq(out, qi->max, "max") < 0) return -1; + if (printq(out, qcordic_circular_gain(-1), "circular-gain") < 0) return -1; + if (printq(out, qcordic_hyperbolic_gain(-1), "hyperbolic-gain") < 0) return -1; + return 0; +} + +static int qconf_print(FILE *out, const qconf_t *qc) { + assert(out); + assert(qc); + if (fprintf(out, "Q Configuration\n") < 0) return -1; + const char *bounds = "unknown"; + qbounds_t bound = qc->bound; + if (bound == qbound_saturate) + bounds = "saturate"; + if (bound == qbound_wrap) + bounds = "wrap"; + if (fprintf(out, "overflow handler: %s\n", bounds) < 0) return -1; + if (fprintf(out, "input/output radix: %u (0 = special case)\n", qc->base) < 0) return -1; + if (fprintf(out, "decimal places: %d\n", qc->dp) < 0) return -1; + return 0; +} + +static FILE *fopen_or_die(const char *file, const char *mode) { + assert(file && mode); + FILE *h = NULL; + errno = 0; + if (!(h = fopen(file, mode))) { + (void)fprintf(stderr, "file open \"%s\" (mode %s) failed: %s\n", file, mode, strerror(errno)); + exit(EXIT_FAILURE); + } + return h; +} + +typedef enum { + EVAL_OK_E, + EVAL_COMMENT_E, + EVAL_ERROR_SCAN_E, + EVAL_ERROR_TYPE_E, + EVAL_ERROR_CONVERT_E, + EVAL_ERROR_OPERATION_E, + EVAL_ERROR_ARG_COUNT_E, + EVAL_ERROR_UNEXPECTED_RESULT_E, + EVAL_ERROR_LIMIT_MODE_E, + + EVAL_ERROR_MAX_ERRORS_E, /**< not an error, but a count of errors */ +} eval_errors_e; + +static const char *eval_error(int e) { + if (e < 0 || e >= EVAL_ERROR_MAX_ERRORS_E) + return "unknown"; + const char *msgs[EVAL_ERROR_MAX_ERRORS_E] = { + [EVAL_OK_E] = "ok", + [EVAL_COMMENT_E] = "(comment)", + [EVAL_ERROR_SCAN_E] = "invalid input line", + [EVAL_ERROR_TYPE_E] = "unknown function type", + [EVAL_ERROR_CONVERT_E] = "numeric conversion failed", + [EVAL_ERROR_OPERATION_E] = "unknown operation", + [EVAL_ERROR_ARG_COUNT_E] = "incorrect argument count", + [EVAL_ERROR_LIMIT_MODE_E] = "unknown limit mode ('|' or '%' allowed)", + [EVAL_ERROR_UNEXPECTED_RESULT_E] = "unexpected result", + }; + return msgs[e] ? msgs[e] : "unknown"; +} + +static int eval_unary_arith(q_t (*m)(q_t), q_t expected, q_t bound, q_t arg1, q_t *result) { + assert(m); + assert(result); + const q_t r = m(arg1); + *result = r; + if (qwithin_interval(r, expected, bound)) + return 0; + return -1; +} + +static int eval_binary_arith(q_t (*f)(q_t, q_t), q_t expected, q_t bound, q_t arg1, q_t arg2, q_t *result) { + assert(f); + assert(result); + const q_t r = f(arg1, arg2); + *result = r; + if (qwithin_interval(r, expected, bound)) + return 0; + return -1; +} + +static int comment(char *line) { + assert(line); + const size_t length = strlen(line); + for (size_t i = 0; i < length; i++) { + if (line[i] == '#') + return 1; + if (!isspace(line[i])) + break; + } + for (size_t i = 0; i < length; i++) { + if (line[i] == '#') { + line[i] = 0; + return 0; + } + } + return 0; +} + +static void bounds_set(char method) { + assert(method == '|' || method == '%'); + if (method == '|') + qconf.bound = qbound_saturate; + qconf.bound = qbound_wrap; +} + +static int eval(char *line, q_t *result) { + assert(line); + assert(result); + *result = qinfo.zero; + if (comment(line)) + return EVAL_COMMENT_E; + char operation[N] = { 0 }, expected[N] = { 0 }, bounds[N], arg1[N] = { 0 }, arg2[N] = { 0 }; + char limit = '|'; + const int count = sscanf(line, "%15s %15s +- %15s %c %15s %15s", operation, expected, bounds, &limit, arg1, arg2); + if (limit != '|' && limit != '%') + return -EVAL_ERROR_LIMIT_MODE_E; + bounds_set(limit); + if (count < 5) + return -EVAL_ERROR_SCAN_E; + const qoperations_t *func = qop(operation); + if (!func) + return -EVAL_ERROR_OPERATION_E; + q_t e = qinfo.zero, b = qinfo.zero, a1 = qinfo.zero, a2 = qinfo.zero; + const int argc = count - 4; + if (func->arity != argc) + return -EVAL_ERROR_ARG_COUNT_E; + if (qconv(&e, expected) < 0) + return -EVAL_ERROR_CONVERT_E; + if (qconv(&b, bounds) < 0) + return -EVAL_ERROR_CONVERT_E; + if (qconv(&a1, arg1) < 0) + return -EVAL_ERROR_CONVERT_E; + switch (func->arity) { + case 1: if (eval_unary_arith(func->eval.unary, e, b, a1, result) < 0) + return -EVAL_ERROR_UNEXPECTED_RESULT_E; + break; + case 2: if (qconv(&a2, arg2) < 0) + return -EVAL_ERROR_CONVERT_E; + if (eval_binary_arith(func->eval.binary, e, b, a1, a2, result) < 0) + return -EVAL_ERROR_UNEXPECTED_RESULT_E; + break; + default: + return -EVAL_ERROR_TYPE_E; + } + return EVAL_OK_E; +} + +static void trim(char *s) { + assert(s); + const int size = strlen(s); + for (int i = size - 1; i >= 0; i--) { + if (!isspace(s[i])) + break; + s[i] = 0; + } +} + +static int eval_file(FILE *input, FILE *output) { + assert(input); + assert(output); + char line[256] = { 0 }; + int rv = 0; + while (fgets(line, sizeof(line) - 1, input)) { + q_t result = 0; + const int r = eval(line, &result); + if (r == EVAL_COMMENT_E) + continue; + if (r < 0) { + const char *msg = eval_error(-r); + trim(line); + if (fprintf(output, "error: %s = ", line) < 0) return -1; + if (qprint(output, result) < 0) return -1; + if (fprintf(output, " : %s/%d\n", msg, r) < 0) return -1; + rv = -1; + continue; + } + trim(line); + char rstring[64 + 1] = { 0, }; + if (qsprint(result, rstring, sizeof(rstring) - 1) < 0) return -1; + if (fprintf(output, "ok: %s | (%s)\n", line, rstring) < 0) return -1; + memset(line, 0, sizeof line); + } + return rv; +} + +/* --- Unit Test Framework --- */ + +typedef struct { + unsigned passed, + run; +} unit_test_t; + +static inline unit_test_t _unit_test_start(const char *file, const char *func, unsigned line) { + assert(file); + assert(func); + unit_test_t t = { .run = 0, .passed = 0 }; + fprintf(stdout, "Start tests: %s in %s:%u\n\n", func, file, line); + return t; +} + +static inline void _unit_test_statement(const char *expr_str) { + assert(expr_str); + if (fprintf(stdout, " STATE: %s\n", expr_str) < 0) + abort(); +} + +static inline void _unit_test(unit_test_t *t, int failed, const char *expr_str, const char *file, const char *func, unsigned line, int die) { + assert(t); + assert(expr_str); + assert(file); + assert(func); + if (failed) { + fprintf(stdout, " FAILED: %s (%s:%s:%u)\n", expr_str, file, func, line); + if (die) { + fputs("VERIFY FAILED - EXITING\n", stdout); + exit(EXIT_FAILURE); + } + } else { + fprintf(stdout, " OK: %s\n", expr_str); + t->passed++; + } + t->run++; +} + +static inline int unit_test_finish(unit_test_t *t) { + assert(t); + fprintf(stdout, "Tests passed/total: %u/%u\n", t->passed, t->run); + if (t->run != t->passed) { + (void)fputs("[FAILED]\n", stdout); + return -1; + } + if (fputs("[SUCCESS]\n", stdout) < 0) + return -1; + return 0; +} + +#define unit_test_statement(T, EXPR) do { (void)(T); EXPR; _unit_test_statement(( #EXPR)); } while (0) +#define unit_test_start() _unit_test_start(__FILE__, __func__, __LINE__) +#define unit_test(T, EXPR) _unit_test((T), 0 == (EXPR), (# EXPR), __FILE__, __func__, __LINE__, 0) +#define unit_test_verify(T, EXPR) _unit_test((T), 0 == (EXPR), (# EXPR), __FILE__, __func__, __LINE__, 1) + +static int test_sanity(void) { + unit_test_t t = unit_test_start(); + q_t t1 = 0, t2 = 0; + unit_test_statement(&t, t1 = qint(1)); + unit_test_statement(&t, t2 = qint(2)); + unit_test(&t, qtoi(qadd(t1, t2)) == 3); + return unit_test_finish(&t); +} + +static int test_pack(void) { + unit_test_t t = unit_test_start(); + q_t p1 = 0, p2 = 0; + char buf[sizeof(p1)] = { 0 }; + unit_test_statement(&t, p1 = qnegate(qinfo.pi)); + unit_test_statement(&t, p2 = 0); + unit_test(&t, qunequal(p1, p2)); + unit_test(&t, qpack(&p1, buf, sizeof p1 - 1) < 0); + unit_test(&t, qpack(&p1, buf, sizeof buf) == sizeof(p1)); + unit_test(&t, qunpack(&p2, buf, sizeof buf) == sizeof(p2)); + unit_test(&t, qequal(p1, p2)); + return unit_test_finish(&t); +} + +static int test_fma(void) { + unit_test_t t = unit_test_start(); + + q_t a, b, c, r; + const q_t one_and_a_half = qdiv(qint(3), qint(2)); + + /* incorrect, but expected, result due to saturation */ + unit_test_statement(&t, a = qinfo.max); + unit_test_statement(&t, b = one_and_a_half); + unit_test_statement(&t, c = qinfo.min); + unit_test_statement(&t, r = qadd(qmul(a, b), c)); + unit_test(&t, qwithin_interval(r, qint(0), qint(1))); + + /* correct result using Fused Multiply Add */ + unit_test_statement(&t, a = qinfo.max); + unit_test_statement(&t, b = one_and_a_half); + unit_test_statement(&t, c = qinfo.min); + unit_test_statement(&t, r = qfma(a, b, c)); + unit_test(&t, qwithin_interval(r, qdiv(qinfo.max, qint(2)), qint(1))); + + return unit_test_finish(&t); +} + +static inline int test_filter(void) { + unit_test_t t = unit_test_start(); + qfilter_t lpf = { .raw = 0 }, hpf = { .raw = 0 }; + const q_t beta = qdiv(qint(1), qint(3)); + qfilter_init(&lpf, qint(0), beta, qint(0)); + qfilter_init(&hpf, qint(0), beta, qint(0)); + for (int i = 0; i < 100; i++) { + char low[64 + 1] = { 0, }, high[64 + 1] = { 0, }; + const q_t step = qdiv(qint(i), qint(100)); + const q_t input = qint(1); + qfilter_low_pass(&lpf, step, input); + qfilter_high_pass(&hpf, step, input); + if (qsprint(qfilter_value(&lpf), low, sizeof(low) - 1ull) < 0) return -1; + if (qsprint(qfilter_value(&hpf), high, sizeof(high) - 1ull) < 0) return -1; + if (fprintf(stdout, "%2d: %s\t%s\n", i, low, high) < 0) return -1; + } + return unit_test_finish(&t); +} + +static int qmatrix_print(FILE *out, const q_t *m) { + assert(out); + assert(m); + const size_t alloc = qmatrix_string_length(m); + char *ms = malloc(alloc + 1); + if (!ms) + return -1; + int r = qmatrix_sprintb(m, ms, alloc + 1, 10); + if (r >= 0) + r = fprintf(out, "%s\n", ms); + free(ms); + return r; +} + +#define QMATRIX(ROW, COLUMN, ...) { 0, ROW * COLUMN, ROW, COLUMN, __VA_ARGS__ } +#define QMATRIXZ(ROW, COLUMN) QMATRIX((ROW), (COLUMN), 0) +#define QMATRIXSZ(ROW, COLUMN) ((((ROW)*(COLUMN)) + 4)*sizeof(q_t)) + +static int test_matrix(void) { + unit_test_t t = unit_test_start(); + FILE *out = stdout; + q_t a[] = QMATRIX(2, 3, + QINT(1), QINT(2), QINT(3), + QINT(4), QINT(5), QINT(6), + ); + q_t b[] = QMATRIX(3, 2, + QINT(2), QINT(3), + QINT(4), QINT(5), + QINT(6), QINT(7), + ); + const q_t abr[] = QMATRIX(2, 2, + QINT(28), QINT(34), + QINT(64), QINT(79), + ); + const q_t abrp[] = QMATRIX(2, 2, + QINT(28), QINT(64), + QINT(34), QINT(79), + ); + q_t ab[QMATRIXSZ(2, 2)] = QMATRIXZ(2, 2); + q_t abp[QMATRIXSZ(2, 2)] = QMATRIXZ(2, 2); + unit_test_verify(&t, 0 == qmatrix_mul(ab, a, b)); + unit_test_verify(&t, 0 == qmatrix_transpose(abp, ab)); + unit_test(&t, qmatrix_equal(ab, abr)); + unit_test(&t, qmatrix_equal(ab, abrp)); + qmatrix_print(out, a); + qmatrix_print(out, b); + qmatrix_print(out, ab); + qmatrix_print(out, abp); + return unit_test_finish(&t); +} + +static int test_matrix_trace(void) { + unit_test_t t = unit_test_start(); + q_t a[] = QMATRIX(2, 2, + QINT(1), QINT(2), + QINT(4), QINT(5), + ); + q_t b[] = QMATRIX(2, 2, + QINT(2), QINT(3), + QINT(4), QINT(5), + ); + q_t ta[sizeof(a) / sizeof(a[0])] = QMATRIX(2, 2, 0); + q_t tb[sizeof(b) / sizeof(b[0])] = QMATRIX(2, 2, 0); + q_t apb[sizeof(a) / sizeof(a[0])] = QMATRIX(2, 2, 0); + BUILD_BUG_ON(sizeof a != sizeof ta); + BUILD_BUG_ON(sizeof a != sizeof tb); + BUILD_BUG_ON(sizeof a != sizeof b); + BUILD_BUG_ON(sizeof a != sizeof apb); + unit_test_verify(&t, 0 == qmatrix_transpose(ta, a)); + unit_test_verify(&t, 0 == qmatrix_transpose(tb, b)); + unit_test_verify(&t, 0 == qmatrix_add(apb, a, b)); + unit_test(&t, qequal(qmatrix_trace(a), QINT(6))); + unit_test(&t, qequal(qmatrix_trace(b), QINT(7))); + unit_test(&t, qequal(qmatrix_trace(a), qmatrix_trace(ta))); + unit_test(&t, qequal(qmatrix_trace(apb), qadd(qmatrix_trace(a), qmatrix_trace(b)))); + printq(stdout, qmatrix_determinant(a), "det(a)"); + return unit_test_finish(&t); +} + +static q_t qid(q_t x) { return x; } +static q_t qsq(q_t x) { return qmul(x, x); } + +static int test_simpson(void) { + unit_test_t t = unit_test_start(); + unit_test(&t, qwithin_interval(qsimpson(qid, QINT(0), QINT(10), 100), QINT(50), QINT(1))); // (x^2)/2 + C + unit_test(&t, qwithin_interval(qsimpson(qsq, qnegate(QINT(2)), QINT(5), 100), QINT(44), QINT(1))); // (x^3)/3 + C + return unit_test_finish(&t); +} + +static int internal_tests(void) { + typedef int (*unit_test_t)(void); + unit_test_t tests[] = { + test_sanity, + test_pack, + test_fma, + // test_filter, + test_matrix, + test_matrix_trace, + test_simpson, + NULL + }; + for (size_t i = 0; tests[i]; i++) + if (tests[i]() < 0) + return -1; + return 0; +} + +static int help(FILE *out, const char *arg0) { + assert(out); + assert(arg0); + const char *h = "\n\ +Program: Q-Number (Q16.16, signed) library testbench\n\ +Author: Richard James Howe (2018)\n\ +License: MIT\n\ +E-mail: howe.r.j.89@gmail.com\n\ +Site: https://github.com/howerj/q\n\n\ +Options:\n\ +\t-s\tprint a sine-cosine table\n\ +\t-h\tprint this help message and exit\n\ +\t-i\tprint library information\n\ +\t-t\trun internal unit tests\n\ +\t-v\tprint version information and exit\n\ +\tfile\texecute commands in file\n\n\ +This program wraps up a Q-Number library and allows it to be tested by\n\ +giving it commands, either from stdin, or from a file. The format is\n\ +strict and the parser primitive, but it is only meant to be used to\n\ +test that the library is doing the correct thing and not as a\n\ +calculator.\n\n\ +Commands evaluated consist of an operator with the require arguments\n\ +(either one or two arguments) and compare the result with an expected\n\ +value, which can within a specified bounds for the test to pass. If\n\ +the test fails the program exits, indicating failure. The format is:\n\ +\n\ +\toperator expected +- allowance | arg1 arg2\n\n\ +operators include '+', '-', '/', 'rem', '<', which require two\n\ +arguments, and unary operators like 'sin', and 'negate', which require\n\ +only one. 'expected', 'allowance', 'arg1' and 'arg2' are all fixed\n\ +numbers of the form '-12.45'. 'expected' is the expected result,\n\ +'allowance' the +/- amount the result is allowed to deviated by, and\n\ +'arg1' and 'arg2' the operator arguments.\n\ +\n\n"; + if (fprintf(out, "usage: %s -h -s -i -v -t -c file\n", arg0) < 0) return -1; + if (fputs(h, out) < 0) return -1; + return 0; +} + +int main(int argc, char **argv) { + bool ran = false; + for (int i = 1; i < argc; i++) { + if (!strcmp("-h", argv[i])) { + if (help(stdout, argv[0]) < 0) + return 1; + return 0; + } else if (!strcmp("-s", argv[i])) { + print_sincos_table(stdout); + ran = true; + } else if (!strcmp("-v", argv[i])) { + fprintf(stdout, "version 1.0\n"); + return 0; + } else if (!strcmp("-t", argv[i])) { + if (internal_tests() < 0) + return 1; + ran = true; + } else if (!strcmp("-i", argv[i])) { + if (qinfo_print(stdout, &qinfo) < 0) + return 1; + if (qconf_print(stdout, &qconf) < 0) + return 1; + ran = true; + } else { + FILE *input = fopen_or_die(argv[i], "rb"); + FILE *output = stdout; + const int r = eval_file(input, output); + ran = true; + fclose(input); + if (r < 0) + return 1; + } + } + if (!ran) + return eval_file(stdin, stdout); + return 0; +} +