Ramanujan

Some Computational Exercises

What follows are four discoveries by Srinivasa Ramanujan that I've confirmed to quad precision, or 32 decimal places. On a Linux system, compile with gfortran -g -std=f2018 -Wall -Wextra -O2 -o [program] [program].f08 and run with ./[program]. Enjoy!

Golden_ratio.f08

PROGRAM Golden_ratio
  USE, INTRINSIC :: ISO_Fortran_env, qp=>REAL128 !modern QUAD PRECISION
  IMPLICIT NONE
  INTEGER :: n, number
  REAL(qp) :: phi
  WRITE(6, '(a6)', advance='no') "  n = "
  READ *, number
  CALL Iterate(1.0_qp, 0, phi, n, number)
  WRITE(6, '(f36.32)') phi
  WRITE(6, '(f36.32)') (1.0_qp + SQRT(5.0_qp)) / 2.0_qp
CONTAINS
  RECURSIVE SUBROUTINE Iterate(phi_0, n_0, phi, n, number)
    INTEGER :: n_1
    INTEGER, VALUE :: n_0 !pass by value
    INTEGER, INTENT(IN) :: number !pass by reference
    INTEGER, INTENT(OUT) :: n !pass by reference
    REAL(qp) :: phi_1
    REAL(qp), VALUE :: phi_0 !pass by value
    REAL(qp), INTENT(OUT) :: phi !pass by reference
    phi_1 = 1.0_qp + (1.0_qp / phi_0)
    n_1 = n_0 + 1
    IF (n_1 >= number) THEN
      phi = phi_1
      n = n_1
    ELSE
      CALL Iterate(phi_1, n_1, phi, n, number)
    END IF
  END SUBROUTINE Iterate
END PROGRAM Golden_ratio

Pi_Ramanujan.f08

PROGRAM Pi_Ramanujan
  USE, INTRINSIC :: ISO_Fortran_env, li=>INT64 !modern LONG INTEGER
  USE, INTRINSIC :: ISO_Fortran_env, qp=>REAL128 !modern QUAD PRECISION
  IMPLICIT NONE
  INTEGER :: k, number
  INTEGER(li) :: Factorial
  REAL(qp) :: Pi
  EXTERNAL Factorial
  WRITE(6, '(a6)', advance='no') "  n = "
  READ *, number
  Pi = 0.0_qp !initialize
  sum: DO k = 0, number
    Pi = Pi + Factorial(4 * k) / Factorial(k)**4 * (26390.0_qp * FLOAT(k) + 1103.0_qp) / 396.0_qp**(4 * k)
  END DO sum
  Pi = ((2.0_qp * SQRT(2.0_qp) / FLOAT(99**2)) * Pi)**(-1)
  WRITE(6, '(f36.32)') Pi
  WRITE(6, '(f36.32)') ATAN2(0.0_qp, -1.0_qp)
END PROGRAM Pi_Ramanujan

INTEGER(li) FUNCTION Factorial(nf)
  USE, INTRINSIC :: ISO_Fortran_env, li=>INT64 !modern LONG INTEGER
  IMPLICIT NONE
  INTEGER :: nf, i
  IF (nf < 0) ERROR STOP " Factorial is singular for negative integers."
  Factorial = 1
  product: DO i = 2, nf
    Factorial = Factorial * i
  END DO product
END FUNCTION Factorial

Ramanujan.f08

PROGRAM Ramanujan
  USE, INTRINSIC :: ISO_Fortran_env, qp=>REAL128 !modern QUAD PRECISION
  IMPLICIT NONE
  REAL(qp) :: R
  REAL(qp), PARAMETER :: PI = ATAN2(0.0_qp, -1.0_qp) !mathematical constant
  R = 1.0_qp / (1.0_qp + EXP(-2.0_qp * PI) &
             / (1.0_qp + EXP(-4.0_qp * PI) &
             / (1.0_qp + EXP(-6.0_qp * PI) &
             / (1.0_qp + EXP(-8.0_qp * PI)))))
  WRITE(6, '(f36.32)') R
  WRITE(6, '(f36.32)') (SQRT((5.0_qp + SQRT(5.0_qp)) / 2.0_qp) - (1.0_qp + SQRT(5.0_qp)) / 2.0_qp) * EXP(2.0_qp * PI / 5.0_qp)
END PROGRAM Ramanujan

Ramanujan_identity.f08

