diff --git a/.github/actions/generate_pyccel_config/action.yml b/.github/actions/generate_pyccel_config/action.yml
index 33ae76dc..fba47e02 100644
--- a/.github/actions/generate_pyccel_config/action.yml
+++ b/.github/actions/generate_pyccel_config/action.yml
@@ -13,11 +13,11 @@ runs:
       run: |
         import json
         config = json.load(open('pyccel_fortran.json'))
-        config['general_flags'].extend(['-O3', '-march=native', '-mtune=native', '-mavx', '-ffast-math'])
+        config['general_flags'].extend(['-O3', '-march=native', '-mtune=native', '-mavx'])
         print(json.dumps(config, indent=4),
                 file=open('pyccel_fortran.json','w'))
         config = json.load(open('pyccel_c.json'))
-        config['general_flags'].extend(['-O3', '-march=native', '-mtune=native', '-mavx', '-ffast-math'])
+        config['general_flags'].extend(['-O3', '-march=native', '-mtune=native', '-mavx'])
         print(json.dumps(config, indent=4),
                 file=open('pyccel_c.json','w'))
       shell: python
diff --git a/benchmarks/pythran.config b/benchmarks/pythran.config
index d5cb0e46..137958eb 100644
--- a/benchmarks/pythran.config
+++ b/benchmarks/pythran.config
@@ -4,7 +4,7 @@ undefs=
 include_dirs=
 libs=
 library_dirs=
-cflags=-std=c++11 -fno-math-errno -fvisibility=hidden -fno-wrapv -Wno-unused-function -Wno-int-in-bool-context -Wno-unknown-warning-option -O3 -march=native -mtune=native -mavx -ffast-math
+cflags=-std=c++11 -fno-math-errno -fvisibility=hidden -fno-wrapv -Wno-unused-function -Wno-int-in-bool-context -Wno-unknown-warning-option -O3 -march=native -mtune=native -mavx
 ldflags=-fvisibility=hidden -Wl,-strip-all
 blas=blas
 CC=
diff --git a/benchmarks/tests/numba_ackermann_mod.py b/benchmarks/tests/numba_ackermann_mod.py
index f5117b59..be3361c4 100644
--- a/benchmarks/tests/numba_ackermann_mod.py
+++ b/benchmarks/tests/numba_ackermann_mod.py
@@ -7,7 +7,7 @@
 """
 from numba import njit
 
-@njit(fastmath=True)
+@njit
 def ackermann(m : int, n : int) -> int:
     """  Total computable function that is not primitive recursive.
     This function is useful for testing recursion
diff --git a/benchmarks/tests/numba_bellman_ford_mod.py b/benchmarks/tests/numba_bellman_ford_mod.py
index feab5003..aef38278 100644
--- a/benchmarks/tests/numba_bellman_ford_mod.py
+++ b/benchmarks/tests/numba_bellman_ford_mod.py
@@ -10,7 +10,7 @@
 import numpy as np
 
 
