Fortran code snippets

[Beliavsky]

The blog post How do calculators compute sine? | Algeo Calculator refers to a paper from Intel,

The Computation of Transcendental Functions on the IA-64 Architecture.

Abstract
The fast and accurate evaluation of transcendental functions (e.g. exp, log, sin, and atan) is vitally important in many fields of scientific computing. Intel provides a software library of these functions that can be called from both the C* and FORTRAN* programming languages. By exploiting some of the key features of the IA-64 floating point architecture, we have been able to provide double precision transcendental functions that are highly accurate yet can typically be evaluated in between 50 and 70 clock cycles. In this paper, we discuss some of the design principles and implementation details of these functions.

Two of the co-authors later wrote New algorithms for improved transcendental functions on IA-64.

[Beliavsky]

The program below simulates the probability that a discretely observed Brownian motion with standard normal increments will hit the barrier x = 1 within N steps, given starting position x=0.

module random_mod
implicit none
integer, parameter :: dp = kind(1.0d0)
contains
function random_normal() result(fn_val)
!   Generate a random normal deviate using the polar method.
!   Reference: Marsaglia,G. & Bray,T.A. 'A convenient method for generating
!              normal variables',Siam Rev.,vol.6,260-264,1964.
implicit none
real(kind=dp)  :: fn_val
! Local variables
real(kind=dp)           :: u,sum
real(kind=dp),save      :: v,sln
logical,save   :: second = .false.
real(kind=dp),parameter :: one = 1.0d0,vsmall = tiny(one)
if (second) then
! If second,use the second random number generated on last call
  second = .false.
  fn_val = v*sln
else
! First call; generate a pair of random normals
  second = .true.
  do
    call random_number(u)
    call random_number(v)
    u = scale(u,1) - one
    v = scale(v,1) - one
    sum = u*u + v*v + vsmall         ! vsmall added to prevent log(zero) / zero
    if (sum < one) exit
  end do
  sln = sqrt(-scale(log(sum), 1) / sum)
  fn_val = u*sln
end if
end function random_normal
end module random_mod

program xxnormal_stopping
use random_mod, only: dp, random_normal
implicit none
integer, parameter :: npaths=10**7, nsteps=10
real(kind=dp) :: zsum
integer :: i, ipath, counts(nsteps)
real(kind=dp), parameter :: thresh = 1.0_dp
counts = 0
do ipath=1,npaths
   zsum = 0.0
   do_i: do i=1,nsteps
      zsum = zsum + random_normal()
      if (zsum > thresh) then
         counts(i:) = counts(i:) + 1
         exit do_i
      end if
   end do do_i
end do
print "(*(i7))", (i, i=1,nsteps)
print "(*(f7.4))", counts / real(npaths, kind=dp)
end program xxnormal_stopping

sample output:

      1      2      3      4      5      6      7      8      9     10
 0.1588 0.2906 0.3785 0.4410 0.4882 0.5254 0.5555 0.5805 0.6018 0.6202

Are there analytical approaches to this problem that are faster than Monte Carlo? The CDF of the normal distribution can be used to obtain the probability that x > 1 after N steps, but for 3 steps, I want the probability that z(1) > 1 or z(1) + z(2) > 1 or z(1) + z(2) + z(3) > 1, where z(i) is a standard normal variate. I think there is an analytic formula for the continuous-time case.

[Beliavsky]

Gnuplot can read and plot binary files produced using unformatted stream in Fortran. Here is an example

program write_data
implicit none
integer, parameter :: dp = kind(1.0d0)
integer, parameter :: n = 10**7
real(kind=dp) :: x(n), y(n)
integer :: i
character(len=*), parameter :: filename = "temp.bin"
integer :: io_unit
open(newunit=io_unit, file=filename, status="replace", access="stream", &
   form="unformatted")
do i = 1, n
    x(i) = real(i, kind=dp)
    y(i) = x(i)**2
    write (io_unit) x(i), y(i)
end do
close(io_unit)
end program write_data

gnuplot script plot.gnu:

set title "Plot of y = x^2"
set xlabel "x"
set ylabel "y"

plot 'temp.bin' binary format='%float64%float64' using 1:2 with lines title "y = x^2"
pause -1

Run it with gnuplot plot.gnu to display the plot. Gnuplot can also read text files, but that would be slower for large files.

[Beliavsky]

John D. Cook discusses the negative binomial distribution, whose discrete PDF f(k, r, p) gives the probability of observing k failures before r successes, given a probability p of success in each trial. It is given by

y = binomial_coeff(k+r-1, k) * ((1-p)**k) * (p**r)

The code

module kind_mod
use iso_fortran_env, only: int64
implicit none
private
public :: i64, dp
integer, parameter :: i64 = int64, dp=kind(1.0d0)
end module kind_mod
!
module binomial_mod
use kind_mod, only: i64, dp
implicit none
private
public :: falling_factorial, factorial, binomial_coeff, negative_binomial, nfail
contains
elemental function falling_factorial(n, k) result(j)
! returns the product n*(n-1)*(n-2)*...*(n-k+1)
integer(kind=i64), intent(in) :: n, k
integer(kind=i64) :: j
integer(kind=i64) :: i
if (n > 0) then
   if (k <= 0) then
      j = 1
   else if (k > n) then
      j = 0
   else
      j = n
      do i=n-1, n-k+1, -1
         j = j*i
      end do
   end if
