diff --git a/algorithms/algorithm.hpp b/algorithms/algorithm.hpp
index 4a5d6fa..3f5615c 100644
--- a/algorithms/algorithm.hpp
+++ b/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);
 	}
diff --git a/algorithms/auto_generated/2_step.hpp b/algorithms/auto_generated/2_step.hpp
index 72f2de4..f901d41 100644
--- a/algorithms/auto_generated/2_step.hpp
+++ b/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),
diff --git a/algorithms/iterative_algorithm.hpp b/algorithms/iterative_algorithm.hpp
index c27e44b..1fa04cf 100644
--- a/algorithms/iterative_algorithm.hpp
+++ b/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);
 			}
diff --git a/algorithms/iterative_riccati_algorithm.hpp b/algorithms/iterative_riccati_algorithm.hpp
index e0d77b3..2cc4e05 100644
--- a/algorithms/iterative_riccati_algorithm.hpp
+++ b/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() );
diff --git a/algorithms/old_algorithms/average_speed_algorithm.hpp b/algorithms/old_algorithms/average_speed_algorithm.hpp
index 35a028b..ab56a40 100644
--- a/algorithms/old_algorithms/average_speed_algorithm.hpp
+++ b/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() );
 
diff --git a/algorithms/old_algorithms/panov_algorithm.hpp b/algorithms/old_algorithms/panov_algorithm.hpp
index 645912b..dd37415 100644
--- a/algorithms/old_algorithms/panov_algorithm.hpp
+++ b/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
diff --git a/algorithms/stuff/input_data.hpp b/algorithms/stuff/input_data.hpp
index 82dbd18..ba39ac7 100644
--- a/algorithms/stuff/input_data.hpp
+++ b/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.");
 	}
 
diff --git a/algorithms/stuff/output_data.hpp b/algorithms/stuff/output_data.hpp
index 97b5943..4f96d56 100644
--- a/algorithms/stuff/output_data.hpp
+++ b/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];
 	}
 
diff --git a/integrator/integrator.hpp b/integrator/integrator.hpp
index 5ed7b78..6a3794e 100644
--- a/integrator/integrator.hpp
+++ b/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;
 
 
 //! Возвращает текущий выбранный алгоритм интегрирования (если никакой не выбран, то возвращает метод Симпсона).
diff --git a/integrator/simpson_integrator.hpp b/integrator/simpson_integrator.hpp
index 92ce0e8..5499f0e 100644
--- a/integrator/simpson_integrator.hpp
+++ b/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
 
diff --git a/math_modelling/artifical_input/artifical_input.hpp b/math_modelling/artifical_input/artifical_input.hpp
index 09c11c2..de57aa5 100644
--- a/math_modelling/artifical_input/artifical_input.hpp
+++ b/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);
 		}
 
diff --git a/math_modelling/artifical_input/plane_angles/artifical_input_plane_angles.hpp b/math_modelling/artifical_input/plane_angles/artifical_input_plane_angles.hpp
index 8f80c9d..aaf5bae 100644
--- a/math_modelling/artifical_input/plane_angles/artifical_input_plane_angles.hpp
+++ b/math_modelling/artifical_input/plane_angles/artifical_input_plane_angles.hpp
@@ -52,7 +52,7 @@ 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_() в заданный момент времени.
@@ -60,25 +60,25 @@ class artifical_input_plane_angles : public artifical_input<quaternion,vector3>
 	 * Классы-потомки должны реализовывать этот метод, задавая тем самым
 	 * конкретный вид колебаний.
 	 */
-	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);
 	}
 
diff --git a/math_modelling/artifical_input/plane_angles/artifical_input_plane_angles_harmonious.hpp b/math_modelling/artifical_input/plane_angles/artifical_input_plane_angles_harmonious.hpp
index c1f8182..66f95ef 100644
--- a/math_modelling/artifical_input/plane_angles/artifical_input_plane_angles_harmonious.hpp
+++ b/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_);
 	}
 
diff --git a/math_modelling/math_modelling.hpp b/math_modelling/math_modelling.hpp
index ac35f07..fafbed1 100644
--- a/math_modelling/math_modelling.hpp
+++ b/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:
diff --git a/types/biquaternion.hpp b/types/biquaternion.hpp
index d2773d4..de86b23 100644
--- a/types/biquaternion.hpp
+++ b/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));
 }
 
