Bill Allombert on Fri, 18 Sep 2020 16:39:51 +0200 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: Fortran compatibility |
On Wed, Sep 16, 2020 at 10:37:41PM +0000, Brereton, Ashley wrote: > Many thanks for your suggestions and your quick responses! > > > I will have a go at using gp2C! > > > Concerning X, this may or not be a problem. > > > X in my case is something like: > > > X=1e36 x log(10000) > > > I hoped I could handle the multiplications inside Pari somehow. Please find the following fortran90 program (you need to link it with -lpari) that should do what you want more or less %gfortran appelc.f90 -Wall -O3 -lpari -o appelc %./appelc log(10000)%(2*Pi) = 2.927155065 Cheers, Bill.
module PARI use ISO_C_BINDING, only : C_LONG, C_DOUBLE, C_PTR interface subroutine pari_init(parisize, maxprime) bind(C,name='pari_init') import C_LONG integer(kind=C_LONG), VALUE :: parisize integer(kind=C_LONG), VALUE :: maxprime end subroutine pari_init ! subroutine pari_close() bind(C,name='pari_close') end subroutine pari_close ! type(C_PTR) function dbltor( r ) bind(C,name='dbltor') import C_DOUBLE, C_PTR real(kind=C_DOUBLE), VALUE :: r end function dbltor ! real(kind=C_DOUBLE) function rtodbl( x ) bind(C,name='rtodbl') import C_DOUBLE, C_PTR type(C_PTR), VALUE :: x end function rtodbl ! type(C_PTR) function gsqr( x ) bind(C,name='gsqr') import C_PTR type(C_PTR), VALUE :: x end function gsqr ! type(C_PTR) function gmul( x , y) bind(C,name='gmul') import C_PTR type(C_PTR), VALUE :: x type(C_PTR), VALUE :: y end function gmul ! type(C_PTR) function gprec( x , d) bind(C,name='gprec') import C_PTR, C_LONG type(C_PTR), VALUE :: x integer(kind=C_LONG), VALUE :: d end function gprec ! type(C_PTR) function gmod( x , y) bind(C,name='gmod') import C_PTR type(C_PTR), VALUE :: x type(C_PTR), VALUE :: y end function gmod ! type(C_PTR) function glog( x , prec) bind(C,name='glog') import C_PTR, C_LONG type(C_PTR), VALUE :: x integer(kind=C_LONG), VALUE :: prec end function glog ! type(C_PTR) function stoi(x) bind(C,name='stoi') import C_PTR, C_LONG integer(kind=C_LONG), VALUE :: x end function stoi ! integer(kind=C_LONG) function itos(x) bind(C,name='itos') import C_PTR, C_LONG type(C_PTR), VALUE :: x end function itos ! type(C_PTR) function Pi2n(x, prec) bind(C,name='Pi2n') import C_PTR, C_LONG integer(kind=C_LONG), VALUE :: x integer(kind=C_LONG), VALUE :: prec end function Pi2n end interface end module PARI ! PROGRAM prog use ISO_C_BINDING, only : C_PTR, C_DOUBLE use PARI implicit none real(kind=C_DOUBLE) :: r = 1e36 type(C_PTR) :: p CALL pari_init(10000000_8,2_8) p = gmul(gprec(dbltor(r),1000_8),glog(stoi(10000_8),20_8)) !p = gmod(p, Pi2n(1_8,20_8)) r = rtodbl(p) CALL pari_close() PRINT '(a,f0.9)','1e36*log(10000)%(2*Pi) = ', r END PROGRAM prog