else if (n == 0) then
   j = merge(1, 0, k < 1)
end if
end function falling_factorial
!
elemental function factorial(n) result(fac)
integer(kind=i64), intent(in) :: n
integer(kind=i64)             :: fac
integer(kind=i64)             :: i
if (n < 2) then
   fac = 1
   return
end if
fac = 2
do i=3, n
   fac = fac*i
end do
end function factorial
!
elemental function binomial_coeff(n, k) result(j)
integer(kind=i64), intent(in) :: n, k
integer(kind=i64) :: j
j = falling_factorial(n, k) / factorial(k)
end function binomial_coeff
!
elemental function negative_binomial(k, r, p) result(y)
! negative binomial probability density at k, given r and p
integer(kind=i64), intent(in)  :: k ! # of failures
integer(kind=i64), intent(in)  :: r ! # of successes
real(kind=dp)    , intent(in)  :: p ! probability of success
real(kind=dp)                  :: y
y = binomial_coeff(k+r-1, k) * ((1-p)**k) * (p**r)
end function negative_binomial
!
function nfail(rmax, p) result(nf)
! simulate the # of failures before 1, 2, ..., rmax successes
integer(kind=i64), intent(in) :: rmax
real(kind=dp)    , intent(in) :: p
integer(kind=i64) :: nf(rmax)
integer(kind=i64) :: i, nsuccess
real(kind=dp) :: x
if (rmax < 1) return
i = 0
nsuccess = 0        
do
   i = i+1
   call random_number(x)
   if (x < p) then
      nsuccess = nsuccess + 1
      nf(nsuccess) = i - nsuccess
      if (nsuccess == rmax) return
   end if   
end do
end function nfail
!
end module binomial_mod
!
program xsim_negative_binomial
use binomial_mod, only: nfail, negative_binomial
use kind_mod, only: i64, dp
integer(kind=i64), parameter :: k=3, nsim=10**6, max_fail=5
integer(kind=i64)            :: isim, ifail, isucc, res(k), counts(0:max_fail, k)
real(kind=dp), parameter     :: p=0.6_dp
real(kind=dp)                :: prob_nb(0:max_fail, k), cumul_prob(0:max_fail)
print "('p =',f7.4)", p
print "('#sim = ',i0)", nsim
counts = 0
do isim=1,nsim
   res = nfail(k, p)
   do isucc=1,k
      ifail = res(isucc)
      if (ifail <= max_fail) counts(ifail, isucc) = counts(ifail, isucc) + 1
   end do
end do
print "(/,'simulated probability',/,a10,*(1x,i9))", "k/r", (ifail, ifail=0, max_fail)
do isucc=1,k
   print "(i10,*(1x,f9.6))", isucc,counts(:,isucc)/real(nsim, kind=dp)
end do
print "(/,'theoretical probability',/,a10,*(1x,i9))", "k/r", (ifail, ifail=0, max_fail)
do isucc=1, k
   prob_nb(:, isucc) = [(negative_binomial(ifail,isucc,p), ifail=0, max_fail)]
   print "(i10,*(1x,f9.6))", isucc,prob_nb(:, isucc)
end do
print "(/,'theoretical cumul probability',/,a10,*(1x,i9))", "k/r", (ifail, ifail=0, max_fail)
do isucc=1, k
   cumul_prob(0) = prob_nb(0, isucc)
   do ifail=1,max_fail
      cumul_prob(ifail) = cumul_prob(ifail-1) + prob_nb(ifail, isucc)
   end do
   print "(i10,*(1x,f9.6))", isucc,cumul_prob
end do
end program xsim_negative_binomial

gives sample output

p = 0.6000
#sim = 1000000

simulated probability
       k/r         0         1         2         3         4         5
         1  0.599892  0.239791  0.096398  0.038449  0.015158  0.006178
         2  0.360540  0.287962  0.172938  0.091871  0.045778  0.022087
         3  0.216419  0.259615  0.207468  0.137676  0.082677  0.046701

theoretical probability
       k/r         0         1         2         3         4         5
         1  0.600000  0.240000  0.096000  0.038400  0.015360  0.006144
         2  0.360000  0.288000  0.172800  0.092160  0.046080  0.022118
         3  0.216000  0.259200  0.207360  0.138240  0.082944  0.046449

theoretical cumul probability
       k/r         0         1         2         3         4         5
         1  0.600000  0.840000  0.936000  0.974400  0.989760  0.995904
         2  0.360000  0.648000  0.820800  0.912960  0.959040  0.981158
         3  0.216000  0.475200  0.682560  0.820800  0.903744  0.950193

Currently the U.S. Open in tennis is underway. Women play best of 3 sets, men best of 5. Looking at the results above, if a player has a 60% chance of winning each set, he has a 0.6480 probability of winning a 3-set match and a 0.6826 probability of winning a 5-set match.

[Beliavsky]