-@njit(fastmath=True)
+@njit
 def bellman_ford ( v_num: int, e_num: int, source: int, e: 'int[:,:]', e_weight: 'real[:]',
                    v_weight: 'real[:]', predecessor: 'int[:]' ):
     """ Calculate the shortest paths from a source vertex to all other
@@ -49,8 +49,7 @@ def bellman_ford ( v_num: int, e_num: int, source: int, e: 'int[:,:]', e_weight:
     return 0
 
 
-@njit(fastmath=True)
-
+@njit
 def bellman_ford_test():
     """ Test bellman ford's algorithm
     """
diff --git a/benchmarks/tests/numba_dijkstra.py b/benchmarks/tests/numba_dijkstra.py
index c3fb7670..b709f37c 100644
--- a/benchmarks/tests/numba_dijkstra.py
+++ b/benchmarks/tests/numba_dijkstra.py
@@ -9,7 +9,7 @@
 import numpy as np
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def find_nearest ( nv: int, mind: 'int[:]', connected: 'bool[:]' ):
     """ Find the nearest node
     """
@@ -26,7 +26,7 @@ def find_nearest ( nv: int, mind: 'int[:]', connected: 'bool[:]' ):
     return d, v
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def update_mind ( nv: int, mv: int, connected: 'bool[:]', ohd: 'int[:,:]', mind: 'int[:]' ):
     """ Update the minimum distance
     """
@@ -39,7 +39,7 @@ def update_mind ( nv: int, mv: int, connected: 'bool[:]', ohd: 'int[:,:]', mind:
                 mind[i] = min ( mind[i], mind[mv] + ohd[mv,i] )
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def dijkstra_distance ( nv: int, ohd: 'int[:,:]', mind: 'int[:]' ):
     """ Find the shortest paths between nodes in a graph
     """
@@ -76,7 +76,7 @@ def dijkstra_distance ( nv: int, ohd: 'int[:,:]', mind: 'int[:]' ):
         update_mind ( nv, mv, connected, ohd, mind )
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def init ( nv: int, ohd: 'int[:,:]' ):
     """ Create a graph
     """
@@ -92,7 +92,7 @@ def init ( nv: int, ohd: 'int[:,:]' ):
     ohd[0,333] = 33
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def dijkstra_distance_test ( ):
     """ Test Dijkstra's algorithm
     """
diff --git a/benchmarks/tests/numba_euler_mod.py b/benchmarks/tests/numba_euler_mod.py
index a873b257..3d79ac22 100644
--- a/benchmarks/tests/numba_euler_mod.py
+++ b/benchmarks/tests/numba_euler_mod.py
@@ -11,7 +11,7 @@
 import numpy as np
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def euler(dydt: '()(real, const real[:], real[:])',
           tspan: 'real[:]', y0: 'real[:]', n: int,
           t: 'real[:]', y: 'real[:,:]'):
@@ -29,7 +29,7 @@ def euler(dydt: '()(real, const real[:], real[:])',
         y[i+1,:] = y[i,:] + dt * y[i+1,:]
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def humps_fun(x: float):
     """
     Humps function
@@ -42,7 +42,7 @@ def humps_fun(x: float):
     return y
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def humps_deriv(x: 'real', y: 'real[:]', out: 'real[:]'):
     """
     Derivative of the humps function
@@ -51,7 +51,7 @@ def humps_deriv(x: 'real', y: 'real[:]', out: 'real[:]'):
     out[0] = - 2.0 * ( x - 0.3 ) / ( ( x - 0.3 )**2 + 0.01 )**2 - 2.0 * ( x - 0.9 ) / ( ( x - 0.9 )**2 + 0.04 )**2
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def euler_humps_test(t0: float, t1: float, n: int):
     """
     Compute an approximate solution y_h(t) ~= y(t) of the initial
diff --git a/benchmarks/tests/numba_laplace_2d_mod.py b/benchmarks/tests/numba_laplace_2d_mod.py
index a475c65b..afd88dc6 100644
--- a/benchmarks/tests/numba_laplace_2d_mod.py
+++ b/benchmarks/tests/numba_laplace_2d_mod.py
@@ -11,7 +11,7 @@
 from numba import njit
 import numpy as np
 
-@njit(fastmath=True)
+@njit
 def laplace_2d(nx: int, ny: int, rtol: float, maxiter: int):
     """
     Solve the 2D Laplace equation for phi(x, y) on the rectangular
diff --git a/benchmarks/tests/numba_linearconv_1d_mod.py b/benchmarks/tests/numba_linearconv_1d_mod.py
index 0baa8e5b..eccf2eb2 100644
--- a/benchmarks/tests/numba_linearconv_1d_mod.py
+++ b/benchmarks/tests/numba_linearconv_1d_mod.py
@@ -11,7 +11,7 @@
 from numba import njit
 import numpy as np
 
-@njit(fastmath=True)
+@njit
 def linearconv_1d(nx: int, dt: float, nt: int):
     """
     Compute an approximation of the solution u(t, x) to the 1D
