Info Guide #08 30 ноября 2005

Программирование

For Coderz - Вычисление тригонометрических и алгебраических функций в языках высокого уровня.

               Смирнов В.В.

      Вычисление тригонометрических
          и алгебраических функций
   (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)