# Effectively calculate the result of geometric series

Given $f(n) = 1 + x + x^2 + x^3 + \ldots + x^n$ and the fact that computers take more time when multiplying two numbers than when adding, how can we work out the result with greater efficiency?

#include <iostream>
using namespace std;

double powerWithIntExponent(double base, int exponent)
{
if(exponent == 0)
return 1;
if (exponent == 1)
return base;

double result = powerWithIntExponent(base, exponent >> 1);
result *= result;
if ((exponent & 0x1) == 1) {
result *= base;
}
return result;
}

double func(double x, int n)
{
if (x == 0) return 1;
if (x == 1) return n+1;

double sum = 0;
double x_n = 0;
x_n = powerWithIntExponent(x, n+1);// the Geometric series start from x^0
sum = (1-x_n)/(1-x);
return sum;

}
int main(int argc, const char * argv[])
{

double result = func(2, 2);
cout<<result<<endl;
return 0;
}


Instead of using a recursive algorithm (which might overflow the stack for large exponents) and calculating the power yourself, use the pow function which is likely pretty optimized.

#include <iostream>
#include <cmath>

using namespace std;

double func(double x, int n)
{
double x_n;

x_n = pow(x, n+1);
return (1-x_n)/(1-x);
}

int main(int argc, const char * argv[])
{

double result = func(2, 2);
cout<<result<<endl;
return 0;
}


Now a few minor remarks about your original code: be sure to be more consistent with your coding style. For example:

if(exponent == 0)
return 1;
if (exponent == 1)
return base;


vs.

if (x == 0) return 1;
if (x == 1) return n+1;


Personally, I'd always either write it on one line as in the second example or if there is a newline, I'd always add braces. For reasons see Apple's goto fail;. But it's actually a personal preference, you don't need to do it that way but I highly recommend that you chose one style and stick to it.

(Same with whitespace: there should be a space after the if in if(exponent == 0))

Consistency would also be nice here:

if ((exponent & 0x1) == 1)


Either use 0x1 both times or 1 both times. I feel that mixing these isn't a good idea.

• You reversed one of the most important points in the OP's code: Not iterating over the $x^i$ but using the finite formula $\frac{1-q^n}{1-q}$. Sep 12, 2014 at 9:36
• Since OP's code uses one stack frame (function call) for each bit of exponent, it will not cause a stack overflow. Sep 12, 2014 at 10:06
• @Nobody: Ouch! You're right, have fixed it. Sep 12, 2014 at 10:09
• @abuzittingillifirca: It can still allocate up to 31 or 63 stack frames when the exponent is positive, depending on sizeof(int) and the value. That's not good but usually not a problem (depending on stack size). But it will cause a stack overflow for n < -1 since it'll recurse endlessly. Of course, n < 0 should not be allowed in the first place. Sep 12, 2014 at 10:19
• I agree with your last points. Now that I've seen your edit+last comment my answer's become somewhat moot. I've got segfault for negative exponents, but it depends on whether the right shift sign-extends. exponent/2 would be better than exponent>>1 anyways. Sep 12, 2014 at 11:23

## Style

func(x, n) is a lousy name. Why not call it geometricSeries(x, n)?

You don't need the x == 0 special case.

## Performance

Your question is based on a faulty premise. The FMUL (floating-point multiplication) and FADD (floating-point addition) instructions take the same number of clock cycles on x87-compatible processors, pretty much universally.

The only time performance would be a concern is when n is large. In such circumstances, you would be better off using the pow() function, which probably works using logarithms. For comparison, I ran the following benchmark:

#include <cmath>
#include <iostream>
#include <random>
#include <sys/time.h>
#include <vector>

double powerWithIntExponent(double base, int exponent)
{
if(exponent == 0)
return 1;
if (exponent == 1)
return base;

double result = powerWithIntExponent(base, exponent >> 1);
result *= result;
if ((exponent & 0x1) == 1) {
result *= base;
}
return result;
}