diff --git a/benchmarks/tests/numba_md_mod.py b/benchmarks/tests/numba_md_mod.py
index f228bbda..013c8f5a 100644
--- a/benchmarks/tests/numba_md_mod.py
+++ b/benchmarks/tests/numba_md_mod.py
@@ -14,7 +14,7 @@
 from numpy import sin
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def compute_kinetic_energy(vel: 'double[:,:]', mass: float):
     """ Compute the kinetic energy associated with the current configuration.
     """
@@ -28,7 +28,7 @@ def compute_kinetic_energy(vel: 'double[:,:]', mass: float):
     return 0.5 * mass * kinetic
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def compute(mass: float, pos: 'double[:,:]', vel: 'double[:,:]',
         force: 'double[:,:]'):
     """
@@ -71,7 +71,7 @@ def compute(mass: float, pos: 'double[:,:]', vel: 'double[:,:]',
     return potential
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def update(dt: float, mass: float, force: 'double[:,:]',
            pos: 'double[:,:]', vel: 'double[:,:]', acc: 'double[:,:]'):
 
@@ -93,7 +93,7 @@ def update(dt: float, mass: float, force: 'double[:,:]',
     acc[:] = force * rmass
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def r8mat_uniform_ab(r: 'double[:, :]', a: float, b: float, seed: int):
     """ Fill r with random numbers with a uniform distribution
     """
@@ -121,7 +121,7 @@ def r8mat_uniform_ab(r: 'double[:, :]', a: float, b: float, seed: int):
     return seed
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def initialize(pos: 'double[:,:]'):
     """ Initialise the positions of the particles
     """
@@ -130,7 +130,7 @@ def initialize(pos: 'double[:,:]'):
     seed = r8mat_uniform_ab(pos, 0.0, 10.0, seed)
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def md(d_num: int, p_num: int, step_num: int, dt: float):
     """
     Run a molecular dynamics simulation. This consists of an N-body
diff --git a/benchmarks/tests/numba_midpoint_explicit_mod.py b/benchmarks/tests/numba_midpoint_explicit_mod.py
index e2c44339..0430406f 100644
--- a/benchmarks/tests/numba_midpoint_explicit_mod.py
+++ b/benchmarks/tests/numba_midpoint_explicit_mod.py
@@ -11,7 +11,7 @@
 import numpy as np
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def midpoint_explicit(dydt: '()(real, const real[:], real[:])',
                       tspan: 'real[:]', y0: 'real[:]', n: int,
                       t: 'real[:]', y: 'real[:,:]'):
