GNU Scientific Library (1000+ math/rand/stats etc funcs)

[Keya]

hello today while i was looking for statistics stuff i found GNU Scientific Library (GSL):
http://www.gnu.org/software/gsl/

The GNU Scientific Library (GSL) is a numerical library for C and C++ programmers. It is free software under the GNU General Public License. The library provides a wide range of mathematical routines such as random number generators, special functions and least-squares fitting. There are over 1000 functions in total with an extensive test suite.

It's originally for Linux, but a quick google suggests there's several OSX builds available, and it also says it supports Cygwin and i found a Win32 port ready to go at http://gnuwin32.sourceforge.net/packages/gsl.htm
(note - thats an old v1.8 build, for v2.1 jack made 32+64 builds available in this post)

Overview of functions - http://www.gnu.org/software/gsl/manual/html_node/

I've only just started with it a few minutes ago so nothing to share just yet, but I did port example.c using the Win32 DLL (libgsl.dll) to test it's ok, and i just wanted to add to the forum to make it available in searches etc, perhaps it'll come in handy to someone else as google suggests it hasn't been used in PB before

Code: Select all

;      int
;      main (void)
;      {
;        double x = 5.0;
;        double y = gsl_sf_bessel_J0 (x);
;        printf ("J0(%g) = %.18e\n", x, y);
;        Return 0;
;      }
; 
; The output is shown below, and should be correct to double-precision accuracy,
;      J0(5) = -1.775967713143382920e-01

Prototype.d proto_gsl_sf_bessel_J0(value.d)

If OpenLibrary(0, "libgsl.dll") 
    gsl_sf_bessel_J0.proto_gsl_sf_bessel_J0 = GetFunction(0, "gsl_sf_bessel_J0")   
    x.d = 5.0
    y.d = gsl_sf_bessel_J0(x)
    Debug("J0(5) = " + StrD(y))    ;-1.775967713
EndIf

[Keya]

Permutations ...

Code: Select all

#GSL_SUCCESS = 0

Structure gsl_permutation
  size.i
  *pdata
EndStructure


Prototype.i proto_gsl_permutation_init(*p.gsl_permutation)
Prototype.i proto_gsl_permutation_free(*p.gsl_permutation)
Prototype.i proto_gsl_permutation_alloc(value.i)
Prototype.i proto_gsl_permutation_next(*p.gsl_permutation)


If OpenLibrary(0, "libgsl.dll") 
  gsl_permutation_init.proto_gsl_permutation_init = GetFunction(0, "gsl_permutation_init")   
  gsl_permutation_free.proto_gsl_permutation_free = GetFunction(0, "gsl_permutation_free")   
  gsl_permutation_alloc.proto_gsl_permutation_alloc = GetFunction(0, "gsl_permutation_alloc")
  gsl_permutation_next.proto_gsl_permutation_next = GetFunction(0, "gsl_permutation_next")
  
  Define *p.gsl_permutation, *pval, numperms, sOut.s
  
  numperms = 3
  *p.gsl_permutation = gsl_permutation_alloc(numperms)
  *plong = *p\pdata
  gsl_permutation_init(*p)
  
  ;Initialize values to permutate, must be lexicographic order 
  ;(order can be reversed if using gsl_permutation_PREV instead of gsl_permutation_next)
  PokeL(*p\pdata+0, 2)
  PokeL(*p\pdata+4, 5)
  PokeL(*p\pdata+8, 8)
  
  Repeat
    sOut = ""
    For i = 0 To numperms-1
      sOut + Str(PeekL(*plong+(i*4)))
      If i < numperms-1: sOut+",": EndIf
    Next
    Debug sOut
  Until gsl_permutation_next(*p) <> #GSL_SUCCESS
  
  gsl_permutation_free(*p)
EndIf


;Output:
;2,5,8
;2,8,5
;5,2,8
;5,8,2
;8,2,5
;8,5,2

[IdeasVacuum]

Wow, you struck oil there Keya!

It would be terrific to have this lib incorporated into PB......

[jack]

the problem is translating all the headers, about 266.

[Keya]

When you read through the helpfile the first impression you might get is that all the functions expect Doubles, but that's just the default version of a function ... each function has a different version for different datatypes, exported with a different name, so there's no need to convert arrays to different datatypes first, phew! (although yes they still all tend to return a Double)
These are:

Code: Select all

      gsl_foo_fn               double
      gsl_foo_long_double_fn   long double
      gsl_foo_float_fn         float
      gsl_foo_long_fn          long
      gsl_foo_ulong_fn         unsigned long
      gsl_foo_int_fn           int
      gsl_foo_uint_fn          unsigned int
      gsl_foo_short_fn         short
      gsl_foo_ushort_fn        unsigned short
      gsl_foo_char_fn          char
      gsl_foo_uchar_fn         unsigned char