To override the random_number intrinsic subroutine, so that the same random numbers are generated across compilers and platforms, one can define impure elemental subroutine random_number in a module such as random_number_mod and use that subroutine, as in the following code. One can also define a suitable random_seed subroutine.

module kind_mod
implicit none
private
public :: dp
integer, parameter :: dp = kind(1.0d0)
end module kind_mod
!
MODULE Ecuyer_random ! from https://jblevins.org/mirror/amiller/taus88.f90
! L'Ecuyer's 1996 random number generator.
! Fortran version by Alan.Miller @ vic.cmis.csiro.au
! N.B. This version is compatible with Lahey's ELF90
! http://www.ozemail.com.au/~milleraj
! Latest revision - 30 March 1999
use kind_mod, only: dp
IMPLICIT NONE
public :: dp, taus88, init_seeds
! INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60)

! These are unsigned integers in the C version
INTEGER, SAVE :: s1 = 1234, s2 = -4567, s3 = 7890

CONTAINS

SUBROUTINE init_seeds(i1, i2, i3)
IMPLICIT NONE

INTEGER, INTENT(IN) :: i1, i2, i3

s1 = i1
s2 = i2
s3 = i3
IF (IAND(s1,-2) == 0) s1 = i1 - 1023
IF (IAND(s2,-8) == 0) s2 = i2 - 1023
IF (IAND(s3,-16) == 0) s3 = i3 - 1023

RETURN
END SUBROUTINE init_seeds

FUNCTION taus88() RESULT(random_numb)
! Generates a random number between 0 and 1.  Translated from C function in:
! Reference:
! L'Ecuyer, P. (1996) `Maximally equidistributed combined Tausworthe
! generators', Math. of Comput., 65, 203-213.

! The cycle length is claimed to be about 2^(88) or about 3E+26.
! Actually - (2^31 - 1).(2^29 - 1).(2^28 - 1).

IMPLICIT NONE
REAL (dp) :: random_numb

INTEGER   :: b

! N.B. ISHFT(i,j) is a bitwise (non-circular) shift operation;
!      to the left if j > 0, otherwise to the right.

b  = ISHFT( IEOR( ISHFT(s1,13), s1), -19)
s1 = IEOR( ISHFT( IAND(s1,-2), 12), b)
b  = ISHFT( IEOR( ISHFT(s2,2), s2), -25)
s2 = IEOR( ISHFT( IAND(s2,-8), 4), b)
b  = ISHFT( IEOR( ISHFT(s3,3), s3), -11)
s3 = IEOR( ISHFT( IAND(s3,-16), 17), b)
random_numb = IEOR( IEOR(s1,s2), s3) * 2.3283064365E-10_dp + 0.5_dp
RETURN
END FUNCTION taus88
END MODULE Ecuyer_random

module random_number_mod
use kind_mod, only: dp
use Ecuyer_random, only: taus88
implicit none
private
public :: dp, random_number
contains
impure elemental subroutine random_number(x)
real(kind=dp), intent(out) :: x
x = taus88()
end subroutine random_number
end module random_number_mod
!
program xmy_random_number
use kind_mod, only: dp
! remove line below to restore intrinsic random_number
use random_number_mod, only: random_number 
implicit none
real(kind=dp) :: x(2)
call random_number(x)
print*,x
end program xmy_random_number

[urbanjost]

Just an FYI that

GitHub - urbanjost/M_sets: basic set functions

was updated to process character strings (as suggested by @Beliavsky ), more integer kinds and with some trepidation floats (comparing floats for equality is problematic, of course) are allowed as set data as well.

It was originally just a one-off. I do not even remember exactly why, but some approximations to Matlab’s set theory functions were needed but just for integer .values; but traffic has been higher than expected for the module so I polished it up a bit more

[jalexiou]

I wanted to show a 3D line plot of the Lorentz attractor results with a gp%lplot() call.

! static 2D multiplot
call gp%multiplot(2,2)    
call gp%plot(x(1,:), x(2,:), 'with lines lt 5 lc rgb "blue"')
call gp%plot(x(1,:), x(3,:), 'with lines lt 5 lc rgb "green"')
call gp%plot(x(2,:), x(3,:), 'with lines lt 5 lc rgb "red"')
! add 3D line plot
call gp%lplot(x(1,:), x(2,:), x(3,:), 'using 1:2:3 with lines')     
call gp%run_script()

But the 4th subplot crashes the program and it quits, only after showing the results for a fraction of a second. I wonder if I am using Gnuplot incorrectly, and if @hkvzjal you can show me the correct way to add a 3D line plot to a multiplot in Fortran and ogpf module.

---

I found a solution.

There is a bug in ogpf.f90 in the lplot3d() subroutine (line 1402 in the file). A space should be added in the 'with lines' string such that it reads like ' with lines'.

The full updated code is

     if ( present(lspec) ) then
         if (hastitle(lspec)) then
             pltstring='splot ' // datablock // ' ' // trim(lspec) // ' with lines'
         else
             pltstring='splot ' // datablock // ' notitle '//trim(lspec) // ' with lines'
         end if
     else
         pltstring='splot ' // datablock // ' notitle with lines'
     end if

and then call the function with use the line color info

