diff --git a/Makefile b/Makefile
index c3351cb..c8287ae 100644
--- a/Makefile
+++ b/Makefile
@@ -1,4 +1,4 @@
-PGSPHERE_VERSION = 1.3.0
+PGSPHERE_VERSION = 1.3.1
 EXTENSION        = pg_sphere
 RELEASE_SQL      = $(EXTENSION)--$(PGSPHERE_VERSION).sql
 USE_PGXS         = 1
@@ -28,7 +28,8 @@ DATA_built  = $(RELEASE_SQL) \
 			  pg_sphere--1.2.0--1.2.1.sql \
 			  pg_sphere--1.2.1--1.2.2.sql \
 			  pg_sphere--1.2.2--1.2.3.sql \
-			  pg_sphere--1.2.3--1.3.0.sql
+			  pg_sphere--1.2.3--1.3.0.sql \
+			  pg_sphere--1.3.0--1.3.1.sql
 
 DOCS        = README.pg_sphere COPYRIGHT.pg_sphere
 REGRESS     = init tables points euler circle line ellipse poly path box index \
@@ -265,6 +266,9 @@ pg_sphere--1.2.2--1.2.3.sql:
 pg_sphere--1.2.3--1.3.0.sql: pgs_brin.sql.in
 	cat upgrade_scripts/$@.in $^ > $@
 
+pg_sphere--1.3.0--1.3.1.sql:
+	cat upgrade_scripts/$@.in > $@
+
 # end of local stuff
 
 src/sscan.o : src/sparse.c
diff --git a/expected/init.out b/expected/init.out
index fc47f40..fb86488 100644
--- a/expected/init.out
+++ b/expected/init.out
@@ -6,6 +6,6 @@ CREATE EXTENSION pg_sphere;
 select pg_sphere_version();
  pg_sphere_version 
 -------------------
- 1.3.0
+ 1.3.1
 (1 row)
 
diff --git a/pg_sphere--1.3.0.sql.in b/pg_sphere--1.3.1.sql.in
similarity index 100%
rename from pg_sphere--1.3.0.sql.in
rename to pg_sphere--1.3.1.sql.in
diff --git a/pg_sphere.control b/pg_sphere.control
index 412cd57..e2b8ad3 100644
--- a/pg_sphere.control
+++ b/pg_sphere.control
@@ -1,5 +1,5 @@
 # pg_sphere extension
 comment = 'spherical objects with useful functions, operators and index support'
-default_version = '1.3.0'
+default_version = '1.3.1'
 module_pathname = '$libdir/pg_sphere'
 relocatable = true
diff --git a/src/pg_sphere.h b/src/pg_sphere.h
index 7c6e24c..24d950a 100644
--- a/src/pg_sphere.h
+++ b/src/pg_sphere.h
@@ -45,10 +45,85 @@
 
 #include "pgs_util.h"
 
+/* On some 32 bit platforms, there is a gcc bug that makes floating point
+ * calculations and comparisons unstable (see the link below). The problem
+ * originates in FPU 80 bits registers where double values are not truncated
+ * to 64 bit values. When gcc compiles some code with enabled optimizations,
+ * the intermediate results may be kept in the FPU registers without truncation
+ * to 64 bit values. Extra bits may produce unstable results when comparing
+ * the numbers.
+ *
+ * The generic solution is to save the intermediate results in the memory where
+ * the values are truncated to 64 bit values. It affects the performance but
+ * makes the tests stable on all platforms.
+ *
+ * PGSPHERE_FLOAT_STORE macro enables storing of intermediate results for FPxx
+ * operations in the memory. It is enabled by default for 32 bit platforms.
+ * It can be explicitly enabled or disabled in CFLAGS. To enable it for all
+ * code the gcc option -ffloat-store may be used as well.
+ *
+ * Link to gcc bug: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=323
+ */
+#if !defined(PGSPHERE_FLOAT_STORE)
+#if _WIN64 || (__GNUC__ && __x86_64__)
+#define PGSPHERE_FLOAT_STORE 0
+#elif _WIN32 || __GNUC__
+#define PGSPHERE_FLOAT_STORE 1
+#else
+#define PGSPHERE_FLOAT_STORE 0
+#endif
+#endif // PGSPHERE_FLOAT_STORE
+
 #define EPSILON					1.0E-09
 
 #define FPzero(A)				(fabs(A) <= EPSILON)
 
+#if PGSPHERE_FLOAT_STORE
+
+static inline bool
+FPeq(double A, double B)
+{
+	const volatile double AB = A - B;
+	return A == B || fabs(AB) <= EPSILON;
+}
+
+static inline bool
+FPne(double A, double B)
+{
+	const volatile double AB = A - B;
+	return A != B && fabs(AB) > EPSILON;
+}
+
+static inline bool
+FPlt(double A, double B)
+{
+	const volatile double AE = A + EPSILON;
+	return AE < B;
+}
+
+static inline bool
+FPle(double A, double B)
+{
+	const volatile double BE = B + EPSILON;
+	return A <= BE;
+}
+
+static inline bool
+FPgt(double A, double B)
+{
+	const volatile double BE = B + EPSILON;
+	return A > BE;
+}
+
+static inline bool
+FPge(double A, double B)
+{
+	const volatile double AE = A + EPSILON;
+	return AE >= B;
+}
+
+#else
+
 static inline bool
 FPeq(double A, double B)
 {
@@ -85,6 +160,8 @@ FPge(double A, double B)
 	return A + EPSILON >= B;
 }
 
+#endif // PGSPHERE_FLOAT_STORE
+
 /*---------------------------------------------------------------------
  * Point - (x,y)
  *-------------------------------------------------------------------*/
diff --git a/upgrade_scripts/pg_sphere--1.3.0--1.3.1.sql.in b/upgrade_scripts/pg_sphere--1.3.0--1.3.1.sql.in
new file mode 100644
index 0000000..8ea3baa
--- /dev/null
+++ b/upgrade_scripts/pg_sphere--1.3.0--1.3.1.sql.in
@@ -0,0 +1 @@
+-- Nothing to upgrade in the schema