LINUX.ORG.RU

вычисление тригонометрии - как?! как вы это делаете?

 ,


1

4

собственно вопрос. Я ни разу не программист, но тут что-то наткнулся на одно старое обсуждение тыц и у меня возник вопрос - собственно как это всё происходит.

Интересует сам механизм, достал свой пыльный конспект по матану всю ночь читал, пытался вспомнить, под утро наиндусил на C разложение в ряд Тейлора(в три итерации) и сравнил выхлоп с системной функцией... и вот вопрос таки возник.

Получилось, что на малых углах разбег начинается после второго знака, но уже к 30° он перескакивает на первый и далее только растёт. Оно понятно что я тот ещё индус, но хотелось бы послушать опыт знающих.(glibc я конечно почитаю потом).

конечный результат системного   косинуса  10:   0.984808
конечный результат вычисленного косинуса  10:   0.984731

конечный результат системного   косинуса  20:   0.939693
конечный результат вычисленного косинуса  20:   0.938460

конечный результат системного   косинуса  30:   0.866025
конечный результат вычисленного косинуса  30:   0.859819

сам код :

#include <stdio.h>
#include <math.h>
#define PI 3.14159265
double  rezult, deg_0,deg,test;
double factor (double iterration);
double stepen (double deg, double step);
double prom( double st, double faktor);
int main()
{
    printf ("введите значение: ");
    scanf("%lf", &deg_0);
    deg=deg_0*(PI/180);
    printf(" \n");
    printf("Processing... \n");
    double q;
    q=2;
    double promt;
    promt=prom(deg,q);
    q=q+2;
    promt=promt+prom(deg,q);
    q=q+2;
    promt=promt-prom(deg,q);
    q=q+2;
    promt=promt+prom(deg,q);
    rezult=1-promt;
    test=cos(deg);
    printf("конечный результат системного   косинуса %10f: %10f\n",deg_0, test);
    printf("конечный результат вычисленного косинуса %10f: %10f\n",deg_0, rezult);
    return 0;
}
double factor (double iterration)
{
    if(iterration < 0)
        return 0;
    if (iterration == 0)
        return 1;
    else
        return (iterration < 2) ? 1 : iterration * factor(iterration - 1);
}

double stepen (double deg, double step)
{
    if(step < 0) step = -step;
    if(step == 0) step = 1;
    double result = 1;
    for (double i = 0; i < step; i++) {
        result *= deg;
    }
    return result;
}
double prom(double st, double faktor)
{
 double znachenie;
  double i,j;
  i=stepen(st,faktor);
  j=factor(faktor);
  znachenie=i/j;
  return znachenie;
}
за код прошу сильно не пинать, я не программер, наиндусил на коленке чтоб проверить сам принцип.

Короче кто как это решает?

PS. интересует именно голый алгоритм,за системные инструменты glibc я в курсе.

★★

Последнее исправление: Gramozeka (всего исправлений: 1)

Дочитай конспект до оценки погрешности вычислительных методов. В частности, что априорная оценка для ряда тейлора, что точно уже не помню, но как минимум погрешность не менее значения последнего вычисленного члена ряда.

anonymous
()
Ответ на: комментарий от anonymous

Дочитай конспект до оценки погрешности вычислительных методов..

хорошо, если найду, почитаю, но по ходу эта часть конспекта была выкурена на каком-то весёлом новом году... но у меня сохранился учебник!(честно купленный примерно тогда же, когда конспект скурили).

Gramozeka ★★
() автор топика

На, дарю:

/*****************************************************************
               sin(x) evaluation routine
 If x >  EPS value, the following recursion is used
    sin(x) = 3*sin(x/3)-4*sin(x/3)^3.
 When x < EPS, a part of Taylor series for sin(x) is employed:
    sin(x) = x - x^3/6 + O(x^4), O<=1/24
 EPS is chosen so that the remainder is  less than.
 macheps = 5.551115e-17 (for double)

*****************************************************************/

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

