# Fast\_math: A collection of functions for fast number crunching using Fortran

**URL:** <https://fortran-lang.discourse.group/t/fast-math-a-collection-of-functions-for-fast-number-crunching-using-fortran/6589>\
**Category:** Announcements\
**Created:** [October 1, 2023, 10:53am UTC](https://fortran-lang.discourse.group/t/fast-math-a-collection-of-functions-for-fast-number-crunching-using-fortran/6589 "2023-10-01T10:53:30Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![hkvzjal](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/hkvzjal/32/3055_2.png) [@hkvzjal](https://fortran-lang.discourse.group/u/hkvzjal)\
**Post date:** [October 1, 2023, 10:53am UTC](https://fortran-lang.discourse.group/t/fast-math-a-collection-of-functions-for-fast-number-crunching-using-fortran/6589/1 "2023-10-01T10:53:30Z")

</div>

Dear all,

Following the steps of [fastmath](https://fortran-lang.discourse.group/t/fast-math-library-for-modern-fortran/4612) by @FedericoPerini, I would like to introduce

> **[GitHub - jalvesz/fast\_math: A fortran library of fast functions](https://github.com/jalvesz/fast_math/tree/main)**
>
> main

This gallery contains several functions, some are own work, others from discussions here in the forum and other forums in the net.

Parts of the original work from @FedericoPerini have been merged here, you will also find implementations for fast sum and fast dot\_product which are not only faster but also more accurate than the intrinsic (at least with gfortran and ifort under my available hardware 🙂 ).

You will find a response file for you to test easily with

```auto
>fpm @testgfortran

```

```auto
>fpm @testifort

```

```auto
>fpm @testifx

```

A teaser:

```auto
================================================================
  SUM on single precision values with a mask
  sum_quad: <time> = 3.9091 ns/eval, speed-up= 0.99X, rel. error= 0.0000E+00
 intrinsic: <time> = 3.8577 ns/eval, speed-up= 1.00X, rel. error= 6.8401E-05
 fsum_pair: <time> = 4.1589 ns/eval, speed-up= 0.93X, rel. error= 1.3716E-05
      fsum: <time> = 0.4318 ns/eval, speed-up= 8.93X, rel. error= 2.1368E-05

================================================================
  dot product on single precision values
  dot_quad: <time> = 1.2195 ns/eval, speed-up= 0.94X, rel. error= 0.0000E+00
 intrinsic: <time> = 1.1446 ns/eval, speed-up= 1.00X, rel. error= 3.6655E-04
     fprod: <time> = 0.1533 ns/eval, speed-up= 7.47X, rel. error= 7.5047E-07

================================================================
  Fast hyperbolic
 ftanh r64: <time> = 2.7130 ns/eval, speed-up= 6.95X, rel. error= 2.5361E-07
 ftanh r64: <time> = 3.7820 ns/eval, speed-up= 5.26X, rel. error= 2.4636E-07
  ferf r32: <time> = 1.8050 ns/eval, speed-up= 8.74X, rel. error= 3.7337E-04
  ferf r64: <time> = 4.2750 ns/eval, speed-up= 3.77X, rel. error= 3.7337E-04

```

Looking forward to your comments and contributions!

---

<div class="post-metadata">

**Author:** ![rwmsu](https://avatars.discourse-cdn.com/v4/letter/r/48db29/32.png) [@rwmsu](https://fortran-lang.discourse.group/u/rwmsu)\
**Post date:** [October 1, 2023, 12:11pm UTC](https://fortran-lang.discourse.group/t/fast-math-a-collection-of-functions-for-fast-number-crunching-using-fortran/6589/2 "2023-10-01T12:11:12Z")

</div>

I guess the question is why are the compiler vendors still using slower and less accurate methods for their intrinsic functions when its been demonstrated that faster and more accurate methods exist. I can only guess that the old “backward compatability” and “would break user codes (ie results)” mantra is the reason.

---

<div class="post-metadata">

**Author:** ![hkvzjal](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/hkvzjal/32/3055_2.png) [@hkvzjal](https://fortran-lang.discourse.group/u/hkvzjal)\
**Post date:** [October 1, 2023, 12:28pm UTC](https://fortran-lang.discourse.group/t/fast-math-a-collection-of-functions-for-fast-number-crunching-using-fortran/6589/3 "2023-10-01T12:28:17Z")

</div>

> [@rwmsu](#):
>
> I can only guess that the old “backward compatibility” and “would break user codes (ie results)” mantra is the reason.

I guess so! for the trigonometric and other functions I’m ok with the compiler giving by default the most accurate version, and giving the user the responsibility of using a “fast and dirty” version if they want to… for the sum and dot\_product, I’m indeed impressed that it was possible to be both faster and more accurate at the same time compared to the compilers, would have expected that only one was possible.

---

<div class="post-metadata">

**Author:** ![FedericoPerini](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/federicoperini/32/1750_2.png) [@FedericoPerini](https://fortran-lang.discourse.group/u/FedericoPerini)\
**Post date:** [October 1, 2023, 12:53pm UTC](https://fortran-lang.discourse.group/t/fast-math-a-collection-of-functions-for-fast-number-crunching-using-fortran/6589/4 "2023-10-01T12:53:09Z")

</div>

> [@rwmsu](#):
>
> I can only guess that the old “backward compatability”

Not an expert, but I would bet that ensuring IEEE arithmetics compliance for all edge cases is also part of the answer.

---

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [October 1, 2023, 1:53pm UTC](https://fortran-lang.discourse.group/t/fast-math-a-collection-of-functions-for-fast-number-crunching-using-fortran/6589/5 "2023-10-01T13:53:37Z")

</div>

For the intrinsic `sum` at least I don’t believe there are many guarantees. There’s been a lengthy discussion before here: [Some Intrinsic SUMS](https://fortran-lang.discourse.group/t/some-intrinsic-sums/5760)

The standard appears to wash it’s hands by stating the result is processor-dependent:

> The result of SUM (ARRAY) has a value equal to a **processor-dependent** approximation to the sum of all the elements of ARRAY or has the value zero if ARRAY has size zero

(Processor being jargon for the Fortran compiler in typical scenarios.)

Especially for global sums across multiple ranks (MPI) or images (coarrays), enhanced accurary is needed to get consistent results. The following paper discusses the issue: [Redirecting](https://doi.org/10.1016/j.parco.2011.02.009)

---

<div class="post-metadata">

**Author:** ![vmagnin](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/vmagnin/32/28_2.png) [@vmagnin](https://fortran-lang.discourse.group/u/vmagnin)\
**Post date:** [October 1, 2023, 2:08pm UTC](https://fortran-lang.discourse.group/t/fast-math-a-collection-of-functions-for-fast-number-crunching-using-fortran/6589/6 "2023-10-01T14:08:29Z")

</div>

> [@ivanpribec](#):
>
> For the intrinsic `sum` at least I don’t believe there are many guarantees.

The floating point addition being not associative contrarily to the mathematical addition, who knows what is the best order to compute a given sum? It depends on the terms in the sum.

We may first sort the numbers then sum them starting from the smallest ones, but it would strongly impede speed.

---

<div class="post-metadata">

**Author:** ![hkvzjal](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/hkvzjal/32/3055_2.png) [@hkvzjal](https://fortran-lang.discourse.group/u/hkvzjal)\
**Post date:** [October 1, 2023, 2:28pm UTC](https://fortran-lang.discourse.group/t/fast-math-a-collection-of-functions-for-fast-number-crunching-using-fortran/6589/7 "2023-10-01T14:28:05Z")

</div>

> [@ivanpribec](#):
>
> There’s been a lengthy discussion before here: [Some Intrinsic SUMS](https://fortran-lang.discourse.group/t/some-intrinsic-sums/5760)

Yes, I participated in that discussion and it motivated me to seek an approach that would give some kind of compromise between the accuracy and speed issues.

> [@vmagnin](#):
>
> We may first sort the numbers then sum them starting from the smallest ones, but it would strongly impede speed.

This, or having a fully random walk before adding up could minimize the error, yet indeed the price in computational time would be too high.

That’s why I came up with the idea of “chunking” the sum/dot\_product within a small vector that would separate sequential values, yet keep crunching incoming batches of memory as they come and avoid needing to jump in memory or sorting. (The idea came in part from [one of the papers](https://epubs.siam.org/doi/10.1137/19M1257780) cited in that thread )

> [@ivanpribec](#):
>
> Especially for global sums across multiple ranks (MPI) or images (coarrays), enhanced accurary is needed to get consistent result

This is extremely important! A key use-case is the dot product required in iterative solvers that has to be synchronized at each iteration. I’ll take a read to the paper!! Thanks 🙂 … I could imagine that a vectorized approach could also be extended to the reduction step for distributed sum/dot product to be as consistent as possible.

---

<div class="post-metadata">

**Author:** ![aerosayan](https://avatars.discourse-cdn.com/v4/letter/a/8c91f0/32.png) [@aerosayan](https://fortran-lang.discourse.group/u/aerosayan)\
**Post date:** [October 2, 2023, 7:47am UTC](https://fortran-lang.discourse.group/t/fast-math-a-collection-of-functions-for-fast-number-crunching-using-fortran/6589/8 "2023-10-02T07:47:30Z")

</div>

> [@rwmsu](#):
>
> I guess the question is why are the compiler vendors still using slower and less accurate methods for their intrinsic functions when its been demonstrated that faster and more accurate methods exist.

I personally don’t mind with some of the math operations being little bit slow, as long as they produce acceptable results.

Most, if not all, numerical codes are heavily memory bound, and the future trend in improving performance of numerical codes, would be heavily biased towards increasing throughput of the codes with parallelization, rather than reducing latency of codes running on each core.

---

<div class="post-metadata">

**Author:** ![hkvzjal](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/hkvzjal/32/3055_2.png) [@hkvzjal](https://fortran-lang.discourse.group/u/hkvzjal)\
**Post date:** [October 2, 2023, 11:23am UTC](https://fortran-lang.discourse.group/t/fast-math-a-collection-of-functions-for-fast-number-crunching-using-fortran/6589/9 "2023-10-02T11:23:36Z")

</div>

Here a challenge if anyone has an idea:

The frsqrt function in the library based on the [Quake III arena algorithm](https://en.wikipedia.org/wiki/Fast_inverse_square_root) is actually not that fast compared to what the compiler can do (for single precision, I see performances between 0.9 and 1.2X 🤷‍♂️) . Looking a bit through the net found this repo: [GitHub - Xerbo/simd\_fastinvsqrt: SIMD (SSE) implementation of the infamous Fast Inverse Square Root algorithm from Quake III Arena.](https://github.com/Xerbo/simd_fastinvsqrt) that compares some versions and shows that there is an intrinsic SSE implementation

```auto
_mm_rsqrt_ps(x)

```

which is the fastest option.

Reading around I do not find ways to call such intrinsic directly in Fortran. Are you aware of any way to achieve that efficiently? Or by any chance are the compilers already detecting a call to 1/sqrt() and replacing for that?

---

<div class="post-metadata">

**Author:** ![nshaffer](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/nshaffer/32/161_2.png) [@nshaffer](https://fortran-lang.discourse.group/u/nshaffer)\
**Post date:** [October 3, 2023, 12:35am UTC](https://fortran-lang.discourse.group/t/fast-math-a-collection-of-functions-for-fast-number-crunching-using-fortran/6589/10 "2023-10-03T00:35:30Z")

</div>

You may be interested in my own little study of the inverse square root. The thread is [A (mis)adventure in optimization](https://fortran-lang.discourse.group/t/a-mis-adventure-in-optimization/551). In that case the Quake III rsqrt was both far less accurate and usually slower than the SSE instruction.

---

<div class="post-metadata">

**Author:** ![FedericoPerini](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/federicoperini/32/1750_2.png) [@FedericoPerini](https://fortran-lang.discourse.group/u/FedericoPerini)\
**Post date:** [October 3, 2023, 6:31am UTC](https://fortran-lang.discourse.group/t/fast-math-a-collection-of-functions-for-fast-number-crunching-using-fortran/6589/11 "2023-10-03T06:31:29Z")

</div>

Nice work @nshaffer - yes the advantage for some of those function is not so clear anymore on modern CPUs. [fast 1/sqrt](https://godbolt.org/z/jcdvn6a8c) does compare pretty favorably(~2x speed-up) if AVX is turned on though. I would bet it is due to using `fma` in the last operation:

```auto
vfnmadd213ss xmm0,xmm1,DWORD PTR [rip+0x555] # 40219c <options.19.2+0x6c>

```

---

<div class="post-metadata">

**Author:** ![hkvzjal](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/hkvzjal/32/3055_2.png) [@hkvzjal](https://fortran-lang.discourse.group/u/hkvzjal)\
**Post date:** [October 3, 2023, 7:24am UTC](https://fortran-lang.discourse.group/t/fast-math-a-collection-of-functions-for-fast-number-crunching-using-fortran/6589/12 "2023-10-03T07:24:47Z")

</div>

Nice work @nshaffer !! So you motivated me to push a bit more and reworked the test and found interesting results:

> **code + testing**
>
> ```auto
> module fast_rsqrt
> use, intrinsic :: iso_fortran_env, only: sp=>real32, dp=>real64
> implicit none
> private
>     
> public :: frsqrt
> 
> interface frsqrt
> !! Retranscript of the original Quake III Arena, see https://en.wikipedia.org/wiki/Fast_inverse_square_root
> !! for pure reference
> module procedure frsqrt_dp
> module procedure frsqrt_sp
> end interface
>     
> contains
> 
> elemental function frsqrt_sp(x) result(y)
> integer, parameter :: wp = sp
> real(wp), intent(in) :: x
> real(wp) :: y
> !-- Internal Variables
> real(wp) :: x2,y2
> integer(wp) :: i
> integer(wp), parameter :: magic = int(Z'5f3759df',kind=wp)
> !-------------------------------------------------
> x2 = 0.5_wp*x
> i = transfer(x,i)
> i = magic - shiftr(i,1)
> y2 = transfer(i,y)
> 
> ! Perform one Newton-Raphson step
> y = y2*(1.5_wp-x2*y2*y2)
> 
> end function frsqrt_sp
> 
> elemental function frsqrt_dp(x) result(y)
> integer, parameter :: wp = dp
> real(wp), intent(in) :: x
> real(wp) :: y
> !-- Internal Variables
> real(wp) :: x2,y2
> integer(wp) :: i
> integer(wp), parameter :: magic = 6910469410427058089_wp
> !-------------------------------------------------
> x2 = 0.5_wp*x
> i = transfer(x,i)
> i = magic - shiftr(i,1)
> y2 = transfer(i,y)
> 
> ! Perform one Newton-Raphson step
> y = y2*(1.5_wp-x2*y2*y2)
> 
> end function frsqrt_dp
>     
> end module fast_rsqrt
> 
> program main
> use, intrinsic :: iso_fortran_env, only: sp=>real32, dp=>real64
> use fast_rsqrt
> !> Internal parameters and variables
> integer, parameter :: n = 1e5, ncalc = 2, niter = 10000
> integer :: iter, i
> real(dp) :: time(0:ncalc), times_tot(ncalc), err
> 1 format(a10,': <time> = ',f9.4,' ns/eval, speed-up=',f5.2,'X, rel. error=',es16.4)
> !====================================================================================
> print *,""
> print *,"================================================================"
> print *," Fast rsqrt = 1/sqrt(x)"
> block
> integer, parameter :: wp=sp
> real(wp), allocatable :: x(:) , y(:), yref(:)
> real(kind=wp) :: tolerance = 1e-2_wp
> !> define a log space
> allocate( x(n) , y(n), yref(n) )
> call random_number(x)
> x = 10._wp**(-10*(1._wp-x) + 10*x)
> times_tot(:) = 0
> do iter=1,niter
> call cpu_time(time(0))
> yref = 1._wp/sqrt(x); call cpu_time(time(1))
> y = frsqrt(x) ; call cpu_time(time(2))
> times_tot(:) = times_tot(:) + time(1:ncalc) - time(0:ncalc-1)
> end do
> err = sqrt( sum(( y - yref )**2) ) / sqrt( sum(( yref )**2) )
> print *, "time default: ",times_tot(1),"time frsqrt: ",times_tot(2)
> print 1, "frsqrt r32" , 1e9*times_tot(2)/(niter*n), times_tot(1)/times_tot(2), err
> end block
> block
> integer, parameter :: wp=dp
> real(wp), allocatable :: x(:) , y(:), yref(:)
> real(kind=wp) :: tolerance = 1e-2_wp
> !> define a log space
> allocate( x(n) , y(n), yref(n) )
> call random_number(x)
> x = 10._wp**(-200*(1._wp-x) + 200*x)
> times_tot(:) = 0
> do iter=1,niter
> call cpu_time(time(0))
> yref = 1._wp/sqrt(x); call cpu_time(time(1))
> y = frsqrt(x) ; call cpu_time(time(2))
> times_tot(:) = times_tot(:) + time(1:ncalc) - time(0:ncalc-1)
> end do
> err = sqrt( sum(( y - yref )**2) ) / sqrt( sum(( yref )**2) )
> print *, "time default: ",times_tot(1),"time frsqrt: ",times_tot(2)
> print 1, "frsqrt r64" , 1e9*times_tot(2)/(niter*n), times_tot(1)/times_tot(2), err
> end block
> 
> end program
> 
> ```

Results from [compiler explorer](https://godbolt.org/z/fMPGYfYWK)

ifort 19 -O3 -xHost

```auto
 time default: 0.112649999999986 time frsqrt: 0.111316999999994     
frsqrt r32: <time> = 0.1113 ns/eval, speed-up= 1.01X, rel. error= 1.1156E-03
 time default: 1.49099300000007 time frsqrt: 0.470805999999973     
frsqrt r64: <time> = 0.4708 ns/eval, speed-up= 3.17X, rel. error= 1.0753E-03

```

ifx2023 -O3 -xHost

```auto
 time default: 0.322242999999995 time frsqrt: 0.317305999999985     
frsqrt r32: <time> = 0.3173 ns/eval, speed-up= 1.02X, rel. error= 1.1156E-03
 time default: 2.53335500000001 time frsqrt: 0.726347999999955     
frsqrt r64: <time> = 0.7263 ns/eval, speed-up= 3.49X, rel. error= 1.0753E-03

```

gfortran 13.2 -Ofast -march=native

```auto
 time default: 0.26977699999999727 time frsqrt: 0.28963799999999806     
frsqrt r32: <time> = 0.2896 ns/eval, speed-up= 0.93X, rel. error= 1.1054E-03
 time default: 1.8799389999999838 time frsqrt: 0.69013799999999925     
frsqrt r64: <time> = 0.6901 ns/eval, speed-up= 2.72X, rel. error= 1.0654E-03

```

gfortran 13.2 -O3 -march=native

```auto
time default: 0.30656400000004769 time frsqrt: 0.13649199999999095     
frsqrt r32: <time> = 0.1365 ns/eval, speed-up= 2.25X, rel. error= 1.1029E-03
 time default: 0.94306099999994464 time frsqrt: 0.27456500000000195     
frsqrt r64: <time> = 0.2746 ns/eval, speed-up= 3.43X, rel. error= 1.0285E-03

```

same tendency with my local machine

So basically:

The double precision frsqrt version is consistently \>2X faster

For the single precision the story is different, with ifort or ifx, they already provide the SSE instructions with -O3

ifort

```auto
        vrsqrtps ymm5, ymm2 

```

ifx:

```auto
        vrsqrt14ps ymm1, ymm0

```

for gfortran one has to compile with -Ofast to also get the fastest intrinsic, otherwise with -O3, frsqrt is still faster

gfortran with -Ofast

```auto
        vrsqrtps ymm1, YMMWORD PTR [rdi+rax]

```

I guess it is instructive to have it there if only just as fun experiment and for learning 🙂

---

<div class="post-metadata">

**Author:** ![FedericoPerini](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/federicoperini/32/1750_2.png) [@FedericoPerini](https://fortran-lang.discourse.group/u/FedericoPerini)\
**Post date:** [October 3, 2023, 1:05pm UTC](https://fortran-lang.discourse.group/t/fast-math-a-collection-of-functions-for-fast-number-crunching-using-fortran/6589/13 "2023-10-03T13:05:45Z")

</div>

Given the speed-up, I think it would make sense to [run at least 2 Newton steps](https://godbolt.org/z/97fb4sxMd) in double precision, there’s minimal penalty but a significant accuracy improvement