GNU Scientific Library (GSL)

[John]

If you have a need for a serious math library, check out the GSL offering. The following is an example (intro.c) I converted to ScriptBasic.

The project was conceived in 1996 by Dr M. Galassi and Dr J. Theiler of Los Alamos National Laboratory. They were joined by other physicists who also felt that the licenses of existing libraries were hindering scientific cooperation. Most of the library has been written by a relatively small number of people with backgrounds in computational physics in order to provide a consistent and reasonably-designed framework.

Regular Cylindrical Bessel Function

Code: Text

1. #include <stdio.h>
2. #include <gsl/gsl_sf_bessel.h>
3.  
4. int
5. main (void)
6. {
7.   double x = 5.0;
8.   double y = gsl_sf_bessel_J0 (x);
9.   printf ("J0(%g) = %.18e\n", x, y);
10.   return 0;
11. }
12.  

Note: I increased the precision to 48 from 18 in the C example.  (SB64 Ubuntu 11.04)

Code: Text

1. DECLARE SUB DLL ALIAS "_idll" LIB "gtk-server"
2. DECLARE SUB REQUIRE ALIAS "_idll_require" LIB "gtk-server"
3. DECLARE SUB DEFINE ALIAS "_idll_define" LIB "gtk-server"
4.  
5. REQUIRE "libgsl.so"
6. DEFINE "gsl_sf_bessel_J0 NONE DOUBLE 1 DOUBLE"
7.  
8. PRINT FORMAT("J0(%g) = %.48e\n", 5.0, DLL("gsl_sf_bessel_J0 " & 5.0))
9.  

J0(5) = -1.775970000000000326156879282279987819492816925049e-01

I may create a SB INCLUDE that wraps the GSL API and provides a suite of SB FUNCTIONs. If you would like to chip in and help create a SB extended math library, indicate what functions you are taking ownership of.

Update

Wolfram Alpha result.

I now need to determine if the SB FORMAT() function is the issue or what GTK-Server is returning as a DOUBLE.  

[John]

Here is the result of the C version.

gcc -Wall intro.c -lgslcblas -lgsl -o intro

J0(5) = -1.775967713143382642471124199801124632358551025391e-01

For grins I DEFINED the printf() and used it instead of the SB FORMAT() function.

DEFINE "printf NONE DOUBLE 3 STRING DOUBLE DOUBLE"

DLL("printf \"J0(%g) = %.48e\n\" 5.0 " & DLL("gsl_sf_bessel_J0 " & 5.0))

J0(5) = -1.775970000000000048601123125990852713584899902344e-01

My guess at this point is the problem is with GTK-Server and returning extended precision values.

Update

Code: [Select]

v = 1.775967713143382642471124199801124632358551025391e-01

PRINT FORMAT("J0(%g) = %.48e\n", 5.0, v)

J0(5) = 1.775967713143382642471124199801124632358551025391e-01

I'm convinced I have a GTK-Server issue with DOUBLE returns.

I thought I would try doing simple math with the extended precision variable to see what the affect might be.

Code: [Select]

v = 1.775967713143382642471124199801124632358551025391e-01

PRINT FORMAT("J0(%g) = %.48e\n", 5.0, v - .01)

J0(5) = 1.675967713143382553653282229788601398468017578125e-01

Doing the same on Wolfram Alpha produces a different result after a certain point in the precision.

[John]

I made a couple changes to GTK-Server to deal with the 64 bit doubles and it didn't seem to help. (changed from 32 to 64 in the following defines)

#define MAX_DIG 64

#define LONG_SIZE 64

I think the best thing in this case is to write a SB extension module in C and forgo the scripting aspect of the API.

[John]

I was able to get the GSL library to work correctly under Wine using the DYC extension module as my FFI rather than GTK-Server native under Ubuntu 64.

Code: Text

1. DECLARE SUB DLL ALIAS "dyc" LIB "dyc"
2.  
3. dv = DLL("ms8,d,libgsl-0.dll,gsl_sf_bessel_J0,d",5.0)
4.  
5. PRINT FORMAT("%.16g",dv)
6.  

C:\scriptbasic\test>scriba testgsl.sb
-0.1775967713143383
C:\scriptbasic\test>


The Bessel functions of the first kind J_n(x) are defined as the solutions to the Bessel differential equation

GSL.zip

The problem with GTK-Server is that it doesn't have a calling method to retrieve it's results for the FPU stack and in a 8 bit word. Maybe Peter can chime in and offer a solution.