Программа была проверена на .вычислительной машине XI математического центра. Для рЫ=4Б°, k=:sinlO°, sm20°,..., sinl80° были вычислены Е и F. Результаты содержали 12 значащих цифр.

Сравнение с 12-разрядными десятичными таблицами Лежандра - Эмде (1931 г.) показало, что 12-я цифра была получена с ошибкой самое большее в четыре единицы. После Юмин счета (т. е. после более чем 100 циклов) не было получено никаких результатов для k=sin89°, phi=\°, так как произошло зацикливание.

Замечание. Поскольку phi не меняется во время вычисления, то мы поставили оператор cos phi: = cos (phi) ,в начало программы, чтобы предотвратить вычисление косинуса 30 или более раз. Кроме того, в выражении для Г[1] и Г(2] значение sqrt{l-sinphi\2) было заменено на cosphi, так как потери значащих цифр не происходит.

Выражение 2Хп было заменено новой переменной для ускорения работы программы.

АЛГОРИТМ 746

Аппрокеимвция с помощью полиномиальной кривой данной степени, проходящей через данные точки (метод наименьших квадратов] [Е2]

Процедура curfit {curve - кривая, fitting - вычерчивание по точкам) находит методом наименьших квадратов полином степени п {knk+m), график которого проходит через точки (ai, bi), (%, bk) и аппроксимирует точки {xi,yi),..., {Хт, Ут), причем Wi - это вес, .приписываемый точке {хиУг)- Процедура curfit использует нелокализованные метки signal74 и signalTAg и стандартную процедуру outreal{\,<A>), выводящую из машины по каналу 1 значение арифметического выражения <Л>.

Детальное описание алгоритма приводится в работе Пека [lOi], где используются такие же обозначения, что и в данной процедуре *.

procedure curfit (k,a,b,m,x,y,w,n,gamma,xO,z,r) result: (alpha,beta,s,sgmsq,c);

value k,m,n,r,xO,gamma; real x0,gamrna; integer k,m,n,r; array a,b,x,y,vir,alpha,beta,s,sgmsq,c,z; begin real p,f,lambda; integer i,j; array virlfhk];

comment Сначала определяются несколько процедур, которые используются в основной программе, начинающейся с метки start;

procedure evalue(x,nu); value x,nu; real x; integer nu;

comment Процедура evalue вычисляет значение f===soPo+

Sipi+ ... + sp, где pi+i(x) = (x-tti)pi(x)-pipi , (x) для

i=0, 1,..., V-1, p i(x)=0, po(x)==l, po==0. Значение p, (x)

остается в p;

begin real pO,t; integer i;

pO:=0; p:=l; f: = slO];

* Назначение основных параметров процедуры curfit раскрывается также и в нижеследующем «Свидетельстве к алгоритму 746». (Прим. ред.)

for i: = 0 step 1 until nu-1 do begin t: = p; p:= (x-alpha[i]) Xp-beta[i]XpO; pO:=t; f: = f+pXs[i + l] end i end evalue; procedure coda (n,c);

value n; integer n; array c;

comment Эта процедура находит коэфф-ициенты с такие, что с0+с1х+ ... +CnX"=sopo(x) + ... +snpn(x); begin real tl,t2; integer i,r; array pm,p[0:n]; for r: = 1 step 1 until n do clrl:=pm[r]: = p[r]:=0; c[0]: = s[0];pm[O]:=0;p[0]:=l; for i:=0 step 1 unti! n-1 do begin t2:=0; for r:=0 step 1 unti! i + l do begin

tl: = (12-alpha[i] X p{r]-beta[ilX pmH) /lambda; t2: -pm[r]:=-pH; p[rj:=tl; cM: = c[r]-btl.Xs[i+,l] end r end i end coda;

procedure gefyt(n,nO,x,y,w,m); [

value n, nO,m; integer n,nO,m; array x,y,w;

comment Это сердце основной процедуры. Процедура gefyt вычисляет аи pi, si, а\ используя метод ортогональных полиномов, описанный в вышеуказанной работе Пека; begin геа! dsq,wpp,wppO,wxpp,wyp,t;

integer i,j,freedom; Boolean exact; array p,pO[l:mj;

if n-nO>m Vn<nO then go to signal74g;

beta[nO]: = dsq:=wpp:=0;

exact: = n-nO m- 1;

for j:=l step 1 unti! m do begin p[j]: = l; pO0]:=0; wpp:=wpp-fw[j];

if-lexact then dsq: = dsq+wlj]XyO]XyDl end расчета начальных условий; for i: = nO step 1 unti! n do

begin freedom: = m-1-(i-nO); wyp:=wxpp:=0; for ]: = 1 step 1 unti! m do begin t:=w[jlxp[jl;

if i<n then wxpp:=wxpp+1X x[j]X p[j]; if freedomO then wyp:==wyp+tXyO] end j;

if freedomO then s[i]:=wyp/wpp; if "~i exact then

begin dsq:==dsq-sn]Xs[i]Xwpp; sgmsq[i]: = dsq/f reedom

end if; if i<n then

begin alpha[i]:=wxpp/wpp; wppO:=wpp; wpp:=0; for j:= 1 step 1 until m do begin

t: = (x[j]-alpha[i]) X pljl-beta[i]x pO[j]; wpp: = wpp+wIj]XtXt;

pO0]: = PJ]; P[j]: = t end j;

beta[i+1]:=wpp/wppO end if end i end gefyt;

comment Здесь начинается основная программа; start: for j: = l step 1 until к do

begin wlli]: = 1; a[j]:= (aO]-xO)/gamma end j; gefyt (k,0,a,b,wl,k);

comment Здесь вычисляется полином степени к-1, график которого проходит через точки (ai,bi),..., (ak,bk), и аь рь Si, где 0<i<k; begin real rho; rho:=0;

for j: = l step 1 until m do

begin rho: = rho+w[j]; x[j]:= (x[j]-xO)/gamma end j; rho: = m/rho;

comment Множитель p используется для нормализации весов.. Теперь для вычисления значения рк(х) и одновременно полинома степени к-1 положим Sk=0; s[k]: = 0;

for j: = l step 1 until m do begin evalue(x[j],k);

if p = 0 then go to signal74; y[i]: = (y[i]-0/p; w[j]: = w[j]XpXpXrho end j . •

end rho;

comment Теперь веса нормализованы, и отрегулированные веса

и ординаты готовы для аппроксимации методом наименьших

квадратов;

gefyt (п,к,х,у,иг,т);

comment Коэффициенты ш, pi (0i<Ti) и si(Oin) теперь готовы. Полином может быть вычислен для x=zi, zi, . .., Zr, но переменная должна быть сначала отрегулирована. Заметим,, что мы можем, уменьшая п, вычислить наилучший полином; меньшей степени; begin real xl;

for j:= 1 step 1 until r do

begin xl: = (z{j]-xO)/gamma; evalue(xl,n); outreaI(l,z[j]); outreal(l,f) end j;

comment Теперь мы можем отрегулировать весовые коэффициенты и затем найти коэффициенты для степенного ряда cn-bcix-b ... -bcnx"=sopo(x)s-b ... +SnPn(x);