CALL gp%xlabel('x-axis, ..')
CALL gp%ylabel('y-axis, ..')
CALL gp%zlabel('z-axis, ..')
call gp%lplot(x(1,:), x(2,:), x(3,:), 'lt 5 lc rgb "black"')     ! show 3D line plot
 !                                     ^
 !                                     no `with lines` needed here
CALL gp%xlabel('x-axis, ..')
CALL gp%ylabel('y-axis, ..')
call gp%plot(x(1,:), x(2,:), 'with lines lt 5 lc rgb "blue"')    ! regular 2D line plot
 !                            ^
 !                            `with lines` is needed here

[hkvzjal]

Nice catch! I’m not not the maintainer of this library though, I just use it occasionally. I think that opening a PR or an issue at the github repo would be appreciated by the author

[Beliavsky]

A relative forwarded a post

Our next calendar year is a mathematical wonder.

Interesting 2025

1. 2025, itself is a square

2. It’s a product of two squares,
   Viz. 9² x 5² = 2025

3. It is the sum of 3-squares,
   viz. 40²+ 20²+5²= 2025

4. It’s the first square after 1936

5. It’s the sum of cubes, of all the single digits, from 1 to 9,
   viz. 1³+2³+3³+4³+5³+6³+7³+8³+9³= 2025.

A Fortran program verifying this is

integer :: i
print*,sqrt(2025.0)
print*,(9*5)**2
print*,44**2
print*,sum([5, 20, 40]**2)
print*,sum([(i**3, i=1, 9)])
end

giving output

   45.0000000    
        2025
        1936
        2025
        2025

[urbanjost]

But is is so evil!

Or why numerologists can always make the facts meet the story …
(Properties of number 2025).

but it is the 45th square, and it’s square root is 45; which seems apocalyptic to me; and sure enough it is one!

[Beliavsky]

Since a Fortran function can return an allocatable array of deferred-length character variables, instead of a verbose character array constructor where one must count the length of the longest word

x = [character (len=5) :: "one", "two", "three"]

one can define a words function that separates a string into words separated by spaces and sets the LEN to that of the longest word and write

x = words("one two three")

Here is the function and a program that calls it.

module mystrings
  implicit none
contains
  pure function words(string) result(word_arr)
    ! Returns an allocatable array of character strings.
    ! Each element’s length equals the length of the longest word.
    character(len=*), intent(in)  :: string
    character(len=:), allocatable :: word_arr(:)
    integer :: lenStr, i, start, endPos, numWords, maxLen
    character (len=1), parameter :: sep = " "
    ! Determine effective length of the input string (ignoring trailing spaces)
    lenStr = len_trim(string)
    numWords = 0
    maxLen   = 0
    i = 1

    ! First pass: count words and determine maximum word length.
    do while (i <= lenStr)
      if (string(i:i) /= sep) then
        start = i
        do while (i <= lenStr .and. string(i:i) /= sep)
          i = i + 1
        end do
        endPos = i - 1
        numWords = numWords + 1
        if (endPos - start + 1 > maxLen) maxLen = endPos - start + 1
      else
        i = i + 1
      end if
    end do

    ! Allocate the output array with each element having length maxLen.
    allocate(character(len=maxLen) :: word_arr(numWords))

    ! Second pass: extract each word into the allocated array.
    i = 1
    numWords = 0
    do while (i <= lenStr)
      if (string(i:i) /= ' ') then
        start = i
        do while (i <= lenStr .and. string(i:i) /= ' ')
          i = i + 1
        end do
        endPos = i - 1
        numWords = numWords + 1
        word_arr(numWords) = string(start:endPos)
      else
        i = i + 1
      end if
    end do
  end function words
end module mystrings

program test_words
  use mystrings, only: words
  implicit none
  character(:), allocatable, dimension(:) :: arr
  integer :: i
  real :: t1, t2
  call cpu_time(t1)
  arr = words("one two three four five six seven")
  call cpu_time(t2)
  print *, "Extracted words (each element has len =", len(arr(1)), "):"
  do i = 1, size(arr)
    print *, "'" // arr(i) // "'"
  end do
  print*,"time taken by words():", t2-t1
end program test_words

Output:

 Extracted words (each element has len =           5 ):
 'one  '
 'two  '
 'three'
 'four '
 'five '
 'six  '
 'seven'
 time taken by words():   0.00000000    

The words function is replacing a character array of literals but is executed in a negligible amount of time. Fortran 2023 has a tokenize subroutine, described here, but a function is more convenient.

[urbanjost]

This is related, but takes up to twenty strings and turns them into an array long enough to accomodate the longest one, avoiding possible truncation when using the mentioned standard syntax.

https://urbanjost.github.io/M_strings/bundle.3m_strings.html

several compilers as an extension do that with a simple declaration which I wish was standard behavior. That is,

 array=['one', 'two', 'three'] 

of course if array is not allocatable you can still get truncation if array has a length too small; and if it is allocatable you may create an unexpectedly long length so not perfect either.

[Harper]

@Beliavsky’s program would not compile with my system using gfortran or ifx because
two inner loops began with things of the form

do while ( A .and. B )

