Repository: xtensor-stack/xtensor (language: python)

Issue:
Precision of the endpoint of linspace
Currently `linspace` generates the array from the starting point step by step. However, it could result in a subtle precision issue with the endpoint:

```
        double at_end = 99.78730976641236;
        auto ls = linspace<double>(0., at_end, 100);
        std::cout << std::setprecision(16) << ls(99) << std::endl;
        // 99.78730976641235
```

This affects some algorithms like `histogram` since the right boundary of the bin edges will be smaller than the largest number in the array. Consequently, the largest number will be excluded. In `numpy`, I do not experience the same issue.

The solution would be to set the last point to the right boundary once the sequence was generated. Could anyone give some hint about how to do that? `linspace` seems to be a generator.


Discussion / hints:
The computation of the value, given the index, is done [here](https://github.com/xtensor-stack/xtensor/blob/master/include/xtensor/xbuilder.hpp#L261). I guess we could pass the number of elements to the constructor of `arange_generator` (or compute it, but it's already done [here](https://github.com/xtensor-stack/xtensor/blob/master/include/xtensor/xbuilder.hpp#L282)) and check in the access method if we ask for the last point?
I did some investigation. Actually, `numpy.arange` has the [accumulated rounding error issue](https://github.com/numpy/numpy/pull/16796) as well. But `numpy.linspace` fixes it
by [explicitly setting the last value](https://github.com/numpy/numpy/blob/fb215c76967739268de71aa4bda55dd1b062bc2e/numpy/core/function_base.py#L164). 

The implementation of `numpy.linspace` is a bit weird which makes use of `numpy.arange`. Currently, the `xtensor` implementations both have
```
detail::make_xgenerator(detail::arange_generator ...
```

Fix patch:
diff --git a/include/xtensor/xbuilder.hpp b/include/xtensor/xbuilder.hpp
index 9cc865f13..5658558b9 100644
--- a/include/xtensor/xbuilder.hpp
+++ b/include/xtensor/xbuilder.hpp
@@ -191,7 +191,7 @@ namespace xt
         // These methods should be private methods of arange_generator, however thi leads
         // to ICE on VS2015
         template <class R, class E, class U, class X, XTL_REQUIRES(xtl::is_integral<X>)>
-        inline void arange_assign_to(xexpression<E>& e, U start, X step) noexcept
+        inline void arange_assign_to(xexpression<E>& e, U start, U, X step, bool) noexcept
         {
             auto& de = e.derived_cast();
             U value = start;
@@ -204,15 +204,20 @@ namespace xt
         }
 
         template <class R, class E, class U, class X, XTL_REQUIRES(xtl::negation<xtl::is_integral<X>>)>
-        inline void arange_assign_to(xexpression<E>& e, U start, X step) noexcept
+        inline void arange_assign_to(xexpression<E>& e, U start, U stop, X step, bool endpoint) noexcept
         {
             auto& buf = e.derived_cast().storage();
             using size_type = decltype(buf.size());
             using mult_type = get_mult_type<U, X>;
-            for(size_type i = 0; i < buf.size(); ++i)
+            size_type num = buf.size();
+            for(size_type i = 0; i < num; ++i)
             {
                 buf[i] = static_cast<R>(start + step * mult_type(i));
             }
+            if (num > 0 && endpoint)
+            {
+                buf[num-1] = static_cast<R>(stop);
+            }
         }
 
         template <class T, class R = T, class S = T>
@@ -223,8 +228,8 @@ namespace xt
             using value_type = R;
             using step_type = S;
 
-            arange_generator(T start, T stop, S step)
-                : m_start(start), m_stop(stop), m_step(step)
+            arange_generator(T start, T stop, S step, size_t num_steps, bool endpoint=false)
+                : m_start(start), m_stop(stop), m_step(step), m_num_steps(num_steps), m_endpoint(endpoint)
             {
             }
 
@@ -237,15 +242,13 @@ namespace xt
             template <class It>
             inline R element(It first, It) const
             {
-                // Avoids warning when T = char (because char + char => int!)
-                using mult_type = get_mult_type<T, S>;
-                return static_cast<R>(m_start + m_step * mult_type(*first));
+                return access_impl(*first);
             }
 
             template <class E>
             inline void assign_to(xexpression<E>& e) const noexcept
             {
-                arange_assign_to<R>(e, m_start, m_step);
+                arange_assign_to<R>(e, m_start, m_stop, m_step, m_endpoint);
             }
 
         private:
@@ -253,10 +256,17 @@ namespace xt
             T m_start;
             T m_stop;
             step_type m_step;
+            size_t m_num_steps;
+            bool m_endpoint; // true for setting the last element to m_stop
 
             template <class T1, class... Args>
             inline R access_impl(T1 t, Args...) const
             {
+                if (m_endpoint && m_num_steps > 1 && size_t(t) == m_num_steps - 1)
+                {
+                    return static_cast<R>(m_stop);
+                }
+                // Avoids warning when T = char (because char + char => int!)
                 using mult_type = get_mult_type<T, S>;
                 return static_cast<R>(m_start + m_step * mult_type(t));
             }
@@ -280,7 +290,7 @@ namespace xt
         inline auto arange_impl(T start, T stop, S step = 1) noexcept
         {
             std::size_t shape = static_cast<std::size_t>(std::ceil((stop - start) / step));
-            return detail::make_xgenerator(detail::arange_generator<T, T, S>(start, stop, step), {shape});
+            return detail::make_xgenerator(detail::arange_generator<T, T, S>(start, stop, step, shape), {shape});
         }
 
         template <class T, class S = T, XTL_REQUIRES(integer_with_signed_integer<T, S>)>
@@ -293,7 +303,7 @@ namespace xt
                 shape = stop > start ? static_cast<std::size_t>((stop - start + step - S(1)) / step)
                                      : static_cast<std::size_t>((start - stop - step - S(1)) / -step);
             }
-            return detail::make_xgenerator(detail::arange_generator<T, T, S>(start, stop, step), {shape});
+            return detail::make_xgenerator(detail::arange_generator<T, T, S>(start, stop, step, shape), {shape});
         }
 
         template <class T, class S = T, XTL_REQUIRES(integer_with_unsigned_integer<T, S>)>
@@ -305,7 +315,7 @@ namespace xt
             {
                 shape = static_cast<std::size_t>((stop - start + step - S(1)) / step);
             }
-            return detail::make_xgenerator(detail::arange_generator<T, T, S>(start, stop, step), { shape });
+            return detail::make_xgenerator(detail::arange_generator<T, T, S>(start, stop, step, shape), { shape });
         }
 
         template <class F>
