Human writable but not human readable

[waspagv]

Понадобился мне алгоритм преобразования полинома в разложение по многочленам Чебышева. Для чего? Объяснял школьникам, как вычисляются элементарные функции. Показал, что тейлоровское разложение плохое, потому что дает большую ошибку на краях интервала, что заставляет удерживать много членов ряда. Что можно экономизировать ряд и распределить ошибку, выбрав в качестве базовых функций многочлены Чебышева. И, конечно, теоретические рассуждения требуется подкрепить примером.

Так вот, возвращаясь к алгоритму. Чтобы долго не думать, взял готовую программу из БЧА НИВЦ МГУ. Но ведь надо объяснить, как она работает! А программа совершенно в духе 60-х - начала 80-х гг.!

      SUBROUTINE ZP42D(A,B,C,N)
      DIMENSION A(N),B(N),C(N)
      INTEGER N,I,K,L,M,O,P,Q,R
      DOUBLE PRECISION A,B,C,E
      L=0
      DO 7 I=1,N
      L=L+I
      P=L
      Q=L-I
      DO 6 K=I,N
      C(P)=0.D0
      IF(I-2) 1,2,3
    1 IF(K.GT.1) GO TO 4
      C(1)=1.D0
      GO TO 5
    2 IF(K.GT.2) GO TO 3
      C(3)=C(1)
      GO TO 5
    3 C(P)=2.D0*C(Q)
    4 IF(K.LE.I+1) GO TO 5
      C(P)=C(P)-C(R)
    5 R=Q+1
      Q=P-1
      P=P+K
    6 CONTINUE
    7 CONTINUE
      I=N
      M=N*(N+1)/2
      B(I)=A(I)
    8 B(I)=B(I)/C(M)
      IF(I-1) 9,11,9
    9 O=M-1
      M=M-I
      E=0.D0
      DO 10 K=I,N
      E=E+C(O)*B(K)
      O=O+K
   10 CONTINUE
      I=I-1
      B(I)=A(I)-E
      GO TO 8
   11 RETURN
      END

Жуть! В ранней молодости я тоже так писал, но потом исправился.

Эти три оператора все объяснили:

    3 C(P)=2.D0*C(Q)
    4 IF(K.LE.I+1) GO TO 5
      C(P)=C(P)-C(R)

Немного усовершенствовал, выбрав вычисление матрицы коэффициентов по столбцам, а не по строкам, как в оригинале. На стандартном Паскале для школьников получилось так:

const
   CConvPolyChebN = N;

type
   TConvPolyChebAr = array [0..CConvPolyChebN] of Real;

{====================================================================================}
{ Вычисление коэффициентов разложения заданного многочлена по многочленам Чебышева   }
{------------------------------------------------------------------------------------}
{ Параметры:                                                                         }
{    A - коэффициенты заданного многочлена A[0]+A[1]*x+A[2]*x^2...;                  }
{    B - вычисленные коэффициенты разложения по многочленам Чебышева                 }
{        C[0]*T_0(X)+C[1]*T_1(x)+...                                                 }
{    N - степень полинома (N >= 0).                                                  }
{ A и B могут быть одним массивом.                                                   }
{====================================================================================}
procedure ConvPolyCheb(var A,C: TConvPolyChebAr; n: Integer);
var
   M: array [0..CConvPolyChebN] of TConvPolyChebAr;
   i,k: Integer;
   s: Real;
begin
   { Формируем матрицу перехода }
   M[0,0] := 1.0;

   if n > 0 then
      begin
         M[0,1] := 0.0;
         M[1,1] := 1.0;
      end;

   for i := 2 to n do
      begin
         M[0,i] := -M[0,i-2];
         for k := 1 to i-2 do
            M[k,i] := 2*M[k-1,i-1] - M[k,i-2];
         M[i-1,i] := 2*M[i-2,i-1];
         M[i,i] := 2*M[i-1,i-1]
      end;

   { Обратный ход для вычислений коэффициентов разложения Чебышёва }
   C[n] := A[n]/M[n,n];
   for k := n-1 downto 0 do
      begin
         s := A[k];
         for i := k+1 to n do
            s := s - M[k,i]*C[i];
         C[k] := s/M[k,k]
      end;
end;