/* Algorithm:
	a=1
	b=1/sqrt(2)
	t=1/4
	x=1
	loop until a-b is tiny
		y=a
		a=(a+b)/2
		b=sqrt(b*y)
		t=t-x*pow(a-y,2)
		x=2*x
	pi is pow(a+b,2)/(4*t)
*/

#include <stdio.h>
#include <strings.h>

#define NDIG 1050

/* implicit radix point after first "digit" */
static unsigned char a[NDIG];
static unsigned char t[NDIG];
static unsigned char b[NDIG];
static unsigned char y[NDIG];
static unsigned char tmp[NDIG];
static unsigned char tmp2[NDIG];
static int lgx;

static __inline__ int imin(int a, int b)
{
 return((a<b)?a:b);
}

static int incdigs(unsigned char *dv, int x)
{
 while (x >= 0)
  { if (dv[x] == 255)
     { dv[x] = 0;
       x --;
     }
    else
     { dv[x] ++;
       return(0);
     }
  }
 return(1);
}

static int subdigs(unsigned char *r, unsigned char *s, unsigned char *m, int n)
{
 int a;

 a = 0;
 for (n--;n>=0;n--)
  { a += (int)s[n] - (int)m[n];
    if (a < 0)
     { r[n] = a + 256;
       a = -1;
     }
    else
     { r[n] = a;
       a = 0;
     }
  }
 return(a);
}

static int adddigs(unsigned char *a, unsigned char *b, unsigned char *c, int n)
{
 int t;

 t = 0;
 for (n--;n>=0;n--)
  { t += b[n] + c[n];
    a[n] = t & 0xff;
    t >>= 8;
  }
 return(t);
}

static int submuldigs(unsigned char *s, unsigned char *m, int f, int n)
{
 int a;

 a = 0;
 for (n--;n>=0;n--)
  { a += (int)s[n] - (f * (int)m[n]);
    if (a < 0)
     { a += 65536;
       s[n] = a & 0xff;
       a = (a >> 8) - 256;
     }
    else
     { s[n] = a & 0xff;
       a >>= 8;
     }
  }
 return(a);
}

static void div3(unsigned char *q, unsigned char *n, unsigned char *d)
{
 unsigned char work[NDIG*2];
 int i;
 int qd;
 int r;
 int a;

 bcopy(n,&work[0],NDIG);
 bzero(&work[NDIG],NDIG);
 r = work[0];
 i = 0;
 while (1)
  { qd = ((r * 0x10000U) + (work[i+1] * 0x100U) + work[i+2]) /
	 ((d[0] * 0x10000U) + (d[1] * 0x100U) + d[2] + 1);
    if (qd > 0xff)
     { if (incdigs(q,i-1)) abort();
       r += subdigs(work+i+1,work+i+1,d+2,NDIG-2) - ((256 * (int)d[0]) + (int)d[1]);
       continue;
     }
    if (i >= NDIG) break;
    q[i] = qd;
    a = submuldigs(work+i+1,d+1,qd,NDIG-1);
    r += a - (d[0] * qd);
    if (r < 0) abort();
    r = (r * 256) + work[i+1];
    i ++;
  }
}

static void sub3(unsigned char *r, unsigned char *s, unsigned char *m)
{
 int i;

 for (i=0;i<NDIG;i++) if (s[i] != m[i]) break;
 if (i >= NDIG)
  { bzero(&r[0],NDIG);
    return;
  }
 if (s[i] > m[i]) subdigs(r,s,m,NDIG); else subdigs(r,m,s,NDIG);
}

static void add3(unsigned char *a, unsigned char *b, unsigned char *c)
{
 if (adddigs(a,b,c,NDIG)) abort();
}

static void mul3h(unsigned char *a, unsigned char *b, unsigned char *c)
{
 int i;
 int j;
 int k;
 unsigned int v;

 v = 0;
 for (i=NDIG-1;i>0;i--)
  { for (j=NDIG-1,k=i;j>=i;j--,k++)
     { v += (unsigned int)b[j] * (unsigned int)c[k];
     }
    v >>= 8;
  }
 for (i=NDIG-1;i>=0;i--)
  { for (j=i,k=0;j>=0;j--,k++)
     { v += (unsigned int)b[j] * (unsigned int)c[k];
     }
    a[i] = v & 0xff;
    v >>= 8;
  }
 if (v) abort();
}

