PR libstdc++/87641 correctly initialize accumulator in valarray::sum()

diff --git a/libstdc++-v3/include/bits/valarray_array.h b/libstdc++-v3/include/bits/valarray_array.h
index 6759d6003e9..2dd1ec836ac 100644
--- a/libstdc++-v3/include/bits/valarray_array.h
+++ b/libstdc++-v3/include/bits/valarray_array.h
@@ -338,33 +338,21 @@  _GLIBCXX_BEGIN_NAMESPACE_VERSION
     }
 
   //
-  // Compute the sum of elements in range [__f, __l)
+  // Compute the sum of elements in range [__f, __l) which must not be empty.
   // This is a naive algorithm.  It suffers from cancelling.
-  // In the future try to specialize
-  // for _Tp = float, double, long double using a more accurate
-  // algorithm.
+  // In the future try to specialize for _Tp = float, double, long double
+  // using a more accurate algorithm.
   //
   template<typename _Tp>
     inline _Tp
     __valarray_sum(const _Tp* __f, const _Tp* __l)
     {
-      _Tp __r = _Tp();
+      _Tp __r = *__f++;
       while (__f != __l)
 	__r += *__f++;
       return __r;
     }
 
-  // Compute the product of all elements in range [__f, __l)
-  template<typename _Tp>
-    inline _Tp
-    __valarray_product(const _Tp* __f, const _Tp* __l)
-    {
-      _Tp __r = _Tp(1);
-      while (__f != __l)
-	__r = __r * *__f++;
-      return __r;
-    }
-
   // Compute the min/max of an array-expression
   template<typename _Ta>
     inline typename _Ta::value_type
diff --git a/libstdc++-v3/src/c++98/valarray.cc b/libstdc++-v3/src/c++98/valarray.cc
index 3cec1843be2..284db21e81c 100644
--- a/libstdc++-v3/src/c++98/valarray.cc
+++ b/libstdc++-v3/src/c++98/valarray.cc
@@ -45,15 +45,21 @@  _GLIBCXX_BEGIN_NAMESPACE_VERSION
   template size_t valarray<size_t>::size() const;
   template size_t& valarray<size_t>::operator[](size_t);
 
+  // Compute the product of all elements in the non-empty range [__f, __l)
+  template<typename _Tp>
+    inline _Tp
+    __valarray_product(const _Tp* __f, const _Tp* __l)
+    {
+      _Tp __r = *__f++;
+      while (__f != __l)
+	__r = __r * *__f++;
+      return __r;
+    }
+
   inline size_t
   __valarray_product(const valarray<size_t>& __a)
   {
-    const size_t __n = __a.size();
-    // XXX: This ugly cast is necessary because
-    //      valarray::operator[]() const return a VALUE!
-    //      Try to get the committee to correct that gross error.
-    valarray<size_t>& __t = const_cast<valarray<size_t>&>(__a);
-    return __valarray_product(&__t[0], &__t[0] + __n);
+    return __valarray_product(&__a[0], &__a[0] + __a.size());
   }
 
   // Map a gslice, described by its multidimensional LENGTHS
diff --git a/libstdc++-v3/testsuite/26_numerics/valarray/87641.cc b/libstdc++-v3/testsuite/26_numerics/valarray/87641.cc
new file mode 100644
index 00000000000..eae5440e60b
--- /dev/null
+++ b/libstdc++-v3/testsuite/26_numerics/valarray/87641.cc
@@ -0,0 +1,75 @@ 
+// Copyright (C) 2018 Free Software Foundation, Inc.
+//
+// This file is part of the GNU ISO C++ Library.  This library is free
+// software; you can redistribute it and/or modify it under the
+// terms of the GNU General Public License as published by the
+// Free Software Foundation; either version 3, or (at your option)
+// any later version.
+
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+
+// You should have received a copy of the GNU General Public License along
+// with this library; see the file COPYING3.  If not see
+// <http://www.gnu.org/licenses/>.
+
+#include <valarray>
+#include <testsuite_hooks.h>
+
+void
+test01()
+{
+  // PR libstdc++/87641
+  std::valarray<int> v1(3);
+  v1[0] = 1;
+  v1[1] = 2;
+  v1[2] = 3;
+  std::valarray< std::valarray<int> > v2(v1, 3);
+  std::valarray<int> v3 = v2.sum();
+  VERIFY( v3.size() == v1.size() );
+  VERIFY( v3[0] == 3 );
+  VERIFY( v3[1] == 6 );
+  VERIFY( v3[2] == 9 );
+}
+
+struct X
+{
+  X() : val(1) { }
+
+  X& operator+=(const X& x) { val += x.val; return *this; }
+  bool operator==(const X& x) { return val == x.val; }
+
+  int val;
+};
+
+void
+test02()
+{
+  std::valarray<X> v1(1);
+  VERIFY( v1.sum() == v1[0] );
+
+  std::valarray<X> v2(2);
+  VERIFY( v2.sum().val == 2 );
+}
+
+struct Y
+{
+  X& operator+=(const Y&) { throw 1; }
+};
+
+void
+test03()
+{
+  std::valarray<Y> v1(1);
+  (void) v1.sum(); // no addition performed for a single element
+}
+
+int
+main()
+{
+  test01();
+  test02();
+  test03();
+}