diff --git a/tests/bessel_precision/compileall.sh b/tests/bessel_precision/compileall.sh
index 50bf761..e16e440 100755
--- a/tests/bessel_precision/compileall.sh
+++ b/tests/bessel_precision/compileall.sh
@@ -1,8 +1,8 @@
 #!/bin/bash
 
-c99 -ggdb -o jtest -DJTEST besseltest.c
-c99 -ggdb -o ytest -DYTEST besseltest.c
-c99 -ggdb -o djtest -DDJTEST besseltest.c
-c99 -ggdb -o dytest -DDYTEST besseltest.c
-c99 -ggdb -o djtest_steed -DJTEST_STEED besseltest.c
+c99 -ggdb -o jtest -DJTEST besseltest.c -lgsl -lblas -lm
+c99 -ggdb -o ytest -DYTEST besseltest.c -lgsl -lblas -lm
+c99 -ggdb -o djtest -DDJTEST besseltest.c -lgsl -lblas -lm
+c99 -ggdb -o dytest -DDYTEST besseltest.c -lgsl -lblas -lm
+c99 -ggdb -o djtest_steed -DJTEST_STEED besseltest.c -lgsl -lblas -lm
 
diff --git a/tests/bessel_precision/sagebesselgen.py b/tests/bessel_precision/sagebesselgen.py
index 6e90aaa..d232f0d 100644
--- a/tests/bessel_precision/sagebesselgen.py
+++ b/tests/bessel_precision/sagebesselgen.py
@@ -31,8 +31,30 @@ def printbesselDJrow(lMax, x, file=sys.stdout):
 def printbesselDYrow(lMax, x, file=sys.stdout):
     print(lMax, N(x), end=' ', file=file);
     for l in range(lMax+1):
-        print(N(-spherical_bessel_J(l+1,x)
-            + l/x*spherical_bessel_J(l,x)), end = ' ', file=file)
+        print(N(-spherical_bessel_Y(l+1,x)
+            + l/x*spherical_bessel_Y(l,x)), end = ' ', file=file)
     print('', file=file)
 
 
+def genall(lMax):
+    f = open('besselDJcases', 'w')
+    for o in IntegerRange(1,100):
+        printbesselDJrow(lMax, o, file=f)
+        printbesselDJrow(lMax, 1/o, file=f)
+        printbesselDJrow(lMax, o/sqrt(3), file=f)
+    f = open('besselDYcases', 'w')
+    for o in IntegerRange(1,100):
+        printbesselDYrow(lMax, o, file=f)
+        printbesselDYrow(lMax, 1/o, file=f)
+        printbesselDYrow(lMax, o/sqrt(3), file=f)
+    f = open('besselJcases', 'w')
+    for o in IntegerRange(1,100):
+        printbesselJrow(lMax, o, file=f)
+        printbesselJrow(lMax, 1/o, file=f)
+        printbesselJrow(lMax, o/sqrt(3), file=f)
+    f = open('besselYcases', 'w')
+    for o in IntegerRange(1,100):
+        printbesselYrow(lMax, o, file=f)
+        printbesselYrow(lMax, 1/o, file=f)
+        printbesselYrow(lMax, o/sqrt(3), file=f)
+