Поиск

Очередной велосипед для вычисления синусов-косинусов


Раздумывая о повышении точности табличного вычисления синусов-косинусов с попутным уменьшением размера используемой таблицами флеш-памяти, я подумал: почему бы не сделать таблицу неравномерной? И вот, что получилось…

Для генерирования таблицы с заданной точностью я сделал вот такой велосипед:

#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 precntwhere 'prec' is desired precisionn", argv[0]);
		return 1;
	}
	char *eptr;
	double prec = strtod(argv[1], &eptr);
	if(*eptr || eptr == argv[1]){
		fprintf(stderr, "Bad number: %sn", 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 %dn", L, MAX_LEN);
		return 3;
	}
	calc_table(prec, angles, sinuses);
	printf("%s angles[%d] = {", DATATYPE, L);
	for(i = 0; i < L; ++i){
		printf("%ld, ", *angles++);
	}
	printf("bb");
	printf("};n%s sinuses[%d] = {", DATATYPE, L);
	for(i = 0; i < L; ++i){
		printf("%ld, ", *sinuses++);
	}
	printf("bb");
	printf("};n");
	return 0;
}


Понятно, что для удобства работы на мелкоконтроллерах лучше точность задавать кратную 0.1, чтобы удобней было работать с фиксированной точкой визуально (да и на экранчики проще выводить).

Полностью все тригонометрические функции я пока не реализовал — отложу это до необходимости, пока только сделал проверку на синусах в первой четверти:

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <stdint.h>

int32_t 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};
int32_t 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};

#define SZ (61)

int32_t calcsin(int32_t rad){
	// find interval for given angle
	int32_t 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]=%dn", iter, ind, angles[ind], angles[ind+1]);
	int32_t ai = angles[ind], si = sinuses[ind];
	if(ai == rad) return si;
	++ind;
	if(angles[ind] == rad) return sinuses[ind];
	return (si + (sinuses[ind] - si)*(rad - ai)/(angles[ind] - ai));
}

int main(int argc, char **argv){
	if(argc != 2) return 1;
	double ang = strtod(argv[1], NULL), drad = ang * M_PI/180.;
	//printf("rad: %gn", drad);
	int32_t rad = (int32_t)(drad * 10000);
	printf("Approximate sin of %gdegr = %d, exact = %gn",
		ang, calcsin(rad), sin(drad));
	return 0;
}


Простые тесты:

for angle in $(seq 0 5 90); do ./calcsin $angle; done
Approximate sin of 0degr = 0, exact = 0
Approximate sin of 5degr = 871, exact = 0.0871557
Approximate sin of 10degr = 1736, exact = 0.173648
Approximate sin of 15degr = 2587, exact = 0.258819
Approximate sin of 20degr = 3419, exact = 0.34202
Approximate sin of 25degr = 4225, exact = 0.422618
Approximate sin of 30degr = 4997, exact = 0.5
Approximate sin of 35degr = 5735, exact = 0.573576
Approximate sin of 40degr = 6427, exact = 0.642788
Approximate sin of 45degr = 7069, exact = 0.707107
Approximate sin of 50degr = 7659, exact = 0.766044
Approximate sin of 55degr = 8191, exact = 0.819152
Approximate sin of 60degr = 8659, exact = 0.866025
Approximate sin of 65degr = 9062, exact = 0.906308
Approximate sin of 70degr = 9396, exact = 0.939693
Approximate sin of 75degr = 9658, exact = 0.965926
Approximate sin of 80degr = 9847, exact = 0.984808
Approximate sin of 85degr = 9961, exact = 0.996195
Approximate sin of 90degr = 9999, exact = 1

На угле 30° виден косяк: разница с нормальным значением составляет целых 0.0003 вместо положенных 0.0001! Где-то в моих вычислениях точки максимального отклонения вышел косяк. Ну, да ладно. Таблицы из 121 значения для точности не хуже трех десятитысячных — вполне приемлемо на мой взгляд.
Кстати, для повышения точности на порядок требуется уже 191 значение, а на два порядка — аж 634!

eddy-em.livejournal.com

Добавить комментарий