solution Rosen(matrix x0, matrix s0, double alfa, double beta, double epsilon, int Nmax, matrix O)
{
    //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();
    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();
            if (Xt.y < X.y)
            {
                X = Xt;
                l(i) += s(i);
                s(i) *= alfa;
            }
            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));
        if (solution::f_calls > Nmax || max_s <= epsilon)
            return X;
    }
}