gmp 6.2.1 and mpfr 4.1.0

General FreeBASIC programming questions.
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: gmp 6.2.1 and mpfr 4.1.0

Post by srvaldez »

that's a quirk of FB, see Global Operators at https://www.freebasic.net/wiki/ProPgOperatorOverloading
jdebord
Posts: 547
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

Re: gmp 6.2.1 and mpfr 4.1.0

Post by jdebord »

So, OPERATOR is for the functions already defined in FB (Sin, Cos ...) and FUNCTION for the new ones (Sinh, Cosh ...).

I should have known !

Thanks again !
jdebord
Posts: 547
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

Re: gmp 6.2.1 and mpfr 4.1.0

Post by jdebord »

It seems that the Beta function is lacking, both in mpfr.bi and mpfr-fb.bi

So, I tried to add it.

In mpfr.bi, add the line :

Code: Select all

declare function mpfr_beta(byval as mpfr_ptr, byval as mpfr_srcptr, byval as mpfr_srcptr, byval as mpfr_rnd_t) as long
In mpfr-fb.bi, add the lines :

Code: Select all

function Beta(byref lhs as mpfr, byref rhs as mpfr) as mpfr
	Dim As mpfr result
	mpfr_beta (@result.num, @lhs.num, @rhs.num, mpfr_rounding_mode)
	Return result
end function
Test program :

Code: Select all

dim as mpfr x, y

x = 1.23
y = 4.56

? Beta(x, y)
? Gamma(x) * Gamma(y) / Gamma(x + y)
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: gmp 6.2.1 and mpfr 4.1.0

Post by srvaldez »

Hi jdebord
my only comment is that floating point literals need to be in quotes otherwise the precision will be that of double, so your example should be

Code: Select all

dim as mpfr x, y

x = "1.23"
y = "4.56"

? Beta(x, y)
? Gamma(x) * Gamma(y) / Gamma(x + y)
jdebord
Posts: 547
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

Re: gmp 6.2.1 and mpfr 4.1.0

Post by jdebord »

Number of primes less than x = 10^i for i = 1..25

See : 'https://en.wikipedia.org/wiki/Prime_number_theorem

Exact value and 2 approximations : x / log(x) and Li(x)

It seems that the logarithm integral function Li(x) is absent from MPFR, so I used Eint(log(x)) instead.

Code: Select all

''=============================================================================
#include once "gmp.bi"
#include once "mpfr.bi"

'we need to set mpfr_digits before including the overloded operators
dim shared as long mpfr_digits_precision = 100 '<- set mpfr digits-precision

dim shared as mpfr_rnd_t mpfr_rounding_mode = MPFR_RNDN '<- set the rounding mode
' rounding options are
'  MPFR_RNDN    round to nearest, with ties to even
'  MPFR_RNDZ    round toward zero
'  MPFR_RNDU    round toward +Inf
'  MPFR_RNDD    round toward -Inf
'  MPFR_RNDA    round away from zero
'  MPFR_RNDF    faithful rounding (not implemented yet)
'  MPFR_RNDNA   round to nearest, with ties away from zero (mpfr_round)

#include once "mpfr-fb.bi" ' mpfr overloaded operators
''=============================================================================

data "4", "25", "168", "1229", "9592", "78498", "664579", "5761455"
data "50847534", "455052511", "4118054813", "37607912018", "346065536839"
data "3204941750802", "29844570422669", "279238341033925", "2623557157654233"
data "24739954287740860", "234057667276344607", "2220819602560918840"
data "21127269486018731928", "201467286689315906290", "1925320391606803968923"
data "18435599767349200867866", "176846309399143769411680"

dim as long i
dim as string sp
dim as mpfr n, x, logx, approx1, approx2

print " i        N primes <= 10^i = x          x / log(x)                    Li(x)"  
print string(95, "-")

for i = 1 to 25
  read sp
  n = sp
  x = 10^i
  logx = log(x)
  approx1 = x / logx
  approx2 = Eint(logx)
  print i; tab(10); n; tab(40); int(approx1); tab(70); int(approx2)
