mirror of
https://github.com/eddyem/eddys_snippets.git
synced 2025-12-06 10:45:12 +03:00
110 lines
3.3 KiB
C
110 lines
3.3 KiB
C
/*
|
|
* gentab.c - generate sin/cos table
|
|
*
|
|
* Copyright 2016 Edward V. Emelianov <eddy@sao.ru, edward.emelianoff@gmail.com>
|
|
*
|
|
* 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 <stdio.h>
|
|
#include <math.h>
|
|
#include <stdlib.h>
|
|
#include <stdint.h>
|
|
|
|
#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;
|
|
}
|