for my current project, I wrote a function in plane Julia. It reads

#Kappa squared or the square of the inverse debye length
Debyelengthinnversesquared(is) = e^2*beta*2*is*n_a/(epsilon*epsilon_r)
#Kappa in SI
Kappa(is) = sqrt(Debyelengthinnversesquared(1000*is))
#Kappa in Angstroms
KappaA(is) = sqrt(Debyelengthinnversesquared(1000*is))*10^(-10)
#Surfacepotential
function surface_potential(kappa::Float64, sigma::Float64)::Float64
log_arg = c*sigma/(epsilon*epsilon_r* kappa) +sqrt((c *sigma/(epsilon*epsilon_r*kappa))^2+1)
return 1/c*log(log_arg)
end
#Defining gamma functions
function gamma(kappa, sigma)
return (0.25*e*surface_potential(kappa, sigma)*beta) |> tanh
end
gamma_daly(is) = gamma(Kappa(is), -0.02)
function Guoy_chapman_potential_AU(is::Float64, z::Float64)::Float64
arg1 = 1 + gamma_daly(is)*exp(-KappaA(is)*z)
arg2 = 1 - gamma_daly(is)*exp(-KappaA(is)*z)
result = 2*log(arg1/arg2)
return result
end

Later since I head a code base in Fortran (I will not repeat the code here). I then used the great Fortran capabilities via Ccall to call the function from Fortran (Just to check that it gives the right values.)

Just for fun I also did a speed test, and I found for the Fortran version

@time values = [Guoy_chapman_fortran(i*1.0, 0.005, -0.02) for i=0:10^6];

the output was

0.101856 seconds (23.92 k allocations: 9.236 MiB, 15.77% compilation time)

However, when I performed the same task in plane Julia, I received for

@time value_GC = [Guoy_chapman_potential_AU(0.005::Float64, i*1.0::Float64)::Float64 for i=0:10^6]

an output of

1.053988 seconds (83.02 M allocations: 1.246 GiB, 2.63% gc time, 0.81% compilation time)

So that is slower by a factor of 10. Does anyone know what I am doing wrong?

note that this includes compilation time as well. Part of the reason for using BenchmarkTools is to avoid things like this, but also to get a more reliable measurement by running the expression multiple times in case your operating system is busy with other things just when you ran your timing test.

Note also that BenchmarkTools also automatically calls cheap functions in a loop as many times as needed to get good timing resolution, so that you can directly do:

There is actually a standard syntax for floating point scientific notation, namely 1.38e-23 instead of 1.38*10^(-23).

At first I thought that there could be a performance difference, since the former is directly parsed as a Float64 while the latter is a Float64 multiplied by a exponentiated Int, but the compiler is smart enough to calculate the value for you.

But there is a difference in accuracy, since your version necessarily must be equal to the product of two floats (of very different magnitudes):

julia> 1.38e-23
1.38e-23
julia> 1.38*10^(-23)
1.3799999999999998e-23 # <- slight loss of precision

So the standard notation is more accurate, as well as more compact and idiomatic.

Note that, despite appearances, this value is also inexact.

Remember that the floating-point printing algorithm used in Julia (the Ryu algorithm) follows the well-established rubric of printing the shortest decimal representation that rounds to the same value. The Float64 constant 1.38e-23 is actually:

Small moment of appreciation for our collective vulnerability to nerdsniping. Some may find it annoying but I love it when people pile up even after a topic has been solved to explore multiple ways to make the code better or more interesting