PROGRAM Ramanujan_identity
  USE, INTRINSIC :: ISO_Fortran_env, qp=>REAL128 !modern QUAD PRECISION
  IMPLICIT NONE
  INTEGER :: n, number
  REAL(qp) :: sum_left, term_left, left_left, sum_right, term_right, left_right, left, middle, right
  REAL(qp), PARAMETER :: theta = 0.0_qp !less than PI
  REAL(qp), PARAMETER :: PI = ATAN2(0.0_qp, -1.0_qp) !mathematical constant
  WRITE(6, '(a6)', advance='no') "  n = "
  READ *, number
  sum_left = 0.0_qp
  sum_right = 0.0_qp
  sum: DO n = 1, number
    term_left = COS(FLOAT(n) * theta) / COSH(FLOAT(n) * PI)
    sum_left = sum_left + term_left
    term_right = COSH(FLOAT(n) * theta) / COSH(FLOAT(n) * PI)
    sum_right = sum_right + term_right
  END DO sum
  left_left = (1.0_qp + 2.0_qp * sum_left)**(-2)
  left_right = (1.0_qp + 2.0_qp * sum_right)**(-2)
  left = left_left + left_right
  middle = 2.0_qp * GAMMA(3.0_qp / 4.0_qp)**4 / PI
  right = 8.0_qp * PI**3 / GAMMA(1.0_qp / 4.0_qp)**4
  WRITE(6, '(f36.32)') left
  WRITE(6, '(f36.32)') middle
  WRITE(6, '(f36.32)') right
END PROGRAM Ramanujan_identity

The Golden Ratio

$\displaystyle{ 1 + \cfrac{1}{1 + {\cfrac{1}{1 + {\cfrac{1}{1 + \ddots}}}}} = \frac{1 + \sqrt{5}}{2} }$

A golden ratio can be identified for quantities $a$ and $b$, with $a$ > $b$ > 0⁠, if $\frac{a + b}{a} = \frac{a}{b}$; examples can be found everywhere. Written as a continued fraction, it converges quite slowly, requiring $n$ = 77 terms to reach quad precision. The best way to tackle continued fractions numerically is through the use of a recursion; in this case, the golden ratio is represented in the program as phi_1 = 1.0_qp + (1.0_qp / phi_0).

joel@Ultrabook:~/Programming/F08/Ramanujan/Golden ratio$ ./Golden_ratio n = 77 1.61803398874989484820458683436564 1.61803398874989484820458683436564

Ramanujan's Formula for $\pi$

$\displaystyle{ \frac{1}{\pi} = \frac{2 \sqrt{2}}{99^2} \sum_{k=0}^{\infty} \frac{(4k)!}{k!^4} \frac{26390 k + 1103}{396^{4k}} }$

Ramanujan's formula for $\pi$ converges very quickly compared to other methods, requiring a mere $n$ = 4 terms for quad precision. This result is checked using the standard way to fetch $\pi$ to machine precision in Fortran, ATAN2(0.0_qp, -1.0_qp).

joel@Ultrabook:~/Programming/F08/Ramanujan/Pi$ ./Pi_Ramanujan n = 4 3.14159265358979323846264338327950 3.14159265358979323846264338327950

0.99813604459850933215002445904707

$\displaystyle{ \cfrac{1}{1 + {\cfrac{\mathrm{e}^{-2 \pi}}{1 + {\cfrac{\mathrm{e}^{-4 \pi}}{1 + {\cfrac{\mathrm{e}^{-6 \pi}}{\ddots}}}}}}} = \left( \sqrt{ \frac{5 + \sqrt{5}}{2} } - \frac{1 + \sqrt{5}}{2} \right) \mathrm{e}^{2 \pi / 5} }$

I don't know what to call this one, except "a Ramanujan number". As you can see from the code, I stopped when the partial numerator reached $\mathrm{e}^{-8 \pi}$; this was enough to give quad precision. As done in the golden ratio, a RECURSIVE SUBROUTINE could have been used here, but given how fast the continued fraction converges, it's fine the way it is.

joel@Ultrabook:~/Programming/F08/Ramanujan$ ./Ramanujan 0.99813604459850933215002445904707 0.99813604459850933215002445904707

A Ramanujan Identity

$\displaystyle{ \begin{align*} & \left( 1 + 2 \sum_{n=1}^{\infty} \frac{\cos(n \theta)}{\cosh(n \pi)} \right)^{-2} + \left( 1 + 2 \sum_{n=1}^{\infty} \frac{\cosh(n \theta)}{\cosh(n \pi)} \right)^{-2} \\[6pt] & = \frac{2 \Gamma^4 \bigl( \frac{3}{4} \bigr) }{\pi} = \frac{8 \pi^3}{\Gamma^4 \bigl( \frac{1}{4} \bigr) } \end{align*} }$

This is more of a triple identity. The expression on the left-hand side converges rather slowly, requiring $n$ = 24 terms to meet quad precision. Interestingly, the middle and right-hand sides include the gamma function, $\Gamma(z) = \int_0^\infty t^{z - 1} \mathrm{e}^{-t} \, \mathrm{d} t$.

joel@Ultrabook:~/Programming/F08/Ramanujan/Identity$ ./Ramanujan_identity n = 24 1.43554002209225999564238644733156 1.43554002209225999564238644733156 1.43554002209225999564238644733156