#include<stdio.h>
#include<time.h>
#include<math.h>

double compute_iterative(int);
double compute_recursive(int n);
double compute_formula(int);
double compute_matrix(int);

double power(double, int);
void power_matrix(double*,double*,double*,double*,int);

int main(){
	int n;
	printf("Write n: ");
	scanf("%d",&n);
	
	clock_t start, end;
	double cpu_time_used;
	
	start=clock();
	printf("I_%d (computed using iterative method) = %lf\n", n, compute_iterative(n));
	end=clock();
	printf("Time taken in iterative method = %lf seconds\n\n",((double) (end - start)) / CLOCKS_PER_SEC);
	
	start=clock();
	printf("I_%d (computed using recursive method) = %lf\n", n, compute_recursive(n));
	end=clock();
	printf("Time taken in recursive method = %lf seconds\n\n",((double) (end - start)) / CLOCKS_PER_SEC);
	
	start=clock();
	printf("I_%d (computed using formula) = %lf\n", n, compute_formula(n));
	end=clock();
	printf("Time taken in formula = %lf seconds\n\n",((double) (end - start)) / CLOCKS_PER_SEC);
	
	start=clock();
	printf("I_%d (computed using matrix multiplication) = %lf\n", n, compute_matrix(n));
	end=clock();
	printf("Time taken in matrix = %lf seconds\n\n",((double) (end - start)) / CLOCKS_PER_SEC);
	
	return 0;
}

double compute_iterative(int n){
	double previous=0, current =1, next;
	if(n==0)	return 0;
	else{
		for(int i=1; i<n; i++){
			next=2*current+previous;
			previous=current;
			current=next;
		}
		return current;
	}
}

double compute_recursive(int n){
	if (n==0) return 0;
	else if (n==1) return 1;
	else return 2*compute_recursive(n-1)+compute_recursive(n-2);	
}


double compute_formula(int n){
	return (power(1+sqrt(2),n) - power(1-sqrt(2),n))/(2*sqrt(2));
}

double power(double x, int n){
	if (n == 0) return 1;
	else if (n%2 == 1) return x*power(x*x, n/2);
	else return power(x*x, n/2);
}

void power_matrix(double *a11, double *a12, double *a21, double *a22, int n){
	double b11, b12, b21, b22; // temporary matrix
	if (n == 0) {
		*a11 = 1.0; *a22=1.0; *a12=0.0; *a21=0.0;
	}
	else{	// matrix b stores a^2
		b11 = (*a11)*(*a11)+(*a12)*(*a21);
		b12 = (*a11)*(*a12)+(*a12)*(*a22);
		b21 = (*a21)*(*a11)+(*a22)*(*a21);
		b22 = (*a21)*(*a12)+(*a22)*(*a22);
		
		// we compute b^{n/2}
		power_matrix(&b11, &b12, &b21, &b22, n/2);
			
		if (n%2 == 0){ 
			*a11=b11, *a12=b12, *a21=b21, *a22=b22;
		}
		else{
			double c11, c12, c21, c22; // temporary matrix, c=a*b
			c11 = (*a11)*(b11)+(*a12)*(b21);
			c12 = (*a11)*(b12)+(*a12)*(b22);
			c21 = (*a21)*(b11)+(*a22)*(b21);
			c22 = (*a21)*(b12)+(*a22)*(b22);
			
			*a11=c11, *a12=c12, *a21=c21, *a22=c22;
		}
	}
}	

double compute_matrix(int n){
	double a11, a12, a21, a22;
	a11=2.0, a12=1.0, a21=1.0, a22=0.0;
	power_matrix(&a11, &a12, &a21, &a22, n);
	return a12;
}