diff --git a/types/dual_number.hpp b/types/dual_number.hpp
index 74320a2..7be77f5 100644
--- a/types/dual_number.hpp
+++ b/types/dual_number.hpp
@@ -16,20 +16,20 @@ class dual_number {
 public:
 
 	//! Действительная компонента числа.
-	double real;
+	long double real;
 	//! Мнимая компонента числа (домножающаяся на мнимую единицу s).
-	double imag;
+	long double imag;
 
 
 	dual_number()
 		: real(0), imag(0)
 	{ }
 
-	explicit dual_number (double real)
+	explicit dual_number (long double real)
 		: real(real), imag(0)
 	{ }
 
-	dual_number (double real, double imag)
+	dual_number (long double real, long double imag)
 		: real(real), imag(imag)
 	{ }
 
@@ -53,11 +53,11 @@ class dual_number {
 		);
 	}
 
-	dual_number operator* (double num) const {
+	dual_number operator* (long double num) const {
 		return dual_number (real * num, imag * num);
 	}
 
-	dual_number operator/ (double num) const {
+	dual_number operator/ (long double num) const {
 		return dual_number (real / num, imag / num);
 	}
 
@@ -71,12 +71,12 @@ class dual_number {
 		return *this;
 	}
 
-	dual_number & operator*= (double num) {
+	dual_number & operator*= (long double num) {
 		*this = *this * num;
 		return *this;
 	}
 
-	dual_number & operator/= (double num) {
+	dual_number & operator/= (long double num) {
 		*this = *this / num;
 		return *this;
 	}
@@ -85,7 +85,7 @@ class dual_number {
 
 
 
-inline dual_number operator* (double num, const dual_number & dual_num) {
+inline dual_number operator* (long double num, const dual_number & dual_num) {
 	return dual_num * num;
 }
 
@@ -118,7 +118,7 @@ inline dual_number log (const dual_number & num) {
 }
 
 inline dual_number exp (const dual_number & num) {
-	double ex = exp (num.real);
+	long double ex = exp (num.real);
 	return dual_number (
 		ex,
 		ex * num.imag
@@ -126,7 +126,7 @@ inline dual_number exp (const dual_number & num) {
 }
 
 inline dual_number sqrt (const dual_number & num) {
-	double sq = sqrt (num.real);
+	long double sq = sqrt (num.real);
 	return dual_number (
 		sq,
 		num.imag / 2 / sq
diff --git a/types/dual_vector.hpp b/types/dual_vector.hpp
index 2a990a1..d820bf5 100644
--- a/types/dual_vector.hpp
+++ b/types/dual_vector.hpp
@@ -72,11 +72,11 @@ class dual_vector {
 		return dual_vector (real - other.real, imag - other.imag);
 	}
 
-	dual_vector operator* (double num) const {
+	dual_vector operator* (long double num) const {
 		return dual_vector (real * num, imag * num);
 	}
 
-	dual_vector operator/ (double num) const {
+	dual_vector operator/ (long double num) const {
 		return dual_vector (real / num, imag / num);
 	}
 
@@ -112,7 +112,7 @@ class dual_vector {
 
 
 
-inline dual_vector operator* (double num, const dual_vector & dual_vec) {
+inline dual_vector operator* (long double num, const dual_vector & dual_vec) {
 	return dual_vec * num;
 }
 
diff --git a/types/matrix33.hpp b/types/matrix33.hpp
index c9eed8e..ba87243 100644
--- a/types/matrix33.hpp
+++ b/types/matrix33.hpp
@@ -26,13 +26,13 @@ class matrix33 {
 		memset (data_, 0, sizeof data_);
 	}
 
-	matrix33 (double data[3][3]) {
+	matrix33 (long double data[3][3]) {
 		for (int i=0; i<3; ++i)
 			for (int j=0; j<3; ++j)
 				data_[i][j] = data[i][j];
 	}
 
-	matrix33 (double a11, double a12, double a13, double a21, double a22, double a23, double a31, double a32, double a33) {
+	matrix33 (long double a11, long double a12, long double a13, long double a21, long double a22, long double a23, long double a31, long double a32, long double a33) {
 		data_[0][0] = a11;
 		data_[0][1] = a12;
 		data_[0][2] = a13;
@@ -49,7 +49,7 @@ class matrix33 {
 	 *
 	 * @throws std::invalid_argument В случае, если передано некорректное значение row или column.
 	 */
-	double operator() (int row, int column) const {
+	long double operator() (int row, int column) const {
 		if (row < 0 || row >= 3)
 			throw std::invalid_argument ("Invalid row value.");
 		if (column < 0 || column >= 3)
@@ -61,7 +61,7 @@ class matrix33 {
 	 *
 	 * @throws std::invalid_argument В случае, если передано некорректное значение row или column.
 	 */
-	double & operator() (int row, int column) {
+	long double & operator() (int row, int column) {
 		if (row < 0 || row >= 3)
 			throw std::invalid_argument ("Invalid row value.");
 		if (column < 0 || column >= 3)
@@ -81,7 +81,7 @@ class matrix33 {
 	}
 
 	//! Умножение на константу.
-	matrix33 operator* (double num) const {
+	matrix33 operator* (long double num) const {
 		matrix33 result = *this;
 		return result *= num;
 	}
@@ -108,7 +108,7 @@ class matrix33 {
 	}
 
 	//! Деление на константу.
-	matrix33 operator/ (double num) const {
+	matrix33 operator/ (long double num) const {
 		matrix33 result = *this;
 		return result /= num;
 	}
@@ -134,7 +134,7 @@ class matrix33 {
 	}
 
 	//! Умножение на константу.
-	matrix33 & operator*= (double num) {
+	matrix33 & operator*= (long double num) {
 		for (int i=0; i<3; ++i)
 			for (int j=0; j<3; ++j)
 				data_[i][j] *= num;
@@ -148,7 +148,7 @@ class matrix33 {
 	}
 	
 	//! Деление на константу.
-	matrix33 & operator/= (double num) {
+	matrix33 & operator/= (long double num) {
 		for (int i=0; i<3; ++i)
 			for (int j=0; j<3; ++j)
 				data_[i][j] /= num;
@@ -160,7 +160,7 @@ class matrix33 {
 
 
 	//! Компоненты матрицы.
-	double data_[3][3];
+	long double data_[3][3];
 
 
 }; // class matrix33
@@ -168,7 +168,7 @@ class matrix33 {
 
 
 //! Умножение на константу
-inline matrix33 operator* (double num, const matrix33 & m) {
+inline matrix33 operator* (long double num, const matrix33 & m) {
 	return m * num;
 }
 
diff --git a/types/plane_angles.hpp b/types/plane_angles.hpp
index 607200b..3d4ccdb 100644
--- a/types/plane_angles.hpp
+++ b/types/plane_angles.hpp
@@ -26,7 +26,7 @@ class plane_angles {
 	/** @name Тройка углов.
 	 */
 	//@{
-	double psi, //!< угол psi - географический курс
+	long double psi, //!< угол psi - географический курс
 		teta, //!< угол teta - угол тангажа
 		gamma; //!< угол gamma - угол крена
 	//@}
@@ -37,7 +37,7 @@ class plane_angles {
 	{ }
 
 	//! Конструктор от тройки углов.
-	plane_angles (double psi, double teta, double gamma)
+	plane_angles (long double psi, long double teta, long double gamma)
 		: psi(psi), teta(teta), gamma(gamma)
 	{ }
 
@@ -75,7 +75,7 @@ class plane_angles {
 		);
 	}
 
-	plane_angles operator* (double p) const {
+	plane_angles operator* (long double p) const {
 		return plane_angles (
 			psi   * p,
 			teta  * p,
diff --git a/types/quaternion.hpp b/types/quaternion.hpp
index 5e40e35..9ad434c 100644
--- a/types/quaternion.hpp
+++ b/types/quaternion.hpp
@@ -27,7 +27,7 @@ class quaternion {
 	
 	/** Компоненты кватерниона. */
 	//@{
-	double w, x, y, z;
+	long double w, x, y, z;
 	//@}
 
 	
@@ -36,17 +36,17 @@ class quaternion {
 		: w(0), x(0), y(0), z(0)
 	{ }
 
-	quaternion (double w, double x, double y, double z)
+	quaternion (long double w, long double x, long double y, long double z)
 		: w(w), x(x), y(y), z(z)
 	{ }
 
 	//! Конструктор кватерниона от скалярной и векторной частей.
-	quaternion (double w, const vector3 & v)
+	quaternion (long double w, const vector3 & v)
 		: w(w), x(v.x), y(v.y), z(v.z)
 	{ }
 
 	//! Конструктор кватерниона с нулевой векторной частью.
-	explicit quaternion (double num)
+	explicit quaternion (long double num)
 		: w(num), x(0), y(0), z(0)
 	{ }
 
@@ -60,7 +60,7 @@ class quaternion {
 	 *
 	 * @throws std::invalid_argument В случае, если передано некорректное значение idx.
 	 */
-	double operator[] (int idx) const {
+	long double operator[] (int idx) const {
 		if (idx == 0)  return w;
 		if (idx == 1)  return x;
 		if (idx == 2)  return y;
@@ -72,7 +72,7 @@ class quaternion {
 	 *
 	 * @throws std::invalid_argument В случае, если передано некорректное значение idx.
 	 */
-	double & operator[] (int idx) {
+	long double & operator[] (int idx) {
 		if (idx == 0)  return w;
 		if (idx == 1)  return x;
 		if (idx == 2)  return y;
@@ -82,7 +82,7 @@ class quaternion {
 
 
 	//! Скалярная часть кватерниона - член w.
-	double get_scalar() const {
+	long double get_scalar() const {
 		return w;
 	}
 
@@ -92,7 +92,7 @@ class quaternion {
 	}
 
 	
-	quaternion operator+ (double num) const {
+	quaternion operator+ (long double num) const {
 		return quaternion (w+num, x, y, z);
 	}
 
@@ -100,7 +100,7 @@ class quaternion {
 		return quaternion (w+q.w, x+q.x, y+q.y, z+q.z);
 	}
 
-	quaternion operator- (double num) const {
+	quaternion operator- (long double num) const {
 		return quaternion (w-num, x, y, z);
 	}
 
@@ -109,12 +109,12 @@ class quaternion {
 	}
 
 	//! Умножение на константу.
-	quaternion operator* (double num) const {
+	quaternion operator* (long double num) const {
 		return quaternion (w*num, x*num, y*num, z*num);
 	}
 
 	//! Деление на константу.
-	quaternion operator/ (double num) const {
+	quaternion operator/ (long double num) const {
 		return quaternion (w/num, x/num, y/num, z/num);
 	}
 
@@ -143,12 +143,12 @@ class quaternion {
 		return *this;
 	}
 
-	quaternion & operator*= (double num) {
+	quaternion & operator*= (long double num) {
 		*this = *this * num;
 		return *this;
 	}
 
-	quaternion & operator/= (double num) {
+	quaternion & operator/= (long double num) {
 		*this = *this / num;
 		return *this;
 	}
@@ -171,12 +171,12 @@ class quaternion {
 
 
 	//! Возвращает норму кватерниона - сумму квадратов компонент.
-	double norm() const {
+	long double norm() const {
 		return w*w + x*x + y*y + z*z;
 	}
 
 	//! Возвращает длину (тензор) кватерниона - квадратный корень из суммы квадратов компонент.
-	double length() const {
+	long double length() const {
 		return sqrt (norm());
 	}
 
@@ -191,22 +191,22 @@ class quaternion {
 
 
 
-inline quaternion operator+ (double a, const quaternion & b) {
+inline quaternion operator+ (long double a, const quaternion & b) {
 	return b + a;
 }
 
-inline quaternion operator- (double a, const quaternion & b) {
+inline quaternion operator- (long double a, const quaternion & b) {
 	return quaternion (a-b.w, -b.x, -b.y, -b.z);
 }
 
-inline quaternion operator* (double a, const quaternion & b) {
+inline quaternion operator* (long double a, const quaternion & b) {
 	return b * a;
 }
 
 
 
 //! Возвращает расстояние между кватернионами - т.е. модуль их разности.
-inline double distance (const quaternion & a, const quaternion & b) {
+inline long double distance (const quaternion & a, const quaternion & b) {
 	return (a-b).length();
 }
 
diff --git a/types/vector3.hpp b/types/vector3.hpp
index c5d9073..837ee4a 100644
--- a/types/vector3.hpp
+++ b/types/vector3.hpp
@@ -21,7 +21,7 @@ class vector3 {
 	
 	/** Компоненты вектора. */
 	//@{
-	double x, y, z;
+	long double x, y, z;
 	//@}
 
 	
@@ -30,7 +30,7 @@ class vector3 {
 		: x(0), y(0), z(0)
 	{ }
 
-	vector3 (double x, double y, double z)
+	vector3 (long double x, long double y, long double z)
 		: x(x), y(y), z(z)
 	{ }
 
@@ -43,7 +43,7 @@ class vector3 {
 	 *
 	 * @throws std::invalid_argument В случае, если передано некорректное значение idx.
 	 */
-	double operator[] (int idx) const {
+	long double operator[] (int idx) const {
 		if (idx == 0)  return x;
 		if (idx == 1)  return y;
 		if (idx == 2)  return z;
@@ -54,7 +54,7 @@ class vector3 {
 	 *
 	 * @throws std::invalid_argument В случае, если передано некорректное значение idx.
 	 */
-	double & operator[] (int idx) {
+	long double & operator[] (int idx) {
 		if (idx == 0)  return x;
 		if (idx == 1)  return y;
 		if (idx == 2)  return z;
@@ -71,12 +71,12 @@ class vector3 {
 	}
 
 	//! Умножение на константу.
-	vector3 operator* (double num) const {
+	vector3 operator* (long double num) const {
 		return vector3 (x*num, y*num, z*num);
 	}
 
 	//! Деление на константу.
-	vector3 operator/ (double num) const {
+	vector3 operator/ (long double num) const {
 		return vector3 (x/num, y/num, z/num);
 	}
 
@@ -101,7 +101,7 @@ class vector3 {
 
 
 	//! Скалярное произведение.
-	double dotProduct (const vector3 & v) const {
+	long double dotProduct (const vector3 & v) const {
 		return x * v.x + y * v.y + z * v.z;
 	}
 
@@ -116,12 +116,12 @@ class vector3 {
 
 
 	//! Возвращает норму вектора - сумму квадратов компонент.
-	double norm() const {
+	long double norm() const {
 		return x*x + y*y + z*z;
 	}
 
 	//! Возвращает длину (тензор) вектора - квадратный корень из суммы квадратов компонент.
-	double length() const {
+	long double length() const {
 		return sqrt (norm());
 	}
 
@@ -131,17 +131,17 @@ class vector3 {
 
 
 //! Умножение на константу
-inline vector3 operator* (double num, const vector3 & v) {
+inline vector3 operator* (long double num, const vector3 & v) {
 	return v * num;
 }
 
 //! Возвращает расстояние между векторами - т.е. модуль их разности.
-inline double distance (const vector3 & a, const vector3 & b) {
+inline long double distance (const vector3 & a, const vector3 & b) {
 	return (a-b).length();
 }
 
 //! Скалярное произведение.
-inline double dotProduct (const vector3 & a, const vector3 & b) {
+inline long double dotProduct (const vector3 & a, const vector3 & b) {
 	return a.dotProduct (b);
 }
 
diff --git a/unit_tests/utility_functions.cpp b/unit_tests/utility_functions.cpp
index 349279e..77373cc 100644
--- a/unit_tests/utility_functions.cpp
+++ b/unit_tests/utility_functions.cpp
@@ -12,11 +12,11 @@ BOOST_AUTO_TEST_CASE( sqr_test )
 {
 	const double tolerance = 1E-7;
 
-    BOOST_CHECK_CLOSE( sqr(0), 0, tolerance );
-    BOOST_CHECK_CLOSE( sqr(2), 4, tolerance );
-    BOOST_CHECK_CLOSE( sqr(-2), 4, tolerance );
-    BOOST_CHECK_CLOSE( sqr(9), 81, tolerance );
-    BOOST_CHECK_CLOSE( sqr(-9), 81, tolerance );
+    BOOST_CHECK_CLOSE( sqr((long double)0), 0, tolerance );
+    BOOST_CHECK_CLOSE( sqr((long double)2), 4, tolerance );
+    BOOST_CHECK_CLOSE( sqr((long double)-2), 4, tolerance );
+    BOOST_CHECK_CLOSE( sqr((long double)9), 81, tolerance );
+    BOOST_CHECK_CLOSE( sqr((long double)-9), 81, tolerance );
 }
 
 
diff --git a/utility_functions.hpp b/utility_functions.hpp
index cd9e2a2..086d7a0 100644
--- a/utility_functions.hpp
+++ b/utility_functions.hpp
@@ -8,7 +8,8 @@
 
 
 //! Вычисляет квадрат числа.
-inline double sqr (double x) {
+template <class T>
+inline T sqr (const T & x) {
 	return x * x;
 }