GNU Scientific Library in Assembler

[jj2007]

that suggests that there is a library init routine someplace
perhaps there are other items that require initialization, as well

Couldn't find any, Dave. But anyway, the next edition of MB will take care of that.

[donkey]

that suggests that there is a library init routine someplace
perhaps there are other items that require initialization, as well

Couldn't find any, Dave. But anyway, the next edition of MB will take care of that.

Maybe gsl_ieee_env_setup

http://www.gnu.org/software/gsl/manual/html_node/Setting-up-your-IEEE-environment.html

[MichaelW]

Looking in the GSL module definition file for exported functions that contain "init", I don't see any that look like some sort of global initialization, but I do see what look like initialization for specific groups of functions, for example gsl_interp_init and gsl_monte_vegas_init. And in my quick look at the source, virtually all of the floating-point variables/parameters/returns are doubles and the configuration is done via defines. The Cygwin library also apparently does not, at least by default, control the FPU precision, and I'm guessing that for both libraries it's a portability thing.

Code Select Expand


;====================================================================
    .486                                ; create 32 bit code
    .model flat, stdcall                ; 32 bit memory model
    option casemap :none                ; case sensitive

    include \masm32\include\windows.inc
    include \masm32\include\masm32.inc
    include \masm32\include\gdi32.inc
    include \masm32\include\user32.inc
    include \masm32\include\kernel32.inc
    include \masm32\include\Comctl32.inc
    include \masm32\include\comdlg32.inc
    include \masm32\include\shell32.inc
    include \masm32\include\oleaut32.inc
    include \masm32\include\msvcrt.inc

    include \masm32\macros\macros.asm

    includelib \masm32\lib\masm32.lib
    includelib \masm32\lib\gdi32.lib
    includelib \masm32\lib\user32.lib
    includelib \masm32\lib\kernel32.lib
    includelib \masm32\lib\Comctl32.lib
    includelib \masm32\lib\comdlg32.lib
    includelib \masm32\lib\shell32.lib
    includelib \masm32\lib\oleaut32.lib

    includelib cygwin1.lib
    includelib libgsl.lib

    includelib \masm32\lib\msvcrt.lib
;====================================================================

cygwin_dll_init PROTO C
dll_dllcrt0 PROTO C
printf PROTO C :VARARG

j0               PROTO C :REAL8
gsl_sf_bessel_J0 PROTO C :REAL8

;====================================================================
    .data
        r85 REAL8 5.0
        r8  REAL8 ?
        r10 REAL10 ?
    .code
;====================================================================

FPC_24BITS equ 0000000000b
FPC_53BITS equ 1000000000b
FPC_64BITS equ 1100000000b

ShowPC proc
    invoke printf, cfm$("PC = ")
    push eax
    fstcw [esp]
    pop eax
    and eax, FPC_64BITS
    .IF eax == FPC_24BITS
        invoke printf, cfm$("24 bits\n\n")
    .ELSEIF eax == FPC_53BITS
        invoke printf, cfm$("53 bits\n\n")
    .ELSE
        invoke printf, cfm$("64 bits\n\n")
    .ENDIF
    ret
ShowPC endp

;====================================================================
start:
;====================================================================
    invoke cygwin_dll_init
    invoke dll_dllcrt0

    invoke ShowPC

    invoke j0, r85
    fstp r10
    invoke printf, cfm$("cygwin j0 to r10       : %.19Lf\n\n"), r10

    invoke j0, r85
    fstp r8
    invoke printf, cfm$("cygwin j0 to r8        : %.19f\n\n"), r8

    invoke gsl_sf_bessel_J0, r85
    fstp r10
    invoke printf, cfm$("gsl_sf_bessel_J0 to r10: %.19Lf\n\n"), r10

    invoke gsl_sf_bessel_J0, r85
    fstp r8
    invoke printf, cfm$("gsl_sf_bessel_J0 to r8 : %.19f\n\n"), r8

    finit

    invoke ShowPC

    invoke j0, r85
    fstp r10
    invoke printf, cfm$("cygwin j0 to r10       : %.19Lf\n\n"), r10

    invoke j0, r85
    fstp r8
    invoke printf, cfm$("cygwin j0 to r8        : %.19f\n\n"), r8

    invoke gsl_sf_bessel_J0, r85
    fstp r10
    invoke printf, cfm$("gsl_sf_bessel_J0 to r10: %.19Lf\n\n"), r10

    invoke gsl_sf_bessel_J0, r85
    fstp r8
    invoke printf, cfm$("gsl_sf_bessel_J0 to r8 : %.19f\n\n"), r8

    inkey "Press any key to exit..."
    exit
;====================================================================
end start


Code Select Expand


PC = 53 bits

cygwin j0 to r10       : -0.1775967713143382920

cygwin j0 to r8        : -0.1775967713143382920

