diff --git a/pipepar/pi_pipe.c b/pipepar/pi_pipe.c
new file mode 100644
index 0000000..578afa7
--- /dev/null
+++ b/pipepar/pi_pipe.c
@@ -0,0 +1,265 @@
+/*
+ * Computation of the n'th decimal digit of pi with very little memory.
+ * Written by Fabrice Bellard on February 26, 1997.
+ * 
+ * We use a slightly modified version of the method described by Simon
+ * Plouffe in "On the Computation of the n'th decimal digit of various
+ * transcendental numbers" (November 1996). We have modified the algorithm
+ * to get a running time of O(n^2) instead of O(n^3log(n)^3).
+ *
+ * This program uses a variation of the formula found by Gosper in 1974 :
+ * 
+ * pi = sum( (25*n-3)/(binomial(3*n,n)*2^(n-1)), n=0..infinity);
+ *
+ * This program uses mostly integer arithmetic. It may be slow on some
+ * hardwares where integer multiplications and divisons must be done by
+ * software. We have supposed that 'int' has a size of at least 32 bits. If
+ * your compiler supports 'long long' integers of 64 bits, you may use the
+ * integer version of 'mul_mod' (see HAS_LONG_LONG).  
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <limits.h>
+#include <errno.h>
+#include <math.h>
+
+
+#define STRTOL_OVERFLOW(n) (((n == LONG_MIN) || (n == LONG_MAX)) && \
+			    (errno == ERANGE))
+
+/* uncomment the following line to use 'long long' integers */
+#define HAS_LONG_LONG
+
+#ifdef HAS_LONG_LONG
+#define mul_mod(a,b,m) (( (long long) (a) * (long long) (b) ) % (m))
+#else
+#define mul_mod(a,b,m) fmod( (double) a * (double) b, m)
+#endif
+
+int nb_cores = 8; /* Bossa */
+//int nb_cores = 12; /* Quad Hexa */
+
+/* General note on variable names: _ denotes a minus sign */
+
+/* return the inverse of x mod y */
+int inv_mod(int x,int y) {
+  int q,u,v,a,c,t;
+
+  u=x;
+  v=y;
+  c=1;
+  a=0;
+  do {
+    q=v/u;
+    
+    t=c;
+    c=a-q*c;
+    a=t;
+    
+    t=u;
+    u=v-q*u;
+    v=t;
+  } while (u!=0);
+  a=a%y;
+  if (a<0) a=y+a;
+  return a;
+}
+
+/* return the inverse of u mod v, if v is odd */
+int inv_mod2(int u,int v) {
+  int u1,u3,v1,v3,t1,t3,y4;
+  
+  y4=0;
+
+  u1=1;
+  u3=u;
+  
+  v1=v;
+  v3=v;
+  
+  if ((u&1)!=0) {
+    t1=0;
+    t3=-v;
+  } else {
+    t1=1;
+    t3=u;
+    y4=1;
+  }
+  
+  do {
+    
+    while (y4 || ((t3&1)==0)) {
+      if ((t1&1)==0) {
+	t1=t1>>1;
+	t3=t3>>1;
+      } else {
+	t1=(t1+v)>>1;
+	t3=t3>>1;
+      }
+      y4=0;
+    };
+    
+    if (t3>=0) {
+      u1=t1;
+      u3=t3;
+    } else {
+      v1=v-t1;
+      v3=-t3;
+    }
+    t1=u1-v1;
+    t3=u3-v3;
+    if (t1<0) {
+      t1=t1+v;
+    }
+  } while (t3 != 0);
+  return u1;
+}
+
+
+/* return (a^b) mod m */
+int pow_mod(int a,int b,int m)
+{
+  int r,aa;
+   
+  r=1;
+  aa=a;
+  while (1) {
+    if (b&1) r=mul_mod(r,aa,m);
+    b=b>>1;
+    if (b == 0) break;
+    aa=mul_mod(aa,aa,m);
+  }
+  return r;
+}
+      
+/* return true if n is prime */
+int is_prime(int n)
+{
+   int r,i;
+   if ((n % 2) == 0) return 0;
+
+   r=(int)(sqrt(n));
+   for(i=3;i<=r;i+=2) if ((n % i) == 0) return 0;
+   return 1;
+}
+
+/* return the prime number immediatly after n */
+int next_prime(int n)
+{
+   do {
+      n++;
+   } while (!is_prime(n));
+   return n;
+}
+
+#define DIVN(t,a,v,vinc,kq,kqinc)		\
+{						\
+  kq+=kqinc;					\
+  if (kq >= a) {				\
+    do { kq-=a; } while (kq>=a);		\
+    if (kq == 0) {				\
+      do {					\
+	t=t/a;					\
+	v+=vinc;				\
+      } while ((t % a) == 0);			\
+    }						\
+  }						\
+}
+
+int nth_digit_pi(int n)
+{
+  int av,a,vmax,N,num,den,k,kq1,kq2,kq3,kq4,t,v,s,i,t1;
+  double sum;
+
+  N=(int)((n+20)*log(10)/log(13.5));
+  sum=0;
+
+  for(a=2;a<=(3*N);a=next_prime(a)) {
+    vmax=(int)(log(3*N)/log(a));
+    if (a==2) {
+      vmax=vmax+(N-n);
+      if (vmax<=0) continue;
+    }
+    av=1;
+    for(i=0;i<vmax;i++) av=av*a;
+
+    s=0;
+    den=1;
+    kq1=0;
+    kq2=-1;
+    kq3=-3;
+    kq4=-2;
+    if (a==2) {
+      num=1;
+      v=-n; 
+    } else {
+      num=pow_mod(2,n,av);
+      v=0;
+    }
+
+    for(k=1;k<=N;k++) {
+
+      t=2*k;
+      DIVN(t,a,v,-1,kq1,2);
+      num=mul_mod(num,t,av);
+      
+      t=2*k-1;
+      DIVN(t,a,v,-1,kq2,2);
+      num=mul_mod(num,t,av);
+
+      t=3*(3*k-1);
+      DIVN(t,a,v,1,kq3,9);
+      den=mul_mod(den,t,av);
+
+      t=(3*k-2);
+      DIVN(t,a,v,1,kq4,3);
+      if (a!=2) t=t*2; else v++;
+      den=mul_mod(den,t,av);
+      
+      if (v > 0) {
+	if (a!=2) t=inv_mod2(den,av);
+	else t=inv_mod(den,av);
+	t=mul_mod(t,num,av);
+	for(i=v;i<vmax;i++) t=mul_mod(t,a,av);
+	t1=(25*k-3);
+	t=mul_mod(t,t1,av);
+	s+=t;
+	if (s>=av) s-=av;
+      }
+    }
+
+    t=pow_mod(5,n-1,av);
+    s=mul_mod(s,t,av);
+    sum=fmod(sum+(double) s/ (double) av,1.0);
+  }
+  printf("Decimal digits of pi at position %d: %09d\n",n,(int)(sum*1e9));
+  return 0;
+}
+   
+int main(int argc, char *argv[])
+{
+	long pi_digit = -1;
+
+	if (argc == 2)
+	{
+		char *endptr;
+
+		pi_digit = strtol(argv[1], &endptr, 10);
+		if ((endptr == argv[1]) || (STRTOL_OVERFLOW(pi_digit)) ||
+			pi_digit <= 0)
+		{
+			fprintf(stderr, "Invalid digit for pi: %s\n", argv[1]);
+			return -1;
+		}
+	}
+	else
+	{
+		printf("This program computes the n'th decimal digit of pi.\n"
+			"Usage: %s n , where n is the digit you want\n", argv[0]
+			);
+		return -1;
+	}
+
+	return nth_digit_pi(pi_digit);
+}