but that does not prescribe the order in which A, B are calculated.
A workaround that compiled and ran correctly with gfortran, ifx and lfortran is

1. In the First pass inner loop replace

do while (i <= lenStr .and. string(i:i) /= sep)

by

do
    if (i > lenStr) exit
    if (string(i:i) == sep) exit

2. In the Second pass inner loop replace

do while (i <= lenStr .and. string(i:i) /= ' ')

by

 do
    if (i>lenStr) exit
    if (string(i:i) == ' ') exit

[Beliavsky]

Thanks. I should have compiled with gfortran -fbounds-check. A module that fixes the errors, incorporating the suggestions, is

module mystrings
  implicit none
  private
  public :: words
contains
  pure function words(string) result(word_arr)
    ! Returns an allocatable array of character strings.
    ! Each element’s length equals the length of the longest word.
    character(len=*), intent(in)  :: string
    character(len=:), allocatable :: word_arr(:)
    integer :: lenStr, i, start, endPos, numWords, maxLen
    character (len=1), parameter :: sep = " "
    ! Determine effective length of the input string (ignoring trailing spaces)
    lenStr = len_trim(string)
    numWords = 0
    maxLen   = 0
    i = 1

    ! First pass: count words and determine maximum word length.
    do while (i <= lenStr)
      if (string(i:i) /= sep) then
        start = i
        do
           if (i > lenStr) exit
           if (string(i:i) == sep) exit
           i = i + 1
        end do
        endPos = i - 1
        numWords = numWords + 1
        if (endPos - start + 1 > maxLen) maxLen = endPos - start + 1
      else
        i = i + 1
      end if
    end do

    ! Allocate the output array with each element having length maxLen.
    allocate(character(len=maxLen) :: word_arr(numWords))

    ! Second pass: extract each word into the allocated array.
    i = 1
    numWords = 0
    do while (i <= lenStr)
      if (string(i:i) /= sep) then
        start = i
        do
           if (i > lenStr) exit
           if (string(i:i) == sep) exit
           i = i + 1
        end do
        endPos = i - 1
        numWords = numWords + 1
        word_arr(numWords) = string(start:endPos)
      else
        i = i + 1
      end if
    end do
  end function words
end module mystrings

[Beliavsky]

John D. Cook’s math blog can often inspire a Fortran program.

Chebyshev’s inequality says that the probability of a random variable being more than k standard deviations away from its mean is less than 1/k 2. In symbols,

This inequality is very general, but also very weak. It assumes very little about the random variable X but it also gives a loose bound. If we assume slightly more, namely that X has a unimodal distribution, then we have a tighter bound, the Vysochanskiĭ-Petunin inequality.

However, the Vysochanskiĭ-Petunin inequality does require that k be larger than √(8/3). In exchange for the assumption of unimodality and the restriction on k we get to reduce our upper bound by more than half.

program verify_vp
  !---------------------------------------------------------------------
  ! Verify the Vysochanskiĭ–Petunin inequality for the standard normal,
  ! and compute corresponding tail probabilities for the standardized
  ! uniform and double‐exponential (Laplace) distributions.
  !
  !   P(|Z| ≥ k) ≤ 4/(9 k^2)   for Z ~ N(0,1)
  !
  ! Tail probabilities for
  !   U ~ Uniform(–√3,√3)  (Var=1):   P(|U| ≥ k) = max(0, 1 – k/√3)
  !   L ~ Laplace(0,b) with b=1/√2 (Var=1): P(|L| ≥ k) = exp(−k/b)
  !
  !---------------------------------------------------------------------
  implicit none
  ! parameters controlling k loop
  integer, parameter :: dp = kind(1.0d0), nk = 3
  real(kind=dp), parameter :: k_min = 2.0_dp, k_inc = 1.0_dp

  ! parameters for the other distributions
  real(kind=dp), parameter :: a = sqrt(3.0_dp)         ! half‐width of Uniform(–√3,√3)
  real(kind=dp), parameter :: b = 1.0_dp/sqrt(2.0_dp)  ! scale of Laplace for Var=1
  real(kind=dp) :: k, phi_val, p_norm, bound, p_unif, p_lap
  integer :: i

  print "(A6,*(2X, A12))", "k", "Bound", " P_uniform", " P_normal", " P_laplace"

  do i=1,nk
    k = k_min + (i-1)*k_inc
    ! --- normal(0,1) tail and VP bound ---
    phi_val = 0.5_dp*(1.0_dp + erf(k/sqrt(2.0_dp)))
    p_norm  = 2.0_dp*(1.0_dp - phi_val)
    bound   = 4.0_dp/(9.0_dp*k**2)

    ! --- uniform(–√3,√3), Var=1 ---
    p_unif = merge(1.0_dp - k/a, 0.0_dp, k < a)

    ! --- Laplace(0,b), Var=1 ---
    p_lap = exp(-k/b)

    print "(F6.3,2X,*(F12.6,2X))", k, bound, p_unif, p_norm, p_lap
  end do

end program verify_vp

output:

     k         Bound     P_uniform      P_normal     P_laplace
 2.000      0.111111      0.000000      0.045500      0.059106
 3.000      0.049383      0.000000      0.002700      0.014370
 4.000      0.027778      0.000000      0.000063      0.003493