const double EPS = 2.0e-4;

double fabs (double x) {return (x>0? x: -x);}

double
      p_sin3x(double x)
      {
       return x*(4*(-x+1)*(x+1)-1);
      }

double
      sin(double x)
      {
       if(fabs(x)>EPS)
        return p_sin3x(sin(x/3));
       else
        return x*((-x+2)*(x+2)+2)/6;
       }

int
    main(int argc, char* argv[])
    {
     if(argc>1)
      printf("%.17e\n", sin(atof(argv[1])));
     else
      printf("Usage: %s x\n", argv[0]);

     return 0;
    }

unanimous ★★★★★
()
Ответ на: комментарий от unanimous

ну окей, умыл(собсно другого то и не стоило ожидать), но простыми словами то можешь объяснить сам механизм? для чугуниевых чайников так сказать.

Gramozeka ★★
() автор топика
double factor (double iterration) {
...
if (iterration == 0)
...

Правильно. Нахера придумали эти все int'ы, когда везде можно ставить double. Не боишься, что препод оторвёт тебе руки в воспитательных целях? Загугли 0.1 + 0.2 (https://0.30000000000000004.com/)

q=q+2;

q+=2;
Вот же делаешь

result *= deg;

double factor (double iterration)
double stepen(double deg, double step)

В math что, совсем нет факториала и возведения в степень? Ну и функция полное говно, потому что степень может быть не целой и отрицательной по твоей же сигнатуре.

#define PI 3.14159265

Сразу бы делал 4, чо выделываться. Какую ты хочешь точность с таким pi вместо нормального кода? Хорошо, что он у тебя вообще, хоть что-то считает, а не рандомно сегфолтится.

crutch_master ★★★★★
()
Последнее исправление: crutch_master (всего исправлений: 5)
Ответ на: комментарий от crutch_master

Не боишься, что препод оторвёт тебе руки в воспитательных целях?

ты сделал мой день!))).

Ты сразу критиковать мой индуизм, а вопрос то не про него, я вообще в Си не умею от слова совсем, а последнего своего препода я видел тридцать с лишним лет назад(бухали на поминках декана).

Я спросил про алгоритм «как это делается», выше человек поделился своим кодом, правда тему не раскрыл, ну да ладно, я уже распарсил чё к чему. Ты если критикуешь - предлагай, но лучше словами опиши(можно без научностей, для чайников), я хоть и старый дурак, но понятливый, человечью речь понимаю ишо.

Gramozeka ★★
() автор топика
Последнее исправление: Gramozeka (всего исправлений: 1)
Ответ на: комментарий от Gramozeka

можешь объяснить сам механизм?

Только интересен не рекурсивный механизм, а использование обратной «лесенки». А так очень даже ничего.

anonymous
()
Ответ на: комментарий от Gramozeka

простыми словами то можешь объяснить сам механизм?

С достаточно маленьким Х тейлоровая серия колапсирует очень быстро (дележ на факториал + умножение на степень от числа достаточно близкого к нулю = двойной мега-дележ). По этому ему хватает вычислить только несколько первых членов серии чтоб получить довольно точный результат. По этому он сначала стремительно уменьшает большой Х, вычисляет очень точное значение, точности которого хватает и при обратном возврате на большой Х. Правда, я бы сначала Х к значению между 0 и 2*pi привел.

anonymous
()
Ответ на: комментарий от invy

спасибо за наколку, почитаю на досуге(но от английского уже воротит после полной сборки CBLFS, да со свежими глюками, да новыми бирюлькамикедами).

Вот немного мозги остынут, почитаю таки свои старые учебники, только маразм своё берёт - тупею активно...

Gramozeka ★★
() автор топика

Не стоит сравнивать double с нулем (можно получить не тот результат, на который рассчитываешь).

anonymous
()

Численные методы самарский в гугл, найдешь одну из 3 моих любимых книжек, должно решить все вопросы. Матан тут не поможет.

