jack
19-08-2011, 02:34
for the last week or so, while messing with the Lanczos approximation to the Gamma function
I was in need of a multiprecision math library that had besides the basic math operators the function exp
I know there's the gmp library but it lacks the exp function and there's the mpfr lib, but I wanted
something small and simple, so I gave MAPM (http://www.tc.umn.edu/~ringx004/mapm-main.html)a try
I added a dllmain, compiled and did some tests, but could not figure out how to declare in ThinBasic the
MAPM built-in constants, I tried declaring the symbols As Sub() And Then for example
Dim mpi As DWord At Function_GetPtr(sub_name)
but it did Not work.
but the following seems to work
Dim MM_PI As DWord At &h6F4110B8
but is it guaranteed to work?
anyway, here are the function declarations
Type mapm
m_apm_data As Byte Ptr
m_apm_id As Long
m_apm_refcount As Long
m_apm_malloclength As Long
m_apm_datalength As Long
m_apm_exponent As Long
m_apm_sign As Long
End Type
Declare Sub m_apm_init CDECL Lib "mapm.dll" Alias "m_apm_init" ()
Declare Sub m_apm_free CDECL Lib "mapm.dll" Alias "m_apm_free" ( ByVal n As mapm Ptr)
Declare Sub m_apm_free_all_mem CDECL Lib "mapm.dll" Alias "m_apm_free_all_mem" ()
Declare Sub m_apm_trim_mem_usage CDECL Lib "mapm.dll" Alias "m_apm_trim_mem_usage" ()
Declare Sub m_apm_set_string CDECL Lib "mapm.dll" Alias "m_apm_set_string" ( ByVal n As mapm Ptr,ByVal sn As String)
Declare Sub m_apm_set_long CDECL Lib "mapm.dll" Alias "m_apm_set_long" ( ByVal mp As mapm Ptr, ByVal Lng As Long)
Declare Sub m_apm_set_double CDECL Lib "mapm.dll" Alias "m_apm_set_double" ( ByVal mp As mapm Ptr, ByVal dbl As Double)
Declare Sub m_apm_to_string CDECL Lib "mapm.dll" Alias "m_apm_to_string" ( ByRef buffer As Asciiz, ByVal dplaces As Long, ByVal n As mapm Ptr)
Declare Sub m_apm_to_fixpt_string CDECL Lib "mapm.dll" Alias "m_apm_to_fixpt_string" ( ByRef buffer As Asciiz, ByVal dplaces As Long, ByVal n As mapm Ptr)
Declare Sub m_apm_to_fixpt_stringex CDECL Lib "mapm.dll" Alias "m_apm_to_fixpt_stringex" ( ByRef buffer As Asciiz, ByVal dplaces As Long, ByVal n As mapm Ptr, ByVal radix As Byte, ByVal separator_char As Byte, ByVal separator_count As Long)
Declare Function m_apm_to_fixpt_stringexp CDECL Lib "mapm.dll" Alias "m_apm_to_fixpt_stringexp" ( ByVal dplaces As Long, ByVal n As mapm Ptr, ByVal radix As Byte, ByVal separator_char As Byte, ByVal separator_count As Long) As Asciiz Ptr
Declare Sub m_apm_to_integer_string CDECL Lib "mapm.dll" Alias "m_apm_to_integer_string" ( ByRef buffer As Asciiz, ByVal n As mapm Ptr)
Declare Sub m_apm_absolute_value CDECL Lib "mapm.dll" Alias "m_apm_absolute_value" ( ByVal r As mapm Ptr, ByVal n As mapm Ptr)
Declare Sub m_apm_negate CDECL Lib "mapm.dll" Alias "m_apm_negate" ( ByVal r As mapm Ptr, ByVal n As mapm Ptr)
Declare Sub m_apm_copy CDECL Lib "mapm.dll" Alias "m_apm_copy" ( ByVal r As mapm Ptr, ByVal n As mapm Ptr)
Declare Sub m_apm_round CDECL Lib "mapm.dll" Alias "m_apm_round" ( ByVal r As mapm Ptr, ByVal decimal_places As Long, ByVal n As mapm Ptr)
Declare Function m_apm_compare CDECL Lib "mapm.dll" Alias "m_apm_compare" ( ByVal x As mapm Ptr, ByVal y As mapm Ptr) As Long
Declare Function m_apm_sign CDECL Lib "mapm.dll" Alias "m_apm_sign" ( ByVal x As mapm Ptr) As Long
Declare Function m_apm_exponent CDECL Lib "mapm.dll" Alias "m_apm_exponent" ( ByVal x As mapm Ptr) As Long
Declare Function m_apm_significant_digits CDECL Lib "mapm.dll" Alias "m_apm_significant_digits" ( ByVal x As mapm Ptr) As Long
Declare Function m_apm_is_integer CDECL Lib "mapm.dll" Alias "m_apm_is_integer" ( ByVal x As mapm Ptr) As Long
Declare Function m_apm_is_even CDECL Lib "mapm.dll" Alias "m_apm_is_even" ( ByVal x As mapm Ptr) As Long
Declare Function m_apm_is_odd CDECL Lib "mapm.dll" Alias "m_apm_is_odd" ( ByVal x As mapm Ptr) As Long
Declare Sub m_apm_gcd CDECL Lib "mapm.dll" Alias "m_apm_gcd" ( ByVal result As mapm Ptr, ByVal x As mapm Ptr, ByVal y As mapm Ptr)
Declare Sub m_apm_lcm CDECL Lib "mapm.dll" Alias "m_apm_lcm" ( ByVal result As mapm Ptr, ByVal x As mapm Ptr, ByVal y As mapm Ptr)
Declare Sub m_apm_add CDECL Lib "mapm.dll" Alias "m_apm_add" ( ByVal result As mapm Ptr, ByVal x As mapm Ptr, ByVal y As mapm Ptr)
Declare Sub m_apm_subtract CDECL Lib "mapm.dll" Alias "m_apm_subtract" ( ByVal result As mapm Ptr, ByVal x As mapm Ptr, ByVal y As mapm Ptr)
Declare Sub m_apm_multiply CDECL Lib "mapm.dll" Alias "m_apm_multiply" ( ByVal result As mapm Ptr, ByVal x As mapm Ptr, ByVal y As mapm Ptr)
Declare Sub m_apm_divide CDECL Lib "mapm.dll" Alias "m_apm_divide" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr, ByVal y As mapm Ptr)
Declare Sub m_apm_integer_divide CDECL Lib "mapm.dll" Alias "m_apm_integer_divide" ( ByVal result As mapm Ptr, ByVal x As mapm Ptr, ByVal y As mapm Ptr)
Declare Sub m_apm_integer_div_rem CDECL Lib "mapm.dll" Alias "m_apm_integer_div_rem" ( ByVal quot As mapm Ptr, ByVal Rem As mapm Ptr, ByVal x As mapm Ptr, ByVal y As mapm Ptr)
Declare Sub m_apm_reciprocal CDECL Lib "mapm.dll" Alias "m_apm_reciprocal" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_factorial CDECL Lib "mapm.dll" Alias "m_apm_factorial" ( ByVal result As mapm Ptr, ByVal x As mapm Ptr)
Declare Sub m_apm_floor CDECL Lib "mapm.dll" Alias "m_apm_floor" ( ByVal result As mapm Ptr, ByVal x As mapm Ptr)
Declare Sub m_apm_ceil CDECL Lib "mapm.dll" Alias "m_apm_ceil" ( ByVal result As mapm Ptr, ByVal x As mapm Ptr)
Declare Sub m_apm_get_random CDECL Lib "mapm.dll" Alias "m_apm_get_random" ( ByVal result As mapm Ptr)
Declare Sub m_apm_set_random_seed CDECL Lib "mapm.dll" Alias "m_apm_set_random_seed" ( ByVal s As String)
Declare Sub m_apm_sqrt CDECL Lib "mapm.dll" Alias "m_apm_sqrt" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_cbrt CDECL Lib "mapm.dll" Alias "m_apm_cbrt" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_log CDECL Lib "mapm.dll" Alias "m_apm_log" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_log10 CDECL Lib "mapm.dll" Alias "m_apm_log10" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_exp CDECL Lib "mapm.dll" Alias "m_apm_exp" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_pow CDECL Lib "mapm.dll" Alias "m_apm_pow" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr, ByVal y As mapm Ptr)
Declare Sub m_apm_integer_pow CDECL Lib "mapm.dll" Alias "m_apm_integer_pow" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr, ByVal y As Long)
Declare Sub m_apm_integer_pow_nr CDECL Lib "mapm.dll" Alias "m_apm_integer_pow_nr" ( ByVal result As mapm Ptr, ByVal x As mapm Ptr, ByVal y As Long)
Declare Sub m_apm_sin_cos CDECL Lib "mapm.dll" Alias "m_apm_sin_cos" ( ByVal Sin As mapm Ptr, ByVal Cos As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_sin CDECL Lib "mapm.dll" Alias "m_apm_sin" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_cos CDECL Lib "mapm.dll" Alias "m_apm_cos" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_tan CDECL Lib "mapm.dll" Alias "m_apm_tan" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_arcsin CDECL Lib "mapm.dll" Alias "m_apm_arcsin" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_arccos CDECL Lib "mapm.dll" Alias "m_apm_arccos" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_arctan CDECL Lib "mapm.dll" Alias "m_apm_arctan" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_arctan2 CDECL Lib "mapm.dll" Alias "m_apm_arctan2" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr, ByVal y As mapm Ptr)
Declare Sub m_apm_sinh CDECL Lib "mapm.dll" Alias "m_apm_sinh" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_cosh CDECL Lib "mapm.dll" Alias "m_apm_cosh" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_tanh CDECL Lib "mapm.dll" Alias "m_apm_tanh" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_arcsinh CDECL Lib "mapm.dll" Alias "m_apm_arcsinh" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_arccosh CDECL Lib "mapm.dll" Alias "m_apm_arccosh" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_arctanh CDECL Lib "mapm.dll" Alias "m_apm_arctanh" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Dim MM_0_5 As DWord At &h6F4110B0
Dim MM_0_85 As DWord At &h6F4110DC
Dim MM_2_PI As DWord At &h6F4110C0
Dim MM_5x_125R As DWord At &h6F4110E0
Dim MM_5x_256R As DWord At &h6F4110E8
Dim MM_5x_64R As DWord At &h6F4110E4
Dim MM_5x_Eight As DWord At &h6F4110EC
Dim MM_5x_Sixteen As DWord At &h6F4110F0
Dim MM_5x_Twenty As DWord At &h6F4110F4
Dim MM_E As DWord At &h6F4110B4
Dim MM_Five As DWord At &h6F4110A8
Dim MM_Four As DWord At &h6F4110A4
Dim MM_HALF_PI As DWord At &h6F4110BC
Dim MM_LOG_10_BASE_E As DWord At &h6F4110FC
Dim MM_LOG_2_BASE_E As DWord At &h6F411100
Dim MM_LOG_3_BASE_E As DWord At &h6F411104
Dim MM_LOG_E_BASE_10 As DWord At &h6F4110F8
Dim MM_One As DWord At &h6F411098
Dim MM_PI As DWord At &h6F4110B8
Dim MM_Ten As DWord At &h6F4110AC
Dim MM_Three As DWord At &h6F4110A0
Dim MM_Two As DWord At &h6F41109C
Dim MM_Zero As DWord At &h6F411094
Sub gauss ( ByVal a() As DWord, ByVal y() As DWord, ByVal coef() As DWord, ByVal ncol As Long, ByVal prec As Long, ByRef error_code As Long)
' matrix solution by Gaussian Elimination
Dim As DWord b(ncol, ncol), w(ncol) ' work array, ncol long
Dim As Long i,j,i1,k,l,n
Dim As DWord hold, sm, t, t2, t3, ab, big
'Const TRUE = -1
'Const FALSE = Not TRUE
error_code=FALSE
n=ncol
hold=m_apm_init()
sm=m_apm_init()
t=m_apm_init()
t2=m_apm_init()
t3=m_apm_init()
ab=m_apm_init()
big=m_apm_init()
For i=1 To n
' copy to work arrays
For j=1 To n
b(i,j)=m_apm_init()
m_apm_copy(b(i,j),a(i,j))
Next j
w(i)=m_apm_init()
m_apm_copy(w(i),y(i))
Next
For i=1 To n-1
m_apm_absolute_value(big,b(i,i))
l=i
i1=i+1
For j=i1 To n
' search for largest element
m_apm_absolute_value(ab,b(j,i))
If m_apm_compare(ab,big)=1 Then
m_apm_copy(big,ab)
l=j
End If
Next
If m_apm_compare(big,MM_Zero)=0 Then
error_code= TRUE
Else
If l<>i Then
' interchange rows to put largest element on diagonal
For j=1 To n
m_apm_copy(hold,b(l,j))
m_apm_copy(b(l,j),b(i,j))
m_apm_copy(b(i,j),hold)
Next
m_apm_copy(hold,w(l))
m_apm_copy(w(l),w(i))
m_apm_copy(w(i),hold)
End If
For j=i1 To n
m_apm_divide( t, prec, b(j,i),b(i,i))
For k=i1 To n
m_apm_multiply(t2,t,b(i,k))
m_apm_copy(t3,b(j,k))
m_apm_subtract(b(j,k),t3,t2)
Next
m_apm_multiply(t2,t,w(i))
m_apm_copy(t3,w(j))
m_apm_subtract(w(j),t3,t2)
Next ' j-loop
End If ' if big
Next ' i-loop
If m_apm_compare(b(n,n),MM_Zero)=0 Then
error_code=TRUE
Else
m_apm_divide(coef(n),prec,w(n),b(n,n))
i=n-1
' back substitution
Do
m_apm_copy(Sm,MM_Zero)
For j=i+1 To n
m_apm_multiply(t2,b(i,j),coef(j))
m_apm_copy(t3,sm)
m_apm_add(Sm,t3,t2)
Next
m_apm_subtract(t2,w(i),sm)
m_apm_divide(coef(i),prec,t2,b(i,i))
i=i-1
Loop Until i=0
End If ' if b(n,n)=0
If error_code Then Print "ERROR: Matrix is singular"
For i=1 To n
For j=1 To n
m_apm_free(b(i,j))
Next
m_apm_free(w(i))
Next
m_apm_free(hold)
m_apm_free(sm)
m_apm_free(t)
m_apm_free(t2)
m_apm_free(t3)
m_apm_free(ab)
m_apm_free(big)
End Sub ' GAUSS
Uses "console"
Dim i, j As Long
Dim outb As Asciiz * 2000
Dim s As String
Dim As DWord sp
'Const n = 14
Const n = 7
Const prec = 56
Const prec_internal = 3*prec+n\2
Dim As Long code
Dim A(n, n) As DWord
Dim b(n) As DWord
Dim x(n) As DWord
Dim As DWord g, tmp1, tmp2, tmp3, tmp4, tmp5, sqr2pi
For i=1 To n
For j=1 To n
a(i,j)=m_apm_init()
Next
b(i)=m_apm_init()
x(i)=m_apm_init()
Next
g=m_apm_init()
tmp1=m_apm_init()
tmp2=m_apm_init()
tmp3=m_apm_init()
tmp4=m_apm_init()
tmp5=m_apm_init()
sqr2pi=m_apm_init()
'm_apm_set_long(g,607)
'm_apm_set_long(tmp1,128)
'm_apm_divide(g,prec_internal,g,tmp1)
m_apm_set_long(g,5)
m_apm_sqrt(sqr2pi, prec_internal, MM_2_PI)
For i=1 To n
m_apm_set_long(a(i,1), 1)
For j=2 To n
m_apm_set_long(tmp1, i+j)
m_apm_subtract(tmp2, tmp1, MM_Two)
m_apm_divide(a(i,j),prec_internal,MM_One,tmp2)
Next
Next
For i=1 To n
m_apm_set_long(tmp1, i)
m_apm_add(tmp2, tmp1, g)
m_apm_subtract(tmp1, tmp2, MM_0_5)
m_apm_set_long(tmp3, i)
m_apm_subtract(tmp2, tmp3, MM_0_5)
m_apm_pow (tmp3, prec_internal, tmp1, tmp2)
m_apm_set_long(tmp4, i)
m_apm_add(tmp5, tmp4, g)
m_apm_subtract(tmp4, tmp5, MM_0_5)
m_apm_exp (tmp5, prec_internal, tmp4)
m_apm_divide(tmp2,prec_internal, tmp3, tmp5)
m_apm_multiply(tmp3, tmp2, sqr2pi)
m_apm_set_long(tmp1, i-1)
m_apm_factorial ( tmp5, tmp1)
m_apm_divide(b(i),prec_internal, tmp5, tmp3)
' sp = m_apm_to_fixpt_stringexp ( prec, b(i), Asc("."), Asc(","), 5)
' Console_WriteLine Peek$(Asciiz, sp)
Next
gauss(a, b, x, n, prec_internal, code)
For i=1 To n
Reset outb
' m_apm_to_string(outb, prec, x(i))
' If LEFT$(outb,1) <>"-" Then outb=" "&outb
' outb="C("&Using$("##", i)&") = "&outb
' Console_WriteLine outb
sp = m_apm_to_fixpt_stringexp ( prec, x(i), Asc("."), Asc(","), 5)
s = Peek$(Asciiz, sp)
If LEFT$(s,1) <>"-" Then s=" "&s
s="C("&Using$("##", i)&") = "&s
Console_WriteLine s
Next
For i=1 To n
For j=1 To n
m_apm_free(a(i,j))
Next
m_apm_free(b(i))
m_apm_free(x(i))
Next
m_apm_free(g)
m_apm_free(tmp1)
m_apm_free(tmp2)
m_apm_free(tmp3)
m_apm_free(tmp4)
m_apm_free(tmp5)
m_apm_free(sqr2pi)
Console_WriteLine "----------------------------------------"
Console_WriteLine "All done. Press any key to finish"
Console_WaitKey
and here the comple package with the dllmain
I was in need of a multiprecision math library that had besides the basic math operators the function exp
I know there's the gmp library but it lacks the exp function and there's the mpfr lib, but I wanted
something small and simple, so I gave MAPM (http://www.tc.umn.edu/~ringx004/mapm-main.html)a try
I added a dllmain, compiled and did some tests, but could not figure out how to declare in ThinBasic the
MAPM built-in constants, I tried declaring the symbols As Sub() And Then for example
Dim mpi As DWord At Function_GetPtr(sub_name)
but it did Not work.
but the following seems to work
Dim MM_PI As DWord At &h6F4110B8
but is it guaranteed to work?
anyway, here are the function declarations
Type mapm
m_apm_data As Byte Ptr
m_apm_id As Long
m_apm_refcount As Long
m_apm_malloclength As Long
m_apm_datalength As Long
m_apm_exponent As Long
m_apm_sign As Long
End Type
Declare Sub m_apm_init CDECL Lib "mapm.dll" Alias "m_apm_init" ()
Declare Sub m_apm_free CDECL Lib "mapm.dll" Alias "m_apm_free" ( ByVal n As mapm Ptr)
Declare Sub m_apm_free_all_mem CDECL Lib "mapm.dll" Alias "m_apm_free_all_mem" ()
Declare Sub m_apm_trim_mem_usage CDECL Lib "mapm.dll" Alias "m_apm_trim_mem_usage" ()
Declare Sub m_apm_set_string CDECL Lib "mapm.dll" Alias "m_apm_set_string" ( ByVal n As mapm Ptr,ByVal sn As String)
Declare Sub m_apm_set_long CDECL Lib "mapm.dll" Alias "m_apm_set_long" ( ByVal mp As mapm Ptr, ByVal Lng As Long)
Declare Sub m_apm_set_double CDECL Lib "mapm.dll" Alias "m_apm_set_double" ( ByVal mp As mapm Ptr, ByVal dbl As Double)
Declare Sub m_apm_to_string CDECL Lib "mapm.dll" Alias "m_apm_to_string" ( ByRef buffer As Asciiz, ByVal dplaces As Long, ByVal n As mapm Ptr)
Declare Sub m_apm_to_fixpt_string CDECL Lib "mapm.dll" Alias "m_apm_to_fixpt_string" ( ByRef buffer As Asciiz, ByVal dplaces As Long, ByVal n As mapm Ptr)
Declare Sub m_apm_to_fixpt_stringex CDECL Lib "mapm.dll" Alias "m_apm_to_fixpt_stringex" ( ByRef buffer As Asciiz, ByVal dplaces As Long, ByVal n As mapm Ptr, ByVal radix As Byte, ByVal separator_char As Byte, ByVal separator_count As Long)
Declare Function m_apm_to_fixpt_stringexp CDECL Lib "mapm.dll" Alias "m_apm_to_fixpt_stringexp" ( ByVal dplaces As Long, ByVal n As mapm Ptr, ByVal radix As Byte, ByVal separator_char As Byte, ByVal separator_count As Long) As Asciiz Ptr
Declare Sub m_apm_to_integer_string CDECL Lib "mapm.dll" Alias "m_apm_to_integer_string" ( ByRef buffer As Asciiz, ByVal n As mapm Ptr)
Declare Sub m_apm_absolute_value CDECL Lib "mapm.dll" Alias "m_apm_absolute_value" ( ByVal r As mapm Ptr, ByVal n As mapm Ptr)
Declare Sub m_apm_negate CDECL Lib "mapm.dll" Alias "m_apm_negate" ( ByVal r As mapm Ptr, ByVal n As mapm Ptr)
Declare Sub m_apm_copy CDECL Lib "mapm.dll" Alias "m_apm_copy" ( ByVal r As mapm Ptr, ByVal n As mapm Ptr)
Declare Sub m_apm_round CDECL Lib "mapm.dll" Alias "m_apm_round" ( ByVal r As mapm Ptr, ByVal decimal_places As Long, ByVal n As mapm Ptr)
Declare Function m_apm_compare CDECL Lib "mapm.dll" Alias "m_apm_compare" ( ByVal x As mapm Ptr, ByVal y As mapm Ptr) As Long
Declare Function m_apm_sign CDECL Lib "mapm.dll" Alias "m_apm_sign" ( ByVal x As mapm Ptr) As Long
Declare Function m_apm_exponent CDECL Lib "mapm.dll" Alias "m_apm_exponent" ( ByVal x As mapm Ptr) As Long
Declare Function m_apm_significant_digits CDECL Lib "mapm.dll" Alias "m_apm_significant_digits" ( ByVal x As mapm Ptr) As Long
Declare Function m_apm_is_integer CDECL Lib "mapm.dll" Alias "m_apm_is_integer" ( ByVal x As mapm Ptr) As Long
Declare Function m_apm_is_even CDECL Lib "mapm.dll" Alias "m_apm_is_even" ( ByVal x As mapm Ptr) As Long
Declare Function m_apm_is_odd CDECL Lib "mapm.dll" Alias "m_apm_is_odd" ( ByVal x As mapm Ptr) As Long
Declare Sub m_apm_gcd CDECL Lib "mapm.dll" Alias "m_apm_gcd" ( ByVal result As mapm Ptr, ByVal x As mapm Ptr, ByVal y As mapm Ptr)
Declare Sub m_apm_lcm CDECL Lib "mapm.dll" Alias "m_apm_lcm" ( ByVal result As mapm Ptr, ByVal x As mapm Ptr, ByVal y As mapm Ptr)
Declare Sub m_apm_add CDECL Lib "mapm.dll" Alias "m_apm_add" ( ByVal result As mapm Ptr, ByVal x As mapm Ptr, ByVal y As mapm Ptr)
Declare Sub m_apm_subtract CDECL Lib "mapm.dll" Alias "m_apm_subtract" ( ByVal result As mapm Ptr, ByVal x As mapm Ptr, ByVal y As mapm Ptr)
Declare Sub m_apm_multiply CDECL Lib "mapm.dll" Alias "m_apm_multiply" ( ByVal result As mapm Ptr, ByVal x As mapm Ptr, ByVal y As mapm Ptr)
Declare Sub m_apm_divide CDECL Lib "mapm.dll" Alias "m_apm_divide" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr, ByVal y As mapm Ptr)
Declare Sub m_apm_integer_divide CDECL Lib "mapm.dll" Alias "m_apm_integer_divide" ( ByVal result As mapm Ptr, ByVal x As mapm Ptr, ByVal y As mapm Ptr)
Declare Sub m_apm_integer_div_rem CDECL Lib "mapm.dll" Alias "m_apm_integer_div_rem" ( ByVal quot As mapm Ptr, ByVal Rem As mapm Ptr, ByVal x As mapm Ptr, ByVal y As mapm Ptr)
Declare Sub m_apm_reciprocal CDECL Lib "mapm.dll" Alias "m_apm_reciprocal" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_factorial CDECL Lib "mapm.dll" Alias "m_apm_factorial" ( ByVal result As mapm Ptr, ByVal x As mapm Ptr)
Declare Sub m_apm_floor CDECL Lib "mapm.dll" Alias "m_apm_floor" ( ByVal result As mapm Ptr, ByVal x As mapm Ptr)
Declare Sub m_apm_ceil CDECL Lib "mapm.dll" Alias "m_apm_ceil" ( ByVal result As mapm Ptr, ByVal x As mapm Ptr)
Declare Sub m_apm_get_random CDECL Lib "mapm.dll" Alias "m_apm_get_random" ( ByVal result As mapm Ptr)
Declare Sub m_apm_set_random_seed CDECL Lib "mapm.dll" Alias "m_apm_set_random_seed" ( ByVal s As String)
Declare Sub m_apm_sqrt CDECL Lib "mapm.dll" Alias "m_apm_sqrt" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_cbrt CDECL Lib "mapm.dll" Alias "m_apm_cbrt" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_log CDECL Lib "mapm.dll" Alias "m_apm_log" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_log10 CDECL Lib "mapm.dll" Alias "m_apm_log10" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_exp CDECL Lib "mapm.dll" Alias "m_apm_exp" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_pow CDECL Lib "mapm.dll" Alias "m_apm_pow" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr, ByVal y As mapm Ptr)
Declare Sub m_apm_integer_pow CDECL Lib "mapm.dll" Alias "m_apm_integer_pow" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr, ByVal y As Long)
Declare Sub m_apm_integer_pow_nr CDECL Lib "mapm.dll" Alias "m_apm_integer_pow_nr" ( ByVal result As mapm Ptr, ByVal x As mapm Ptr, ByVal y As Long)
Declare Sub m_apm_sin_cos CDECL Lib "mapm.dll" Alias "m_apm_sin_cos" ( ByVal Sin As mapm Ptr, ByVal Cos As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_sin CDECL Lib "mapm.dll" Alias "m_apm_sin" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_cos CDECL Lib "mapm.dll" Alias "m_apm_cos" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_tan CDECL Lib "mapm.dll" Alias "m_apm_tan" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_arcsin CDECL Lib "mapm.dll" Alias "m_apm_arcsin" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_arccos CDECL Lib "mapm.dll" Alias "m_apm_arccos" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_arctan CDECL Lib "mapm.dll" Alias "m_apm_arctan" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_arctan2 CDECL Lib "mapm.dll" Alias "m_apm_arctan2" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr, ByVal y As mapm Ptr)
Declare Sub m_apm_sinh CDECL Lib "mapm.dll" Alias "m_apm_sinh" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_cosh CDECL Lib "mapm.dll" Alias "m_apm_cosh" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_tanh CDECL Lib "mapm.dll" Alias "m_apm_tanh" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_arcsinh CDECL Lib "mapm.dll" Alias "m_apm_arcsinh" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_arccosh CDECL Lib "mapm.dll" Alias "m_apm_arccosh" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Declare Sub m_apm_arctanh CDECL Lib "mapm.dll" Alias "m_apm_arctanh" ( ByVal result As mapm Ptr, ByVal dplaces As Long, ByVal x As mapm Ptr)
Dim MM_0_5 As DWord At &h6F4110B0
Dim MM_0_85 As DWord At &h6F4110DC
Dim MM_2_PI As DWord At &h6F4110C0
Dim MM_5x_125R As DWord At &h6F4110E0
Dim MM_5x_256R As DWord At &h6F4110E8
Dim MM_5x_64R As DWord At &h6F4110E4
Dim MM_5x_Eight As DWord At &h6F4110EC
Dim MM_5x_Sixteen As DWord At &h6F4110F0
Dim MM_5x_Twenty As DWord At &h6F4110F4
Dim MM_E As DWord At &h6F4110B4
Dim MM_Five As DWord At &h6F4110A8
Dim MM_Four As DWord At &h6F4110A4
Dim MM_HALF_PI As DWord At &h6F4110BC
Dim MM_LOG_10_BASE_E As DWord At &h6F4110FC
Dim MM_LOG_2_BASE_E As DWord At &h6F411100
Dim MM_LOG_3_BASE_E As DWord At &h6F411104
Dim MM_LOG_E_BASE_10 As DWord At &h6F4110F8
Dim MM_One As DWord At &h6F411098
Dim MM_PI As DWord At &h6F4110B8
Dim MM_Ten As DWord At &h6F4110AC
Dim MM_Three As DWord At &h6F4110A0
Dim MM_Two As DWord At &h6F41109C
Dim MM_Zero As DWord At &h6F411094
Sub gauss ( ByVal a() As DWord, ByVal y() As DWord, ByVal coef() As DWord, ByVal ncol As Long, ByVal prec As Long, ByRef error_code As Long)
' matrix solution by Gaussian Elimination
Dim As DWord b(ncol, ncol), w(ncol) ' work array, ncol long
Dim As Long i,j,i1,k,l,n
Dim As DWord hold, sm, t, t2, t3, ab, big
'Const TRUE = -1
'Const FALSE = Not TRUE
error_code=FALSE
n=ncol
hold=m_apm_init()
sm=m_apm_init()
t=m_apm_init()
t2=m_apm_init()
t3=m_apm_init()
ab=m_apm_init()
big=m_apm_init()
For i=1 To n
' copy to work arrays
For j=1 To n
b(i,j)=m_apm_init()
m_apm_copy(b(i,j),a(i,j))
Next j
w(i)=m_apm_init()
m_apm_copy(w(i),y(i))
Next
For i=1 To n-1
m_apm_absolute_value(big,b(i,i))
l=i
i1=i+1
For j=i1 To n
' search for largest element
m_apm_absolute_value(ab,b(j,i))
If m_apm_compare(ab,big)=1 Then
m_apm_copy(big,ab)
l=j
End If
Next
If m_apm_compare(big,MM_Zero)=0 Then
error_code= TRUE
Else
If l<>i Then
' interchange rows to put largest element on diagonal
For j=1 To n
m_apm_copy(hold,b(l,j))
m_apm_copy(b(l,j),b(i,j))
m_apm_copy(b(i,j),hold)
Next
m_apm_copy(hold,w(l))
m_apm_copy(w(l),w(i))
m_apm_copy(w(i),hold)
End If
For j=i1 To n
m_apm_divide( t, prec, b(j,i),b(i,i))
For k=i1 To n
m_apm_multiply(t2,t,b(i,k))
m_apm_copy(t3,b(j,k))
m_apm_subtract(b(j,k),t3,t2)
Next
m_apm_multiply(t2,t,w(i))
m_apm_copy(t3,w(j))
m_apm_subtract(w(j),t3,t2)
Next ' j-loop
End If ' if big
Next ' i-loop
If m_apm_compare(b(n,n),MM_Zero)=0 Then
error_code=TRUE
Else
m_apm_divide(coef(n),prec,w(n),b(n,n))
i=n-1
' back substitution
Do
m_apm_copy(Sm,MM_Zero)
For j=i+1 To n
m_apm_multiply(t2,b(i,j),coef(j))
m_apm_copy(t3,sm)
m_apm_add(Sm,t3,t2)
Next
m_apm_subtract(t2,w(i),sm)
m_apm_divide(coef(i),prec,t2,b(i,i))
i=i-1
Loop Until i=0
End If ' if b(n,n)=0
If error_code Then Print "ERROR: Matrix is singular"
For i=1 To n
For j=1 To n
m_apm_free(b(i,j))
Next
m_apm_free(w(i))
Next
m_apm_free(hold)
m_apm_free(sm)
m_apm_free(t)
m_apm_free(t2)
m_apm_free(t3)
m_apm_free(ab)
m_apm_free(big)
End Sub ' GAUSS
Uses "console"
Dim i, j As Long
Dim outb As Asciiz * 2000
Dim s As String
Dim As DWord sp
'Const n = 14
Const n = 7
Const prec = 56
Const prec_internal = 3*prec+n\2
Dim As Long code
Dim A(n, n) As DWord
Dim b(n) As DWord
Dim x(n) As DWord
Dim As DWord g, tmp1, tmp2, tmp3, tmp4, tmp5, sqr2pi
For i=1 To n
For j=1 To n
a(i,j)=m_apm_init()
Next
b(i)=m_apm_init()
x(i)=m_apm_init()
Next
g=m_apm_init()
tmp1=m_apm_init()
tmp2=m_apm_init()
tmp3=m_apm_init()
tmp4=m_apm_init()
tmp5=m_apm_init()
sqr2pi=m_apm_init()
'm_apm_set_long(g,607)
'm_apm_set_long(tmp1,128)
'm_apm_divide(g,prec_internal,g,tmp1)
m_apm_set_long(g,5)
m_apm_sqrt(sqr2pi, prec_internal, MM_2_PI)
For i=1 To n
m_apm_set_long(a(i,1), 1)
For j=2 To n
m_apm_set_long(tmp1, i+j)
m_apm_subtract(tmp2, tmp1, MM_Two)
m_apm_divide(a(i,j),prec_internal,MM_One,tmp2)
Next
Next
For i=1 To n
m_apm_set_long(tmp1, i)
m_apm_add(tmp2, tmp1, g)
m_apm_subtract(tmp1, tmp2, MM_0_5)
m_apm_set_long(tmp3, i)
m_apm_subtract(tmp2, tmp3, MM_0_5)
m_apm_pow (tmp3, prec_internal, tmp1, tmp2)
m_apm_set_long(tmp4, i)
m_apm_add(tmp5, tmp4, g)
m_apm_subtract(tmp4, tmp5, MM_0_5)
m_apm_exp (tmp5, prec_internal, tmp4)
m_apm_divide(tmp2,prec_internal, tmp3, tmp5)
m_apm_multiply(tmp3, tmp2, sqr2pi)
m_apm_set_long(tmp1, i-1)
m_apm_factorial ( tmp5, tmp1)
m_apm_divide(b(i),prec_internal, tmp5, tmp3)
' sp = m_apm_to_fixpt_stringexp ( prec, b(i), Asc("."), Asc(","), 5)
' Console_WriteLine Peek$(Asciiz, sp)
Next
gauss(a, b, x, n, prec_internal, code)
For i=1 To n
Reset outb
' m_apm_to_string(outb, prec, x(i))
' If LEFT$(outb,1) <>"-" Then outb=" "&outb
' outb="C("&Using$("##", i)&") = "&outb
' Console_WriteLine outb
sp = m_apm_to_fixpt_stringexp ( prec, x(i), Asc("."), Asc(","), 5)
s = Peek$(Asciiz, sp)
If LEFT$(s,1) <>"-" Then s=" "&s
s="C("&Using$("##", i)&") = "&s
Console_WriteLine s
Next
For i=1 To n
For j=1 To n
m_apm_free(a(i,j))
Next
m_apm_free(b(i))
m_apm_free(x(i))
Next
m_apm_free(g)
m_apm_free(tmp1)
m_apm_free(tmp2)
m_apm_free(tmp3)
m_apm_free(tmp4)
m_apm_free(tmp5)
m_apm_free(sqr2pi)
Console_WriteLine "----------------------------------------"
Console_WriteLine "All done. Press any key to finish"
Console_WaitKey
and here the comple package with the dllmain