next i 

sleep
jdebord
Posts: 547
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

Re: gmp 6.2.1 and mpfr 4.1.0

Post by jdebord »

Compute Pi by the method of Chudnovsky

See : https://en.wikipedia.org/wiki/Chudnovsky_algorithm

Code: Select all

#include once "gmp.bi"
#include once "mpfr.bi"

const Digits = 100

'we need to set mpfr_digits before including the overloded operators
dim shared as long mpfr_digits_precision = Digits '<- set mpfr digits-precision

dim shared as mpfr_rnd_t mpfr_rounding_mode = MPFR_RNDN '<- set the rounding mode
' rounding options are
'  MPFR_RNDN    round to nearest, with ties to even
'  MPFR_RNDZ    round toward zero
'  MPFR_RNDU    round toward +Inf
'  MPFR_RNDD    round toward -Inf
'  MPFR_RNDA    round away from zero
'  MPFR_RNDF    faithful rounding (not implemented yet)
'  MPFR_RNDNA   round to nearest, with ties away from zero (mpfr_round)

#include once "mpfr-fb.bi" ' mpfr overloaded operators
''=============================================================================

' The formula adds 14 digits at each iteration
 
const MaxIter = Digits \ 14 + 1

dim as mpfr C, L, DL, X, FX, M, S
dim as long iter, iter1
    
C   = mpfr(426880) * sqr(mpfr(10005))
L   = mpfr(13591409)
DL  = mpfr(545140134)
X   = 1
FX  = mpfr(-262537412640768000)
S   = L
    
for iter = 1 to MaxIter
  M = facui(6 * iter) / (facui(3 * iter) * facui(iter) ^ 3)
  L = L + DL
  X = X * FX
  S = S + M * L / X
next iter
    
print C / S
print
print pi_const

sleep
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: gmp 6.2.1 and mpfr 4.1.0

Post by srvaldez »

hi jdebord
nice example, feel free to improve the includes and post them on this forum :-)
jdebord
Posts: 547
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

Re: gmp 6.2.1 and mpfr 4.1.0

Post by jdebord »

Thank you srvaldez !

Here is another example : Pi by the AGM method.

Code: Select all

' -------------------------------------------------------------------
' Compute Pi by Arithmetic-Geometric Mean (AGM) 
' https://en.wikipedia.org/wiki/Arithmetic–geometric_mean
' -------------------------------------------------------------------

#include once "gmp.bi"
#include once "mpfr.bi"

const Digits = 100

'we need to set mpfr_digits before including the overloded operators
dim shared as long mpfr_digits_precision = Digits '<- set mpfr digits-precision

dim shared as mpfr_rnd_t mpfr_rounding_mode = MPFR_RNDN '<- set the rounding mode
' rounding options are
'  MPFR_RNDN    round to nearest, with ties to even
'  MPFR_RNDZ    round toward zero
'  MPFR_RNDU    round toward +Inf
'  MPFR_RNDD    round toward -Inf
'  MPFR_RNDA    round away from zero
'  MPFR_RNDF    faithful rounding (not implemented yet)
'  MPFR_RNDNA   round to nearest, with ties away from zero (mpfr_round)

#include once "mpfr-fb.bi" ' mpfr overloaded operators
''=============================================================================

' The formula doubles the number of digits at each iteration
const MaxIter = int(log(Digits) / log(2)) + 1

dim as mpfr a, a1, b, b1, d, s
dim as long iter
    
a = 1
b = 1 / sqr(mpfr(2))
d = 2
s = mpfr(0.5)
  
for iter = 1 to MaxIter
  a1 = (a + b) / 2
  b1 = sqr(a * b)
  s = s - d * (square(a1) - square(b1))
  a = a1
  b = b1
  d = 2 * d
next iter
    
print 2 * square(a) / s
print
print pi_const

sleep
Post Reply