Смирнов В.В.
Вычисление тригонометрических
и алгебраических функций
(SIN, COS, TAN, ATAN, LOG, EXP, ...)
в языках высокого уровня
Все эти функции в БЕЙСИКЕ именуют "ста─
ндартными" или "встроенными". В ФОРТРАНЕ,
ПАСКАЛЕ и в др. языках высокого уровня -
"стандартными" или "библиотечными".
Для разработки алгоритмов вычисления
этих функций всюду,в т.ч. в микрокалькуля─
торах,применена тема т.н. "рядов".Ряд есть
сумма бесконечного числа слагаемых. Напри─
мер, для функции sin (x) ряд является зна─
копеременным, по нечётным степеням аргуме─
нта (аргумент в радианах):
sin(x)=x-x^3/3!+x^5/5!-x^7/7!+x^9/9!- ...
... +(-1)^n*x^(2*n+1)/(2*n+1)!+ ...
Для некоторых других функций ряды пред─
ставляются такими:
cos(x)=1-x^2/2!+x^4/4!-x^6/6!+ ...
... +(-1)^n*x^(2*n)/(2*n)!+ ...
т.е.знакопеременный ряд по чётным степеням
аргумента;
e^x = 1+x/1!+x^2/2!+x^3/3!+x^4/4!+ ...
... +x^n/n!+ ...
обозначив (x-1)/(x+1)=z,
запишем следующий ряд:
ln(x)=2(z+z^3/3+z^5/5+z^7/7+z^9/9+ ...
... +z^(2*n+1)/(2*n+1)+ ...
arcsin(x)=x+x^3/(2*3)+1*3*x^5/(2*4*5)+
+1*3*5*x^7/(2*4*6*7)+ ...
... +(1*3*5*...*(2*n-1)*x^(2*n+1)/
(2*4*6*...*(2*n)*(2*n+1))+...
(Ряды некоторых функций,например,танге─
нса и секанса, содержат в качестве коэффи─
циентов специальные числа (здесь - числа
Бернулли и числа Эйлера соответственно).
Для вычисления этих чисел потребуется со─
ответствующий алгоритм,который станет,сле─
довательно,частью алгоритма вычисления ря─
да соответствующей функции.)
Эти представления рядов различных функ─
ций приводятся в математических справочни─
ках, в математических словарях,в математи─
ческой энциклопедии.
Исходно ряды получают и исследуют в
вузовских курсах математического анализа
(иначе - курс дифференциального и интегра─
льного исчисления), в состав которых (кур─
сов) тематика "рядов" входит в обязатель─
ном порядке - можете поискать в вузовских
многочисленных томах "матана".В "рядах" из
справочников следует иметь в виду для да─
льнейшего, что приводимое в них выражение
т.н."общего члена ряда" (позволяющее полу─
чить по номеру слагаемого выражение этого
слагаемого) может начинаться как с нуля,
так и с единицы - по выбору составителей
справочника. Этот начальный номер следует
уточнить,в конкретном случае подставив n=1
и оценив, получается ли в полном представ─
лении первый член ряда, или же второй. (В
последнем случае, следовательно,номера на─
чинаются от нуля.)
При вычислении этих "стандартных" функ─
ций широко используется функция "фактори─
ал",обозначаемая "восклицательным знаком".
Напомним здесь:целочисленная функция цело─
численного аргумента, "эн факториал" (n!)
означает произведение всех целых чисел от
1 до n.
(Ред.: 0! принято равным 1. Для дробно─
го аргумента факториал можно доопределить
через гамма-функцию.)
Итак, о вычислении тригонометрических и
алгебраических функций.
Для практической реализации алгоритма,
очевидно, следует вычислить и просуммиро─
вать некоторое КОНЕЧНОЕ КОЛИЧЕСТВО СЛАГАЕ─
МЫХ - ЧЛЕНОВ РЯДА.Для построения алгоритма
вычисления, скажем,синуса следует задаться
точностью вычислений - допустим, 0.00001.
Это значение будет использовано в единст─
венном в программе операторе условного пе─
рехода:когда очередное слагаемое по модулю
станет меньше заданной точности (перемен─
ная Т ),следует прекратить накопление сум─
мы и перейти к печати результата.
В алгоритме сначала назначаем перемен─
ные, используемые в программе: X - аргу─
мент; N - текущее значение показателя сте─
пени (и аргумента факториала); S - пере─
менная, в которой накапливается сумма и, в
конечном итоге,в которой оказывается вычи─
сленное значение синуса (соответственно, С
- переменная суммирования косинуса); H -
очередное слагаемое для суммирования.
Алгоритм может быть построен с вычисле─
нием очередного члена ряда (слагаемого):
как по полной формуле "общего члена ряда",
так и с использованием для вычислений оче─
редного слагаемого, результата вычислений
предыдущего слагаемого (т.е. с переопреде─
лением переменной, где новое её значение
определяется через старое значение). Этот
второй путь лучше,проще,короче,т.к.не при─
ходится отдельно вычислять функцию "факто─
риал". Кроме того, второй способ, видимо,
короче по времени. А самое главное, прямой
метод может отказать из-за переполнения,
т.к. факториал растёт даже быстрее,чем по─
казательная функция.
Вот полный алгоритм вычисления функции
"синус":
(А)
10 LET X=1.2345
20 LET T=0.00001
30 LET N=1
40 LET S=X
50 LET H=X
60 LET H=H*(-1)*X^2/(N+1)/(N+2)
70 LET S=S+H
80 LET N=N+2
90 IF ABS(H)>T THEN GOTO 60
100 PRINT "S=";S
В строках 10 - 50 задаём начальные зна─
чения соответствующих переменных. (Вообще
накопление суммы чаще начинают с нуля, но
здесь удобнее с первого члена ряда,с исхо─
дного X ).
В строке 60 происходит вычисление оче─
редного члена ряда с использованием его
"старого" значения. "Старое" - это выраже─
ние предыдущего члена ряда. Формула строки
60 получена за счет сравнения ДВУХ ЛЮБЫХ
СОСЕДНИХ членов ряда. Скажем, 2-го и 3-го.
Чем они отличаются? На "икс в квадрате" -
в числителе, на два новых следующих поряд─
ковых числа - в знаменателе,а также знаком
"плюс-минус". Всё это и отражено в строке
60.
В строке 70 происходит накопление сум─
мы, т.е. "главное вычисление". В строке 80
вычисляется новое значение параметра N для
последующих членов ряда,если они будут вы─
числяться.
В строке 90 как раз определяется: "пре─
кратить вычисления", или "продолжать",если
точность ещё не достигнута.
Программа вычислений по первому вариан─
ту (с использованием выражения общего чле─
на ряда) дана ниже:
(В)
10 LET X=1.23456
20 LET N=1
30 LET T=0.00001
40 LET S=0
50 LET A=1
60 LET F=1
70 FOR K=1 TO (2*N+1):LET F=F*K:NEXT K
80 LET H=A*X^(2*N+1)/F
90 LET S=S+H
100 LET A=-A
110 LET N=N+1
120 IF ABS(H)>T THEN GOTO 60
130 PRINT "S=";S
Здесь в строках 60, 70 вычисляется фун─
кция "факториал": в цикле накапливается
произведение чисел от единицы до текущего
параметра.Использование переменной A - это
"ухищрения" в связи с тем, что нужное по
формуле общего члена ряда вычисление "зна─
коизменяющего" сомножителя ПРИ ПЕРВОМ ПРО─
ХОДЕ, для "минус единицы в степени нуль",
может происходить НЕВЕРНО (или с сообщени─
ем об ошибке) на разных машинах и в разных
версиях Бейсика.Нужно отметить,что показа─
тельная функция с отрицательным основанием
не определена. Видим, что эта программа
более громоздка,чем предыдущая. К тому же,
имеет недостатки, названные выше. Поэтому
с этим типом программ (В) далее не работа─
ем.
* * *
Если бы в программе (А) потребовалась
выдача значений функций для аргумента,вво─
димого с клавиатуры ("справочник Бради─
са"), то в строке 10 вместо назначения ар─
гумента будет оператор INPUT.
Или, скажем,вместо вывода на печать ре─
зультат может быть передан другой програм─
ме, а эта станет, т.о.,подпрограммой - для
тех версий Бейсика или другого языка, где
исходно нет тригонометрии, например, в
"Бейсик-микро" (ж. "Радио", 1980-е гг.),
либо в каком-нибудь /разработанном вами!/
простейшем языке высокого уровня. В "само─
пальном" ЯВУ этот путь позволяет получить
тригонометрию более простыми средствами
самого же языка, а не более трудоёмкой ма─
шиннокодовой процедурой. Это экономит па─
мять (ПЗУ), но,например,ведёт к потерям во
времени.
Для вычисления значений функции "коси─
нус" нетрудно составить аналогичную прог─
рамму. Если же нужно иметь и синус,и коси─
нус, можно доработать исходную программу,
добавив в (А) следующие строки:
45 LET C=1
75 LET C=C+H/X*(N+2)
120 PRINT "C=";C
Здесь в строке 45 назначается перемен─
ная С, в которой будет накапливаться зна─
чение ряда косинуса. Этой переменной прис─
воено, соответственно, начальное значение
ряда косинуса.
Единственная строка вычислений здесь
(75) использует вычисленное к этому момен─
ту значение очередного (с тем же порядко─
вым номером) члена ряда синуса. В этой
строке происходит накопление ряда косину─
са, с откорректированным значением H. (Пе─
ременная H этой строкой не меняется.) Кор─
рекция заключается в домножении на множи─
тель (N+2)/X. Вид множителя получен из
сравнения "одноимённых" членов ряда синуса
и косинуса.
Обычно при практической реализации вы─
числений, используя периодичность тригоно─
метрических функций, пересчитывают аргу─
мент,вычитая целое число периодов (кратное
2*пи=6.283185307179586476925286766559...).
Иначе ряд "хуже сходится", т.е. требуется
гораздо большее число итераций (см. ниже).
Кроме того,в таких же целях аргумент после
этого ещё уменьшают: используя "школьные"
формулы приведения,приводят в первый квад─
рант. Тогда вычисления очень коротки.
* * *
Далее,развивая и изменяя программу (А),
проведём небольшое исследование. (Дальней─
шее имеет общих характер для численных ме─
тодов.Т.е.приложимо не только к вычислению
стандартных функций ЯВУ.)
Любая задача, после её решения, требует
ПРОВЕРКИ.
Здесь можно выполнить проверку двумя
способами:
а) можно вычислять известные функции от
известных аргументов, например, синус от
пи/6;
б) вычисляя "встроенными" функциями SIN,
COS, ... которые, конечно, имеются во всех
современных версиях ЯВУ.
Для проверки по второму способу надо
добавить строку:
110 PRINT "SIN=";SIN(X)
При этом полезно подкорректировать фор─
мат вывода предыдущей строки (100), чтобы
удобно было сравнивать одноимённые разряды
(там - добавить 2 пробела внутри кавычек).
Далее. Добавим в (А) строку:
78 PRINT S;" ";C;" ";H
В результате того,что эта печать задана
НА ПОВТОРЯЮЩЕМСЯ УЧАСТКЕ ПРОГРАММЫ, в про─
цессе решения машина напечатает в три сто─
лбца текущие значения названых переменных.
По этим данным можно видеть, как происхо─
дило приближение к настоящему значению си─
нуса (косинуса). (Настоящее - вычисленное
встроенной соответствующей функцией.) Так─
же видно уменьшение модуля H до значения,
ниже заданного в Т. Причём меньшее Т -
единственная строка перед печатью итогово─
го результата.Из этих трёх столбцов следу─
ет и ещё кое-какая информация: СКОЛЬКО РАЗ
машина прошла на этом повторяющемся участ─
ке программы.А именно:сколько строк - сто─
лько и повторений.
Для оценки количества проходов по пов─
торяющемуся участку программы обычно на
циклическом её участке добавляют счётчик:
23 LET K=0
65 LET K=K+1
133 PRINT "K=";K
Изменяя требуемую точность - задавая
0.0001 или 0.000001 - будем получать раз─
ное число строк в столбцах.Чем больше точ─
ность, тем, очевидно,больше число проходов
(итераций), и соответственно, больше пара─
метр К и больше число строк вспомогатель─
ной печати (здесь - переменных S,C,H внут─
ри цикла).
Далее. Скажем, при точности "десять в
минус пятой" в итоговых результатах совпа─
дут (со встроенной функцией) по пять цифр
после запятой (Ред.: вернее сказать, резу─
льтаты будут отличаться не больше чем на
единицу 5-го разряда после запятой. Ниже
аналогично).
При точности "десять в минус шестой" -
совпадут по шесть цифр.
Если назначить точность "десять в минус
восьмой", то совпадут все цифры, и, таким
образом, результат по этому алгоритму ока─
зывается в точности тем же,что и по встро─
енным функциям SIN(X), COS(X).
Стало быть, своей программой вы можете
превзойти точность вашей версии языка: до─
пустим,если в нем не заложена "двойная то─
чность" (16 значащих цифр) - задайте дос─
таточно малое значение Т.
Если задать значения аргумента более
"пи/2" или ещё более - "Х+2*пи", - то мож─
но видеть,как резко замедляется сходимость
ряда.Причём сначала появляются сильные ко─
лебания модуля суммы.Из всего этого и сле─
дует целесообразность пересчёта аргумента
в первый квадрант.
SWW (A.E.-2)
