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