Главная Промышленная автоматика.

if (abs(kX;l) Д (abs(phi) < 1.57079633) then go to large; ee: = ff: = io8; go to fin; small: al:=phi;

stepl: n: = n+2; ee:=(n-l)/n; tl:=kt2Xhl; h2:=eeXtl;

a2: = eeXal-snf (n-1) Xcs/n; dl:=:h2Xa2; d2:=-tlXa2/n; sl:=sH-dl; s2: = s2-bd2; hl:=h2; al: = a2;

if abs((sl + dl)-sl)>0Aabs(phiXsntn)abs(a2) then

go to stepl; ff: = phi + sl; ee: = phi+s2; go to fin;

kp: = l-к2; al: = l;

pl:=ml: = nl:=0; tl: = abs(k) Xcs; t2:= 1-(kXsn) 2; n:=n + 2; ee:=(n-l)/n; ff: = tlXt2t((n-l)/2)/n; h2: = eeXhl; а2: = ее2ХкрХа1; p2: = pl + l/((n-l)X(n/2)); m2:= (ml-ffXh2) X ((n+ l)/(n+2))\2Xkp; n2:= (nl-ffXhl) XeeX (n+1) Xkp/(n+2); dl:=m2-a2Xp2;

d2:=n2+kpx (al/nt2-eeXalXp2); d3: = a2; d4:= (n+I) Xa2/(n+2);

sl: = sl+dl; s2: = s2 + d2; s3: = s3+d3; s4:=s4+d4-; hl: = h2; al: = a2; pl: = p2; ml:=m2; nl:=n2;

if abs((sl+dl)-sl)>0 then go to step2; nl: = sqrt(t2); n2: = abs(kxsn); t2: = tl/nl; tl: = ln(4/(nl+tl)); ff:=tlX(I+s3)+t2Xln(0.5X(l + n2))-bsl; ee:= (0.5+s4) XkpXtl +1-t2X (1-n2) +s2; if phiO then go to fin; ee:=-ее; ff:=-ff; fin: end ellint;

Свидетельство к алгоритму 736

Алгоритм 7Эб получен в результате дальнейшей оптимизации алгоритма 73а путем замены некоторых неоднократно повторяющихся выражений простыми переменными. Например, в строке, начинающейся с метки «stepl:», добавлен оператор й:=2Хп1; что дает возможность заменить в последующих строках выражение 2Хл1 на простую переменную tl. Аналогично простыми переменными были заменены выражения abs(k)Xcs и 1-(Xsft)f,2 (начиная со второй строки после метки «large:») и выражения sqrt{\-(Xsft)f2) и abs(kXsn) (начиная со строки, следующей за оператором if G6s((sl-bdl)-sl)>0 then go to step2;). Кроме того, после метки «large:» оператор



kp:-sqrt{\-k\2) был заменен оператором kp\ = \-fef2, что дало возможность всюду в процедуре заменить выражение kp\2 на kp.

Алгоритм 736 был транслирован на машине М.-220 в системе ТА-IM, и были вычислены значения Р(ф, sin а), £({р, sin а) для ф=0°, Г, 10°, 20°,..., 90° и а=0°, 1°, 10°, 20°,..., 80°, 89°, 90°. Результаты трансляции совпали с контрольными значениями взятыми из семизначных таблиц В. П. Ветчинкина [39], за исключением F(60°, sin60°) == 1.212597 вместо 1.212537 у Ветчинкина, что объясняется, по-видимому, опечаткой. Некоторые из полученных значений приведены в табл. 8. Следует заметить, что при ф=1° и а=89° не наблюдалось никакого зацикливания, о котором говорится в нижеследующем «Подтверждении к алгоритму 73» Р. П. Ван де Рита, и были получены правильные результаты.

Таблица 8

Е (<f, sin а)

F {if, sin а)

Е (<f, sin а)

F {а, sin а)

0.52359877

0.52359877

0.017452406

0.017454178

0.52359187

0.52360567

0.017453071

0.017453514

0.52291511

0.52428401

0.17431249

0.17475385

0.51408861

0.53342745

0.68506023

0.71164727

0.50286804

0.54593191

1.1631768

1.2853005

0.50000000

0.54930614

1.4674622

1.6857503

Свидетельство к алгоритму 73а

Алгоритм 73а получен в результате исправления, оптимизации и ординарной переработки алгоритма 73 (J е f f е г s о п D. К. «САСМ», 1962, № 12). Перевод замечания Д. Джефферсона («САСМ», 1962, № 10) не приводится, поскольку оно не добавляет новой информации к остальным материалам по алгоритму 73.

Кроме ошибок, указанных в приводимых ниже «Подтверждениях» Н. А. Мейера и Р. П. Ван де Рита, в алгоритме 73 исправлено следующее.

1. Выражение я/2 заменено его значением 1.57079633, а tanh{\) заменен на g"(l) = 1.557393 согласно справочнику [19].

2. Переменная п отнесена к типу integer, иначе для рЫ<.0 выражения, подобные sn\{n-1), были бы не определены.

Переменные алгоритма 73а соответствуют переменным алгоритма 73 согласно следующему:

- delll]

- sinphi

- slgma[\\

- del [2]

- cosphi

- Sigma [2]

-т\\\

- del[3]

~H{\\

- sigvm [3]

-г [2]

- deim

-H[2\

- Sigma [4]

-А\\\

-L[l]

~M[\\

- £

-L12]



6, град

<Р. грая

-1X10-S

-1X10-

-IXIO-

-3X10-S

1X10-

2X10-

-3X10-S

2X10-

6Х10-«

Таблица 10

Ф, град

е, град

-1X10-

-1X10-

-1X10-

-1X10-

1X10-S

1X10-"

-7X10-S

3X10-"

1XI0-S

Подтверждение к алгоритму 73

Р. П. Ван де Рит (Van de Riet R. P. «САСМ», 1963, № 4) Алгоритм содержал три опечатки ... **

* Далее указываются три ошибки и боз1МОЖ,вость усовершенстБОвания алгоритма 73, учтенные в алгоритме 73а. (Прим. ред.)

** Далее указываются три опечатки в алгоритме 73, учтенные в алгоритме 73а. (Прим. ред.)

Д. Крайбл (Kriebel D. С. «САСМ», 1961, № 12)

Этот алгоритм был первоначально запрограммирован на языке машины NORC, и были получены таблицы неполного эллиптического интеграла первого и второго рода [2i].

Алгоритм был закодирован для транслятора MAD точно так, как он написан на языке АЛГОЛ и проверен на машине IBM 7090. Было сосчитано 40 случаев для k и рЫ, меняюшихся от 0° до 90°. Результаты выдавались с восемью значащими цифрами и совпали с таблицами Дидонато и Хершея с точностью от нуля до двух единиц восьмой цифры. (Это могло произойти из-за переводов из десятичной в двоичную и из двоичной в десятичную систему счисления при вводе и выводе, используемых в двоичной вычислительной машине. Результаты сравнивались с результатами прямого десятичного счета на машине NORC.)

Подтверждение к алгоритму 73

Н. А. Мейер (Meyer N. А. «САСМ», 1963, № 2)

Процедура elUnt была вручную закодирована на языке ФОРТРАН для машины IBM 7070. Были сделаны следующие исправления ... *

Исправленный алгоритм дал удовлетворительные результаты при сравнении с таблицами Дидонато и Хершея. Отличия встречались в восьмой значащей цифре, как показано в табл. 9 (интеграл f) и 10 (интеграл Е).

Таблица 9





0 1 2 3 4 5 6 7 8 9 10 11 12 [13] 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43

0.0037