std::vector<double> randomNumbers(size_t len, double min, double max) {
std::default_random_engine generator;
std::uniform_real_distribution<double> distribution(min, max);
std::vector<double> results(len);
for (int i = 0; i < len; ++i) {
results[i] = distribution(generator);
}
return results;
}

std::vector<int> randomNumbers(size_t len, int min, int max) {
std::default_random_engine generator;
std::uniform_int_distribution<int> distribution(min, max);
std::vector<int> results(len);
for (int i = 0; i < len; ++i) {
results[i] = distribution(generator);
}
return results;
}

int main() {
const size_t len = 10000000;
std::vector<double> xs = randomNumbers(len, 1e-5, 500.0);
std::vector<int>    ns = randomNumbers(len, 0, 50);
std::vector<double> r1(len), r2(len);

timeval t1, t2;
double elapsedMillis;

gettimeofday(&t1, NULL);
for (int i = 0; i < len; ++i) {
r2[i] = pow(xs[i], ns[i]);
}
gettimeofday(&t2, NULL);
elapsedMillis = (t2.tv_sec  - t1.tv_sec)  * 1000 +
(t2.tv_usec - t1.tv_usec) / 1000;
std::cout << "pow(): " << elapsedMillis << std::endl;
gettimeofday(&t1, NULL);

for (int i = 0; i < len; ++i) {
r1[i] = powerWithIntExponent(xs[i], ns[i]);
}
gettimeofday(&t2, NULL);
elapsedMillis = (t2.tv_sec  - t1.tv_sec)  * 1000 +
(t2.tv_usec - t1.tv_usec) / 1000;
std::cout << "powerWithIntExponent(): " << elapsedMillis << std::endl;

}


On my machine, your powerWithIntExponent() runs about 20% slower. Of course, the results will differ widely depending on what range you choose for n.

Note that for extremely large results, the two implementations can produce answers that differ significantly. I don't know which one is more accurate.

• As the implementation via logarithm uses approximations it is highly probable that the "explicit" integer power via multiplication is more accurate for large n. Sep 12, 2014 at 9:51

You need not use recursion in your function:

double powerWithIntExponent(double base, int exponent)
{
if(exponent == 0)
return 1;
if (exponent == 1)
return base;

int k = sizeof(int) * CHAR_BIT - 1;

double result = 1.0;

for (int i = k-1; i >= 0; --i) {
result *= result;
if (exponent & (1 << i)) {
result *= base;
}
}
return result;
}

result *= result;
if (exponent & ...) {
result *= base;
}


Also you can declare exponent unsigned, or support negative exponents by adding another guard clause:

if (exponent < 0)
return 1 / powerWithIntExponent(base, -exponent);


Currently passing in a negative exponents causes runtime error.

I am confused by your question. In the text you seem to imply you actually want to sum the powers of x as part of the program, but the program just uses the formula for a finite part of the geometric series. I would assume the latter (i.e. using the formula) will almost always be faster.

If you did actually want to perform the addition the method I'm familiar with is using the identity

$$a_0+a_1x+a_2x^2\cdots +a_nx^n=((\cdots(a_nx+a_{n-1})x+\cdots)x+a_0)$$

for evaluating general polynomials. The gain in this case is using only $n$ multiplications and $n+1$ additions compared to the usual $2n$ multiplications and $n$ additions.

• You are overgeneralizing the problem here. The OP talks about a geometric series where the usage of Horner's method is not needed as there exists a constant time solution $a\frac{1 - x^n}{1 - x}$. Sep 12, 2014 at 11:21
• Thanks, the reason why I use the formula is that it is the best answer I can figure out. I look through the Internet and find your solution may be the better answer to above. It's a coding interview question, and I don't know the proper solution yet. Sep 12, 2014 at 11:23
• @Nobody While I understand your point, if that is the question it seems particularly trivial. Especially if you consider $x^n$ to be constant time. Asking someone to write a computer program to compute a formula seems hardly worth the question. I suppose for C you actually have to use a function. Still the function double geom(double x,int n){ Sep 12, 2014 at 11:55