1 files changed, 74 insertions, 0 deletions
diff --git a/libm/double/bernum.c b/libm/double/bernum.c
new file mode 100644
index 000000000..e401ff5df
--- /dev/null
+++ b/libm/double/bernum.c
@@ -0,0 +1,74 @@
+/* This program computes the Bernoulli numbers.
+ * See radd.c for rational arithmetic.
+ */
+
+typedef struct{
+	double n;
+	double d;
+	}fract;
+
+#define PD 44
+fract x[PD+1] = {0.0};
+fract p[PD+1] = {0.0};
+#include <math.h>
+#ifdef ANSIPROT
+extern double fabs ( double );
+extern double log10 ( double );
+#else
+double fabs(), log10();
+#endif
+extern double MACHEP;
+
+main()
+{
+int nx, np, nu;
+int i, j, k, n, sign;
+fract r, s, t;
+
+
+for(i=0; i<=PD; i++ )
+	{
+	x[i].n = 0.0;
+	x[i].d = 1.0;
+	p[i].n = 0.0;
+	p[i].d = 1.0;
+	}
+p[0].n = 1.0;
+p[0].d = 1.0;
+p[1].n = 1.0;
+p[1].d = 1.0;
+np = 1;
+x[0].n = 1.0;
+x[0].d = 1.0;
+
+for( n=1; n<PD-2; n++ )
+{
+
+/* Create line of Pascal's triangle */
+/* multiply p = u * p */
+for( k=0; k<=np; k++ )
+	{
+	radd( &p[np-k+1], &p[np-k], &p[np-k+1] );
+	}
+np += 1;
+
+/* B0 + nC1 B1 + ... + nCn-1 Bn-1 = 0 */
+s.n = 0.0;
+s.d = 1.0;
+
+for( i=0; i<n; i++ )
+	{
+	rmul( &p[i], &x[i], &t );
+	radd( &s, &t, &s );
+	}
+
+
+rdiv( &p[n], &s, &x[n] );	/* x[n] = -s/p[n] */
+x[n].n = -x[n].n;
+nx += 1;
+printf( "%2d %.15e / %.15e\n", n, x[n].n, x[n].d );
+}
+
+
+}
+