В математике под функцией понимают отображение множества A на множество B, при котором каждому элементу из A однозначно соответствует некоторый элемент из B. A называют областью определения, B - областью значений.
Функция в программировании - это конструкция языка, представляющая собой алгоритмическую процедуру, описанную отдельно от тела основной программы. В бейсике "МК 85" нет ни функций, ни каких-либо других программных модулей. Это обстоятельство создало ряд проблем на пути унификации бейсик-программ, заставило включить в их тексты служебные обозначения. Так, определяемая пользователем функция обозначается C=f{A,...,B}, где A,...,B - параметры, C - результат.
Среди способов задания функции, с которыми приходится сталкиваться пользователю, можно выделить выражения и таблицы. Выражение дает более полную информацию о функциональной зависимости и поддается обработке алгоритмами математического анализа, использующими понятие предела. Таблично заданную функцию во многих задачах приходится заменять подходящим выражением (см. раздел "Функции в форме таблиц").
Арифметические:
-X - инверсия знака,
"минус";
SGN X
- "знак": 1 при x>0, -1 при x<0 и 0 при x=0;
ABS X - абсолютное
значение, "модуль";
INT X - целая часть;
FRAC X - дробная часть;
RND(X,N) - корректное
округление, оставляет N знаков после десятичной точки.
Тригонометрические и обратные тригонометрические:
SIN X - синус,
ASN X - арксинус (|x|≤1);
COS X - косинус,
ACS X - арккосинус
(|x|≤1);
TAN X - тангенс,
ATN X - арктангенс.
Логарифмические, показательные, степенные:
LOG X - десятичный
логарифм x>0;
LN X - натуральный
логарифм x>0;
EXP X - экспонента
x;
X-P
- x в степени p (если p нецелое, x≥0);
SQR X - квадратный
корень из x≥0.
Вычисление функций последних двух групп производится подпрограммами, использующими разложение в степенной ряд.
Обращение к математическим функциям разрешено везде, где может стоять числовая константа или переменная.
В простейшем случае функция задается одним оператором присвоения. Например, гиперболический тангенс может быть явно выражен через экспоненты: Y=(EXP X-EXP -X)/(EXP X+EXP -X). Он будет вычисляться вдвое быстрее, если экспоненты X и -X вычислить заранее: A=EXP X:B=EXP -X:Y=(A-B)/(A+B). Если вычислять гиперболический тангенс через его собственный степенной ряд, быстродействие увеличится еще вдвое.
Иногда оказывается оправданным хранение функции в виде таблицы значений или аппроксимация ее более простой функцией, однако здесь выигрыш в быстродействии оплачивается потерей точности и сужением области определения. Вычисление функций через степенные ряды описано в разделе "Степенные ряды", обработка таблиц - в разделе "Функции в форме таблиц".
Часто встречаются функции, для которых задано несколько выражений, справедливых на отдельных участках области определения. Например, экспоненциальная плотность вероятности:
j(x)= 0 при x<0;
j(x)= e-lx
при x≥0
Она может быть описана строкой F=0:IF X>0;F=EXP (-L*X). Первое присвоение выполняется всегда. В данном случае это оправдано, так как делает компактнее описание процедуры, не увеличивая ее трудоемкости. В общем случае следовало бы записать:
50 IF X≤0;F=0
51 IF X>0;F=EXP (-L*X)
Аналогично описываются особые точки. Например, для f(x)=sin(x)/x:
50 F=1:IF X¦0;F=SIN X/X
Любой вычислительный алгоритм может быть описан на бейсике набором операторов присвоения и условных арифметических операторов ( IF ... THEN ...), однако бейсик располагает и более мощными средствами управления, облегчающими процесс программирования.
Цикл FOR - NEXT удобно применять при табулировании функции. Пример:
FOR X=A TO B STEP S:Y=SIN X+COS X:PRINT "Y(";X;")=";Y:NEXT X.
Здесь (A,B) - интервал табулирования, S - его шаг.
Если функция зависит от нескольких параметров, их изменение описывается вложенными циклами:
50 FOR X=A TO B STEP S
60 FOR V=C TO D STEP R
70 PRINT X,V,SQR (X-1/V)
80 NEXT V
90 NEXT X
Здесь для каждого значения X вычисляется ряд значений функции SQR (X-1/V) в интервале (C,D) с шагом R. Таблица результатов будет двумерной.
При использовании циклов степень их вложенности не должна превышать 4. Если алгоритм требует большей вложенности, цикл может быть заменен конструкцией вида:
50 X=A
60 ...
...
90 X=X+S:IF SGN S*X¦SGN
S*B THEN 60
Обычно SGN S определен в алгоритме жестко и может быть исключен из условия (X≤B при S>0 и X≥B при S<0).
Чтобы вольное использование цикла не приводило к ошибкам, детально исследуйте его действие, придумав соответствующие тесты.
Отрабатывая оператор FOR - TO - STEP - NEXT, система выполняет операции: установить или переустановить цикл, продолжить цикл, сбросить цикл. Цикл устанавливается, когда выполнение программы доходит до оператора FOR. Если управляющая переменная цикла уже есть в списке циклов, новый цикл не устанавливается, а старый восстанавливается в исходное положение. (Именно поэтому разрешен выход из цикла минуя оператор NEXT.) В системную область памяти записывается имя управляющей переменной, выполняются все операторы присвоения, встречающиеся в заголовке цикла, запоминаются конечное значение и шаг. Обработка заголовка программы завершается запоминанием адреса следующего оператора, чтобы NEXT I знал, куда передать управление. Выполнение системой оператора NEXT может закончиться либо сбросом цикла, либо его продолжением. По параметру (I) система определяет, к какому циклу относится NEXT. После этого изменяется значение управляющей переменной (I=I+S) и проверяется, не превысило ли оно заданное конечное значение. Если превысило, цикл сбрасывается и управление передается следующему за NEXT оператору. (Физически сброс выражается в освобождении соответствующей строки списка циклов в системной области памяти.) Если условие завершения не выполнилось, управление передается первому оператору тела цикла.
Подпрограмма GOSUB - RETURN часто используется для того, чтобы вынести функцию за пределы алгоритма ее обработки. Это делает программу нагляднее. Если значение функции должно вычисляться более, чем в одной точке алгоритма, выделение ее в подпрограмму экономит память.
Определенный интерес представляет рекурсивный вызов подпрограммы. Вот, например, один из способов вычисления факториала:
10 INPUT N:M=1:GOSUB 20:PRINT
"N!=";M:END
20 IF N>1;M=M*N:N=N-1:GOSUB 20
30 RETURN
Бейсик "МК 85" выдаст ошибку при N>8, так как уровень вложения подпрограмм не может превышать 8. Как и цикл, подпрограмма хранит свои параметры в системной области памяти. Оператор GOSUB #P,L запоминает адрес следующего за ним оператора и передает управление в программную область P в начало строки L. Оператор RETURN передает управление оператору, адрес которого записан в последней строке адресов возврата из подпрограмм и "вычеркивает" эту строку. Следовательно, нельзя выходить из подпрограммы, минуя RETURN.
В разделе представлены быстрые алгоритмы вычисления специальных функций.
En(x) = т1гe-xt/tn dt.
Рекуррентная формула и разложение в степенной ряд:
En+1(x)=(e-x-xEn(x))/n
E0(x)=e-x/x
E1(x)= -g-ln x -
е(i)(-1)ixi/i!/i,
i=1...г
10 INPUT "n",A,"x",B:G=EXP -B:IF
A=0;C=G/B:GOTO 60
20 E=1:C=0:D=0
30 D=D+1:F=C:E=-E*B/D:C=C-E/D:IF C¦F THEN 30
40 C=C-.5772156649-LN B:D=1
50 IF D¦A;C=(G-B*C)/D:D=D+1:GOTO 50
60 PRINT C:GOTO 10
Размер: 150, G
Пример: n=2; x=0,1.
Ответ: E2(0,1)=0,722545036.
Ei(x)= тxгe-t/t dt, x>0
Разложение в ряд:
Ei(x)= g+ln x + е(i)xi/i!/i, i=1...г
10 INPUT A:B=0:C=.5772156649+LN
A:D=1
20 E=C:B=B+1:D=D*A/B:C=C+D/B:IF C¦E THEN 20
30 PRINT C:GOTO 10
Размер: 75, F
Пример: x=0,1.
Ответ: Ei(0,1)=-1,62281282.
an(x) = т1гtne-xt dt.
Рекурсивная формула:
an(x)=(e-x+nan-1(x))/x
a0(x)=e-x/x
10 INPUT "n",A,"x",B:C=0:E=EXP -B:D=E/B
20 IF C=A;PRINT D:GOTO 10
30 C=C+1:D=(E+C*D)/B:GOTO 20
Размер: 69, E
Пример: n=3; x=2.
Ответ: a3(2)=0,321421302.
bn(x) = т-11tne-xt dt.
Рекурсивная формула:
bn(x)=( (-1)nex-e-x+nbn-1(x)
)/x
b0(x)= (ex-e-x)/x
10 INPUT "n",A,"x",B:C=0:F=1:E=EXP
-B:D=(1/E-E)/B
20 IF C=A;PRINT D:GOTO 10
30 C=C+1:F=-F:D=(F/E-E+C*D)/B:GOTO 20
Размер: 88, F
Пример: n=2; x=0,1.
Ответ: b2(0,1)=0,6686087.
x!=G(x+1)= т0гe-tt x dt
Программа 53. Гамма-функция. Аппроксимация в интервале (0;1) многочленом восьмой степени (максимальная погрешность 3Ч10-7).
Расчетная формула:
G(x+1)= е(i)cixi; i=0...8
где ci - числовые коэффициенты (см. текст программы).
При |x|>1 действует рекуррентная формула:
G(x+1)=x G(x).
При x<0 действует формула
G(x+1)= |x|p / G|x+1| / sin(px+p)
100 C=ABS A:B=FRAC C:D=(35868343*B-193527818)*B+482199394
110 D=(((D*B-756704078)*B+918206857)*B-897056937)*B+988205891
120 D=(D*B-577191652)*B*1E-9+1
130 B=B+1:IF B≤C;D=D*B:GOTO 130
140 IF A≥0;RETURN
150 IF FRAC A¦0;MODE 5:D=p*A/D/SIN
(p*A):RETURN
160 D=1E4000:RETURN
Размер: 211, D
Пример: x=-0,5.
Обращение к подпрограмме:
10 INPUT A:GOSUB 100:PRINT D:GOTO 10
Ответ: G(0,5)=1,77245369.
Программа 54. Гамма-функция. Аппроксимация формулой Стирлинга (дает 5 верных знаков в интервале (-20;50). Расчетная формула:
G(x+1) ¦ Ц{(2p)/(x+21)} exp{(x+21)(ln(x+21)-1) +1/12(x+21)}/Х(i)(x+i), i=0...20
10 INPUT A:B=1:FOR C=A+1 TO A+20:B=B*C:NEXT
C
20 A=EXP (C*LN C-C+1/12/C)*SQR (2*p/C)/B:PRINT
A:GOTO 10
Размер: 69, C
Пример: x=-0,5.
Ответ: G(0,5)=1,77245371.
Программа 55. Дигамма-функция.
Y(x) = d ln G(x)/dx
Усеченный ряд
Y(x)=ln x+1/2/x + 1/12/x2+1/120/x4
обеспечивает хорошую точность только при x>20. При x<20 используется рекуррентная формула
Y(x-1) = Y(x) - 1/x
10 INPUT A:B=20+FRAC A:IF B<A;B=A
20 C=LN B+1/2/B-1/12/B/B+1/120/B-4
30 IF A>B;PRINT C:GOTO 10
40 C=C-1/B:B=B-1:GOTO 30
Размер: 88, C
Пример: x=1.
Ответ: Y(1)=-0,57721562.
Программа 56. Тригамма-функция.
Y'(x)=dY(x) / dx
Усеченный ряд для x>20:
Y'(x)=1/x - 1/2/x2+1/6/x3 - 1/30/x5
Рекуррентная формула для x<20:
Y'(x+1)= Y'(x)+1/x2.
10 INPUT A:B=20-FRAC A:IF B<A;B=A
20 C=1/B-1/2/B/B+1/6/B-3-1/30/B-5
30 IF A>B;PRINT C:GOTO 10
40 C=C+1/B/B:B=B-1:GOTO 30
Размер: 91, C
Пример: x=0.
Ответ: Y'(0)=1,64493406.
Программа 57. Неполная Гамма-функция. Разложение в ряд:
g(a, x)=e-xxa/a е(i)xia!/(a+i)!, i=0...г
10 INPUT "a",A,"x",B:C=1:D=A:E=1
20 F=C:D=D+1:E=E*B/D:C=C+E:IF F¦C THEN 20
30 PRINT C/EXP B*B-A/A:GOTO 10
Размер: 79, F
Пример: a =0,1; x=10.
Ответ: g(0,1; 10)=9,51350196.
Программа 58. Функции Бесселя дробного порядка. Расчетные формулы:
Jv(x) = (x/2)v е(i)(-x2/4)v/(v+i)!,
i=0...г
Yv(x) = (Jv(x)cos(vp)
- J-v(x))/sin(vp)
В окрестности целых значений n погрешность вычисления Jv(x) не ограничена. Если порядок отличается от целого менее, чем на 10-4, пользуйтесь программой 59.
10 INPUT "n",A,"x",E:E=E/2:GOSUB
100:GOSUB 40
20 PRINT "J";F:G=F:A=-A:GOSUB 100:GOSUB 40:B=-A*p
30 PRINT "Y";(G*COS B-F)/SIN B:GOTO 10
40 D=E-A/D:C=0:F=D
50 B=F:C=C+1:D=-D*E*E/C/(A+C):F=F+D:IF B¦F THEN 50
60 RETURN
100 Программа 53
Размер: 363, G
Пример: n=2,1; x=3.
Ответ: J2.1(3)=0,4761627444; Y2.1(3)=-0,205851878.
Программа 59. Функции Бесселя натурального порядка.
10 INPUT "n",A,"x",E:E=E/2
20 C=E-A:G=0:D=2*LN E+1.1544313298:IF A=0
THEN 50
30 F=1/A:B=0
40 B=B+1:D=D-1/B:C=C/B:G=G-F:IF B<A;F=F*E*E/B/(A-B):GOTO 40
50 B=0:G=G/C:F=0
60 F=F+C:G=G+D*C:IF F+C=F;PRINT "J";F,"Y";G/p:GOTO
10
70 B=B+1:C=-C*E*E/B/(B+A):D=D-1/B-1/(B+A):GOTO 60
Размер: 229, G
Пример: n=2; x=3.
Ответ: J2(3)=0,4860912606; Y2(3)=-0,16040039.
Программа 60. Модифицированные функции Бесселя дробного порядка. Расчетные формулы:
В окрестности целых значений n погрешность вычисления Kv(x) не ограничена. Если порядок отличается от целого менее, чем на 10-4, пользуйтесь программой 61.
Iv(x) = (x/2)v е(i)(x2/4)v/(v+i)!/i!,
i=0...г
Kv(x) = p/2 (I-v(x)
- Iv(x))/sin(vp)
10 INPUT "n",A,"x",E:E=E/2:GOSUB
100:GOSUB 40
20 PRINT "I";F:G=F:A=-A:GOSUB 100:GOSUB 40
30 PRINT "K";p/2/SIN
(A*p)*(G-F):GOTO
10
40 D=E-A/D:C=0:F=D
50 B=F:C=C+1:D=D*E*E/C/(A+C):F=F+D:IF B¦F THEN 50
60 RETURN
100 Программа 53
Размер: 360, G
Пример: n=2,1; x=3.
Ответ: I2.1(3)=2,086691962; K2.1(3)=0,06513768076.
Программа 61. Модифицированные функции Бесселя натурального порядка.
10 INPUT "n",A,"x",E:E=E/2
20 C=E-A:G=0:D=2*LN E+1.1544313298:IF A=0
THEN 50
30 F=1/A:B=0
40 B=B+1:D=D-1/B:C=C/B:G=G+F:IF B<A;F=-F*E*E/B/(A-B):GOTO 40
50 B=0:G=G/C:F=0
60 F=F+C:G=G-D*C*-1-A:IF F+C=F;PRINT "I";F,"K";G/2:GOTO
10
70 B=B+1:C=C*E*E/B/(B+A):D=D-1/B-1/(B+A):GOTO 60
Размер: 234, G
Пример: n=2; x=3.
Ответ: I2(3)=2,245212441; K2(3)=0,0615104577.
Программа 62. Сферические функции Бесселя.
jn(x)= Ц(p/2/x)
Jn+0.5(x);
yn(x)= Ц(p/2/x)
Jn+0.5(x);
10 INPUT "n",A,"x",E:E=E/2:A=A+.5:MODE
5:GOSUB 70:GOSUB 40
20 G=F:A=-A:GOSUB 70:GOSUB 40:C=SQR (p/E)/2
30 PRINT "j";G*C,"y";-1-(.5-A)*F*C:GOTO
10
40 D=E-A/D:C=0:F=D
50 B=F:C=C+1:D=-D*E*E/C/(A+C):F=F+D:IF B¦F THEN 50
60 RETURN
70 D=SQR p:C=ABS
A:FOR B=.5 TO C:D=D*B:NEXT B
75 IF A<0;D=-1-(A+.5)*p*C/D
80 RETURN
Размер: 232, G
Пример: n=2; x=3.
Ответ: j2(3)=0,298637497; y2(3)=-0,267038335.
Программа 63. Функции Кельвина bervx и beivx.
bervx = е(i)cos(3pv/4+pi/2)/(v+i)!/i!(x/2)v+2i,
i=0...г
beivx = е(i)sin(3pv/4+pi/2)/(v+i)!/i!(x/2)v+2i,
i=0...г
10 INPUT "n",A,"x",E:E=E/2:MODE
5
20 GOSUB 100:D=E-A/D:C=0:F=0:G=0:H=D
30 IF ABS D≤1E-11*ABS H;PRINT "ber";G,"bei";F:GOTO
10
40 FOR B=0 TO 1:F(B)=F(B)-D*SIN (p*(A/4+(C+B)/2)):NEXT
B
50 C=C+1:D=D*E*E/C/(C+A):GOTO 30
100 Программа 53
Размер: 372, H
Пример: n=2; x=3.
Ответ: ber23=0,808368452; bei23=-0,891022351.
Программа 64. Функции Кельвина bernx и beinx, kernx и keinx.
Область P0:
10 INPUT "n",A,"x",E:E=E/2:MODE
5:D=E-A:FOR B=1 TO A
20 D=D/B:NEXT B:C=0:F=0:G=0:H=D
30 IF ABS D≤1E-11*ABS H;PRINT "ber";G,"bei";F:GOTO
#1
40 FOR B=0 TO 1:F(B)=F(B)-D*SIN (p*(A/4+(C+B)/2)):NEXT
B
50 C=C+1:D=D*E*E/C/(C+A):GOTO 30
Область P1:
10 H=p/4*G+LN
E*F:F=p/4*F-LN
E*G:I=0:J=0:K=0:L=0
20 C=1:G=1:IF A>1;FOR D=2 TO A:C=C*D-C:G=G*D:NEXT D
30 D=1:M=.75*A*p:B=0:GOSUB
120
40 R=COS M:S=SIN M:T=(Q+P)/G*D:K=K+R*T:L=L+S*T:IF B=0;U=T
50 IF B<A;I=I+R*C*D:J=J+S*C*D:GOTO 70
60 IF U-T=U THEN 100
70 B=B+1:D=D*E*E:G=G*B*(A+B):GOSUB 120:M=M+p/2
80 IF B<A;C=C/B:IF B<A-1;C=C/(A-B-1)
90 GOTO 40
100 G=E-A:F=F+I/G/2+K*G/2:H=L*G/2-J/G/2-H
110 PRINT "ker";F,"kei";H:GOTO #0
120 N=B+1:GOSUB 130:Q=P:N=N+A:GOSUB 130:RETURN
130 P=.577215664:IF N=1;RETURN
140 FOR O=1 TO N-1:P=P+1/O:NEXT O:RETURN
Размер: 589, U
Пример: n=3; x=1
Ответ: ber31=0,0156287683;
bei31=-0,0137879837;
ker31=4,0858681054;
kei31=-6,8922228687.
Программа 65. Гипергеометрическая функция.
F(a,b,c,x)=1+ е(i)Х(j)x(a+j)(b+j)/(1+j)/(c+j), |x|<1< i=0...г, j=0...i-1.
Ряд сходится медленно. При плохой сходимости становится значительной вычислительная погрешность. Для точных вычислений используется рекурсивная формула, позволяющая скомпенсировать ошибку:
F(a,b,c,x)=y0.
yj=yj+1 x(a+j)(b+j)/(1+j)/(c+j)+1,
j=0...n-1;
yn=lg n;
Расчет начинается с последнего члена, поэтому для достижения требуемой точности следует произвести вычисления при различных n, наблюдая за устойчивостью результата.
10 INPUT "a",A,"b",B,"c",C,"x",D,"n",F:E=LOG
F
20 FOR F=F-1 TO 0 STEP -1:E=E*(A+F)/(C+F)*(B+F)/(1+F)*D+1
30 NEXT F:PRINT E:GOTO 10
Размер: 96, F
Пример: a=1; b=1; c=2; x=0,5; n=30.
Ответ: F(1; 1; 2; 0,5)=1,386294361.
Программа 66. Вырожденная гипергеометрическая функция. Разложение в ряд:
F(a,c,x)=1+ е(i)Х(j)x(a+j)/(1+j)/(c+j), |x|<1< i=0...г, j=0...i.
10 INPUT "a",A,"c",B,"x",C:D=1:E=1:F=0
20 G=D:E=E*(A+F)/(B+F)/(1+F)*C:D=D+E:IF D¦G;F=F+1:GOTO
20
30 PRINT D:GOTO 10
Размер: 93, G
Пример: a=1; c=1; x=0,5.
Ответ: F(1; 1; 0,5)=1,64872127.
Программа 67. Полные эллиптические интегралы.
K(m) = т0p/21/Ц(1-sin2q)
dq;
E(m) = т0p/2Ц(1-sin2q)
dq
Разложение в ряды:
K(m) = p/2{1+(1/2)2m+(1/2Ч3/4)2m2+(1/2Ч3/4Ч5/6)2m3+...};
E(m) = p/2{1 - (1/2)2m
- (1/2Ч3/4)2m2/3 - (1/2Ч3/4Ч5/6)2m3/5
-...};
10 INPUT A:B=-1:C=1:D=1:E=1:F=1
20 G=F:B=B+2:C=C*B/(B+1):D=D*A
30 E=E+C*C*D:F=F-C*C*D/B:IF G¦F THEN 20
40 PRINT "K";E*p/2,"E";F*p/2:GOTO
10
Размер: 115, G
Пример: m=0,5.
Ответ: K(0,5)=1,85407452; E(0,5)=1,35064387.
Программа 68. Эллиптические функции Якоби sn(u,m), cn(u,m), dn(u,m). Для вычисления используется процесс арифметикогеометрического среднего.
Начальные условия: a0=1; b0=Ц(1-m); c0=Цm.
Пока ai¦bi,
повторять: i=i+1;
ai=(ai-1+bi-1)/2, bi=Ц(ai-1bi-1),
ci=(ai-1 - bi-1)/2,
Обратный ход: n=i; jn=2nanu;
ji-1=(arcsin(ci/ai
sin ji)+ji)/2;
i=n...1.
Результат: sn(u,m)=sin j0; cn(u,m)=cos j0; dn(u,m)= cos j0/ cos (j1 - j0)
10 MODE 5:A=1:INPUT "m",B:B=SQR
(1-B):D=0
20 D=D+1:C=A+B:E(D)=(A-B)/C:B=SQR (A*B):A=C/2:IF A>B THEN 20
30 E=D:INPUT "u",C:C=2-D*A*C
40 B=C:C=(ASN (E(E)*SIN B)+B)/2:E=E-1:IF E>0 THEN 40
50 PRINT "sn";SIN C,"cn";COS C,"dn";COS C/COS (C-B):GOTO 30
Размер: 179, L
Пример: m=1-10-9; u=2.
Ответ: sn(2;1)=0,964027615; cn(2;1)=0,265802135; dn(2;1)=0,265802152.
Программа 69. Тэта-функция Якоби. Для вычисления используется процесс арифметико-геометрического среднего (см. описание программы 68). Результат:
X(u,m)=(cos(j1 - j0)/an/cos j0)1/2(1-m)1/4 Х(j)cos(2ji-1- ji)-(-1/2i+1), i=1...n
10 MODE 5:A=1:INPUT "m",B:B=SQR
(1-B):F=B:D=0
20 D=D+1:C=A+B:G(D)=(A-B)/C:B=SQR (A*B):A=C/2:IF A>B THEN 20
30 E=D:G=1:INPUT "u",C:C=2-D*A*C
40 B=C:C=(ASN (G(E)*SIN B)+B)/2
50 G=SQR (COS (C+C-B)*G):E=E-1:IF E>0 THEN 40
60 G=COS (B-C)/A/COS C*F/G:PRINT SQR G:GOTO 30
Размер: 195, N
Пример: m=1-10-9; u=2.
Ответ: X(2;1)=0,04880075718.
Программа 70. Дзэта-функция Якоби. Для вычисления используется процесс арифметико-геометрического среднего (см. описание программы 68). Результат:
Z(u,m) = е(i)ci sin ji, i=1...n
10 MODE 5:A=1:INPUT "m",B:B=SQR
(1-B):D=0
20 D=D+2:C=A+B:E(D)=(A-B)/C:F(D)=(A-B)/2
30 B=SQR (A*B):A=C/2:IF A>B THEN 20
40 E=D:INPUT "u",C:C=2-(D/2)*A*C:F=0
50 B=C:C=SIN B:F=F+F(E)*C
60 C=(ASN (E(E)*C)+B)/2:E=E-2:IF E>0 THEN 50
70 PRINT F:GOTO 40
Размер: 190, T
Пример: m=1-10-9; u=2.
Ответ: Z(2;1)=0,793784783.
Программа 71. Неполные эллиптические интегралы.
F(j,a) = т0j1/Ц(1-
sin2a sin2q)
dq;
E(j,a) = т0jЦ(1-
sin2a sin2q)
dq
Процесс арифметико-геометрического среднего.
Начальные условия: a0=1; b0=cos a; c0=sin a; j0=j.
Пока ai¦bi,
повторять:
ji+1= Arctg(tg ji
bi/ai)+ji;
ai+1=(ai+bi)/2, bi+1=Ц(aibi),
ci+1=(ai - bi)/2, i=i+1
Результат: n=i;
K(a)= p/2/an;
E(a)=K(a)(1-
е(i) 2i-1ci2);
i=0...n
F(j,a)= jn/2n/an;
Z(j,a) = е(i)ci
sin ji, i=1...n
E(j,a)= Z(j,a)+
F(j,a)E(a)/K(a)
10 MODE 5:INPUT "a",B,"ф",C:E=SIN
B-2:A=1:B=COS B:G=1:F=0
20 C=ATN (B/A*TAN C)+C+p*INT
(C/p+.5):G=G+G:D=(A-B)/2
30 E=E+G*D*D:F=F+D*SIN C:D=A+B:B=SQR (A*B):A=D/2
40 IF A>B THEN 20
50 C=C/G/A:A=p/2/A:B=A-A*E/2:F=F+B/A*C
60 PRINT "K";A,"E";B,"Fv";C,"Ev";F:GOTO 10
Размер: 209, G
Пример: a=80*p/180; j =85*p/180.