gsl_sf_bessel_J0 to r10: -0.1775967713143382642

gsl_sf_bessel_J0 to r8 : -0.1775967713143382642

PC = 64 bits

cygwin j0 to r10       : -0.1775967713143382642

cygwin j0 to r8        : -0.1775967713143382642

gsl_sf_bessel_J0 to r10: -0.1775967713143382920

gsl_sf_bessel_J0 to r8 : -0.1775967713143382920


I have no idea how to explain the opposite effect that the FPU PC has on the two libraries.

[jj2007]

Michael,

Edgar found the FPU init function: gsl_ieee_env_setup:
"When [no]* options are specified using this method, the resulting mode is based on a default setting of the highest available precision (double precision or extended precision, depending on the platform) in round-to-nearest mode, with all exceptions enabled apart from the inexact exception. The inexact  exception is generated whenever rounding occurs"

When my proggie called the Bessel routine, the FPU happened to be set to 64 bits (i.e. full precision) but rounding down instead of nearest. Which triggered the wrong results.

I have fixed that in the meantime, the MasmBasic package of today 23 May calls finit after the initialisation routine. The gsl macro has also become more powerful and accurate, but exotic cases may still choke - I haven't seen one, but with over 1,000 functions there surely will be some problematic cases.

*: Original says "When options are specified" but they obviously mean NO options.

[MichaelW]

New results, this time controlling the RC (which apparently defaults to Nearest):

Code Select Expand


PC = 53 bits

RC = Nearest

cygwin j0 to r10       : -0.1775967713143382920

cygwin j0 to r8        : -0.1775967713143382920

gsl_sf_bessel_J0 to r10: -0.1775967713143382642

gsl_sf_bessel_J0 to r8 : -0.1775967713143382642

PC = 64 bits

cygwin j0 to r10       : -0.1775967713143382642

cygwin j0 to r8        : -0.1775967713143382642

gsl_sf_bessel_J0 to r10: -0.1775967713143382920

gsl_sf_bessel_J0 to r8 : -0.1775967713143382920

RC = Down

cygwin j0 to r10       : -0.1775967713143383198

cygwin j0 to r8        : -0.1775967713143383198

gsl_sf_bessel_J0 to r10: -0.1775967713143383198

gsl_sf_bessel_J0 to r8 : -0.1775967713143383198


[jj2007]

Really strange, that "exchange of results" for rising precision. By the way, the correct value seems to be

-0.177596771314338304347397013075 - according to Wolfram Mathematica

There is another issue to observe: Those gsl functions that return a double leave it on the FPU, in ST(0). So you need to cleanup yourself, either by a fstp MyReal8MemVar^*), or with fstp st. Otherwise strange things happen - the Bessel algo uses all 8 FPU regs!

include \masm32\MasmBasic\MasmBasic.inc   ; Download
WantGoodResult = 1   ; test it with ... = 0
.data
MyReal8 REAL8 5.0

   Init
   gsl   double gsl_sf_bessel_J0 (double)
   gsl_INIT
   gsl_sf_bessel_J0(MyReal8)      ; method A: use the macro
   Print Str$("BesselJ(0)=%Jf\n", ST(0))      ; display what is still on the FPU
   if WantGoodResult
      fstp st      ; the Bessel routine uses all eight FPU regs, so you better clear the result
   endif
   invoke hgsl_sf_bessel_J0, MyReal8      ; method B: invoke hXXX
   Print Str$("BesselJ(0)=%Jf\n", ST(0))
   gsl_EXIT
   Exit
end start

Output:
BesselJ(0)=-0.1775967713143382920
BesselJ(0)=-0.1775967713143382920

^*) Another option is this macro (works for Real8 and Real4 destinations):

Code Select Expand

include \masm32\MasmBasic\MasmBasic.inc

PopReal MACRO dest, arg
  arg
  fstp dest
ENDM

.data
MyData REAL8 5.0, 6.0, 3.2, 1.8, 9.0 ; define a double precision array
MyReal8 REAL8 ?

Init
gsl double gsl_stats_mean (const double data[], size_t stride, size_t n)
gsl_INIT
PopReal MyReal8, gsl_stats_mean(offset MyData, 1, 5)
PrintLine Str$("Mean = %7f", MyReal8)
gsl_EXIT
Exit
end start

[jj2007]

Just in case you are using the gsl_INIT macro inside a procedure:
There was a bug in line 6234 of \Masm32\MasmBasic\MasmBasic.inc - a plain ret instead of the correct retn. Strange things could happen if you used it inside a proc, such as (through the eyes of Olly):

Code Select Expand


  pop ebx
  pop edi
  pop esi
  mov esp, ebp
  pop ebp
  retn 8

The plain ret is not as plain as it looks - it is a actually a macro taking care of adjusting the stack and the regs that were declared in xx proc use esi edi ebx :bg

Version 10 Feb 2012 is ok.