For example the gsl_stats_mean() function calculates the Mean of a Double array, but in my example i just have a single array of signed Bytes (PB's .b), not Doubles, so instead i call the gsl_stats_char_mean() function (there's also a uchar version for unsigned Ascii .a bytes):

Code: Select all

PrototypeC.d proto_gsl_stats_char_mean (*pdata, stride.i, nsize.i)

If OpenLibrary(0, "libgsl.dll") 
  gsl_stats_char_mean.proto_gsl_stats_char_mean  = GetFunction(0, "gsl_stats_char_mean")   
  
  Dim B.b(3)
  B(0)=80
  B(1)=88
  B(2)=82
  B(3)=85
  
  Debug("Mean=" + StrD( gsl_stats_char_mean(@B(), 1, 4)) )  
EndIf

There are several other functions that fit the same prototype... do a simple string replace of "mean" with either of the following: average, absdev, kurtosis, skew (there are many more but im still just scratching the surface myself!)

also it seems cdecl PrototypeC is required, as the above demo will crash with regular Prototype

also, Stride and Size i find a little bit confusing to begin with but i believe they always refer to the number of elements (so not bytes, unless using char/uchar functions)

[jack]

hi Keya, just in case someone needs/wants the windows dll's I compiled and bundle them here gsl-2.1 32-bit and gsl-2.1 64-bit

[davido]

@Keya
@jack,

Thank you, very much.

[jack]

I am trying to parse the gsl single page html doc file found here https://www.gnu.org/software/gsl/manual/
and generate prototypes, I think it's easier than translating the numerous headers, it may be next week before I am finished.
I will be happy if someone wants to beat me to it.

[Keya]

excellent jack, thanks! the win32 port didnt seem to have a 64bit version, so your build is a relief

[edit] Oooh, also i just noticed that the Win32 port that i was using is still only v1.8 (2006)!? whereas your builds are the latest v2.1 (Nov 2015), nice! what did you use to compile them and were there many challenges?

[jack]

you are welcome Keya
I used MinGW, not too challenging, the 64-bit version compiled OK with msys2 ( https://msys2.github.io ) but for some reason the 32-bit version would not pass all the tests so had to use an older version of MinGW.

[Lunasole]

Serious library Added it to my list of libs, maybe it might be useful to optimize calculations in something like game, thanks

[Keya]

maybe it might be useful to optimize calculations

or in my case confirm my calculations are wrong

[jack]

Keya, this library has many functions that return a udt for example gsl_complex, how to handle that?
since PB does not allow udt as return value.

[Keya]

it should just be returning a pointer to the structure though, i assume/hope! so just a .i or *ptr.StructureName should do the trick? im looking for example now

so for the function "gsl_complex_add", which returns a "gsl_complex" structure:

Code: Select all

Function: gsl_complex gsl_complex_add (gsl_complex a, gsl_complex b)

i would assume it's something like this?

Code: Select all

Structure gsl_complex
  ;whatever it is
EndStructure
 
ImportC "gsllib"
  gsl_complex_add.i (*a.gsl_complex, *b.gsl_complex) As "gsl_complex_add"
EndImport

i dont know if gsl_complex_add.i can be changed to be a *ptr to gsl_complex like the parameters are but the end result is the same when calling it:

Code: Select all

Define *a.gsl_complex, *b.gsl_complex, *ret.gsl_complex
*ret = gsl_complex_add(*a,*b)

but don't ask me what complex addition is .... numbers greater than 9 maybe?

I tried this, but its not working... perhaps it is actually passing full structures not pointers? i think this is where in VB6 it would by ByVal Structure

Code: Select all

;typedef struct
;  {
;    double dat[2];
;  }
;gsl_complex;
Structure gsl_complex
  dat.d[2]
EndStructure


Prototype.i proto_gsl_complex_add (*a.gsl_complex, *b.gsl_complex)

Define a.gsl_complex, b.gsl_complex, *ret.gsl_complex

If OpenLibrary(0, "libgsl.dll") 
  gsl_complex_add.proto_gsl_complex_add  = GetFunction(0, "gsl_complex_add")   
  
  a\dat[0] = 3
  b\dat[0] = 5
  
  *ret = gsl_complex_add(@a,@b)
  
  Debug("Result0 = " + StrD(*ret\dat[0]))
  Debug("Result1 = " + StrD(*ret\dat[1]))
EndIf 

[Dude]

It is free software under the GNU General Public License.

Which means if you use it, your app must also be open-source for anyone to download and modify. Otherwise, you can't legally use it.

From Wikipedia: "The GPL is a copyleft license, which means that derived works can only be distributed under the same license terms."