@@ -446,7 +456,7 @@ namespace xt
     {
         using fp_type = std::common_type_t<T, double>;
         fp_type step = fp_type(stop - start) / std::fmax(fp_type(1), fp_type(num_samples - (endpoint ? 1 : 0)));
-        return detail::make_xgenerator(detail::arange_generator<fp_type, T>(fp_type(start), fp_type(stop), step), {num_samples});
+        return detail::make_xgenerator(detail::arange_generator<fp_type, T>(fp_type(start), fp_type(stop), step, num_samples, endpoint), {num_samples});
     }
 
     /**


Test patch:
diff --git a/test/test_xbuilder.cpp b/test/test_xbuilder.cpp
index 042e89793..d9e415eb1 100644
--- a/test/test_xbuilder.cpp
+++ b/test/test_xbuilder.cpp
@@ -271,6 +271,13 @@ namespace xt
         ASSERT_EQ(m_assigned(3), at_3);
     }
 
+    TEST(xbuilder, linspace_endpoint)
+    {
+        double at_end = 99.78730976641236;
+        xt::xtensor<double, 1> ls = linspace<double>(0., at_end, 100);
+        ASSERT_EQ(ls(99), at_end);
+    }
+
     TEST(xbuilder, linspace_integer)
     {
         xarray<int> ls = linspace<int>(0, 10, 13);
@@ -686,8 +693,10 @@ namespace xt
 
     TEST(xbuilder, linspace_double)
     {
-        xt::xarray<double> a = xt::linspace(0., 100.);
-        auto b = xt::linspace(0., 100.);
+        // access_imp will also be tested here
+        double at_end = 99.78730976641236;
+        xt::xarray<double> a = xt::linspace(0., at_end, 100);
+        auto b = xt::linspace(0., at_end, 100);
         EXPECT_EQ(a, b);
     }
 }