The results match Cook’s for the three distributions implemented.

[jkd2022]

Fortran for the Common Man

I modeled a trumpet with the function

A(t)=\sum_n A_n\sum_i g_i \exp[i f_n \omega_i(t-t_n)]\exp[-(t-t_n)^2/d_n^2]({\rm erf}[(t-t_n)/s]+1)

using the dominant frequencies \omega_i and complex amplitudes g_i obtained from the Fourier transform of a real trumpet, and added this together for the notes n in Copland’s Fanfare for the Common Man with the frequency formula

f_n=f_02^{n/12}

where f_0=466.16 Hz – the frequency at which the trumpet was sampled.

Here’s the result:

program Fortran_for_the_Common_Man
implicit none
real, allocatable :: t(:),at(:),notes(:,:),w(:)
complex, allocatable :: g(:)
integer nw,nn,sr,nt,i,j
real tmax,atm,dt
! dominant frequencies and amplitudes of a sampled trumpet in A#4 (466.16 Hz)
w=[5852.79, 8780.75, 2927.96, 5874.78, 8809.03, 5892.06, 14635.11, 11743.27, &
 2945.24, 11711.86, 17564.64, 8763.47, 8831.02, 20491.04, 14708.94, 11760.55, &
 14668.10, 11802.96, 14734.07, 8741.48, 5909.34, 14693.23, 5833.94, 8719.49, &
 11821.81, 17581.92, 14751.35, 17547.37, 20467.48, 11677.30, 5934.47, 11787.26,&
 2912.26, 8848.30, 5816.66, 11694.58, 14768.63, 17597.63, 11647.45, 17652.61, &
 17616.48, 5797.81, 11837.52, 20525.60, 8689.65, 11625.46, 23420.57, 32215.46, &
 17635.33, 14609.98, 5973.74, 17522.23, 2968.81, 5956.46, 8879.71, 20437.63, &
 17677.74, 9401.22, 26398.80, 98111.94, 3026.92, 20582.14, 20508.32, 8909.56, &
 20610.42, 29285.93, 17495.53, 3044.20, 23510.11, 20649.69, 20550.73, 35143.43,&
 95188.69, 2998.65, 29260.79, 14785.91, 3091.33, 32191.90, 11548.49, 26348.54, &
 8648.80, 26312.41, 14555.00, 11873.65, 9418.49, 8926.84, 11854.80, 2883.98, &
 11573.63, 26331.26]
g=[(-12.847,73.162), (12.656,66.056), (-23.249,-36.155), (15.392,-29.912), &
 (-9.734,29.441), (-10.744,-22.579), (-15.942,18.515), (1.342,-21.755), &
 (-6.791,-19.484), (-7.380,-19.161), (-5.987,-19.299), (-16.223,6.601), &
 (-2.903,14.915), (-12.222,-6.706), (13.797,1.916), (4.005,-12.722), &
 (-11.090,-6.728), (5.215,11.284), (-6.107,10.629), (-9.961,-6.950), &
 (-10.454,0.239), (3.417,9.764), (2.810,9.689), (-1.275,-9.433), &
 (6.637,5.847), (6.224,-5.807), (-2.296,7.834), (-1.440,7.980), &
 (-6.721,-4.521), (4.822,6.308), (7.258,-1.395), (1.987,6.964), &
 (-1.086,7.157), (6.035,3.417), (-6.105,-2.844), (-2.598,-6.208), &
 (1.848,6.094), (2.791,-5.262), (5.642,-1.187), (5.346,-1.735), &
 (2.827,4.786), (1.816,-5.042), (3.186,-4.157), (2.638,-4.399), &
 (-4.582,-1.908), (3.117,3.826), (-3.858,-2.813), (-3.107,-3.526), &
 (4.307,1.661), (-4.192,-1.810), (-4.278,-0.335), (-3.431,2.212), &
 (-2.222,-3.244), (2.699,-2.672), (2.997,-1.720), (-2.918,1.811), &
 (2.379,2.108), (-3.148,-0.344), (1.851,-2.560), (-0.129,-3.147), &
 (-2.952,0.815), (1.090,2.855), (2.997,0.341), (2.929,0.608), &
 (-2.493,-1.624), (-2.771,0.838), (-1.860,2.179), (2.589,1.162), &
 (0.492,2.790), (-0.055,-2.769), (-2.245,-1.594), (0.044,-2.660), &
 (-2.604,0.409), (-2.566,-0.530), (-2.241,1.320), (2.418,0.954), &
 (-2.043,-1.496), (-0.340,-2.474), (-2.477,-0.181), (-1.400,-1.984), &
 (0.551,2.364), (-2.192,0.925), (-2.196,0.807), (-1.703,-1.589), &
 (0.695,2.207), (-0.914,-2.125), (-0.725,-2.189), (-0.246,2.268), &
 (-1.769,1.429), (-0.640,2.177)]
