1    	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2    	/*                                                                           */
3    	/*                  This file is part of the class library                   */
4    	/*       SoPlex --- the Sequential object-oriented simPlex.                  */
5    	/*                                                                           */
6    	/*  Copyright (c) 1996-2023 Zuse Institute Berlin (ZIB)                      */
7    	/*                                                                           */
8    	/*  Licensed under the Apache License, Version 2.0 (the "License");          */
9    	/*  you may not use this file except in compliance with the License.         */
10   	/*  You may obtain a copy of the License at                                  */
11   	/*                                                                           */
12   	/*      http://www.apache.org/licenses/LICENSE-2.0                           */
13   	/*                                                                           */
14   	/*  Unless required by applicable law or agreed to in writing, software      */
15   	/*  distributed under the License is distributed on an "AS IS" BASIS,        */
16   	/*  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17   	/*  See the License for the specific language governing permissions and      */
18   	/*  limitations under the License.                                           */
19   	/*                                                                           */
20   	/*  You should have received a copy of the Apache-2.0 license                */
21   	/*  along with SoPlex; see the file LICENSE. If not email to soplex@zib.de.  */
22   	/*                                                                           */
23   	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24   	
25   	#ifndef _SOPLEX_STABLE_SUM_H_
26   	#define _SOPLEX_STABLE_SUM_H_
27   	#include <type_traits>
28   	
29   	// #define SOPLEX_CHECK_STABLESUM  // double check the stable sum computation
30   	
31   	namespace soplex
32   	{
33   	
34   	template <typename T>
35   	class StableSum
36   	{
37   	   typename std::remove_const<T>::type sum;
38   	
39   	public:
40   	   StableSum() : sum(0) {}
41   	   StableSum(const T& init) : sum(init) {}
42   	
43   	   void operator+=(const T& input)
44   	   {
45   	      sum += input;
46   	   }
47   	
48   	   void operator-=(const T& input)
49   	   {
50   	      sum -= input;
51   	   }
52   	
53   	   operator typename std::remove_const<T>::type() const
54   	   {
55   	      return sum;
56   	   }
57   	
58   	};
59   	
60   	template <>
61   	class StableSum<double>
62   	{
63   	   double sum = 0;
64   	   double c = 0;
65   	#ifdef SOPLEX_CHECK_STABLESUM
66   	   double checksum = 0;
67   	#endif
68   	
69   	public:
70   	   StableSum() = default;
71   	   StableSum(double init) : sum(init), c(0) {}
72   	
73   	   void operator+=(double input)
74   	   {
75   	#if defined(_MSC_VER) || defined(__INTEL_COMPILER)
76   	#pragma float_control( precise, on )
77   	#endif
78   	#ifdef SOPLEX_CHECK_STABLESUM
79   	      checksum += input;
80   	#endif
81   	      double t = sum + input;
82   	      double z = t - sum;
83   	      double y = (sum - (t - z)) + (input - z);
84   	      c += y;
85   	
86   	      sum = t;
87   	   }
88   	
89   	   void operator-=(double input)
90   	   {
91   	      (*this) += -input;
92   	   }
93   	
94   	   operator double() const
95   	   {
96   	#ifdef SOPLEX_CHECK_STABLESUM
97   	
98   	      if(spxAbs(checksum - (sum + c)) >= 1e-6)
99   	         printf("stablesum viol: %g, rel: %g, checksum: %g\n", spxAbs(checksum - (sum + c)),
100  	                spxAbs(checksum - (sum + c)) / SOPLEX_MAX(1.0, SOPLEX_MAX(spxAbs(checksum), spxAbs(sum + c))),
101  	                checksum);
102  	
103  	      assert(spxAbs(checksum - (sum + c)) < 1e-6);
104  	#endif
105  	      return sum + c;
106  	   }
107  	};
108  	
109  	/// Output operator.
110  	template < class T >
111  	std::ostream& operator<<(std::ostream& s, const StableSum<T>& sum)
112  	{
113  	   s << static_cast<T>(sum);
114  	
115  	   return s;
116  	}
117  	
118  	}
119  	
120  	#endif
121