peregrine ★★★★★
()
Последнее исправление: peregrine (всего исправлений: 1)
Ответ на: комментарий от anonymous

Не стоит сравнивать double с нулем (можно получить не тот результат, на который рассчитываешь).

[зануда on]

Более корректная формулировка - проверять на равенство, неважно с нулем или другой величиной, а сравнивать можно всегда и хоть с чем. Но это уже серьёзная стадия, когда человек начинает понимать что такое число с плавающей запятой. Чуть-чуть осталось до следующей стадии - когда человек понимает, что в некоторых случаях сравнивать на равенство все-таки можно и это корректно.

[зануда off]

yetanother ★★
()
Последнее исправление: yetanother (всего исправлений: 1)
Ответ на: комментарий от Gramozeka

По таблице из справочника sin(0)=0 sin(30)=0.5 … До 90. Все остальное включая косинус и тангенс выводится.

invy ★★★★★
()

Вот из стандартной библиотеки go (это собирается при условии отсутствии реализации на ассемблере)

// sin coefficients
var _sin = [...]float64{
	1.58962301576546568060e-10, // 0x3de5d8fd1fd19ccd
	-2.50507477628578072866e-8, // 0xbe5ae5e5a9291f5d
	2.75573136213857245213e-6,  // 0x3ec71de3567d48a1
	-1.98412698295895385996e-4, // 0xbf2a01a019bfdf03
	8.33333333332211858878e-3,  // 0x3f8111111110f7d0
	-1.66666666666666307295e-1, // 0xbfc5555555555548
}

// cos coefficients
var _cos = [...]float64{
	-1.13585365213876817300e-11, // 0xbda8fa49a0861a9b
	2.08757008419747316778e-9,   // 0x3e21ee9d7b4e3f05
	-2.75573141792967388112e-7,  // 0xbe927e4f7eac4bc6
	2.48015872888517045348e-5,   // 0x3efa01a019c844f5
	-1.38888888888730564116e-3,  // 0xbf56c16c16c14f91
	4.16666666666665929218e-2,   // 0x3fa555555555554b
}


func sin(x float64) float64 {
	const (
		PI4A = 7.85398125648498535156e-1  // 0x3fe921fb40000000, Pi/4 split into three parts
		PI4B = 3.77489470793079817668e-8  // 0x3e64442d00000000,
		PI4C = 2.69515142907905952645e-15 // 0x3ce8469898cc5170,
	)
	// special cases
	switch {
	case x == 0 || IsNaN(x):
		return x // return ±0 || NaN()
	case IsInf(x, 0):
		return NaN()
	}

	// make argument positive but save the sign
	sign := false
	if x < 0 {
		x = -x
		sign = true
	}

	var j uint64
	var y, z float64
	if x >= reduceThreshold {
		j, z = trigReduce(x)
	} else {
		j = uint64(x * (4 / Pi)) // integer part of x/(Pi/4), as integer for tests on the phase angle
		y = float64(j)           // integer part of x/(Pi/4), as float

		// map zeros to origin
		if j&1 == 1 {
			j++
			y++
		}
		j &= 7                               // octant modulo 2Pi radians (360 degrees)
		z = ((x - y*PI4A) - y*PI4B) - y*PI4C // Extended precision modular arithmetic
	}
	// reflect in x axis
	if j > 3 {
		sign = !sign
		j -= 4
	}
	zz := z * z
	if j == 1 || j == 2 {
		y = 1.0 - 0.5*zz + zz*zz*((((((_cos[0]*zz)+_cos[1])*zz+_cos[2])*zz+_cos[3])*zz+_cos[4])*zz+_cos[5])
	} else {
		y = z + z*zz*((((((_sin[0]*zz)+_sin[1])*zz+_sin[2])*zz+_sin[3])*zz+_sin[4])*zz+_sin[5])
	}
	if sign {
		y = -y
	}
	return y
}
nikolnik ★★★
()

