diff --git a/lib/sfrsd2/sfrsd_paper/Makefile b/lib/sfrsd2/sfrsd_paper/Makefile
new file mode 100644
index 000000000..172adc275
--- /dev/null
+++ b/lib/sfrsd2/sfrsd_paper/Makefile
@@ -0,0 +1,31 @@
+CC = gcc
+FC = gfortran
+
+# Default rules
+%.o: %.c
+	${CC} ${CFLAGS} -c $<
+%.o: %.f
+	${FC} ${FFLAGS} -c $<
+%.o: %.F
+	${FC} ${FFLAGS} -c $<
+%.o: %.f90
+	${FC} ${FFLAGS} -c $<
+%.o: %.F90
+	${FC} ${FFLAGS} -c $<
+
+all:    probs.out
+
+OBJS1 = prob.o binomial_subs.o
+prob:   $(OBJS1)
+	$(FC) -o prob $(OBJS1)
+
+probs.out: prob
+#             x  N  X  s
+	prob 35 63 40 40 > probs.out
+	prob 37 63 40 45 >> probs.out
+	prob 37 53 40 45 >> probs.out
+	prob 38 53 40 47 >> probs.out
+
+clean:
+	rm -rf *.o prob probs.out
+
diff --git a/lib/sfrsd2/sfrsd_paper/binomial_subs.c b/lib/sfrsd2/sfrsd_paper/binomial_subs.c
new file mode 100644
index 000000000..32b9e94ae
--- /dev/null
+++ b/lib/sfrsd2/sfrsd_paper/binomial_subs.c
@@ -0,0 +1,55 @@
+#include <stdio.h>
+#include <limits.h>
+
+/* Original code copied from 
+   http://rosettacode.org/wiki/Evaluate_binomial_coefficients
+*/
+ 
+/* We go to some effort to handle overflow situations */
+ 
+static unsigned long gcd_ui(unsigned long x, unsigned long y) {
+  unsigned long t;
+  if (y < x) { t = x; x = y; y = t; }
+  while (y > 0) {
+    t = y;  y = x % y;  x = t;  /* y1 <- x0 % y0 ; x1 <- y0 */
+  }
+  return x;
+}
+ 
+unsigned long binomial(unsigned long n, unsigned long k) {
+  unsigned long d, g, r = 1;
+  if (k == 0) return 1;
+  if (k == 1) return n;
+  if (k >= n) return (k == n);
+  if (k > n/2) k = n-k;
+  for (d = 1; d <= k; d++) {
+    if (r >= ULONG_MAX/n) {  /* Possible overflow */
+      unsigned long nr, dr;  /* reduced numerator / denominator */
+      g = gcd_ui(n, d);  nr = n/g;  dr = d/g;
+      g = gcd_ui(r, dr);  r = r/g;  dr = dr/g;
+      if (r >= ULONG_MAX/nr) return 0;  /* Unavoidable overflow */
+      r *= nr;
+      r /= dr;
+      n--;
+    } else {
+      r *= n--;
+      r /= d;
+    }
+  }
+  return r;
+}
+
+unsigned long binomial_(int *n, int *k)
+{
+  //  printf("n=%d  k=%d  %lu\n",*n,*k,binomial(*n,*k));
+  return binomial(*n,*k);
+}
+
+double hypergeo_(int *x, int *NN, int *XX, int *s)
+{    
+  double a,b,c;
+  a=(double)binomial(*XX, *x);
+  b=(double)binomial(*NN-*XX, *s-*x);
+  c=(double)binomial(*NN, *s);
+  return a*b/c;
+}
diff --git a/lib/sfrsd2/sfrsd_paper/prob.f90 b/lib/sfrsd2/sfrsd_paper/prob.f90
new file mode 100644
index 000000000..c11a579d8
--- /dev/null
+++ b/lib/sfrsd2/sfrsd_paper/prob.f90
@@ -0,0 +1,41 @@
+program prob
+
+  implicit real*8 (a-h,o-z)
+  integer*8 binomial
+  integer x,s,XX,NN
+  real*8 hypergeo
+  character arg*8
+
+  nargs=iargc()
+  if(nargs.ne.4) then
+     print*,'Usage:    prob  x  N  X  s'
+     print*,'Example:  prob 35 63 40 40'
+     go to 999
+  endif
+  call getarg(1,arg)
+  read(arg,*) x
+  call getarg(2,arg)
+  read(arg,*) NN
+  call getarg(3,arg)
+  read(arg,*) XX
+  call getarg(4,arg)
+  read(arg,*) s
+
+!  print*,binomial(5, 3)             ! 10
+!  print*,binomial(40, 19)           ! 131282408400
+
+  write(*,1010) x,NN,XX,s
+1010 format(//' x=',i2,'   N=',i2,'   X=',i2,'   s=',i2)
+  write(*,1012) 
+1012 format(/' x    P(x|N,X,s)  P(>=x|N,X,s) '/  &
+             '-------------------------------')
+
+  hsum=0.d0
+  do ix=x,XX
+     h=hypergeo(ix,NN,XX,s)
+     hsum=hsum + h
+     write(*,1020) ix,h,hsum
+1020 format(i3,2d13.4)
+  enddo
+
+999 end program prob