/* pi.c #2
 Original code from http://www.newdream.net/~sage/old/numbers/approxpi.htm
 Modifications (stern-brocot method code) by Robert Chin
 */

#define pi   "3.14159265358979323846264338327950288419716939937510"
#define ZERO " 0.0000000000000000000000000000000000000000000000000"
#define LEN 100


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

void sdiv(long num, long den, char *res, int prec);
void ssub(char *a, char *b, char *res);
void sadd(char *a, char *b, char *res);
int scomp(char *a, char *b);


main(int argc, char *argv[])
{
	long num[3],den[3];
	char acc[3][LEN],bacc[LEN],res[3][LEN], tmp[LEN];
	char *acc2;
    int swi;
	char tr;
	long cnt = 0;

	num[0]=3;
    den[0]=1;
    num[1]=4;
	den[1]=1;

	if (argc > 1) {
		num[0] = atol(argv[1]);
		num[1] = atol(argv[2]);
		den[0] = atol(argv[3]);
		den[1] = atol(argv[4]);
		if (den[0] == 0) den[0] = 1;
		if (den[1] == 0) den[1] = 1;
	}

	strcpy(bacc," 9.000000000000000000000000000000000000000000");

	printf("<pre>pi = %-30s\n\n",pi);
	printf("%13s %11s = %-32s (%-23s)\n","Num.","Den.","Result","Accuracy");
	printf("%13s %11s = %-32s (%-33s)\n","-------------","-----------","--------------------------------","---------------------------------");

    /* Uses stern-brocot method of numerical approximation. We only need to
    keep track of the two entries to be worked on, and a working entry. After
    calculating the working entry, we can then figure out whether it should
    replace the top or bottom entry to calculate a new entry. */
    swi = 0;
    sdiv(num[0],den[0],res[0],30);
    sdiv(num[1],den[1],res[1],30);
    ssub(pi,res[0],acc[0]);
    ssub(pi,res[1],acc[1]);
	while (1) {
        /* Add the two entries together to create a new entry. */
        num[2] = num[0] + num[1];
        den[2] = den[0] + den[1];
        sdiv(num[2], den[2], res[2], 30);
        ssub(pi, res[2], acc[2]);
        /* Figure out whether this is an entry should replace the top or
           bottom entry. */
        if(acc[2][0] == '-'){
            acc2 = acc[2] + sizeof(char);
            swi = 1;
            tr = scomp(acc[1], acc2);
        } else {
            acc2 = acc[2];
            swi = -1;
            tr = scomp(acc[0], acc2);
        }

/*		if (*kb & 1) cprintf("%7ld/%7ld = %.20s (%.20s)\r",num,den,res,acc);
*/
 		if (tr == 0) {
				cnt++;
/*				cprintf(" (%d)\n\r",cnt); */
		} else {
            /* Copy the number to either the top or boottom of the working
               area */
            if(swi == 1){
                strcpy(acc[1], acc2);
                strcpy(res[1], res[2]);
                num[1] = num[2];
                den[1] = den[2];
            } else {
                strcpy(acc[0], acc2);
                strcpy(res[0], res[2]);
                num[0] = num[2];
                den[0] = den[2];
            }
/*			cprintf("\n\r"); */
            tr = scomp(bacc, acc2);
            /* Uncomment the below line to get a trace of the computed values.*/
            //printf("\n%12ld/%12ld = %-30s (%33s)",num[2],den[2],res[2],acc[2]);
            if(tr > 0){
                strcpy(bacc, acc2);
                if (cnt) printf(" (x %ld)", cnt + 1);
            printf("\n%12ld/%12ld = %-30s (%33s)",num[2],den[2],res[2],acc[2]);
                cnt = 0;
                for (cnt=0; acc2[cnt+2] == '0'; cnt++);
                printf(" [%2d]", cnt);
                fflush(stdout);
                cnt = 0;
            }
        }
    }
	printf("\n</pre>\n");
}


/*  ssub - subtracts two strings
 *   MUST BE DECIMAL ALIGNED!
 */

int ssub_borrow(char *a, char p)
{
	if (p == -1) return 0;
	if (a[p] == '.') {
		return ssub_borrow(a,p-1);
	} else if (a[p] == '0') {
		a[p] = '9';
		return ssub_borrow(a,p-1);
	} else {
		a[p]--;
	}
	return 1;
}

void ssub(char *a, char *b, char *res)
{
	int p;
	int rp;
	int tr;

	p = 0;
	rp = 0;

	while (a[p] && b[p]) {
		if (a[p] == '.') {
			res[rp++] = '.';
			p++;
		}
		tr = a[p] - b[p];
		if (tr < 0) {
			if (ssub_borrow(res,p-1)) {
				tr += 10;
			} else {
				res[0] = '-';
				ssub(b,a,res+1);
				return;
			}
		}
		res[rp++] = '0' + tr;
		p++;
	}
	res[rp] = 0;
}

/* scomp - compares two strings
 *   MUST BE DECIMAL ALIGNED!
 *  return 0 - a < b
 *  return 1 - a > b
 */

int scomp(char *a, char *b)
{
	int p;

	p = 1;

    if (a[0] == '-' && b[0] != '-')
        return -1;
    if(a[0] != '-' && b[0] == '-')
        return 1;
	while (a[p] && b[p]) {
		if (a[p] < b[p]) return -1;
		if (a[p] > b[p]) return 1;
		p++;
	}
	return 0;
}

/* sdiv - divide two ints - result put in a string
 *
 */

void sdiv(long num, long den, char *res, int prec)
{
	int p,pp;
	unsigned long car,ocar,t;
	int rp;

	char nums[20];
	char unums;

	sprintf(nums,"%ld",num);

	rp = 0;
	unums = 1;
	car = 0;
	p = 0;
	pp = 0;
	res[0] = 0;

	while (pp < prec) {
		car *= 10;
		if (unums) {
			if (nums[p]) {
				car += nums[p] - '0';
			} else {
				unums = 0;
				res[rp++] = '.';
				pp++;
			}
		} else pp++;
		ocar = car;
		t = car / den;
		if (t || res[0]) res[rp++] = '0' + t;
		t *= den;
		car = ocar - t;
		p++;
		if (car == 0 && unums == 0) break;
	}
	while (pp < prec) {
		res[rp++] = '0';
		pp++;
	}
	res[rp] = 0;
}


/*  		001
	234 | 2354908123
		  02
			0
			23
		  = 0
			235
			234
		  =  14


*/

/*


3.1297321
3.1278412
0.0018909

 3.1297321
 3.1278412
-0.001
*/