static void mul3hshift(unsigned char *a, unsigned char *b, unsigned char *c, int sc)
{
 int i;
 int j;
 int k;
 int s;
 int o;
 unsigned int val;
 unsigned int buf;

 o = sc / 8;
 s = sc & 7;
 val = 0;
 buf = 0;
 for (i=NDIG-1;i>0;i--)
  { for (j=NDIG-1,k=i;j>=i;j--,k++)
     { val += (unsigned int)b[j] * (unsigned int)c[k];
     }
    buf |= (val & 0xff) << s;
    if (o >= i) a[NDIG-1+i-o] = buf & 0xff;
    val >>= 8;
    buf >>= 8;
  }
 for (i=NDIG-1;i>=0;i--)
  { for (j=i,k=0;j>=0;j--,k++)
     { val += (unsigned int)b[j] * (unsigned int)c[k];
     }
    buf |= (val & 0xff) << s;
    if (o > i)
     { if (buf) abort();
     }
    else
     { a[i-o] = buf & 0xff;
     }
    val >>= 8;
    buf >>= 8;
  }
 if (val) abort();
}

static void halve1(unsigned char *a)
{
 int i;
 int r;

 r = 0;
 for (i=0;i<NDIG;i++)
  { r = (256 * r) + a[i];
    a[i] = r >> 1;
    r &= 1;
  }
}

static void sqrtstep(unsigned char *n, unsigned char *a, unsigned char *b)
{
 unsigned char tmp[NDIG];

 div3(&tmp[0],n,a);
 add3(b,a,&tmp[0]);
 halve1(b);
}

static void printnum(unsigned char *n, int len)
{
 int i;

 for (i=0;i<len;i++) printf("%02x",n[i]);
}

static void dosqrt(unsigned char *n, unsigned char *r)
{
 unsigned char t[3][NDIG];
 int i3;
 int i2;
 int i;

 sqrtstep(&n[0],&r[0],&t[0][0]);
 sqrtstep(&n[0],&t[0][0],&t[1][0]);
 sqrtstep(&n[0],&t[1][0],&t[2][0]);
 i3 = 0;
 i2 = 1;
 i = 2;
 while (1)
  { i3 = i2;
    i2 = i;
    i ++;
    if (i > 2) i -= 3;
    sqrtstep(&n[0],&t[i2][0],&t[i][0]);
    if (!bcmp(&t[i][0],&t[i3][0],NDIG))
     { bcopy(&t[i][0],r,NDIG);
       return;
     }
  }
}

int main(void);
int main(void)
{
setlinebuf(stdout);
 b[0] = 0x00;
 b[1] = 0xb5;
 bzero(&b[2],NDIG-2);
 t[0] = 0;
 t[1] = 0x80;
 bzero(&t[2],NDIG-2);
 dosqrt(&t[0],&b[0]);
printf("sqrt(1/2) = ");
printnum(&b[0],NDIG);
printf("\n");
 a[0] = 1;
 bzero(&a[1],NDIG-1);
 t[1] = 0x40;
 lgx = 0;
 while (1)
  { bcopy(&a[0],&y[0],NDIG);
    add3(&a[0],&a[0],&b[0]);
    halve1(&a[0]);
printf("a = ");
printnum(&a[0],NDIG);
printf("\n");
    mul3h(&tmp[0],&b[0],&y[0]);
    dosqrt(&tmp[0],&b[0]);
printf("b = ");
printnum(&b[0],NDIG);
printf("\n");
    sub3(&tmp[0],&a[0],&y[0]);
    mul3hshift(&tmp2[0],&tmp[0],&tmp[0],lgx);
    sub3(&t[0],&t[0],&tmp2[0]);
printf("t = ");
printnum(&t[0],NDIG);
printf("\n");
    lgx ++;
    if (!bcmp(&a[0],&b[0],NDIG-1)) break;
  }
 mul3h(&tmp[0],&a[0],&b[0]);
 div3(&a[0],&tmp[0],&t[0]);
printf("pi = ");
printnum(&a[0],NDIG);
printf("\n");
 exit(0);
}