HiveBrain v1.2.0
Get Started
← Back to all entries
patterncppMinor

Series sum calculator

Submitted by: @import:stackexchange-codereview··
0
Viewed 0 times
seriescalculatorsum

Problem

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 

using namespace std;
int main() {
    int fac = 1;
    float x, sum = 0, p = 1;

    cout > x;
    p =  x;
    for(int i = 2; i > " << 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?

Solution

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 
#include 

int main() {
    float x;
    std::cout > x;

    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.

Code Snippets

#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";
}
ENTER x: 1.2
 2: 0.720000029
 4: 0.806400061
 6: 0.810547233
 8: 0.810653865
10: 0.810655594
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
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

Context

StackExchange Code Review Q#101194, answer score: 4

Revisions (0)

No revisions yet.