From d82692516259ce539d841640555633b1e1b440a9 Mon Sep 17 00:00:00 2001 From: eddyem Date: Mon, 20 Jun 2016 14:40:26 +0300 Subject: [PATCH] Add trigonometry functions for calculations on MCU --- getopt/cmdlnopts/parseargs.c | 4 +- table_sincos/Makefile | 22 +++++++ table_sincos/README | 13 +++++ table_sincos/calcsin.c | 87 ++++++++++++++++++++++++++++ table_sincos/gentab.c | 109 +++++++++++++++++++++++++++++++++++ table_sincos/sqrt.c | 43 ++++++++++++++ 6 files changed, 276 insertions(+), 2 deletions(-) create mode 100644 table_sincos/Makefile create mode 100644 table_sincos/README create mode 100644 table_sincos/calcsin.c create mode 100644 table_sincos/gentab.c create mode 100644 table_sincos/sqrt.c diff --git a/getopt/cmdlnopts/parseargs.c b/getopt/cmdlnopts/parseargs.c index 2f57be9..64bfab1 100644 --- a/getopt/cmdlnopts/parseargs.c +++ b/getopt/cmdlnopts/parseargs.c @@ -233,7 +233,7 @@ void parseargs(int *argc, char ***argv, myoption *options){ loptr->flag = opts->flag; loptr->val = opts->val; // fill short options if they are: - if(!opts->flag){ + if(!opts->flag && opts->val){ #ifdef EBUG shortlist[i] = (char) opts->val; #endif @@ -396,7 +396,7 @@ void showhelp(int oindex, myoption *options){ qsort(opts, N, sizeof(myoption), argsort); do{ int p = sprintf(buf, " "); // a little indent - if(!opts->flag) // .val is short argument + if(!opts->flag && opts->val) // .val is short argument p += snprintf(buf+p, bufsz-p, "-%c, ", opts->val); p += snprintf(buf+p, bufsz-p, "--%s", opts->name); if(opts->has_arg == 1) // required argument diff --git a/table_sincos/Makefile b/table_sincos/Makefile new file mode 100644 index 0000000..1ce273e --- /dev/null +++ b/table_sincos/Makefile @@ -0,0 +1,22 @@ +PROGRAM = gentab +LDFLAGS = -lm +SRCS = gentab.c +CC = gcc +DEFINES = -D_XOPEN_SOURCE=1111 +CXX = gcc +CFLAGS = -Wall -Werror -Wextra $(DEFINES) +OBJS = $(SRCS:.c=.o) +all : $(PROGRAM) +$(PROGRAM) : $(OBJS) + $(CC) $(CFLAGS) $(OBJS) $(LDFLAGS) -o $(PROGRAM) + +# some addition dependencies +# %.o: %.c +# $(CC) $(LDFLAGS) $(CFLAGS) $< -o $@ +#$(SRCS) : %.c : %.h $(INDEPENDENT_HEADERS) +# @touch $@ + +clean: + /bin/rm -f *.o *~ +depend: + $(CXX) -MM $(CXX.SRCS) diff --git a/table_sincos/README b/table_sincos/README new file mode 100644 index 0000000..43c01e8 --- /dev/null +++ b/table_sincos/README @@ -0,0 +1,13 @@ +Generate table and functions for calculation of sin and cos on MCU (in integers) + +compile: make +defines: +-DMAX_LEN=val - maximal length of generated array (128 by default) +-DDATATYPE=\"type\" - type of data (\"int32_t\" by default) + + +run: +gentab number + where 'number' is desired precision - step of sin/cos (for example, 0.001) + sin is integer, for "1" if precision is 0.001 would be 1000 + diff --git a/table_sincos/calcsin.c b/table_sincos/calcsin.c new file mode 100644 index 0000000..44bb24a --- /dev/null +++ b/table_sincos/calcsin.c @@ -0,0 +1,87 @@ +/* + * calcsin.c + * + * Copyright 2016 Edward V. Emelianov + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#include +#include +#include +#include +#include + +#define ONE_FIXPT (10000) +#define SZ (61) +typedef int32_t TYPE; +#define _PI (31415) +#define _2PI (62831) +#define _PI_2 (15707) +TYPE angles[61] = {0, 966, 1700, 2025, 2331, 2622, 2898, 3163, 3418, 3663, 3900, 4129, 4351, 4566, 4982, 5375, 5747, 6102, 6441, 6764, 7073, 7369, 7652, 7925, 8186, 8438, 8680, 8912, 9136, 9352, 9560, 9761, 9954, 10141, 10321, 10495, 10663, 10825, 10982, 11133, 11425, 11698, 11953, 12191, 12414, 12622, 12817, 12999, 13170, 13329, 13479, 13759, 14003, 14217, 14404, 14567, 14710, 14959, 15147, 15427, 15708}; +TYPE sinuses[61] = {0, 965, 1692, 2011, 2310, 2592, 2858, 3111, 3352, 3582, 3802, 4013, 4215, 4409, 4778, 5119, 5436, 5731, 6004, 6260, 6498, 6720, 6927, 7121, 7302, 7472, 7630, 7778, 7917, 8047, 8169, 8283, 8390, 8490, 8584, 8672, 8754, 8831, 8904, 8972, 9097, 9207, 9303, 9388, 9462, 9528, 9585, 9635, 9680, 9718, 9753, 9811, 9855, 9889, 9915, 9935, 9950, 9972, 9984, 9996, 10000}; + +TYPE calcsin(TYPE rad){ + rad %= _2PI; + if(rad < 0) rad += _2PI; + TYPE sign = 1; + if(rad > _PI){ + rad -= _PI; + sign = -1; + } + if(rad > _PI_2) rad = _PI - rad; + // find interval for given angle + int ind = 0, half = SZ/2, end = SZ - 1, iter=0; + do{ + ++iter; + if(angles[half] > rad){ // left half + end = half; + }else{ // right half + ind = half; + } + half = ind + (end - ind) / 2; + }while(angles[ind] > rad || angles[ind+1] < rad); + //printf("%d iterations, ind=%d, ang[ind]=%d, ang[ind+1]=%d\n", iter, ind, angles[ind], angles[ind+1]); + TYPE ai = angles[ind], si = sinuses[ind]; + if(ai == rad) return si; + ++ind; + if(angles[ind] == rad) return sinuses[ind]; + return sign*(si + (sinuses[ind] - si)*(rad - ai)/(angles[ind] - ai)); +} + +TYPE calccos(TYPE rad){ + return calcsin(_PI_2 - rad); +} + +TYPE calctan(TYPE rad){ + TYPE s = calcsin(rad) * ONE_FIXPT, c = calccos(rad); + if(c) return s/calccos(rad); + else return INT_MAX; +} + +int main(int argc, char **argv){ + if(argc != 2) return 1; + double ang = strtod(argv[1], NULL), drad = ang * M_PI/180.; + //printf("rad: %g\n", drad); + TYPE rad = (TYPE)(drad * ONE_FIXPT); + printf("Approximate sin of %g (%d) = %g, exact = %g\n", + drad, rad, calcsin(rad)/((double)ONE_FIXPT), sin(drad)); + printf("Approximate cos of %g (%d) = %g, exact = %g\n", + drad, rad, calccos(rad)/((double)ONE_FIXPT), cos(drad)); + printf("Approximate tan of %g (%d) = %g, exact = %g\n", + drad, rad, calctan(rad)/((double)ONE_FIXPT), tan(drad)); + return 0; +} diff --git a/table_sincos/gentab.c b/table_sincos/gentab.c new file mode 100644 index 0000000..7c8d2f6 --- /dev/null +++ b/table_sincos/gentab.c @@ -0,0 +1,109 @@ +/* + * gentab.c - generate sin/cos table + * + * Copyright 2016 Edward V. Emelianov + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ +#include +#include +#include +#include + +#define MAX_LEN (256) +#define DATATYPE "int32_t" + +/** + * Calculate angle \alpha_{i+1} for given precision + * @param ai - \alpha_i + * @param prec - precision + * @return calculated angle + */ +double get_next_angle(double ai, double prec){ + double anxt, amid = M_PI_2, diff = 1.; +int i = 0; + do{ + anxt = amid; + double si = sin(ai), snxt = sin(anxt); + /* + * We have: sin(a[i]+a) \approx sin(a[i]) + (sin(a[i+1])-sin(a[i]))/(a[i+1]-a[i])*a + * so maximal deviation from real sin(a) would be at point + * amax = acos((sin(a[i+1]) - sin(a[i]))/(a[i+1] - a[i])) + * to calculate a[i+1] we always start from pi/2 + */ + amid = acos((snxt - si)/(anxt - ai)); + double sinapp = si + (snxt - si)/(anxt - ai)*(amid - ai); + diff = fabs(sinapp - sin(amid)); + ++i; + }while(diff > prec); + return anxt; +} + +/** + * calculate sin/cos table + * sin(a) = s1 + (s2 - s1)*(a - a1)/(a2 - a1) + * @param ang - NULL (to calculate tabular length only) or array of angles + * @param sin - NULL (-//-) or array of sinuses from 0 to 1 + */ +int calc_table(double prec, uint64_t *ang, uint64_t *sinuses){ + int L = 1; // points 0 and pi/2 included! + double ai = 0.; + if(ang) *ang++ = 0; + if(sinuses) *sinuses++ = 0; + while(M_PI_2 > ai){ + ai = get_next_angle(ai, prec); + if(ang) *ang++ = (uint64_t)round(ai/prec); + if(sinuses) *sinuses++ = (uint64_t)round(sin(ai)/prec); + ++L; + } + return L; +} + +int main(int argc, char **argv){ + if(argc != 2){ + fprintf(stderr, "Usage: %s prec\n\twhere 'prec' is desired precision\n", argv[0]); + return 1; + } + char *eptr; + double prec = strtod(argv[1], &eptr); + if(*eptr || eptr == argv[1]){ + fprintf(stderr, "Bad number: %s\n", argv[1]); + return 2; + } + uint64_t *angles = malloc(sizeof(uint64_t)*MAX_LEN); + uint64_t *sinuses = malloc(sizeof(uint64_t)*MAX_LEN); + int L = calc_table(prec, NULL, NULL), i; + if(L > MAX_LEN){ + fprintf(stderr, "Error! Get vector of length %d instead of %d\n", L, MAX_LEN); + return 3; + } + calc_table(prec, angles, sinuses); + printf("#define ONE_FIXPT (%d)\n#define SZ (%d)\ntypedef %s TYPE;\n", (int)(1./prec), L, DATATYPE); + printf("#define _PI (%ld)\n#define 2_PI (%ld)\n#define _PI_2 (%ld)\n\n", (int64_t)(M_PI/prec), (int64_t)(2*M_PI/prec), + (int64_t)(M_PI/2/prec)); + printf("TYPE angles[%d] = {", L); + for(i = 0; i < L; ++i){ + printf("%ld, ", *angles++); + } + printf("\b\b"); + printf("};\nTYPE sinuses[%d] = {", L); + for(i = 0; i < L; ++i){ + printf("%ld, ", *sinuses++); + } + printf("\b\b"); + printf("};\n"); + return 0; +} diff --git a/table_sincos/sqrt.c b/table_sincos/sqrt.c new file mode 100644 index 0000000..d741cf3 --- /dev/null +++ b/table_sincos/sqrt.c @@ -0,0 +1,43 @@ +/* + * sqrt.c - Newton sqrt + * + * Copyright 2016 Edward V. Emelianov + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ +#define _XOPEN_SOURCE 1111 +#include +#include +#include + +float SQRT(float x, int iter){ + float Result; + int i = 0; + Result = x; + for (i; i < iter; i++) + Result = (Result + x / Result) / 2; + return Result; +} + +int main(int argc, char **argv){ + if(argc != 2) return 1; + float x = strtof(argv[1], NULL); + if(x < 0) return 2; + printf("sqr(%f) = %f, by 6iter: %f, by 8iter: %f, by 16iter: %f\n", x, sqrtf(x), + SQRT(x, 6), SQRT(x, 8), SQRT(x, 16)); + return 0; +} +