Series sum calculator

Question

This is my program currently to calculate sum of following series:

$$\frac{x^2}{2!} + \frac{x^4}{4!} + \frac{x^6}{6!} \dots$$

#include <iostream>
using namespace std;
int main() {
 int fac = 1;
 float x, sum = 0, p = 1;
 cout << "\nENTER x: ";
 cin >> x;
 p = x;
 for(int i = 2; i <= 10; i++) {
 fac *= i;
 p *= x;
 cout << p << " ";
 if(i%2 == 0) {
 sum += p/fac;
 }
 }
 cout << " >> " << sum;
 return 0;
}

Is it well written? Can it be used without any optimisation for larger n values (currently only tried with small values)? How can I optimize these calculations?

• Rob Audenaerde: For large n you might introduce a big error, as you do multiple multiplications. See: docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html
• Rob Audenaerde: Also, it looks like a taylor series (rusty knowledge).. Care to elaborate what you are trying to calculate (it is not e^x, but looks similar :))
• Max Payne: @RobAu well not exactly trying to calculate a function, but this should be cosh - 1 ;)
• Rob Audenaerde: You mean this one: mathworld.wolfram.com/InverseHyperbolicCosine.html ?
• Juho: In general, it is a good idea to precompute as much as possible, especially if that something won't take much space. For instance, one can precompute the factorial values and avoid computing them altogether.

Answer 1

Here is a rewritten version of your sum which incorporates some techniques which improve performance, numeric stability and the accuracy of the final result.

• The two multiplicative terms (the power of x and the factorial) are combined into a single accumulated value (the variable named "term"). Obviously, for larger x, the two separate values grow rapidly, quickly surpassing the capacity of their types. But the combined term (power / factorial) is well-behaved for a much broader range of input values, usually going to zero as iterations accumulate.

• Each iteration of the loop corresponds to one term in the summation. There is no fundamental reason for each iteration to correspond to each integer power of x. For this particular series, the term's power of x goes up by 2, so I precompute the square of x in order to do this with a single multiplication. This also gets rid of that awkward conditional inside the loop (if(i%2==0)) (probably optimized away by the compiler, but still best to be rid of it).

• The last trick automatically determines the iteration count by comparing this iteration's updated sum (the variable next_sum) with the last iteration's sum (the running sum variable). If the sum hasn't changed, it stops iterating and returns that value. This technique causes problems in several situations, however and cannot be universally relied upon. For example, some series don't sum monotonically to a limit -- the sum will bounce around it, and the summation's state can end up oscillatory forever. Also, this relies on a fixed-precision data type (like float or double). In general, more work should probably be put into deciding when the sum's precision is acceptable and no further iterations need to be calculated. But, this trick is so super-cheap and easy to implement and it is usually so good at detecting the best final result that I include it here, but use it at your peril.

  #include <iostream>
  #include <iomanip>
  int main() {
   float x;
   std::cout << "\nENTER x: ";
   std::cin >> x;
   std::cout << std::setprecision(9);
   const float x_squared = x * x;
   float sum = 0;
   float term = 1;
   for(int i = 2;; i += 2) {
   term *= x_squared / (i * (i - 1));
   const float next_sum = sum + term;
   if(next_sum == sum) break;
   sum = next_sum;
   std::cout << std::setw(2) << i << ": " << sum << "\n";
   }
   // std::cout << ">> " << sum << "\n";
  }
  

Here are some results. I changed the code to dump the term's x exponent and the running sum at each iteration, rather than dumping the exponentiated x.

The original version and my version dump the same figures when given an input value of 1.2:

ENTER x: 1.2
 2: 0.720000029
 4: 0.806400061
 6: 0.810547233
 8: 0.810653865
10: 0.810655594

For an input value of 10, Wolfram says we should get about 11012.232920... But with 10, the original code has problems (I went ahead and set the max iterated power to 30, just to see where things end up).

ENTER x: 10
 2: 50
 4: 466.666656
 6: 1855.55554
 8: 4335.71436
10: 7091.44629
12: 9179.12207
14: 87368.5547
16: 5076917.5
18: -1.1079721e+009
20: -4.86787072e+010
22: -1.91795581e+013
24: -1.30792867e+015
26: -5.52487843e+016
28: -7.32410864e+018
30: 7.02255016e+020

My version (For an input value of 10, it stops at a max x power of 30):

ENTER x: 10
 2: 50
 4: 466.666656
 6: 1855.55542
 8: 4335.71387
10: 7091.44531
12: 9179.12109
14: 10326.1953
16: 10804.1426
18: 10960.335
20: 11001.4385
22: 11010.335
24: 11011.9463
26: 11012.1943
28: 11012.2275
30: 11012.2314

And this is very close to Wolfram's cosh(10)-1, given the precision of a 32-bit float.

• Christopher Oicles: Thanks for the formatting help Quill, I'm not very good at getting the nuances right.
• Max Payne: THanks! i too thought those extra iterations were useless! recursive approach is better.

Answer 2

Growth rate of factorial is a well-known problem. A standard technique to overcome it is to not calculate factorials at all: rewrite your series as

\$\frac{x^2}{1\cdot2}(1 + \frac{x^2}{3\cdot4}(1 + \frac{x^2}{5\cdot6}(\cdot\cdot\cdot(1 + \frac{x^2}{(2n-1)\cdot2n})))\$

and work from inside.