diff --git a/.gitattributes b/.gitattributes
index 2aa291d5261..d70c3437bbc 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -1255,6 +1255,7 @@ Polyhedron/examples/Polyhedron/corner_with_hole.off -text
 Polyhedron/examples/Polyhedron/corner_with_sharp_edge.off -text
 Polyhedron/examples/Polyhedron/cross.off -text
 Polytope_distance_d/examples/Polytope_distance_d/polytope_distance_d.C -text
+Polytope_distance_d/examples/Polytope_distance_d/polytope_distance_d_fast_exact.C -text
 Polytope_distance_d/test/Polytope_distance_d/create_test_PD_cin -text
 Polytope_distance_d/test/Polytope_distance_d/test_PD.C -text
 Polytope_distance_d/test/Polytope_distance_d/test_PD.cin -text
diff --git a/Polytope_distance_d/examples/Polytope_distance_d/polytope_distance_d.C b/Polytope_distance_d/examples/Polytope_distance_d/polytope_distance_d.C
index a9ddc049076..aaa088416b7 100644
--- a/Polytope_distance_d/examples/Polytope_distance_d/polytope_distance_d.C
+++ b/Polytope_distance_d/examples/Polytope_distance_d/polytope_distance_d.C
@@ -1,7 +1,8 @@
 // file: examples/Polytope_distance_d/polytope_distance.C
 
-// computes the distance between two cubes in R^3
-
+// computes the distance between two cubes in R^3 using double
+// as input type and internal type (prone to roundoff errors)
+#include <iostream>
 #include <CGAL/Simple_cartesian.h>
 #include <CGAL/Polytope_distance_d.h>
 #include <CGAL/Optimisation_d_traits_3.h>
@@ -25,8 +26,8 @@ int main()
   Polytope_distance pd(P, P+8, Q, Q+8); 
 
   // get squared distance (2,2,2)-(1,1,1))^2 = 3
-  std::cout << "Distance (as rational number): " <<
-    pd.squared_distance_numerator() << " / " <<
+  std::cout << "Squared distance: " <<
+    pd.squared_distance_numerator() /
     pd.squared_distance_denominator() << std::endl;
 
   // get points that realize the distance
diff --git a/Polytope_distance_d/examples/Polytope_distance_d/polytope_distance_d_fast_exact.C b/Polytope_distance_d/examples/Polytope_distance_d/polytope_distance_d_fast_exact.C
new file mode 100644
index 00000000000..71a47a587a9
--- /dev/null
+++ b/Polytope_distance_d/examples/Polytope_distance_d/polytope_distance_d_fast_exact.C
@@ -0,0 +1,57 @@
+// file: examples/Polytope_distance_d/polytope_distance_fast_exact.C
+
+// computes the distance between two cubes in R^3 using double
+// as input type and CGAL::Double as exact internal type; this
+// is guaranteed to have no roundoff errors. Note: CGAL::Double
+// is based on GMP but not yet an official CGAL number type; in 
+// this respect, the example represents experimental code
+#include <iostream>
+#include <CGAL/QP_solver/gmp_double.h> // will become CGAL number type
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/QP_solver/Double.h>     // will become CGAL number type
+#include <CGAL/Polytope_distance_d.h>
+#include <CGAL/Optimisation_d_traits_3.h>
+
+typedef CGAL::Simple_cartesian<double>    K;  
+typedef K::Point_3                        Point;
+typedef CGAL::Optimisation_d_traits_3<K, CGAL::Double, double>  
+                                          Traits;
+typedef CGAL::Polytope_distance_d<Traits> Polytope_distance;
+
+int main()
+{
+  // the cube [0,1]^3
+  Point P[8] = { Point(0,0,0), Point(0,0,1), Point(0,1,0), Point(0,1,1),
+                 Point(1,0,0), Point(1,0,1), Point(1,1,0), Point(1,1,1)};
+
+  // the cube [2,3]^3
+  Point Q[8] = { Point(2,2,2), Point(2,2,3), Point(2,3,2), Point(2,3,3),
+                 Point(3,2,2), Point(3,2,3), Point(3,3,2), Point(3,3,3)};
+
+  Polytope_distance pd(P, P+8, Q, Q+8); 
+
+  // get squared distance (2,2,2)-(1,1,1))^2 = 3
+  std::cout << "Squared Distance: " <<
+    CGAL::to_double(pd.squared_distance_numerator()) /
+    CGAL::to_double(pd.squared_distance_denominator()) << std::endl;
+
+  // get points that realize the distance
+  Polytope_distance::Coordinate_iterator  coord_it;
+
+  std::cout << "p:"; // homogeneous point from first cube, (1,1,1,1)
+  for (coord_it = pd.realizing_point_p_coordinates_begin();
+       coord_it != pd.realizing_point_p_coordinates_end();
+       ++coord_it)
+    std::cout << " " <<  CGAL::to_double(*coord_it);
+  std::cout << std::endl;
+
+  std::cout << "q:"; // homogeneous point from second cube, (2,2,2,1)
+  for (coord_it = pd.realizing_point_q_coordinates_begin();
+       coord_it != pd.realizing_point_q_coordinates_end();
+       ++coord_it)
+    std::cout << " " <<  CGAL::to_double(*coord_it);
+  std::cout << std::endl;
+
+  return 0;
+}
+