| 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