Метод Пауэлла

#include <math.h>
double Powell(double *a, double *b);
// Программа поиска минимума на интервале [a, b]
// методом квадратичной аппроксимации Пауэлла.
// Значение функции f(x) вычисляется в подпрограмме Fx(x)
// В результате интервал неопределённости уменьшается и
// его границы доступны в вызывающей программе.
// Также «наверх» передаётся минимальное значение функции.
extern double Fx(double);
double dx, x[4], f[4], xt, ft, dn, nm, tolerance = 0.001;
int j, k, s1, s2, s3;
double dx=0.05; // шаг
// Начать процесс с первых трёх точек
x[0] = *a;        f[0]=Fx(x[0]);
x[1]=*a + dx; f[1] = Fx(x[1]);
if (f[0] < f[1]) x[2] = *a - dx;
 else x[2] = *a + 2.0 * dx;
f[2] = Fx(x[2]);
//Первый интерполированный минимум
dn=(x[1] - x[2]) * f[0];
dn+=(x[2] - x[0]) * f[1] + (x[0] - x[1]) * f[2];
nm=(x[1] * x[1] - x[2] * x[2]) * f[0];
nm+=(x[2] * x[2] - x[0] * x[0]) * f[1];
nm+=(x[0] * x[0] - x[1] * x[1]) * f[2];
×[3]=nm / (2.0 * dn); f[3] = Fx(x[3]);
while(1) {//Бесконечный цикл интерполяции
    for(j = 0; j < 3; j++)
        for (k = 0; k < 4; k++) {
            if (f[j] > f[k]{//Упорядочить x и f
            xt = x[j]; x[j] = x[k]; x[k] = xt;
            ft = f[j]; f[j] = f[k]; f[k] = ft;
                }//Конец оператора if( )
        } // Конец циклов по ј и k
    if (fabs(x[0] - x[1]) < tolerance) {
        // Точность по величине интервала достигнута.
        // Выход из цикла интерполяции.
        goto EXIT;
            // Точность по интервалу не достигнута,
            // создаём новую точку
            s1 = sgn(x[1] - x[0]); //sgn() = 1, если скобки > 0
        s2 = sin(x[2] - x[0]); //sgn() = 0, если скобки =
        s3 = sgh(x[3] - x[0]); //sgn() = -1, если скобки < 0
        if (s1 == s2 && s1 == (-s3)) {
            x[2] = x[3];
            f[2] = f[3];
            // Делаем вторую интерполяцию и
            // переходим на упорядочивание
            dn = (x[1] - x[2]) * f[0] + (x[2] - x[0]) * f[1] + (x[0] -
                ×[1]) * f[2];
            ft = (f[0] - f[1]) / (2.0 * dn);
            ft *= (x[1] - x[2]) * (x[2] - x[0]);
            x[3] = (x[0] + x[1]) / 2.0 + ft;
            f[3] = Fx(x[3]);
        }//Конец бесконечного цикла интерполяции
    EXIT:
        *a = x[0];
        *b = x[1];
        return (f[0]);
    }