Из журнала Info Guide #8, Рязань, 12.2005 Синус и компания Смирнов В.В. Вычисление тригонометрических и алгебраических функций (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)