1 files changed, 67 insertions, 0 deletions
diff --git a/libm/ldouble/wronkl.c b/libm/ldouble/wronkl.c
new file mode 100644
index 000000000..bec958f01
--- /dev/null
+++ b/libm/ldouble/wronkl.c
@@ -0,0 +1,67 @@
+/* Wronksian test for Bessel functions. */
+
+long double jnl (), ynl (), floorl ();
+#define PI 3.14159265358979323846L
+
+long double y, Jn, Jnp1, Jmn, Jmnp1, Yn, Ynp1;
+long double w1, w2, err1, max1, err2, max2;
+void wronk ();
+
+main ()
+{
+  long double x, delta;
+  int n, i, j;
+
+  max1 = 0.0L;
+  max2 = 0.0L;
+  delta = 0.6 / PI;
+  for (n = -30; n <= 30; n++)
+    {
+      x = -30.0;
+      while (x < 30.0)
+	{
+	  wronk (n, x);
+	  x += delta;
+	}
+      delta += .00123456;
+    }
+}
+
+void 
+wronk (n, x)
+     int n;
+     long double x;
+{
+
+  Jnp1 = jnl (n + 1, x);
+  Jmn = jnl (-n, x);
+  Jn = jnl (n, x);
+  Jmnp1 = jnl (-(n + 1), x);
+  /* This should be trivially zero.  */
+  err1 = Jnp1 * Jmn + Jn * Jmnp1;
+  if (err1 < 0.0)
+    err1 = -err1;
+  if (err1 > max1)
+    {
+      max1 = err1;
+      printf ("1 %3d %.5Le %.3Le\n", n, x, max1);
+    }
+  if (x < 0.0)
+    {
+      x = -x;
+      Jn = jnl (n, x);
+      Jnp1 = jnl (n + 1, x);
+    }
+  Yn = ynl (n, x);
+  Ynp1 = ynl (n + 1, x);
+  /* The Wronksian.  */
+  w1 = Jnp1 * Yn - Jn * Ynp1;
+  /* What the Wronksian should be. */
+  w2 = 2.0 / (PI * x);
+  err2 = w1 - w2;
+  if (err2 > max2)
+    {
+      max2 = err2;
+      printf ("2 %3d %.5Le %.3Le\n", n, x, max2);
+    }
+}