Есть алгоритм вычислений весов точек для интерполяции Лагранжа, вот такой феньки: prod=(i-pos[j1])(i-pos[j2])(i-pos[j3])... i<>pos[j] for (i=0; i<n; i++) { t=1; for (j=0; j<k; j++) { if (i<>pos[j]) t*=i-pos[j]; } prod=t; } когда наполнение pos[] тривиально, т.е. pos[j]==j, то можно вычислить все веса итеративно за линейное время вместо квадратичного, а именно, prod=prod[i-1]*i/(i-k) Как вычислить за линейное время в общем случае, если наполнение pos[] это наполнение некоторыми произвольными k числами из диапазона 0..n-1? Весь моск уже обламал себе...
prod=(i-pos[j1])(i-pos[j2])(i-pos[j3])... i<>pos[j] - это функция k-того порядка. Любую функцию k-ого порядка можно расписать на k сложений и таким образом избавиться от умножения (на английском это называется "forward differencing"). В прикреплённом файле наглядный пример. Как вычислить за линейное время не знаю.
похоже, что так, только не просто складывать, а еще и прогнать через Быстрое Преобразованиe Волша... Но для этого нужно быть доктором философии в математике =(((. Мне остается только тыбрить чужой код... Что, простых каких нить ходов нет? Яж придумал для pos[j]==j
Ты всегда говоришь какие-то сложные вещи: Рассмотрим функцию из xls-файла в посте №2 f(i)=(i-87)*(i-45)*(i-89) ты мог бы вычислять её так: Код (Text): for i:=1 to 10 do prod[i]:=(i-87)*(i-45)*(i-89); Но лучше сделать так: Код (Text): C:=-442; B:=15443; A:=-332992; for i:=1 to 10 do begin prod[i]:=A; C:=C+6; B:=B+C; A:=A+B; end; При этом инкремент по C равен n! где n - порядок функции.
Сейчас для вычисления B и C я вычисляю n первых шагов функции напрямую (n - порядок), и делаю (n-1)*(n-2)/2+n-1 сложений.
Прикольно! Осталось только проспаться и понять, как там через производные взялись эти А,B и C. То бишь ты предлагаешь для сложности k*O(n^2) уменьшить k, а квадрат оставить? В принципе, в моей проге мне вообщем то все-равно, складывать или умножать, а вот квадрат-это принципиально, поскоку миллион^2 сам понимаешь скока получается... А сам по себе миллион для компа - просто тьфу на палочке!
Обычная арифметика. В прикреплённом файле попытался расписать поподробнее. Как я уже писал погугли "forward differencing". И вот ещё статья русского автора (но на английском) с исходниками http://antigrain.com/research/bezier_interpolation/index.html
murder самую главную изюминку лагранжевой интерполяции я просек однако! Для того чтобы не сравнивать (i<>pos[j]) дофига и дофига раз за квадратичное время, нужно выписать полином (x-pos[1])(x-pos[2])... и взять от него производную! Потом подставлять икс уже в производную.
Остаются две задачки, вычисление полинома вне узлов и вычисление его производной в узлах (поскоку сам полином в узлах даст просто ноль)
murder Задачка наконец успещно рещена =))) Код (Text): pos=[2 3 8]; %Это за квадратичное время for i=1:8, t(i)=1; for j=1:3. if i~=pos(j), t(i)=t(i)*(i-pos(j)); end; end; end; t %А это за n log n =))) x=[0 1 2 3 4 5 6 7 0 -7 -6 -5 -4 -3 -2 -1]; m=[0 1 1 0 0 0 0 1 0 0 0 0 0 0 0 0]; x=log(x); x(1)=0; x(9)=0; t_fast=exp(ifft(fft(x).*fft(m))) То бишь вместо умножения берется сложение логарифмов, а именно вектор логарифмов сворачивается с маской через БПФ. Но логарифмы отрицательных чисел - комплексные, дискретные тоже не всегда доступны - вспоминаем open key криптографию... Можешь тут вместо логарифмов применить свои волшебные суммы и показать, что это работает со сверткой?
Х.з. ты в математике разбираешся гораздо лучше меня. БПФ я пока не изучал. Не знаю погогут ли эти формулы ln(x)=x/1-x^2/2+x^3/3-... +(-1)^(n+1)*x^n/n exp(x)=1+x/1+x^2/2+x^3/6+...+x^n/n!
