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

Свидетельство к алгоритму 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 согласно следующему:

6, град

Подтверждение к алгоритму 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