Получилось, что на малых углах разбег начинается после второго знака, но уже к 30° он перескакивает на первый и далее только растёт.

В ряд Тейлора вокруг нуля раскладывал? Маладец.

yvv ★★☆
()

Так кто тут в треде главный «тригонометр»?

anonymous
()
Ответ на: комментарий от hobbit

Ой, а расскажешь подробнее, как оно?

намана! Две недели неспешного канпеляния и напилинга - и у тебю «идеальная система» которая может всё в обе архитектуры))... а уж сколько говнокода перелопатишь пока соберёшь,.. на пол-жизни хватит впечатлений, да вон в профиле можешь глянуть, предыдущий вариант(ещё с 3.8 питоном и 10 шлангом) я зафиксировал. Но уже пересобрал под 3.9 и 11 шланг, да с новыми кедами(кстати таки торт), но там конечно не всё, в этот раз собрал себе криту и блендер на потыкать и ffmpeg наконец укомплектовал на 100%(скринкасты теперь пишутся без глюков, артефактов и тормозов, даже под взрослой нагрузкой).

...но это такая себе развлекуха, не каждому зайдёт, это даже не гента, это ещё упоротей.

Gramozeka ★★
() автор топика
Ответ на: комментарий от AdonisSexyBoy

Для гцц нашёл -fassociative-math, но результат такой себе:

	vmovsd	xmm2, QWORD PTR y[rip]
	vmulsd	xmm1, xmm2, QWORD PTR PI4C[rip]
	vmovsd	xmm0, QWORD PTR x[rip]
	vsubsd	xmm0, xmm0, xmm1
	vmovsd	xmm1, QWORD PTR PI4A[rip]
	vaddsd	xmm1, xmm1, QWORD PTR PI4B[rip]
	vmulsd	xmm1, xmm1, xmm2
	vsubsd	xmm0, xmm0, xmm1

против ручной оптимизации

	vmovsd	xmm0, QWORD PTR PI4A[rip]
	vaddsd	xmm0, xmm0, QWORD PTR PI4B[rip]
	vaddsd	xmm0, xmm0, QWORD PTR PI4C[rip]
	vmulsd	xmm0, xmm0, QWORD PTR y[rip]
	vmovsd	xmm1, QWORD PTR x[rip]
	vsubsd	xmm0, xmm1, xmm0

Но вообще, да, ты прав. Там даже есть упоминание, что оно violates the ISO C and C++ language standard by possibly changing computation result. Так что выходит, что мат. оптимизации компилятор сам по себе проводит.

anonymous
()

Вариантов масса:

  1. Ряд тейлора
  2. Интерполяции различными функциями (кто сказал полиномы Чебышева и кубические сплайны?) по минимаксному критерию (минимизируем максимальное отклонение от эталона на заданном интервале)
  3. CORDIC в конце концов!
shkolnick-kun ★★★★★
()

Оценка погрешности |x|^n/n!

Все.

Сходимость хоть какая-то начинается от n>2x (в радианах). С формулой Стирлинга можно получить оценку типа:

n lnx - n ln n + n = ln epsilon

anonymous
()

Бля, в радианах углы считай

anonymous
()
Ответ на: комментарий от peregrine

Компилятору об этом сказали? А то он эту портяночку то упростить может.

Написали же: если явно через -fassociative-math не просить, то не может.

monk ★★★★★
()
Ответ на: комментарий от anonymous

мат. оптимизации компилятор сам по себе НЕ проводит

да.

anonymous
()

Если интересует именно математика - то рядом Тейлора никто не считает, используется аппроксимация полиномами Чебышева или аппроксимация Чебышева-Паде (даёт меньшую относительную погрешность на малых значениях). Вот тут есть пример расчёта коэффициентов полинома в Maple.

Есть справочник посчитанных коэффициентов полиномов разных функций для разных степеней полиномов и, соответственно, разной точности, но название я не помню. Всегда можно посчитать коэффициенты для нужной функцию в Maple.