! number of frequencies and amplitudes
nw=size(w)
! excerpt from Fanfare for the Common Man by Aaron Copland
! note, amplitude, duration, time
notes=reshape([ &
 -2., 0.78, 0.078,  0.500,  5., 0.86, 0.093,  0.694, 10., 0.94, 1.548,  0.886, &
  3., 0.77, 0.090,  2.032, 10., 0.72, 0.102,  2.224,  8., 0.89, 2.244,  2.416, &
 12., 0.89, 1.098,  3.956,  8., 0.91, 1.137,  4.718,  3., 0.94, 1.128,  5.488, &
 -2., 0.86, 2.691,  6.260, -4., 0.80, 0.102,  8.176,  0., 0.83, 0.099,  8.374, &
  3., 0.88, 1.347,  8.580,  3., 0.88, 0.099,  9.708,  7., 0.80, 0.108,  9.904, &
 10., 0.95, 1.338, 10.096, -4., 0.80, 0.096, 11.240,  0., 0.84, 0.108, 11.428, &
  3., 0.89, 1.311, 11.632,  7., 0.89, 0.255, 12.782, -2., 0.86, 0.282, 13.164, &
 -9., 0.90, 1.539, 13.572 ], [4,22])
! number of notes
nn=size(notes,2)
! maximum time in seconds
tmax=maxval(notes(4,1:nn))+5.
! sample rate
sr=44100
! number of time steps
nt=tmax*sr
allocate(t(nt),at(nt))
! time steps
do i=1,nt
  t(i)=tmax*real(i-1)/nt
end do
! sum the notes to get the total audio amplitude as a function of time
write(*,'("Generating time-dependent audio amplitude...")')
at(:)=0.0
do i=1,nt
  do j=1,nn
    if (abs(t(i)-notes(4,j)) > 8.*notes(3,j)) cycle
    at(i)=at(i)+atn(j,t(i))
  end do
end do
! find the maximum amplitude
atm=maxval(abs(at(:)))
! normalize to 1
at(:)=at(:)/atm
! write to file in PCM 32-bit floating-point little-endian format
write(*,'("Writing to file...")')
open(50,file='fftcm.f32',access='STREAM',status='REPLACE')
do i=1,nt
  write(50) at(i)
end do
close(50)
! uncomment either line to play from the Fortran executable
!call execute_command_line('sox -r 44100 -c 1 -e floating-point fftcm.f32 -d')
!call execute_command_line('ffplay -ar 44100 -f f32le fftcm.f32')
contains

pure real function atn(j,t)
implicit none
integer, intent(in) :: j
real, intent(in) :: t
real, parameter :: s=0.05  ! skewness
real n,a,d,t0,dt,f
n=notes(1,j); a=notes(2,j); d=notes(3,j); t0=notes(4,j)
dt=t-t0
! frequency in terms of number of semitones away from A#4
f=2.**(n/12.)
! model trumpet note
atn=a*real(sum(g(:)*exp((0.,1.)*f*w(:)*dt)))*exp(-(dt/d)**2)*(erf(dt/s)+1.)
end function

end program

It compiles and runs with both Intel Fortran and GFortran. The audio file can be played with either

sox -r 44100 -c 1 -e floating-point fftcm.f32 -d

or

ffplay -ar 44100 -f f32le fftcm.f32

Alternatively, uncomment a ‘call execute_command_line’ line in the code and the Fortran executable will play it for you.

Here’s a challenge for someone: Fortran the Brave on Bagpipes.

[certik]

@jkd2022

I modeled a trumpet with the function

Nice! 8 years ago I implemented a piano string simulation: Ondřej Čertík / stringsim · GitLab, it solves a non-linear wave equation from the paper in the README using finite differences (130 lines of Fortran: src/tests/test_wave.f90 · 26a0ad4612e3f2e1ac8d5e0eb42b96ae49e7b3e2 · Ondřej Čertík / stringsim · GitLab), and you can then play the sound. It’s not bad, but not as good as a real piano.

[jeff.irwin]

@Beliavsky your post reminds me of the topic of rational interpolation. Any function can be interpolated by a ratio of polynomials using Neville’s algorithm. Plugging in a little code from my numerical-analysis repo, I can compare Kellog’s and Neville’s approximations. Accuracy can be improved by adding more support points xs.

The code:

module approx_mod
implicit none
integer, parameter :: dp = kind(1.0d0)

contains

elemental function log_approx(p, q) result(y)
	real(kind=dp), intent(in) :: p, q
	real(kind=dp)             :: y
	real(kind=dp)             :: s, d
	s = p + q
	d = p - q
	y = 6*d*s/(3*s**2 - d**2)
end function log_approx

elemental function log_neville(x) result(fx)
	real(kind=dp), parameter :: xs(*) = [1.05d0, 0.95d0, 0.85d0, 0.75d0, 0.65d0]
	real(kind=dp), parameter :: fi(*) = log(xs)
	real(kind=dp), intent(in) :: x
	real(kind=dp) :: fx
	fx = neville_rational_interpolator_vals(xs, fi, x)
end function log_neville

