Емельянов Эдуард Владимирович (eddy_em) wrote,
Емельянов Эдуард Владимирович
eddy_em

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

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

Для генерирования таблицы с заданной точностью я сделал вот такой велосипед:
#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;
}

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

А вот этот блок вычисляет синусы-косинусы-тангенсы заданного угла:
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <stdint.h>
#include <limits.h>

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

Простые тесты:
for angle in $(seq 0 20 80); do ./calcsin $angle; echo ""; done
Approximate sin of 0 (0) = 0, exact = 0
Approximate cos of 0 (0) = 0.9999, exact = 1
Approximate tan of 0 (0) = 0, exact = 0

Approximate sin of 0.349066 (3490) = 0.3419, exact = 0.34202
Approximate cos of 0.349066 (3490) = 0.9396, exact = 0.939693
Approximate tan of 0.349066 (3490) = 0.3638, exact = 0.36397

Approximate sin of 0.698132 (6981) = 0.6427, exact = 0.642788
Approximate cos of 0.698132 (6981) = 0.7659, exact = 0.766044
Approximate tan of 0.698132 (6981) = 0.8391, exact = 0.8391

Approximate sin of 1.0472 (10471) = 0.8659, exact = 0.866025
Approximate cos of 1.0472 (10471) = 0.4998, exact = 0.5
Approximate tan of 1.0472 (10471) = 1.7324, exact = 1.73205

Approximate sin of 1.39626 (13962) = 0.9847, exact = 0.984808
Approximate cos of 1.39626 (13962) = 0.1736, exact = 0.173648
Approximate tan of 1.39626 (13962) = 5.6722, exact = 5.67128

for angle in $(seq 80 90); do ./calcsin $angle; echo ""; done
Approximate sin of 1.39626 (13962) = 0.9847, exact = 0.984808
Approximate cos of 1.39626 (13962) = 0.1736, exact = 0.173648
Approximate tan of 1.39626 (13962) = 5.6722, exact = 5.67128

Approximate sin of 1.41372 (14137) = 0.9876, exact = 0.987688
Approximate cos of 1.41372 (14137) = 0.1563, exact = 0.156434
Approximate tan of 1.41372 (14137) = 6.3186, exact = 6.31375

Approximate sin of 1.43117 (14311) = 0.9902, exact = 0.990268
Approximate cos of 1.43117 (14311) = 0.139, exact = 0.139173
Approximate tan of 1.43117 (14311) = 7.1237, exact = 7.11537

Approximate sin of 1.44862 (14486) = 0.9925, exact = 0.992546
Approximate cos of 1.44862 (14486) = 0.1217, exact = 0.121869
Approximate tan of 1.44862 (14486) = 8.1552, exact = 8.14435

Approximate sin of 1.46608 (14660) = 0.9944, exact = 0.994522
Approximate cos of 1.46608 (14660) = 0.1045, exact = 0.104528
Approximate tan of 1.46608 (14660) = 9.5157, exact = 9.51436

Approximate sin of 1.48353 (14835) = 0.9961, exact = 0.996195
Approximate cos of 1.48353 (14835) = 0.0871, exact = 0.0871557
Approximate tan of 1.48353 (14835) = 11.4362, exact = 11.4301

Approximate sin of 1.50098 (15009) = 0.9975, exact = 0.997564
Approximate cos of 1.50098 (15009) = 0.0697, exact = 0.0697565
Approximate tan of 1.50098 (15009) = 14.3113, exact = 14.3007

Approximate sin of 1.51844 (15184) = 0.9985, exact = 0.99863
Approximate cos of 1.51844 (15184) = 0.0522, exact = 0.052336
Approximate tan of 1.51844 (15184) = 19.1283, exact = 19.0811

Approximate sin of 1.53589 (15358) = 0.9993, exact = 0.999391
Approximate cos of 1.53589 (15358) = 0.0348, exact = 0.0348995
Approximate tan of 1.53589 (15358) = 28.7155, exact = 28.6363

Approximate sin of 1.55334 (15533) = 0.9997, exact = 0.999848
Approximate cos of 1.55334 (15533) = 0.0173, exact = 0.0174524
Approximate tan of 1.55334 (15533) = 57.7861, exact = 57.29

Approximate sin of 1.5708 (15707) = 0.9999, exact = 1
Approximate cos of 1.5708 (15707) = 0, exact = 6.12323e-17
Approximate tan of 1.5708 (15707) = 214748, exact = 1.63312e+16

Понятно, что точность хуже заданной 0.0001, т.к. вычисления — целочисленные, но три знака у синусов-косинусов вполне честные. Котангенс для углов ±90° возвращает INT_MAX как индикатор ошибки. По понятной причине точность котангенсов на углах вблизи 90° сильно ухудшается, если понадобится более точный тангенс, можно для него отдельные таблицы сделать (но значений уже будет намного больше).
Две таблицы по 61 значению для точности не хуже ~0.0004 — вполне приемлемо на мой взгляд. Кстати, для повышения точности на порядок требуется уже 191 значение, а на два порядка — аж 634!
Tags: c, всячина
Subscribe

  • Контроллер вентиляторов

    Потихоньку допиливаю контроллер. Забыл сегодня на работе взять платиновый терморезистор для аттестации терморезисторов, которые там используются…

  • "Звезды" на STM8

    Я уже писал о том, что на новогодний концерт сделал управление "звездами" на основе STM8. Сам концерт можете посмотреть на youtube (на всякий…

  • Андроидоуправлялка

    Хоть ардуина - редкостная дрянь, но благодаря ей появилась уйма приложений под андроид, позволяющих через bluetooth-uart модуль контролировать разную…

promo eddy_em august 17, 2019 12:33 3
Buy for 10 tokens
Юра намедни напечатал корпус для хронометра. Для первого блина получилось неплохо: И еще немного фотографий:
  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic
  • 0 comments