для о < si,...,s„ < 1, и этот метод можно рассматривать, как одновременное вычисление 2" линейных полиномов от 2" переменных F{ti,..., t„). Общеизвестный технический прием, предложенный Ф. Ятсом (F. Yates) [The Design and Analysis of Factorial Experiments (Harpenden: Imperial Bureau of Soil Sciences, 1937)], можно использовать для уменьшения числа сложений в (38) от 2"(2" - 1) до п2". Метод Ятса можно понять, рассмотрев случай, когда п = 3: пусть Xtit2t3 - F{ti,t2,t3).

Для вычисления от колонки "Задано" к колонке "Первый шаг" требуется четыре сложения и четыре вычитания; интересной особенностью метода Ятса является то, что точно такое же преобразование, которое действует от колонки "Задано" к колонке "Первый шаг", будет действовать от колонки "Первый шаг" к колонке "Второй шаг" и от колонки "Второй шаг" к колонке "Третий шаг". В каждом случае вьшолняется четыре сложения и затем четыре вычитания, а после трех шагов магически получается требуемое преобразование Фурье /(si,52,83) на месте, первоначально занимаемом F{si,S2, s3).

Этот особый случай часто называют преобразованием Уолша 2" данных элементов, поскольку соответствующая модель знаков изучена Дж. Л. Уолшем (J. L. Walsh) [Amer. J. Math. 45 (1923), 5-24]. Заметим, что число перемен знаков слева направо в колонке "Третий шаг" принимает соответствующие значения

О, 7, 3, 4, 1, 6, 2, 5.

Это перестановки чисел {0,1,2,3,4,5,6,7}. Уолш заметил, что будет точно О, 1, ..., 2" -1 изменений знаков в общем случае, если соответствующим образом переставить преобразованные элементы. При этом коэффициенты обеспечивают дискретную аппроксимацию синусоидами с различными частотами. (См. работу Н. F. Harmuth, IEEE Spectrum 6,11 (November, 1969), 82-91, в которой речь идет об использовании этого свойства; дополнительно коэффициенты Уолша обсуждаются в разделе 7.2.1.)

Метод Уолша можно обобщить для оценки любого дискретного преобразования Фурье и фактически для оценки любого числа сумм, которые можно записать в общем виде

/(Si,S2,...,S„) =

53 9l{Sl,S2,-.-, Sn, h) 92(82, ,Sn,t2) gn{Sn,tn)F{ti,t2, ,tn) (39)

0<«i<mi 0<«„<m„

для о < Sj < rrij и заданных функций gj{sj,... ,Sn,tj). Мы поступим следующим образом:

/о(*1,2,3, •itn) = F{ti,t2,t3, ... ,tn);

fl{Sn,tl,t2,...,tn-l) = 9niSn,tn)foitl,t2,---,tn);

0<t„<m„

/2(s„-l,S„,ti, . . .,tn-2) = 9n-liSn-l,Sn,tn-l)fl{Sn,tl, • ,*n-l);

0«„ i<m„ i

/n(Sl,S2,S3,...,S„) = ffl(Sl,...,S„,ti)/„ l(S2,S3,...,S„,fi); 0<«i<mi

/(SbS2,S3, •••,««) = /n(si,S2,S3,...,S„). (40)

Для метода Ятса, как показано выше, gj{sj,..., Sn,tj) = (-l)**; /о(*1,*2)*з) представляет колонку "Задано", /1(83,1,2)-колонку "Первый шаг" и т. д. Всякий раз, когда множество сумм можно представить в виде (39) для умеренно простых функций gj{sj,... ,Sn,tj), схема (40) будет уменьшать количество вычислений от порядка Л до порядка NlogN или около того, где N = mi ...m„-количество данных точек. К тому же данная схема идеально подходит для параллельного вычисления. Важный особый случай одномерного преобразования Фурье обсуждается в упр. 14 и 53; одномерный случай также приводится в разделе 4.3.3С.

Рассмотрим более частный случай вычисления полинома. Интерполяционный полином Лагранжа порядка п, который мы запишем в виде

и, ,(х)=и ix-Xi){x-X2)...ix-Xn) {х-Хо)(х-Х2)...(х-Хп)

