/*
 
 Demo code for the geometric series to show
 - round-off errors
 - properties of floating point numbers
 - termination errors
 
 use by changing the parameters q, N, and version
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

int main(int argc, char* argv[])
{
	double q,lq,x,sum,invq;
	int i,N,version;
    
    //---------------------------------------
    
    //change q, e.g., 0.9,0.5,0.2
    //change the number of terms added, N
    //version=1 : calculate q^i each time
    //version=2 : calculate q^i=q*q^{i-1}
    
    version=1;
	q=0.9;
	N=2000;
    
    //---------------------------------------
    
    printf("geometric series for q=%le and %d terms:\n",q,N);
    lq=log(q); //used to calculate q^i
    
    //forward summation
    x=1.0;
	sum=0.0;
    
	for(i=0;i<=N;i++) {
		if(version==1) x=exp(i*lq);
		sum+=x;
		if(version==2) x*=q;
	}
	printf("  forward summation: %.17le\n",sum);
    
    
    //reverse summation
	invq=1/q;
	x*=invq; //this is the last term added to the "forward sum" in version 2
	sum=0.0;
    
    for(i=N;i>=0;i--) {
        if(version==1) x=exp(i*lq);
        sum+=x;
        if(version==2) x*=invq;
    }
	printf("  reverse summation: %.17le\n",sum);
    
    //print exact result
	sum=1/(1-q);
	printf("  exact:             %.17le\n",sum);
    
	return 0;
};