#if LAB_NO==2 && LAB_PART==2
		(*ud).add_row((B.x - A.x)());
#endif
			if (B.x - A.x < epsilon || abs(D.x()-D_old.x()<gamma))
			return C;
		D_old = D;
	}
}
#endif
#if LAB_NO>2
//sprawdzam w jednej płaszczyźnie poźniej w drugiej i przesuwam się tam gdzie jest lepiej
solution HJ(matrix x0, double s, double alpha, double epsilon, int Nmax, matrix *ud, matrix *ad)
{
	solution XB(x0), XB_old, X;
	XB.fit_fun(ud, ad);
	while (true)
	{
		X = HJ_trial(XB, s, ud, ad);
		if (X.y <XB.y)
		{
			while (true)
			{
				XB_old = XB;
				XB = X;
#if LAB_NO==3 && LAB_PART==2
				? ? ?
#endif
				X.x = 2 * XB.x - XB_old.x;
				X.fit_fun(ud, ad);
				X = HJ_trial(X, s, ud, ad);
				if (X.y >= XB.y)
					break;
				if (solution::f_calls>Nmax)
					return XB;
			}
		}
		else
			s=s*alpha;
		if (s < epsilon || solution::f_calls > Nmax)
			return XB;
	}
}

solution HJ_trial(solution XB, double s, matrix *ud, matrix *ad)
{
	int n=get_dim(XB);
	matrix D = ident_mat(n);
	solution X;
	for (int i = 0; i < n; ++i)
	{
		X.x = XB.x + s * D[i];
		X.fit_fun(ud, ad);
		if (X.y < XB.y)
			XB = X;
		else
		{
			X.x = XB.x - s * D[i];
			X.fit_fun(ud, ad);
			if (X.y < XB.y)
				XB = X;
		}
	}
	return XB;
}

// różne długości dla różnych kierunków;
//alpha >1
//jeśli nie poszło lepiej zmniejszam krok i zmieniam na przeciwny (s=-1*s*beta)[beta<1]
//zliczam również sumaryczne przesunięcia w danym kierunku
solution Rosen(matrix x0, matrix s0, double alpha, double beta, double epsilon, int Nmax, matrix *ud, matrix *ad)
{
	solution X(x0), Xt;
	int n = get_dim(X);//liczba wymiarów
	matrix l(n, 1),//odległoś od orginalnego punktu
		p(n, 1),//porażki w kierunu
		s(s0),//długości kroku
		D = ident_mat(n);
	X.fit_fun(ud, ad);
	while (true)
	{
		for (int i = 0; i < n; ++i)
		{
			Xt.x = X.x + s(i) * D[i];
			Xt.fit_fun(ud, ad);
			if (Xt.y < X.y)
			{
				X = Xt;
				l(i) += s(i);
				s(i) *= alpha;
			}
			else
			{
				++p(i);
				s(i) *= -beta;
			}
		}
#if LAB_NO==3 && LAB_PART==2
		???
#endif
		bool change = true;
		for (int i = 0; i < n; ++i)
			if (p(i)==0 || l(i)==0)
			{
				change = false;
				break;
			}
		if (change)
		{
			matrix Q(n,n), v(n,1);
			for (int i = 0; i<n; ++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; ++i)
			{
				matrix temp(n,1);
				for (int j = 0; j < i; ++j)
					temp = temp + trans(Q[i]) * D[j] * D[j];
				v = (Q[i] - temp) / norm(Q[i] - temp);
				D.set_col(v,i);
			}
			s = s0;
			l = matrix(n, 1);
			p = matrix(n, 1);
		}
		double max_s = abs(s(0));
		for (int i = 1; i < n; ++i)
			if (max_s < abs(s(i)))
				max_s = abs(s(i));
		if (max_s<epsilon || solution::f_calls>Nmax)
			return X;
	}
}