°{Хо-Хг){Хо-Х2)...{Хо-Хп)(Х1-Хо){Х1-Х2)...{Х1-Хп)

. „ ix-XQ)ix-xi)...ix-Xn-i) , .

"{Хп-Хо){Хп-Хг)...{п-Хп-1Г

является только полиномом от х степени < п, который принимает соответствующие значения уо, 2/1, • • • i Уп в п -Н1 различной точке х = Xq,Xi,. .. ,Хп- (Из (41) очевидно, что U[n]{Xk) = Ук для О < fc < п. Если fix) - любой такой полином степени < п, то д{х) = f(x) - U[n](x) имеет степень < п и д{х) равно нулю для х = Xo,Xi,... ...,Хп\ следовательно, д{х) должен быть кратен полиному {x-xo){x-Xi)... {х-Хп)-Степень последнего полинома больше п, поэтому д{х) = 0.) Если предположить, что значения функции табулированы и хорошо аппроксимируются полиномом, то формулу (41) можно использовать для "интерполирования" значений функции в точках X, не занесенных в таблицу. Лагранж предложил (41) своим ученикам в Paris Ecole Normale в 1795 году [см. CEuvres 7 (Paris, 1877), 286], однако Эдвард Уоринг (Edward Waring) из Кембриджского университета к тому времени уже представил такую же формулу, совершенно точно и ясно сформулированную в Philosophical Tr&nsactions 69 (1779), 59-67.

На первый взгляд кажется, что в формулах Уоринга и Лагранжа совсем немного сложений, вычитаний, умножений и делений; на самом деле - точно п сложений.

2n + 2n вычитаний, 2n + n - 1 умножений и n + l делений. Но, к счастью (как мы подозревали), возможно улучшение.

Основная идея упрошения (41) учитывает тот факт, что Щп]{х) - W[„ i](a;) = О для X = Хо, ,Xn-i; таким образом, Щп]{х) - U[n-i]{x)-полином степени п или .меньше, кратный {х - Xq) .. .{х - x„i). Можно сделать вывод, что

U[„]{x) = ап{Х -Хо)...{х - Xn-l) + Щп-1]{Х),

где а„ - константа. Это приводит к интерполяционной формуле Ньютона

Щп]{х) = ап{х - хо){х -xi)...ix - Xn-i) +

+ а2{х - Хо){х - xi) +ai{x - xq) + ао, (42)

где ak-коэффициенты, которые необходимо определить по заданным числам хо, xi, Хп, Уо, У1, •) Уп- Заметим, что эта формула имеет место для всех п; коэффициент ак не зависит от Xk+i, ..., Хп или от ук+г, ,Уп- Когда ак известны, интерполяционная формула Ньютона удобна для вычисления, так как можно снова обобщить правило Горнера и записать

U[n]{x) = ((... {an{x-Xn-i) + an-i){x-Xn-2) + ){х-хо) + ао). (43)

Это потребует п умножений и 2п сложений. Иначе можно вычислить каждый отдельный член (42) справа налево; выполнив 2п - 1 умножений и 2п сложений, мы вычислим все значения ищ{х), ищ{х), ..., U[n]{x), что во всяком случае покажет, будет ли сходиться интерполяционный процесс.

Коэффициенты ак в формуле Ньютона можно найти, вычисляя отношения разностей по следующей таблице (показано для п = 3).

( )/

L-vi)iZ-Z)к -i;/;=(у-ужг-хо)=уг

У2 . ... . , {Уз-У2)/{хз-Х1) =Уз

у, Уг-У2)1{хг-Х2) = уз (44

Можно доказать, что «о = уо, = yJ, аг = Уг и т. д., и показать, что отношения разностей имеют тесную связь с производными интерполируемой функции (см. упр. 15). Значит, следующие вычисления (соответствующие (44)) можно использовать, чтобы получить ак-

Начать с (ао,а1,...,а„)-ь-(уо,у1,...,у„); затем для А; = 1,2,...,п(в таком порядке)

присвоить aj{aj-aj-i)l{xj-Xj-k) для j = п,п - 1,...,fc (в таком порядке).

Для этого процесса требуется \{п + п) делений ип +п вычитаний, и экономится около трех четвертей работы по сравнению с (41).