@@ -38,7 +38,7 @@ def midpoint_explicit(dydt: '()(real, const real[:], real[:])',
         y[i+1,:] = y[i,:] + dt * y[i+1,:]
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def humps_fun(x: float):
     """
     Humps function
@@ -51,7 +51,7 @@ def humps_fun(x: float):
     return y
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def humps_deriv(x: 'real', y: 'real[:]', out: 'real[:]'):
     """
     Derivative of the humps function
@@ -60,7 +60,7 @@ def humps_deriv(x: 'real', y: 'real[:]', out: 'real[:]'):
     out[0] = - 2.0 * ( x - 0.3 ) / ( ( x - 0.3 )**2 + 0.01 )**2 - 2.0 * ( x - 0.9 ) / ( ( x - 0.9 )**2 + 0.04 )**2
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def midpoint_explicit_humps_test(t0: float, t1: float, n: int):
     """
     Compute an approximate solution y_h(t) ~= y(t) of the initial
diff --git a/benchmarks/tests/numba_midpoint_fixed_mod.py b/benchmarks/tests/numba_midpoint_fixed_mod.py
index fd4cd3f3..474c87c4 100644
--- a/benchmarks/tests/numba_midpoint_fixed_mod.py
+++ b/benchmarks/tests/numba_midpoint_fixed_mod.py
@@ -11,7 +11,7 @@
 import numpy as np
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def midpoint_fixed(dydt: '()(real, const real[:], real[:])',
                    tspan: 'real[:]', y0: 'real[:]', n: int,
                    t: 'real[:]', y: 'real[:,:]'):
@@ -45,7 +45,7 @@ def midpoint_fixed(dydt: '()(real, const real[:], real[:])',
                  + (1.0 - 1.0 / theta) * y[i,:]
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def humps_fun(x: float):
     """
     Humps function
@@ -58,7 +58,7 @@ def humps_fun(x: float):
     return y
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def humps_deriv(x: 'real', y: 'real[:]', out: 'real[:]'):
     """
     Derivative of the humps function
@@ -67,7 +67,7 @@ def humps_deriv(x: 'real', y: 'real[:]', out: 'real[:]'):
     out[0] = - 2.0 * ( x - 0.3 ) / ( ( x - 0.3 )**2 + 0.01 )**2 - 2.0 * ( x - 0.9 ) / ( ( x - 0.9 )**2 + 0.04 )**2
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def midpoint_fixed_humps_test(t0: float, t1: float, n: int):
     """
     Compute an approximate solution y_h(t) ~= y(t) of the initial
diff --git a/benchmarks/tests/numba_nonlinearconv_1d_mod.py b/benchmarks/tests/numba_nonlinearconv_1d_mod.py
index 7be6934c..c20d9ad9 100644
--- a/benchmarks/tests/numba_nonlinearconv_1d_mod.py
+++ b/benchmarks/tests/numba_nonlinearconv_1d_mod.py
@@ -12,7 +12,7 @@
 import numpy as np
 
 
-@njit(fastmath=True)
+@njit
 def nonlinearconv_1d(nx: int, dt: float, nt: int):
     """
     Compute an approximation of the solution u(t, x) to the 1D
diff --git a/benchmarks/tests/numba_poisson_2d_mod.py b/benchmarks/tests/numba_poisson_2d_mod.py
index 8bec8f5a..8112a437 100644
--- a/benchmarks/tests/numba_poisson_2d_mod.py
+++ b/benchmarks/tests/numba_poisson_2d_mod.py
@@ -12,7 +12,7 @@
 import numpy as np
 
 
-@njit(fastmath=True)
+@njit
 def poisson_2d(nx: int, ny: int, nt: int):
     """
     Solve the 2D poisson equation for phi(x, y) on the rectangular
diff --git a/benchmarks/tests/numba_rk4_mod.py b/benchmarks/tests/numba_rk4_mod.py
index 28621371..ce520b0e 100644
--- a/benchmarks/tests/numba_rk4_mod.py
+++ b/benchmarks/tests/numba_rk4_mod.py
@@ -11,7 +11,7 @@
 import numpy as np
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def rk4(dydt: '()(real, const real[:], real[:])',
         tspan: 'real[:]', y0: 'real[:]', n: int,
         t: 'real[:]', y: 'real[:,:]'):
@@ -43,7 +43,7 @@ def rk4(dydt: '()(real, const real[:], real[:])',
         y[i+1,:] = y[i,:] + dt * (f1[:] + 2.0 * f2[:] + 2.0 * f3[:] + f4[:]) / 6.0
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def humps_fun(x: float):
     """
     Humps function
@@ -56,7 +56,7 @@ def humps_fun(x: float):
     return y
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def humps_deriv(x: 'real', y: 'real[:]', out: 'real[:]'):
     """
     Derivative of the humps function
@@ -65,7 +65,7 @@ def humps_deriv(x: 'real', y: 'real[:]', out: 'real[:]'):
     out[0] = - 2.0 * ( x - 0.3 ) / ( ( x - 0.3 )**2 + 0.01 )**2 - 2.0 * ( x - 0.9 ) / ( ( x - 0.9 )**2 + 0.04 )**2
 
 # ================================================================
-@njit(fastmath=True)
+@njit
 def rk4_humps_test(t0: float, t1: float, n: int):
     """
     Compute an approximate solution y_h(t) ~= y(t) of the initial