atan(y/x) на двух умножениях!
Чего-то никак меня не отпустит эта тема, всё кажется, что есть очень простой и эффективный метод, надо только его найти!
Сейчас вот такое приближение "придумал". Напомню: у меня на входе два числа, x и y, причём x2+y2=1.
Допустим, что мы перешли к квадрантам I и IV, т.е x≥0, а потом ещё обеспечили условие x≥|y| (поменяв местами аргументы, если необходимо), т.е угол будет в диапазоне от -45 до +45 градусов. Тогда:
}{\approx}y(1{,}38396806-0{,}387864815x))
Максимальная ошибка составляет 0,041°, или 2,5 угловые минуты, и, возможно, ещё можно коэффициенты чуть оптимизировать.
Многовато вообще для моей задачи: мне ж этот угол ещё на 2 умножать (т.к в кватернионах угол был половинный), это будет 0,082°, при допустимой ошибке в 0,1°. Но может, кому пригодится такое приближение.
Под катом упоротые выкладки, которые привели к этому выражению... Ну как упоротые - вспоминали школьную математику, ну максимум 1-й курс института.
По сути, у нас есть синус и косинус искомого угла, скажем, z, и надо найти угол. Причём sin(z) "загибается вниз" относительно прямой линии, а tg(z) - наоборот, устремляется вверх, то, может быть, нам взять их сумму с некоторыми весами, и будет нам счастье??
То есть, ввести приближение
Правда, чтобы получить тангенс, нам нужно деление, а без него хотелось бы обойтись. Зная, что косинус не шибко далеко убежит от единицы, можно взять приближение:
}\approx{1+(1-\cos{z})}=2-\cos{z})
Вот и приходим к выражению
)
Причём здесь есть "произвол" сразу в обоих коэффициентах! Как их ни выбирай - выражение будет получаться нечётным, как и положено.
Для того, чтобы подобрать коэффициенты, преобразуем выражение следующим образом:
(раскрыли скобки, воспользовались синусом двойного угла)
Чтобы понять, "есть ли здесь какие-то перспективы", можно выбрать коэффициенты из разложения в ряд Тэйлора вблизи единицы:
+\frac{\beta}{2}(2z-\frac{8z^3}{6}))
Хотим, чтобы коэффициент при z в правой части был единичным, а при z3: нулевой. Выходит:
Откуда α=4/3, β=-1/3. Так получается приближение с "красивыми коэффициентами":
)
И как водится для ряда Тэйлора, это выражение шикарно работает вблизи нуля, но к краям ошибка растёт довольно резко. На 45° (когда x=y≈0,707) ошибка составляет 0,53°. И это не так уж плохо: примитивное приближение atan(y/x)≈y давало ошибку 4,49°, а более хитрое линейное приближение atan(y/x)≈1,0811y: 1,2°. В этом что-то есть :)
Метод наименьших квадратов (МНК)
В прошлый раз я настаивал на минимаксе, но сейчас решил применить старый добрый МНК, тупо потому что он проще :) По крайней мере привычнее. Минимакс от одного параметра - ещё куда не шло, а от двух параметров - чего-то не вполне понимаю, как подступиться.
А тут вводим интеграл от квадрата ошибки на всём интересующем нас диапазоне углов [0;A]:
^2dx)
Надо подобрать такие α, β, чтобы его минимизировать. Но проще всего сначала этот интеграл взять :)
Раскрываем скобки:
dx)
И проинтегрируем каждое слагаемое. Первое - мега-халява:
и в то же время абсолютно бесполезное слагаемое, т.к не зависит от α и β, можно было его сразу вычеркнуть.
Следующее заставляет вспомнить интегрирование по частям. Я его так и не запомнил, выписываю сначала производную произведения:
%27=u%27v+uv%27)
потом соображаю, что здесь нам подойдёт выражение x·cos(x):
%27=\cos{x}-x\sin{x})
Переношу в левую часть нужный нам кусочек -x·sin(x), а в правую всё остальное:
%27-\cos{x})
И только тогда вставляю знаки интеграла:
И наконец, добавляю недостающие 2α и границы интегрирования:
dx}=2\alpha{x}\cos{x}\bigg|_0^A-\int_0^A{\cos{x}dx}=2\alpha{A}\cos{A}-2\alpha\sin{A})
Да, мне очень стыдно :)
Следующее слагаемое очень похоже на это, но с двойкой в аргументе косинуса. Не буду торопиться, а то потом замучаюсь ошибку искать. Сначала производную выписываем:
%27=\cos{2x}-2x\sin{2x})
Выражаем -x·sin(2x) и сразу всовываем неопределённый интеграл:
dx}=\frac{1}{2}x\cos{2x}-\int{\frac{1}{2}\cos{2x}dx})
Добавляем множитель β и границы интегрирования:
dx}=\beta(\frac{1}{2}x\cos{2x}\bigg|_0^A-\int_0^A{\frac{1}{2}\cos{2x}dx})=)
)
В следующем слагаемом нужно использовать формулу произведения синусов, но её я тоже никогда не мог запомнить. Поэтому выписываю косинус от суммы, вот её почему-то помню:
}=\cos{\alpha}\cos{\beta}-\sin{\alpha}\sin{\beta})
(может потому, что если подсунуть α=β, отсюда придём к косинусу двойного угла, который "косинус квадрат минус синус квадрат". Был бы справа знак "+" - вышел бы полный бред, единичка)
Ну и косинус разности, который отсюда следует, т.к косинус чётная функция, а синус нечётная:
}=\cos{\alpha}\cos{\beta}+\sin{\alpha}\sin{\beta})
(а тут если α=β единица в правой части будет вполне логична :))
Если из верхней формулы вычесть нижнюю и поделить на два, практически устно получается то что нам надо:
}-\cos{(\alpha+\beta)}}{2})
И подставим наши x и 2x:
Ну а интегрируется оно в пол-пинка:
=)
)
Далее идёт слагаемое с sin2x, и тут "рука тянется" сразу помножить на 1/2 и на длину отрезка (как это в рядах Фурье и при вычислении мощности в цепи переменного тока, и много где ещё), но здесь так нельзя, т.к у нас ни полного периода, ни половины периода не набирается! Надо всё честно:
И ещё не забыть помножить на α2.
И последнее слагаемое очень похожее:
И ещё помножить на β2/4.
Всё, собираем вместе:
+\alpha^2(\frac{A}{2}-\frac{1}{4}\sin{2A})+\beta(\frac{1}{2}A\cos{2A}-\frac{1}{4}\sin{2A})+\beta^2(\frac{A}{8}-\frac{\sin{4A}}{32})+\alpha\beta(\frac{1}{2}\sin{A}-\frac{1}{6}\sin{3A}))
Вот теперь пора брать частные производные по α и β и приравнивать их к нулю:
+\beta(\frac{1}{2}\sin{A}-\frac{1}{6}\sin{3A})+2A\cos{A}-2\sin{A}=0,)
+\beta(\frac{A}{4}-\frac{\sin{4A}}{16})+\frac{1}{2}A\cos{2A}-\frac{1}{4}\sin{2A}=0)
Это система линейных уравнений относительно α и β. Выписывать в явном виде выражения для α и β я не решился - завёл эти 6 коэффициентов в эксель (4 элемента матрицы и 2 элемента вектора), потом решение по правилу Крамера.
В итоге, подставив A=π/4, я получил α=1,379206316, β=-0,382389721267396. Они относительно близки к 4/3 и -1/3, полученные из ряда Тэйлора. Но другие, причём их сумма равна 0,997, вот не единица :)
С такими значениями максимальная ошибка получается на угле π/4, в 0,077° График ошибки выглядит вот так:
(по оси X угол в градусах, косинус и синус которого у нас в качестве входных данных. По оси Y ошибка данного метода, также в градусах)
Замечаем характерную особенность метода наименьших квадратов: ему важнее снизить ошибку на длинном участке графика (там ошибка не более 0,032°), а высокий, но узкий пик в конце диапазона не вносит большого вклада в нашу "интегральную" ошибку.
Тэйлор нервно курит в сторонке: уменьшили ошибку в 16 раз! Но если нас интересует минимакс (минимизация максимальной ошибки), то можно ещё поиграться...
Я пока не стал "всё начинать заново" и теперь подходить к этой проблеме со стороны минимакса. Вместо этого просто стал увеличивать значение A, то есть "как бы расширять диапазон". В итоге, до значения A=0,825 максимальная ошибка на диапазоне от 0 до π/4 УМЕНЬШАЛАСЬ, а потом снова начала расти. На A=0,825 мы получили ровно те значения, что приведены в начале поста, а график ошибки выглядит так:
Да, здесь всё как мы любим: ошибка "в плюс" сравнялась с ошибкой "в минус", но главное, уменьшилась по величине почти что вдвое!
Но я не уверен, что "выжали" из этой формулы всё, что смогли, ведь от значений, полученных по "методу наименьших квадратов", мы вели поиск "по одной оси". Возможно, если поставить целью минимакс конкретно на [0;π/4], получится ещё лучше, хотя сомневаюсь, что существенно.
Кстати, я ещё опробовал МНК на простейшем линейном приближении. Скажем, на [0;π/4] мы предлагали формулу atan(y/x)≈1,0811y и получали максимальную ошибку в 1,2 градуса. А МНК предлагает 1,0634y, что даёт максимальную ошибку в 1,92 градуса, но, понятное дело, меньшую среднеквадратичную.
И примерно такое же соотношение для других диапазонов.
Наверное, эта тема меня так увлекла, потому что всегда интересно было узнать, а как же компьютер вообще считает сложные математические функции? Как в целом работает, складывает и умножает - уже разобрался, когда АЛУ ковырял, а вот эту математику обходили стороной как могли. В math.h не заглянешь, там одни заголовки, а функции уже откомпилированы под конкретную архитектуру. И в процессор не заглянешь (мой-исключение :), как там FPU вкалывает, а там же явно микрокода выше крыши! Кнут дальше арифметики не пошёл, во многих учебниках по программированию не мудрствуя лукаво говорят: ряд Тэйлора. А тут вроде как и "по делу" пришлось :)
Сейчас вот такое приближение "придумал". Напомню: у меня на входе два числа, x и y, причём x2+y2=1.
Допустим, что мы перешли к квадрантам I и IV, т.е x≥0, а потом ещё обеспечили условие x≥|y| (поменяв местами аргументы, если необходимо), т.е угол будет в диапазоне от -45 до +45 градусов. Тогда:
Максимальная ошибка составляет 0,041°, или 2,5 угловые минуты, и, возможно, ещё можно коэффициенты чуть оптимизировать.
Многовато вообще для моей задачи: мне ж этот угол ещё на 2 умножать (т.к в кватернионах угол был половинный), это будет 0,082°, при допустимой ошибке в 0,1°. Но может, кому пригодится такое приближение.
Под катом упоротые выкладки, которые привели к этому выражению... Ну как упоротые - вспоминали школьную математику, ну максимум 1-й курс института.
По сути, у нас есть синус и косинус искомого угла, скажем, z, и надо найти угол. Причём sin(z) "загибается вниз" относительно прямой линии, а tg(z) - наоборот, устремляется вверх, то, может быть, нам взять их сумму с некоторыми весами, и будет нам счастье??
То есть, ввести приближение
Правда, чтобы получить тангенс, нам нужно деление, а без него хотелось бы обойтись. Зная, что косинус не шибко далеко убежит от единицы, можно взять приближение:
Вот и приходим к выражению
Причём здесь есть "произвол" сразу в обоих коэффициентах! Как их ни выбирай - выражение будет получаться нечётным, как и положено.
Для того, чтобы подобрать коэффициенты, преобразуем выражение следующим образом:
(раскрыли скобки, воспользовались синусом двойного угла)
Чтобы понять, "есть ли здесь какие-то перспективы", можно выбрать коэффициенты из разложения в ряд Тэйлора вблизи единицы:
Хотим, чтобы коэффициент при z в правой части был единичным, а при z3: нулевой. Выходит:
Откуда α=4/3, β=-1/3. Так получается приближение с "красивыми коэффициентами":
И как водится для ряда Тэйлора, это выражение шикарно работает вблизи нуля, но к краям ошибка растёт довольно резко. На 45° (когда x=y≈0,707) ошибка составляет 0,53°. И это не так уж плохо: примитивное приближение atan(y/x)≈y давало ошибку 4,49°, а более хитрое линейное приближение atan(y/x)≈1,0811y: 1,2°. В этом что-то есть :)
Метод наименьших квадратов (МНК)
В прошлый раз я настаивал на минимаксе, но сейчас решил применить старый добрый МНК, тупо потому что он проще :) По крайней мере привычнее. Минимакс от одного параметра - ещё куда не шло, а от двух параметров - чего-то не вполне понимаю, как подступиться.
А тут вводим интеграл от квадрата ошибки на всём интересующем нас диапазоне углов [0;A]:
Надо подобрать такие α, β, чтобы его минимизировать. Но проще всего сначала этот интеграл взять :)
Раскрываем скобки:
И проинтегрируем каждое слагаемое. Первое - мега-халява:
и в то же время абсолютно бесполезное слагаемое, т.к не зависит от α и β, можно было его сразу вычеркнуть.
Следующее заставляет вспомнить интегрирование по частям. Я его так и не запомнил, выписываю сначала производную произведения:
потом соображаю, что здесь нам подойдёт выражение x·cos(x):
Переношу в левую часть нужный нам кусочек -x·sin(x), а в правую всё остальное:
И только тогда вставляю знаки интеграла:
И наконец, добавляю недостающие 2α и границы интегрирования:
Да, мне очень стыдно :)
Следующее слагаемое очень похоже на это, но с двойкой в аргументе косинуса. Не буду торопиться, а то потом замучаюсь ошибку искать. Сначала производную выписываем:
Выражаем -x·sin(2x) и сразу всовываем неопределённый интеграл:
Добавляем множитель β и границы интегрирования:
В следующем слагаемом нужно использовать формулу произведения синусов, но её я тоже никогда не мог запомнить. Поэтому выписываю косинус от суммы, вот её почему-то помню:
(может потому, что если подсунуть α=β, отсюда придём к косинусу двойного угла, который "косинус квадрат минус синус квадрат". Был бы справа знак "+" - вышел бы полный бред, единичка)
Ну и косинус разности, который отсюда следует, т.к косинус чётная функция, а синус нечётная:
(а тут если α=β единица в правой части будет вполне логична :))
Если из верхней формулы вычесть нижнюю и поделить на два, практически устно получается то что нам надо:
И подставим наши x и 2x:
Ну а интегрируется оно в пол-пинка:
Далее идёт слагаемое с sin2x, и тут "рука тянется" сразу помножить на 1/2 и на длину отрезка (как это в рядах Фурье и при вычислении мощности в цепи переменного тока, и много где ещё), но здесь так нельзя, т.к у нас ни полного периода, ни половины периода не набирается! Надо всё честно:
И ещё не забыть помножить на α2.
И последнее слагаемое очень похожее:
И ещё помножить на β2/4.
Всё, собираем вместе:
Вот теперь пора брать частные производные по α и β и приравнивать их к нулю:
Это система линейных уравнений относительно α и β. Выписывать в явном виде выражения для α и β я не решился - завёл эти 6 коэффициентов в эксель (4 элемента матрицы и 2 элемента вектора), потом решение по правилу Крамера.
В итоге, подставив A=π/4, я получил α=1,379206316, β=-0,382389721267396. Они относительно близки к 4/3 и -1/3, полученные из ряда Тэйлора. Но другие, причём их сумма равна 0,997, вот не единица :)
С такими значениями максимальная ошибка получается на угле π/4, в 0,077° График ошибки выглядит вот так:
(по оси X угол в градусах, косинус и синус которого у нас в качестве входных данных. По оси Y ошибка данного метода, также в градусах)
Замечаем характерную особенность метода наименьших квадратов: ему важнее снизить ошибку на длинном участке графика (там ошибка не более 0,032°), а высокий, но узкий пик в конце диапазона не вносит большого вклада в нашу "интегральную" ошибку.
Тэйлор нервно курит в сторонке: уменьшили ошибку в 16 раз! Но если нас интересует минимакс (минимизация максимальной ошибки), то можно ещё поиграться...
Я пока не стал "всё начинать заново" и теперь подходить к этой проблеме со стороны минимакса. Вместо этого просто стал увеличивать значение A, то есть "как бы расширять диапазон". В итоге, до значения A=0,825 максимальная ошибка на диапазоне от 0 до π/4 УМЕНЬШАЛАСЬ, а потом снова начала расти. На A=0,825 мы получили ровно те значения, что приведены в начале поста, а график ошибки выглядит так:
Да, здесь всё как мы любим: ошибка "в плюс" сравнялась с ошибкой "в минус", но главное, уменьшилась по величине почти что вдвое!
Но я не уверен, что "выжали" из этой формулы всё, что смогли, ведь от значений, полученных по "методу наименьших квадратов", мы вели поиск "по одной оси". Возможно, если поставить целью минимакс конкретно на [0;π/4], получится ещё лучше, хотя сомневаюсь, что существенно.
Кстати, я ещё опробовал МНК на простейшем линейном приближении. Скажем, на [0;π/4] мы предлагали формулу atan(y/x)≈1,0811y и получали максимальную ошибку в 1,2 градуса. А МНК предлагает 1,0634y, что даёт максимальную ошибку в 1,92 градуса, но, понятное дело, меньшую среднеквадратичную.
И примерно такое же соотношение для других диапазонов.
Наверное, эта тема меня так увлекла, потому что всегда интересно было узнать, а как же компьютер вообще считает сложные математические функции? Как в целом работает, складывает и умножает - уже разобрался, когда АЛУ ковырял, а вот эту математику обходили стороной как могли. В math.h не заглянешь, там одни заголовки, а функции уже откомпилированы под конкретную архитектуру. И в процессор не заглянешь (мой-исключение :), как там FPU вкалывает, а там же явно микрокода выше крыши! Кнут дальше арифметики не пошёл, во многих учебниках по программированию не мудрствуя лукаво говорят: ряд Тэйлора. А тут вроде как и "по делу" пришлось :)