pure function neville_rational_interpolator_vals(xi, fi, x) result(fx)
	! Interpolate data at `x` using Neville's algorithm with rational polynomials
	real(kind=dp), intent(in) :: xi(:), fi(:)
	real(kind=dp), intent(in) :: x
	real(kind=dp) :: fx
	real(kind=dp) :: tk2
	real(kind=dp), allocatable :: t(:), t0(:)
	integer :: i, k, n1
	n1 = size(xi)  ! Degree of polynomial + 1
	allocate(t(n1), t0(n1))
	t0 = 0.d0
	do i = 1, n1
		t(1) = fi(i)
		tk2 = 0.d0
		do k = 2, i
			if (k > 2) tk2 = t0(k-2)
			t(k) = t(k-1) + (t(k-1) - t0(k-1)) &
				/ ((x - xi(i-k+1)) / (x - xi(i)) * &
				( &
					1.d0 - (t(k-1) - t0(k-1)) / (t(k-1) - tk2) &
				) - 1.d0)
		end do
		t0 = t
	end do
	fx = t(n1)
end function neville_rational_interpolator_vals

end module approx_mod

program xapprox
	use approx_mod, only: dp, log_approx, log_neville
	implicit none
	integer :: i
	real(kind=dp), parameter :: x(*) = 0.1_dp*[(i,i=11,1,-1)]
	character (len=*), parameter :: fmt_r = "(a8, *(1x,f10.7))"
	print fmt_r,"x", x
	print fmt_r,"log(x)", log(x)
	print fmt_r,"approx", log_approx(x,1.0_dp)
	print fmt_r,"neville", log_neville(x)
	print fmt_r,"diff", log(x) - log_approx(x,1.0_dp)
	print fmt_r,"diff nev", log(x) - log_neville(x)
end program xapprox

gives:

       x  1.1000000  1.0000000  0.9000000  0.8000000  0.7000000  0.6000000  0.5000000  0.4000000  0.3000000  0.2000000  0.1000000
  log(x)  0.0953102  0.0000000 -0.1053605 -0.2231436 -0.3566749 -0.5108256 -0.6931472 -0.9162907 -1.2039728 -1.6094379 -2.3025851
  approx  0.0953101  0.0000000 -0.1053604 -0.2231405 -0.3566434 -0.5106383 -0.6923077 -0.9130435 -1.1921397 -1.5652174 -2.1063830
 neville  0.0953082  0.0000003 -0.1053607 -0.2231433 -0.3566756 -0.5108169 -0.6930029 -0.9152934 -1.1987629 -1.5841029 -2.1618745
    diff  0.0000000  0.0000000 -0.0000001 -0.0000031 -0.0000316 -0.0001873 -0.0008395 -0.0032473 -0.0118331 -0.0442205 -0.1962021
diff nev  0.0000019 -0.0000003  0.0000002 -0.0000002  0.0000007 -0.0000087 -0.0001443 -0.0009973 -0.0052099 -0.0253350 -0.1407106

[vmagnin]

@certik

Nice! 8 years ago I implemented a piano string simulation: Ondřej Čertík / stringsim · GitLab, it solves a non-linear wave equation from the paper in the README using finite differences (130 lines of Fortran: src/tests/test_wave.f90 · 26a0ad4612e3f2e1ac8d5e0eb42b96ae49e7b3e2 · Ondřej Čertík / stringsim · GitLab), and you can then play the sound. It’s not bad, but not as good as a real piano.

I have added the following loop in your src/tests/test_wave.f90 file to obtain an output file that can be imported in my ForSynth project:

open(newunit=u, file="sound.txt", status="replace")
do m = 1, Nt, Nskip
    write(u, *) yt(m)
end do
close(u)

The ForSynth sonify command can then be used directly to generate a WAV:

$ sonify -i sound.txt
 Nb of tracks, excluding track 0:           1
 Used RAM:                 1765768 bytes
 File size ~               1765768 bytes
 You can now play the file sonification.wav

[Beliavsky]

I like the f0.d format for printing reals, except that it does not print leading zeros, and for very large numbers the digits after the decimal point may be noise. A program with a subroutine that prints with the f0.6 format but adds leading zeros and switches to scientific notation for large numbers is

module util_mod
implicit none
private
public :: dp, print_real
integer, parameter :: dp = kind(1.0d0)
contains
impure elemental subroutine print_real(x)
! print real x with leading zero for abs(x) < 1
real(kind=dp), intent(in) :: x
if (abs(x) < 1.0_dp) then
   if (x >= 0) then
      print "(F8.6)", x
   else
      print "(F9.6)", x ! space for leading negative sign
   end if
else if (abs(x) > 1.0e22_dp) then ! use scientific notation
   if (x >= 0) then
      print "(ES12.6)", x
   else
      print "(ES13.6)", x ! space for leading negative sign
   end if
else
   print "(F0.6)", x
end if
end subroutine print_real
end module util_mod

program xprint_real
use util_mod, only: dp, print_real
implicit none
call print_real([1.1_dp, 0.9_dp, -0.9_dp, &
   1.0e22_dp, -1.0e22_dp, 1.0e23_dp, -1.0e23_dp])
end program xprint_real

with output

1.100000
0.900000
-0.900000
10000000000000000000000.000000
-10000000000000000000000.000000
1.000000E+23
-1.000000E+23

Fortran has added leading zero format control. How would you use that to print with the f0.6 format but with a leading zero? Do any compilers implement it yet?