tim239 ★★
()
Ответ на: комментарий от Minona

остромысл типа?

Вот тебе мозг размять задачка(коли вумный такой):

дано - КПВТ установлен на платформе и пристрелян на заданное расстояние. Обстрел ведётся только в горизонтальной плоскости.

платформа управляется сервоприводом.

Задача - расчитать команду сервоприводу(на какой угол нужно повернуть червячный вал, при получении целеуказания(цель появляется в секторе ±60°). Считать что один полный оборот поворачивает платформу на 1°.

Можешь на гумажке разрисовать только алгоритм, можешь пользоваться брадисом..)))

Gramozeka ★★
() автор топика
Последнее исправление: Gramozeka (всего исправлений: 2)
Ответ на: комментарий от tim239

аппроксимация полиномами Чебышева

да, видел это, всё хорошо, но надо исключить готовые таблицы, вопрос в «горячем» расчёте прямо на месте.(у меня чисто академический интерес).

Gramozeka ★★
() автор топика
Ответ на: комментарий от Gramozeka

Ну вот вариант с таблицами поменьше через комплексную экспоненту. Там тоже таблицы, но, в принципе, их можно и на лету рассчитывать (там написано как).

SkyMaverick ★★★★★
()
Ответ на: комментарий от Gramozeka

на какой угол нужно повернуть червячный вал, при получении целеуказания(цель появляется в секторе ±60°). Считать что один полный оборот поворачивает платформу на 1°

Не пойму в чём подвох. ОБОРОТЫ_ВАЛА = СТАРЫЙ_АЗИМУТ - НОВЫЙ_АЗИМУТ?

anonymous
()
Ответ на: комментарий от anonymous

против ручной оптимизации

Результат-то бенчили? В ручной оптимизации у инструкций 2-4 последовательная зависимость по данным - выполнение следующей инструкции не может начаться до окончания предыдущей.

red75prim ★★★
()
Ответ на: комментарий от anonymous

Не пойму в чём подвох.

товарищ предложил:

таблицы Брадиса предлагали?

ну и на расстоянии в 1км 1° на горизонте на самом валу будет ...(ушёл за счётами)

зы. цель движется, линейная скорость известна..

Gramozeka ★★
() автор топика
Последнее исправление: Gramozeka (всего исправлений: 2)
Ответ на: комментарий от mazdai

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

Harald ★★★★★
()

Ряд тейлора работает, у тебя просто константы говно. Используй нормальные константы:

program sinus;

{$mode objfpc}{$H+}
uses
    SysUtils
  , LazUtf8
  ;

var
  x, y, s, h, eps, t: double;
  n: int32;
begin
  x := 0;
  h := pi/10;
  eps := 1e-7;
  repeat
    s:= x;
    t := x;
    n := 2;
    repeat
      t := -t *sqr(x)/(n*(n+1));
      s := s + t;
      inc(n,2);
    until abs(t) < eps;
    y := sin(x);
    write('x = ', x:0:2);
    write(UTF8ToConsole(' Приближенное значение синуса = '), s:0:7);
    writeln(UTF8ToConsole(' Точное значение синуса = '), y:0:7);
    x := x + h;
  until x > pi;
  writeln(UTF8ToConsole('Для выхода нажмите ENTER'));
  readln;
end.

https://imgur.com/a/kxtPnF9

fsb4000 ★★★★★
()
Ответ на: комментарий от fsb4000

сто лет паскаля не видал, уже и забыл как он выглядит. Да, про константы ещё в самом начале сказали.

Gramozeka ★★
() автор топика
Ответ на: комментарий от red75prim

Результат-то бенчили? В ручной оптимизации у инструкций 2-4 последовательная зависимость по данным

Не, никто не бенчил. А действительно, может быть гцц и нормально оптимизирует.

anonymous
()
Ответ на: комментарий от Gramozeka

зы. цель движется, линейная скорость известна..

Ясно. А пи брать равным четырём или лучше пяти?

anonymous
()
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.