solution HJ(matrix(*ff)(matrix, matrix, matrix), matrix x0, double s, double alpha, double epsilon, int Nmax, matrix ud1, matrix ud2)
{
    try
    {
        solution XB, XB_old, X;
        XB.x = x0;
        XB.fit_fun(ff, ud1, ud2);
        while (true)
        {

            X = HJ_trial(ff, XB, s, ud1, ud2);
            cout << X.x(0) << " " << X.x(1) << endl;
            if (X.y < XB.y)
            {
                while (true)
                {
                    XB_old = XB;
                    XB = X;
                    X.x = XB.x + XB.x - XB_old.x;

                    X.fit_fun(ff, ud1, ud2);
                    X = HJ_trial(ff, X, s, ud1, ud2);
                    if (X.y >= XB.y)
                        break;
                    if (XB.f_calls > Nmax)
                        return XB;
                }
            }
            else
                s *= alpha;
            if (XB.f_calls > Nmax || s < epsilon)
                return XB;
        }
    }
    catch (string ex_info)
    {
        throw ("solution HJ(...):\n" + ex_info);
    }
}

solution HJ_trial(matrix(*ff)(matrix, matrix, matrix), solution XB, double s, matrix ud1, matrix ud2)
{
    try
    {
        int* n = get_size(XB.x); //rozmiar problemu
        matrix D(n[0], n[0]); //macierz kierunku
        for (int i = 0; i < n[0]; i++) {
            D(i, i) = 1;
        }
        solution X;
        for (int i = 0; i < n[0]; ++i)
        {
            X.x = XB.x + s * D[i];
            X.fit_fun(ff, ud1, ud2);
            if (X.y < XB.y)
                XB = X;
            else
            {
                X.x = XB.x - s * D[i];
                X.fit_fun(ff, ud1, ud2);
                if (X.y < XB.y)
                    XB = X;
            }
        }
        return XB;
    }
    catch (string ex_info)
    {
        throw ("solution HJ_trial(...):\n" + ex_info);
    }
}

solution Rosen(matrix(*ff)(matrix, matrix, matrix), matrix x0, matrix s0, double alpha, double beta, double epsilon, int Nmax, matrix ud1 , matrix ud2)// throw (string);
{
    //alfa >0
    //beta<0
    int* n = get_size(x0);
    matrix l(n[0], 1), p(n[0], 1), s(s0);
    matrix D(n[0], n[0]); //macierz kierunku
    for (int i = 0; i < n[0]; i++) {
        D(i, i) = 1;
    }
    solution X, Xt;
    X.x = x0;
    X.fit_fun(ff, ud1, ud2);
    while (true)
    {
        cout << X.x(0) << "\n" << X.x(1) << endl;
        for (int i = 0; i < n[0]; ++i)
        {
            Xt.x = X.x + s(i) * D[i];
            Xt.fit_fun(ff, ud1, ud2);
            if (Xt.y < X.y)
            {
                X = Xt;
                l(i) += s(i);
                s(i) *= alpha;
            }
            else
            {
                p(i) += 1;
                s(i) *= (-beta);
            }
        }
        bool change = true;
        for (int i = 0; i < n[0]; ++i)
            if (l(i) != 0 || p(i) != 0) // l || p
            {
                change = false;
                break;
            }
        if (change)
        {
            matrix Q(n[0], n[0]), v(n[0], 1);
            for (int i = 0; i < n[0]; ++i)
                for (int j = 0; j <= i; ++j)
                    Q(i, j) = l(i);
            Q = D * Q;
            v = Q[0] / norm(Q[0]);
            D.set_col(v, 0);
            for (int i = 1; i < n[0]; ++i)
            {
                matrix temp(n[0], 1);
                for (int j = 0; j < i; ++j)
                    temp = trans(Q[j]) * D[j] * D[j];
                v = Q[i] - temp / norm(Q[i] - temp);
                D.set_col(v, i);
            }
            s = s0;
            l = 0;
            p = 0;
        }
        double max_s = abs(s(0));
        for (int i = 1; i < n[0]; ++i)
            if (max_s < abs(s(i)))
                max_s = abs(s(i));
            cout<<X<<endl;
        if (solution::f_calls > Nmax || max_s <= epsilon)
            return X;
    }
}