Repository: pyccel/pyccel (python)

Issue:
[Feature Request] Implement `expm1` function
### Is there an existing issue for this?

- [x] I have searched the existing issues

### Feature Description

`expm1` is a function to accurately compute values for `exp(x) - 1` which is important for small values of `x` (round-off errors). According to ChatGPT, both Fortran and C support this function.

### Use Case

I could compute values of `exp(x) - 1` more accurately for small `x`

### Benefits

_No response_

### Add Screenshots

_No response_

### Priority

High

### Record

- [x] I have read the [Contributing Guidelines](https://github.com/pyccel/pyccel/blob/devel/docs/CONTRIBUTING.md).
- [x] I have starred the repository

Pull request:
Add support for `numpy.expm1`
Add support for `numpy.expm1`. Fixes #2390.
NumPy's algorithm is used for complex support. The math library `-lm` is used for Fortran support.

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index df10e7ef..f6a0c5ab 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,6 +10,7 @@ All notable changes to this project will be documented in this file.
 -   #1742 : Add support for set method `difference_update()`.
 -   #1849 : Add support for lambda functions in assign statements by treating them as inline functions.
 -   #1585 : Add support for `np.divide` and its alias `np.true_divide`.
+-   #2390 : Add support for `np.expm1`.
 
 ### Fixed
 
diff --git a/docs/numpy-functions.md b/docs/numpy-functions.md
index b0cea8b7..84cf5d21 100644
--- a/docs/numpy-functions.md
+++ b/docs/numpy-functions.md
@@ -432,8 +432,8 @@ In Pyccel we try to support the NumPy functions which developers use the most..
 
 -   Supported [math functions](https://numpy.org/doc/stable/reference/routines.math.html) (optional parameters are not supported):
 
-    `sqrt`, `abs`, `sin`, `cos`, `exp`, `log`, `tan`, `arcsin`, `arccos`, `arctan`, `arctan2`, `sinh`, `cosh`, `tanh`, `arcsinh`, `arccosh` and
-    `arctanh`, `true_divide`, `divide`.
+    `sqrt`, `abs`, `sin`, `cos`, `exp`, `expm1`, `log`, `tan`, `arcsin`, `arccos`, `arctan`, `arctan2`, `sinh`, `cosh`, `tanh`,
+    `arcsinh`, `arccosh`, `arctanh`, `true_divide`, `divide`.
 
 -   Supported [logic functions](https://numpy.org/doc/stable/reference/routines.logic.html) (optional parameters are not supported):
 
diff --git a/pyccel/ast/numpyext.py b/pyccel/ast/numpyext.py
index a581fbce..791b4b09 100644
--- a/pyccel/ast/numpyext.py
+++ b/pyccel/ast/numpyext.py
@@ -79,6 +79,7 @@ __all__ = (
     'NumpyCosh',
     'NumpyDivide',
     'NumpyExp',
+    'NumpyExpm1',
     'NumpyFabs',
     'NumpyFloor',
     'NumpyHypot',
@@ -1953,6 +1954,19 @@ class NumpyExp     (NumpyUfuncUnary):
     """Represent a call to the exp function in the Numpy library"""
     __slots__ = ()
     name = 'exp'
+class NumpyExpm1   (NumpyUfuncUnary):
+    """
+    Represent a call to the np.expm1 function in the Numpy library.
+
+    Represent a call to the np.expm1 function in the Numpy library.
+
+    Parameters
+    ----------
+    x : PyccelAstType
+        The argument of the unary function.
+    """
+    __slots__ = ()
+    name = 'expm1'
 class NumpyLog     (NumpyUfuncUnary):
     """Represent a call to the log function in the Numpy library"""
     __slots__ = ()
@@ -2878,6 +2892,7 @@ numpy_funcs = {
     'absolute'  : PyccelFunctionDef('absolute'  , NumpyAbs),
     'fabs'      : PyccelFunctionDef('fabs'      , NumpyFabs),
     'exp'       : PyccelFunctionDef('exp'       , NumpyExp),
+    'expm1'     : PyccelFunctionDef('expm1'     , NumpyExpm1),
     'log'       : PyccelFunctionDef('log'       , NumpyLog),
     'sqrt'      : PyccelFunctionDef('sqrt'      , NumpySqrt),
     # ---
diff --git a/pyccel/codegen/printing/ccode.py b/pyccel/codegen/printing/ccode.py
index ff6659ef..1454007a 100644
--- a/pyccel/codegen/printing/ccode.py
+++ b/pyccel/codegen/printing/ccode.py
@@ -2141,6 +2141,20 @@ class CCodePrinter(CodePrinter):
         code_arg = self._print(expr.arg)
         return f"isnan({code_arg})"
 
+    def _print_NumpyExpm1(self, expr):
+        arg, = expr.args
+        if expr.dtype.primitive_type is PrimitiveComplexType():
+            self.add_import(c_imports['pyc_math_c'])
+            arg_code = self._print(arg)
+            return f'cpyc_expm1({arg_code})'
+        else:
+            self.add_import(c_imports['math'])
+            if arg.dtype.primitive_type is not PrimitiveFloatingPointType():
+                arg_code = self._print(NumpyFloat(arg))
+            else :
+                arg_code = self._print(arg)
+            return f'expm1({arg_code})'
+
     def _print_MathFunctionBase(self, expr):
         """ Convert a Python expression with a math function call to C
         function call
diff --git a/pyccel/codegen/printing/fcode.py b/pyccel/codegen/printing/fcode.py
index d436e3f1..5b7ad779 100644
--- a/pyccel/codegen/printing/fcode.py
+++ b/pyccel/codegen/printing/fcode.py
@@ -74,7 +74,7 @@ from pyccel.ast.numpyext import NumpySign
 from pyccel.ast.numpyext import NumpyIsFinite, NumpyIsNan
 from pyccel.ast.numpyext import get_shape_of_multi_level_container
 
-from pyccel.ast.numpytypes import NumpyNDArrayType, NumpyInt64Type
+from pyccel.ast.numpytypes import NumpyNDArrayType, NumpyInt64Type, NumpyFloat64Type, NumpyComplex128Type
 
 from pyccel.ast.operators import PyccelAdd, PyccelMul, PyccelMinus, PyccelAnd, PyccelEq
 from pyccel.ast.operators import PyccelMod, PyccelNot, PyccelAssociativeParenthesis
@@ -3291,6 +3291,21 @@ class FCodePrinter(CodePrinter):
             abs_result = self._print(self._apply_cast(expr.dtype, PythonBool(arg)))
             return f'sign({abs_result}, {arg_code})'
 
+    def _print_NumpyExpm1(self, expr):
+        arg = expr.arg
+        if isinstance(expr.dtype.primitive_type, PrimitiveComplexType):
+            if not isinstance(arg.dtype, NumpyComplex128Type):
+                arg = self._apply_cast(NumpyComplex128Type(), arg)
+        else:
+            if not isinstance(arg.dtype, NumpyFloat64Type):
+                arg = self._apply_cast(NumpyFloat64Type(), arg)
+
+        self.add_import(Import('pyc_math_f90', Module('pyc_math_f90',(),())))
+
+        arg_code = self._print(arg)
+
+        return f'expm1({arg_code})'
+
     def _print_NumpyTranspose(self, expr):
         var = expr.internal_var
         arg = self._print(var)
diff --git a/pyccel/codegen/utilities.py b/pyccel/codegen/utilities.py
index a8f504f8..6aec7b3f 100644
--- a/pyccel/codegen/utilities.py
+++ b/pyccel/codegen/utilities.py
@@ -43,7 +43,7 @@ language_extension = {'fortran':'f90', 'c':'c', 'python':'py'}
 # map internal libraries to their folders inside pyccel/stdlib and their compile objects
 # The compile object folder will be in the pyccel dirpath
 internal_libs = {
-    "pyc_math_f90"     : (stdlib_path / "math", "math", CompileObj("pyc_math_f90.f90",folder="math")),
+    "pyc_math_f90"     : (stdlib_path / "math", "math", CompileObj("pyc_math_f90.f90",folder="math", libs = ('m',))),
     "pyc_math_c"       : (stdlib_path / "math", "math", CompileObj("pyc_math_c.c",folder="math")),
     "pyc_tools_f90"    : (stdlib_path / "tools", "tools", CompileObj("pyc_tools_f90.f90",folder="tools")),
     "cwrapper"         : (stdlib_path / "cwrapper", "cwrapper", CompileObj("cwrapper.c",folder="cwrapper", accelerators=('python',))),
diff --git a/pyccel/naming/cnameclashchecker.py b/pyccel/naming/cnameclashchecker.py
index 22725f8b..fbd02984 100644
--- a/pyccel/naming/cnameclashchecker.py
+++ b/pyccel/naming/cnameclashchecker.py
@@ -34,7 +34,7 @@ class CNameClashChecker(LanguageNameClashChecker):
         'array_bool_3d', 'array_float_complex_1d', 'array_float_complex_2d',
         'array_float_complex_3d', 'array_double_complex_1d', 'array_double_complex_2d',
         'array_double_complex_3d', 'c_ALL', 'c_END', 'cspan_slice', 'cspan_transpose',
-        'complex_max', 'complex_min'])
+        'complex_max', 'complex_min', 'expm1', 'complex_expm1'])
 
     def has_clash(self, name, symbols):
         """
diff --git a/pyccel/stdlib/math/pyc_math_c.c b/pyccel/stdlib/math/pyc_math_c.c
index 9094facf..498f7f4c 100644
--- a/pyccel/stdlib/math/pyc_math_c.c
+++ b/pyccel/stdlib/math/pyc_math_c.c
@@ -112,6 +112,11 @@ int64_t ipyc_bankers_round(int64_t arg, int64_t ndigits)
     }
 }
 
+double complex cpyc_expm1(double complex x) {
+    double a = sin(cimag(x) * 0.5);
+    return (expm1(creal(x)) * cos(cimag(x)) - 2. * a * a) + exp(creal(x)) * sin(cimag(x)) * _Complex_I;
+}
+
 extern inline double       pyc_radians(double degrees);
 extern inline double       pyc_degrees(double radians);
 extern inline int64_t      pyc_modulo(int64_t a, int64_t b);
diff --git a/pyccel/stdlib/math/pyc_math_c.h b/pyccel/stdlib/math/pyc_math_c.h
index 62e7db90..195efa69 100644
--- a/pyccel/stdlib/math/pyc_math_c.h
+++ b/pyccel/stdlib/math/pyc_math_c.h
@@ -52,6 +52,8 @@ double complex csign(double complex x);
 double fpyc_bankers_round(double arg, int64_t ndigits);
 int64_t ipyc_bankers_round(int64_t arg, int64_t ndigits);
 
+double complex cpyc_expm1(double complex x);
+
 #define PY_FLOOR_DIV_TYPE(TYPE)                         \
     static inline TYPE py_floor_div_##TYPE(TYPE x, TYPE y) { \
         return (TYPE)(x / y - ((x % y != 0) && ((x < 0) ^ (y < 0)))); \
diff --git a/pyccel/stdlib/math/pyc_math_f90.f90 b/pyccel/stdlib/math/pyc_math_f90.f90
index 6b6c1ae1..308b14a6 100644
--- a/pyccel/stdlib/math/pyc_math_f90.f90
+++ b/pyccel/stdlib/math/pyc_math_f90.f90
@@ -23,7 +23,8 @@ public :: pyc_gcd, &
           csgn, &
           csign, &
           pyc_bankers_round, &
-          pyc_floor_div
+          pyc_floor_div, &
+          expm1
 
 private
 
@@ -76,6 +77,18 @@ interface pyc_floor_div
     module procedure pyc_floor_div_i64
 end interface pyc_floor_div
 
+interface
+   real(f64) pure function c_expm1(x) bind(c, name='expm1')
+     use, intrinsic :: ISO_C_Binding, only : f64 => C_DOUBLE
+     real(f64), intent(in), value :: x
+   end function c_expm1
+end interface
+
+interface expm1
+    module procedure pyc_expm1_f64
+    module procedure pyc_expm1_c64
+end interface expm1
+
 contains
 
 ! Implementation of math factorial function
@@ -454,5 +467,27 @@ elemental pure integer(kind=8) function pyc_floor_div_i64(x, y) result(res)
   res = x / y - merge(1, 0, mod(x, y) /= 0 .and. ((x < 0) .neqv. (y < 0)))
 end function pyc_floor_div_i64
 
+elemental pure function pyc_expm1_f64(x) result(Out_0001)
+    implicit none
+    real(f64) :: Out_0001
+    real(f64), value :: x
+
+    Out_0001 = c_expm1(x)
+
+end function pyc_expm1_f64
+
+elemental pure function pyc_expm1_c64(x) result(Out_0001)
+    implicit none
+
+    complex(c64) :: Out_0001
+    complex(c64), value :: x
+
+    real(f64) :: a
+
+    a = sin(aimag(x) * 0.5_f64)
+    Out_0001 = (c_expm1(real(x)) * cos(aimag(x)) - 2._f64 * a * a) + &
+               (exp(real(x)) * sin(aimag(x))) * cmplx(0,1, kind = c64)
+
+end function pyc_expm1_c64
 
 end module pyc_math_f90


Test patch:
diff --git a/tests/epyccel/recognised_functions/test_numpy_funcs.py b/tests/epyccel/recognised_functions/test_numpy_funcs.py
index 15e4e081..16f98c4b 100644
--- a/tests/epyccel/recognised_functions/test_numpy_funcs.py
+++ b/tests/epyccel/recognised_functions/test_numpy_funcs.py
@@ -771,6 +771,150 @@ def test_exp_phrase_r_i(language):
     assert isclose(f2(x,y), exp_phrase_r_i(x,y), rtol=RTOL, atol=ATOL)
     assert isclose(f2(x,-y), exp_phrase_r_i(x,-y), rtol=RTOL, atol=ATOL)
 
+#--------------------------------- expm1 function ------------------------------#
+def test_expm1_call_i(language):
+    def expm1_call_i(x : 'int'):
+        from numpy import expm1
+        return expm1(x)
+
+    f1 = epyccel(expm1_call_i, language = language)
+    x = randint(100)
+    assert isclose(f1(x), expm1_call_i(x), rtol=RTOL, atol=ATOL)
+    assert isclose(f1(-x), expm1_call_i(-x), rtol=RTOL, atol=ATOL)
+    assert matching_types(f1(x), expm1_call_i(x))
+
+def test_expm1_call_f(language):
+    def expm1_call_f(x : 'float'):
+        from numpy import expm1
+        return expm1(x)
+
+    f1 = epyccel(expm1_call_f, language = language)
+    x = uniform(high=100)
+    assert isclose(f1(x), expm1_call_f(x), rtol=RTOL, atol=ATOL)
+    assert isclose(f1(-x), expm1_call_f(-x), rtol=RTOL, atol=ATOL)
+    assert matching_types(f1(x), expm1_call_f(x))
+
+def test_expm1_call_c(language):
+    def expm1_call_c(x : complex):
+        from numpy import expm1
+        return expm1(x)
+
+    f1 = epyccel(expm1_call_c, language = language)
+    x = uniform(high=100) + uniform(high=100)*1j
+    assert isclose(f1(x), expm1_call_c(x), rtol=RTOL, atol=ATOL)
+    assert isclose(f1(-x), expm1_call_c(-x), rtol=RTOL, atol=ATOL)
+    assert matching_types(f1(x), expm1_call_c(x))
+
+def test_expm1_call_f_array(language):
+    def expm1_call_f(x : 'float[:]'):
+        from numpy import expm1
+        return expm1(x)
+
+    f1 = epyccel(expm1_call_f, language = language)
+    x = uniform(high=100, size=5)
+    assert np.allclose(f1(x), expm1_call_f(x), rtol=RTOL, atol=ATOL)
+    assert np.allclose(f1(-x), expm1_call_f(-x), rtol=RTOL, atol=ATOL)
+
+def test_expm1_call_c_array(language):
+    def expm1_call_c(x : 'complex[:]'):
+        from numpy import expm1
+        return expm1(x)
+
+    f1 = epyccel(expm1_call_c, language = language)
+    x = uniform(high=100, size=5) + uniform(high=100, size=5)*1j
+    assert np.allclose(f1(x), expm1_call_c(x), rtol=RTOL, atol=ATOL)
+    assert np.allclose(f1(-x), expm1_call_c(-x), rtol=RTOL, atol=ATOL)
+
+def test_expm1_call_cast_f(language):
+    def expm1_call_f(x : 'float32'):
+        from numpy import expm1
+        return expm1(x)
+
+    f1 = epyccel(expm1_call_f, language = language)
+    x = np.float32(uniform(high=30))
+    assert isclose(f1(x), expm1_call_f(x), rtol=RTOL32, atol=ATOL32)
+    assert matching_types(f1(x), expm1_call_f(x))
+
+def test_expm1_call_cast_c(language):
+    def expm1_call_c(x : 'complex64'):
+        from numpy import expm1
+        return expm1(x)
+
+    f1 = epyccel(expm1_call_c, language = language)
+    x = np.complex64(uniform(high=15) + uniform(high=15)*1j)
+    assert isclose(f1(x), expm1_call_c(x), rtol=RTOL32, atol=ATOL32)
+    assert matching_types(f1(x), expm1_call_c(x))
+
+def test_expm1_phrase_i_i(language):
+    def expm1_phrase_i_i(x : 'int', y : 'int'):
+        from numpy import expm1
+        a = expm1(x)+expm1(y)
+        return a
+
+    f2 = epyccel(expm1_phrase_i_i, language = language)
+    x = randint(100)
+    y = randint(100)
+    assert isclose(f2(x,y), expm1_phrase_i_i(x,y), rtol=RTOL, atol=ATOL)
+    assert isclose(f2(-x,-y), expm1_phrase_i_i(-x,-y), rtol=RTOL, atol=ATOL)
+    assert isclose(f2(-x,y), expm1_phrase_i_i(-x,y), rtol=RTOL, atol=ATOL)
+    assert isclose(f2(x,-y), expm1_phrase_i_i(x,-y), rtol=RTOL, atol=ATOL)
+
+def test_expm1_phrase_f_f(language):
+    def expm1_phrase_f_f(x : 'float', y : 'float'):
+        from numpy import expm1
+        a = expm1(x)+expm1(y)
+        return a
+
+    f2 = epyccel(expm1_phrase_f_f, language = language)
+    x = uniform(high=100)
+    y = uniform(high=100)
+    assert isclose(f2(x,y), expm1_phrase_f_f(x,y), rtol=RTOL, atol=ATOL)
+    assert isclose(f2(-x,-y), expm1_phrase_f_f(-x,-y), rtol=RTOL, atol=ATOL)
+    assert isclose(f2(-x,y), expm1_phrase_f_f(-x,y), rtol=RTOL, atol=ATOL)
+    assert isclose(f2(x,-y), expm1_phrase_f_f(x,-y), rtol=RTOL, atol=ATOL)
+
+def test_expm1_phrase_i_f(language):
+    def expm1_phrase_i_f(x : 'int', y : 'float'):
+        from numpy import expm1
+        a = expm1(x)+expm1(y)
+        return a
+
+    f2 = epyccel(expm1_phrase_i_f, language = language)
+    x = randint(100)
+    y = uniform(high=100)
+    assert isclose(f2(x,y), expm1_phrase_i_f(x,y), rtol=RTOL, atol=ATOL)
+    assert isclose(f2(-x,-y), expm1_phrase_i_f(-x,-y), rtol=RTOL, atol=ATOL)
+    assert isclose(f2(-x,y), expm1_phrase_i_f(-x,y), rtol=RTOL, atol=ATOL)
+    assert isclose(f2(x,-y), expm1_phrase_i_f(x,-y), rtol=RTOL, atol=ATOL)
+
+def test_expm1_phrase_f_i(language):
+    def expm1_phrase_f_i(x : 'float', y : 'int'):
+        from numpy import expm1
+        a = expm1(x)+expm1(y)
+        return a
+
+    f2 = epyccel(expm1_phrase_f_i, language = language)
+    x = uniform(high=100)
+    y = randint(100)
+    assert isclose(f2(x,y), expm1_phrase_f_i(x,y), rtol=RTOL, atol=ATOL)
+    assert isclose(f2(x,y), expm1_phrase_f_i(x,y), rtol=RTOL, atol=ATOL)
+    assert isclose(f2(x,y), expm1_phrase_f_i(x,y), rtol=RTOL, atol=ATOL)
+    assert isclose(f2(x,-y), expm1_phrase_f_i(x,-y), rtol=RTOL, atol=ATOL)
+
+def test_expm1_phrase_i_c(language):
+    def expm1_phrase_i_c(x : int, y : complex):
+        from numpy import expm1
+        a = expm1(x)+expm1(y)
+        return a
+
+    f2 = epyccel(expm1_phrase_i_c, language = language)
+    x = randint(100)
+    y = uniform(high=100) + uniform(high=100)*1j
+    assert isclose(f2(x,y), expm1_phrase_i_c(x,y), rtol=RTOL, atol=ATOL)
+    assert isclose(f2(x,y), expm1_phrase_i_c(x,y), rtol=RTOL, atol=ATOL)
+    assert isclose(f2(x,y), expm1_phrase_i_c(x,y), rtol=RTOL, atol=ATOL)
+    assert isclose(f2(x,-y), expm1_phrase_i_c(x,-y), rtol=RTOL, atol=ATOL)
+
 #--------------------------------- log function ------------------------------#
 def test_log_call_i(language):
     def log_call_i(x : 'int'):