Перевёл весь проект на long double для повышения точности вычислений …

…под g++.

algorithms/algorithm.hpp

@@ -67,12 +67,12 @@ class algorithm {
 	 *
 	 * Этот метод должен быть вызван хотя бы один раз до запуска алгоритма.
 	 */
-	void set_step (double step) {
+	void set_step (long double step) {
 		step_ = step;
 	}
 
 	//! Возвращает временной шаг алгоритма (т.е. шаг, с которым алгоритм будет выдавать решение).
-	double get_step() const {
+	long double get_step() const {
 		return step_;
 	}
 
@@ -81,12 +81,12 @@ class algorithm {
 	 *
 	 * Этот метод должен быть вызван хотя бы один раз до запуска алгоритма.
 	 */
-	void set_last_time (double last_time) {
+	void set_last_time (long double last_time) {
 		last_time_ = last_time;
 	}
 
 	//! Возвращает время, до которого должен работать алгоритм.
-	double get_last_time() const {
+	long double get_last_time() const {
 		return last_time_;
 	}
 
@@ -105,10 +105,10 @@ class algorithm {
 	t_input_data_ptr input_data_;
 
 	//! Временной шаг алгоритма (т.е. шаг, с которым алгоритм будет выдавать решение).
-	double step_;
+	long double step_;
 
 	//! Время, до которого должен работать алгоритм.
-	double last_time_;
+	long double last_time_;
 
 	/** Возвращает созданную и подготовленную структуру output_data.
 	 *
@@ -128,7 +128,7 @@ class algorithm {
 
 		t_output_data_ptr result (new t_output_data);
 
-		for (double t=0; t<=last_time_+EPS; t+=step_)
+		for (long double t=0; t<=last_time_+EPS; t+=step_)
 			result->ts.push_back (t);
 		result->qs.resize (result->ts.size());
 
@@ -138,18 +138,18 @@ class algorithm {
 	}
 
 	//! Возвращает мгновенные входные данные в данной точке.
-	I get_instanteous_data (double t) {
+	I get_instanteous_data (long double t) {
 		return input_data_->get_instanteous (t);
 	}
 
 	//! Возвращает интегральные входные данные на указанном отрезке.
-	I get_integrated_data (double t1, double t2) {
+	I get_integrated_data (long double t1, long double t2) {
 		return input_data_->get_integrated (t1, t2);
 	}
 
 	//! Возвращает интегральные входные данные на idx-ом по счёту отрезке (начиная с единицы).
 	I get_integrated_data (size_t idx) {
-		double t1 = step_ * (idx - 1),
+		long double t1 = step_ * (idx - 1),
 			t2 = t1 + step_;
 		return input_data_->get_integrated (t1, t2);
 	}

algorithms/auto_generated/2_step.hpp

@@ -30,7 +30,7 @@ class auto_generated_2_step_algorithm : public algorithm<Q,I> {
 			quaternion lambda;
 
 			if (i == 1) {
-				double phi_m = phi.length();
+				long double phi_m = phi.length();
 
 				lambda = quaternion (
 					cos (phi_m / 2),

algorithms/iterative_algorithm.hpp

@@ -70,7 +70,7 @@ class iterative_algorithm : public algorithm<Q,I> {
 	 * @param t Время, в которое требуется найти решение.
 	 * @param gamma Входные данные на временном отрезке (это вектор, поскольку для многошаговых алгоритмов (см. get_algorithm_steps_count_()) передаётся соответствующее число входных данных: на нескольких подотрезках, в порядке их следования во времени).
 	 */
-	virtual Q get_local_solution_ (double t, const std::vector<I> & gamma) = 0;
+	virtual Q get_local_solution_ (long double t, const std::vector<I> & gamma) = 0;
 
 
 private:
@@ -89,14 +89,14 @@ class iterative_algorithm : public algorithm<Q,I> {
 		int steps_count = this->get_algorithm_steps_count_();
 		std::vector<I> gamma (steps_count);
 
-		double step = this->step_;
-		double delta_t = step / steps_count;
+		long double step = this->step_;
+		long double delta_t = step / steps_count;
 
 		for (size_t i=1; i<result->get_count(); ++i) {
 			// вычисляем входные данные
-			double t = result->ts[i];
+			long double t = result->ts[i];
 			for (int j=0; j<steps_count; ++j) {
-				double t1 = t - step + j * delta_t,
+				long double t1 = t - step + j * delta_t,
 					t2 = t1 + delta_t;
 				gamma[j] = this->input_data_->get_integrated (t1, t2);
 			}

algorithms/iterative_riccati_algorithm.hpp

@@ -55,7 +55,7 @@ class iterative_riccati_algorithm : public iterative_algorithm<Q,I> {
 	 * @param t Время, в которое требуется найти решение.
 	 * @param gamma Входные данные на временном отрезке (это вектор, поскольку для многошаговых алгоритмов (см. get_algorithm_steps_count_()) передаётся соответствующее число входных данных: на нескольких подотрезках, в порядке их следования во времени).
 	 */
-	virtual Q get_local_riccati_solution_ (double t, const std::vector<I> & gamma) = 0;
+	virtual Q get_local_riccati_solution_ (long double t, const std::vector<I> & gamma) = 0;
 
 
 private:
@@ -66,7 +66,7 @@ class iterative_riccati_algorithm : public iterative_algorithm<Q,I> {
 	 * @param t Время, в которое требуется найти решение.
 	 * @param gamma Входные данные на временном отрезке (это вектор, поскольку для многошаговых алгоритмов (см. get_algorithm_steps_count_()) передаётся соответствующее число входных данных: на нескольких подотрезках, в порядке их следования во времени).
 	 */
-	virtual Q get_local_solution_ (double t, const std::vector<I> & gamma) {
+	virtual Q get_local_solution_ (long double t, const std::vector<I> & gamma) {
 		Q x = this->get_local_riccati_solution_ (t, gamma);
 
 		BOOST_AUTO( norm, x.norm() );

algorithms/old_algorithms/average_speed_algorithm.hpp

@@ -36,7 +36,7 @@ class average_speed_algorithm : public iterative_algorithm<Q,I> {
 	 * @param t Время, в которое требуется найти решение.
 	 * @param gamma Вектор, состоящий из единственного элемента - входных данных на текущем временном отрезке.
 	 */
-	virtual Q get_local_solution_ (double t, const std::vector<I> & gamma) {
+	virtual Q get_local_solution_ (long double t, const std::vector<I> & gamma) {
 		const I & phi = gamma[0];
 		BOOST_AUTO( phi_m, phi.length() );
 

algorithms/old_algorithms/panov_algorithm.hpp

@@ -44,7 +44,7 @@ class panov_algorithm : public iterative_algorithm<quaternion,vector3> {
 	 * @param t Время, в которое требуется найти решение.
 	 * @param gamma Вектор, состоящий из единственного элемента - входных данных на текущем временном отрезке.
 	 */
-	virtual quaternion get_local_solution_ (double t, const std::vector<vector3> & gamma) {
+	virtual quaternion get_local_solution_ (long double t, const std::vector<vector3> & gamma) {
 		matrix33 Gamma[4];
 		for (int j=0; j<4; ++j) {
 			const vector3 & g = gamma[j];
@@ -69,7 +69,7 @@ class panov_algorithm : public iterative_algorithm<quaternion,vector3> {
 		phi += delta_phi;
 
 
-		double phi_m = phi.length();
+		long double phi_m = phi.length();
 		return quaternion (
 			cos (phi_m / 2),
 			phi * sin (phi_m / 2) / phi_m

algorithms/stuff/input_data.hpp

@@ -51,7 +51,7 @@ class input_data {
 	 *
 	 * @throws std::runtime_exception Всегда кидает это исключение, если вызов дошёл до данной реализации.
 	 */
-	virtual I get_instanteous (double t) {
+	virtual I get_instanteous (long double t) {
 		throw std::runtime_error ("Not implemented: instanteous input data was not calculated.");
 	}
 
@@ -63,7 +63,7 @@ class input_data {
 	 *
 	 * @throws std::runtime_exception Всегда кидает это исключение, если вызов дошёл до данной реализации.
 	 */
-	virtual I get_integrated (double t1, double t2) {
+	virtual I get_integrated (long double t1, long double t2) {
 		throw std::runtime_error ("Not implemented: integrated input data was not calculated.");
 	}
 

algorithms/stuff/output_data.hpp

@@ -30,7 +30,7 @@ class output_data {
 	 *
 	 * Времена указаны в порядке возврастания.
 	 */
-	std::vector<double> ts;
+	std::vector<long double> ts;
 	//! Список решений - по одному объекту класса Q для каждого момента времени.
 	std::vector<Q> qs;
 
@@ -39,7 +39,7 @@ class output_data {
 	 *
 	 * Предполагается, что этот метод вызывается в порядке увеличения времён.
 	 */
-	void add (double t, const Q & q) {
+	void add (long double t, const Q & q) {
 		ts.push_back (t);
 		qs.push_back (q);
 	}
@@ -50,7 +50,7 @@ class output_data {
 	}
 
 	//! Метод для удобного доступа к списку времён.
-	double get_time (int idx) const {
+	long double get_time (int idx) const {
 		return ts[idx];
 	}
 

integrator/integrator.hpp

@@ -34,7 +34,7 @@ class integrator {
 	 * @param x0 Левая граница интегрирования (откуда начинаем интегрировать).
 	 * @param x1 Правая граница интегрирования (до которой интегрируем).
 	 */
-	virtual T integrate (boost::function < T(double) > f, double x0, double x1) = 0;
+	virtual T integrate (boost::function < T(long double) > f, long double x0, long double x1) = 0;
 
 };
 
@@ -44,7 +44,7 @@ class integrator {
 
 
 //! Погрешность метода Симпсона по умолчанию.
-static double default_integrator_step = 1E-4;
+static long double default_integrator_step = 1E-5;
 
 
 //! Возвращает текущий выбранный алгоритм интегрирования (если никакой не выбран, то возвращает метод Симпсона).

integrator/simpson_integrator.hpp

@@ -26,7 +26,7 @@ class simpson_integrator : public integrator<T> {
 	//! Конструктор, принимает в качестве параметра число N точек разбиения
 	/** \param h шаг разбиения - параметр, определяющий точность алгоритма (меньше h - выше точность)
 	*/
-	simpson_integrator (double h)
+	simpson_integrator (long double h)
 		: h(h)
 	{ }
 
@@ -38,15 +38,15 @@ class simpson_integrator : public integrator<T> {
 	}
 
 	//! Выполняет численное интегрирование и возвращает результат
-	T integrate (boost::function < T(double) > f, double x0, double x1) {
+	T integrate (boost::function < T(long double) > f, long double x0, long double x1) {
 		int N = int ((x1 - x0) / h);
 		N = std::max (N, 10);
 		if (N % 2)  ++N;
 
-		double h = (x1 - x0) / N;
+		long double h = (x1 - x0) / N;
 		T res = T();
 		for (int i=0; i<=N; ++i) {
-			double x = x0 + h * i;
+			long double x = x0 + h * i;
 			res += f(x) * ((i==0 || i==N) ? 1 : (i%2==0) ? 2 : 4);
 		}
 		return res * (x1 - x0) / N / 3.0;
@@ -55,7 +55,7 @@ class simpson_integrator : public integrator<T> {
 protected:
 
 	//! параметр метода - h - шаг интегрирования
-	double h;
+	long double h;
 
 }; // class simpson_integrator
 

math_modelling/artifical_input/artifical_input.hpp

@@ -64,9 +64,9 @@ class artifical_input {
 	}
 
 	//! Вычисляет и возвращает точное решение во все требуемые моменты времени.
-	t_output_data_ptr get_exact_solution (double step, double last_time) {
+	t_output_data_ptr get_exact_solution (long double step, long double last_time) {
 		t_output_data_ptr result (new t_output_data);
-		for (double t=0; t<=last_time+EPS; t+=step)
+		for (long double t=0; t<=last_time+EPS; t+=step)
 			result->add (t, this->internal_get_exact_solution_ (t));
 		return result;
 	}
@@ -83,7 +83,7 @@ class artifical_input {
 	 *
 	 * @throws std::runtime_exception Всегда кидает это исключение, если вызов дошёл до данной реализации.
 	 */
-	virtual I internal_get_instanteous_ (double t) {
+	virtual I internal_get_instanteous_ (long double t) {
 		throw std::runtime_error ("Not implemented: instanteous input data was not calculated.");
 	}
 
@@ -95,12 +95,12 @@ class artifical_input {
 	 *
 	 * @throws std::runtime_exception Всегда кидает это исключение, если вызов дошёл до данной реализации.
 	 */
-	virtual I internal_get_integrated_ (double t1, double t2) {
+	virtual I internal_get_integrated_ (long double t1, long double t2) {
 		throw std::runtime_error ("Not implemented: integrated input data was not calculated.");
 	}
 
 	//! Возвращает точное решение в указанный момент времени.
-	virtual Q internal_get_exact_solution_ (double t) = 0;
+	virtual Q internal_get_exact_solution_ (long double t) = 0;
 
 
 private:
@@ -119,11 +119,11 @@ class artifical_input {
 			return that->internal_get_exact_solution_ (0);
 		}
 
-		virtual I get_instanteous (double t) {
+		virtual I get_instanteous (long double t) {
 			return that->internal_get_instanteous_ (t);
 		}
 
-		virtual I get_integrated (double t1, double t2) {
+		virtual I get_integrated (long double t1, long double t2) {
 			return that->internal_get_integrated_ (t1, t2);
 		}
 

math_modelling/artifical_input/plane_angles/artifical_input_plane_angles.hpp

@@ -52,33 +52,33 @@ class artifical_input_plane_angles : public artifical_input<quaternion,vector3>
 	 * Классы-потомки должны реализовывать этот метод, задавая тем самым
 	 * конкретный вид колебаний.
 	 */
-	virtual plane_angles get_angs_ (double t) = 0;
+	virtual plane_angles get_angs_ (long double t) = 0;
 
 
 	/** Возвращает производную от get_angs_() в заданный момент времени.
 	 * 
 	 * Классы-потомки должны реализовывать этот метод, задавая тем самым
 	 * конкретный вид колебаний.
 	 */
-	virtual plane_angles get_angs_diff_ (double t) = 0;
+	virtual plane_angles get_angs_diff_ (long double t) = 0;
 
 
 private:
 
 
 	//! Возвращает мгновенные входные данные в данный момент времени.
-	virtual vector3 internal_get_instanteous_ (double t) {
+	virtual vector3 internal_get_instanteous_ (long double t) {
 		return calc_omega_ (this->get_angs_ (t), this->get_angs_diff_ (t));
 	}
 
 	//! Возвращает интегральные входные данные за указанный промежуток времени.
-	virtual vector3 internal_get_integrated_ (double t1, double t2) {
+	virtual vector3 internal_get_integrated_ (long double t1, long double t2) {
 		BOOST_AUTO( func, boost::bind (&artifical_input_plane_angles::internal_get_instanteous_, this, _1) );
 		return default_integrator<vector3>()->integrate (func, t1, t2);
 	}
 
 	//! Возвращает точное решение в указанный момент времени.
-	virtual quaternion internal_get_exact_solution_ (double t) {
+	virtual quaternion internal_get_exact_solution_ (long double t) {
 		return (quaternion) get_angs_ (t);
 	}
 

math_modelling/artifical_input/plane_angles/artifical_input_plane_angles_harmonious.hpp

@@ -66,13 +66,13 @@ class artifical_input_plane_angles_harmonious : public artifical_input_plane_ang
 
 
 	//! Возвращает ориентацию в заданный момент времени.
-	virtual plane_angles get_angs_ (double t) {
+	virtual plane_angles get_angs_ (long double t) {
 		return amp_ * sin (freq_ * t + shift_);
 	}
 
 
 	//! Возвращает производную от get_angs_() в заданный момент времени.
-	virtual plane_angles get_angs_diff_ (double t) {
+	virtual plane_angles get_angs_diff_ (long double t) {
 		return amp_ * freq_ * cos (freq_ * t + shift_);
 	}
 

math_modelling/math_modelling.hpp

@@ -80,9 +80,9 @@ class math_modelling {
 		//! Точное решение.
 		t_output_data_ptr exact_solution;
 		//! Массив погрешностей - погрешность подсчитана для каждого момента времени.
-		std::vector<double> differences;
+		std::vector<long double> differences;
 		//! Наибольшая погрешность.
-		double max_difference;
+		long double max_difference;
 
 
 	private:

types/biquaternion.hpp

@@ -75,12 +75,12 @@ class biquaternion {
 	}
 
 	//! Умножение на константу.
-	biquaternion operator* (double num) const {
+	biquaternion operator* (long double num) const {
 		return biquaternion (a*num, b*num);
 	}
 
 	//! Деление на константу.
-	biquaternion operator/ (double num) const {
+	biquaternion operator/ (long double num) const {
 		return biquaternion (a/num, b/num);
 	}
 
@@ -102,7 +102,7 @@ class biquaternion {
 
 
 //! Возвращает расстояние между бикватернионами - максимум из расстояний между частями a и b.
-inline double distance (const biquaternion & p, const biquaternion & q) {
+inline long double distance (const biquaternion & p, const biquaternion & q) {
 	return std::max (distance (p.a, q.a), distance (p.b, q.b));
 }