Repository: pyccel/pyccel
Pull request #2142 (2025-01-10 10:48:34)

improve floordiv translation
Floor division translation to `C` and `Fortran` relies on float division then rounding when dealing with integers. This introduces rounding errors that create discrepancies between `Python` and `Pyccel` translated code.
This PR aims to keep operations over the integers whenever possible.

pyccel-bot[bot]: Hello again! Thank you for this new pull request :star_struck:.

Please begin by requesting your checklist using the command `/bot checklist`

jalalium: /bot checklist
pyccel-bot[bot]: It seems like you have forgotten to complete your checklist for this pull request. Please check the list [here](https://github.com/pyccel/pyccel/pull/2142#issuecomment-2582523497), fix the missing items and let me know when everything is ticked and fixed with `/bot mark as ready`.

pyccel-bot[bot]: Unfortunately your PR is not passing the tests so it is not quite ready for review yet. Let me know when it is fixed with `/bot mark as ready`.

pyccel-bot[bot]: Unfortunately your PR is not passing the tests so it is not quite ready for review yet. Let me know when it is fixed with `/bot mark as ready`.

pyccel-bot[bot]: Unfortunately your PR is not passing the tests so it is not quite ready for review yet. Let me know when it is fixed with `/bot mark as ready`.

pyccel-bot[bot]: Unfortunately your PR is not passing the tests so it is not quite ready for review yet. Let me know when it is fixed with `/bot mark as ready`.

pyccel-bot[bot]: Unfortunately your PR is not passing the tests so it is not quite ready for review yet. Let me know when it is fixed with `/bot mark as ready`.

pyccel-bot[bot]: Hey @pyccel/pyccel-dev ! @jalalium has just created this great new pull request! Check it out and let me know what you think!

pyccel-bot[bot]: @jalalium, @EmilyBourne has a few questions/comments about your code. Can you go through and see if you agree with them. If not go ahead and explain why. Once you've adressed all the comments let me know with `/bot mark as ready` and we will see if we can get approval.


Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 90a013a0..dbdf4a2d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -122,6 +122,7 @@ All notable changes to this project will be documented in this file.
 -   #2175 : Fix the shape of multi-level containers.
 -   Catch all internal errors arising in the syntactic, semantic, printing or code generation stages.
 -   #2206 : Fix returning an array of unknown literal size.
+-   #2112 : Improve floor division.
 
 ### Changed
 
diff --git a/pyccel/codegen/printing/ccode.py b/pyccel/codegen/printing/ccode.py
index 45417c56..cbc9b5f2 100644
--- a/pyccel/codegen/printing/ccode.py
+++ b/pyccel/codegen/printing/ccode.py
@@ -2346,17 +2346,19 @@ class CCodePrinter(CodePrinter):
         return  ' / '.join(self._print(a) for a in args)
 
     def _print_PyccelFloorDiv(self, expr):
-        self.add_import(c_imports['math'])
         # the result type of the floor division is dependent on the arguments
-        # type, if all arguments are integers the result is integer otherwise
-        # the result type is float
-        need_to_cast = all(a.dtype.primitive_type is PrimitiveIntegerType() for a in expr.args)
+        # type, if all arguments are integers or booleans the result is integer
+        # otherwise the result type is float
+        need_to_cast = all(a.dtype.primitive_type in (PrimitiveIntegerType(), PrimitiveBooleanType()) for a in expr.args)
+        if need_to_cast:
+            self.add_import(c_imports['pyc_math_c'])
+            cast_type = self.get_c_type(expr.dtype)
+            return f'py_floor_div_{cast_type}({self._print(expr.args[0])}, {self._print(expr.args[1])})'
+
+        self.add_import(c_imports['math'])
         code = ' / '.join(self._print(a if a.dtype.primitive_type is PrimitiveFloatingPointType()
                                         else NumpyFloat(a)) for a in expr.args)
-        if (need_to_cast):
-            cast_type = self.get_c_type(expr.dtype)
-            return "({})floor({})".format(cast_type, code)
-        return "floor({})".format(code)
+        return f"floor({code})"
 
     def _print_PyccelRShift(self, expr):
         return ' >> '.join(self._print(a) for a in expr.args)
diff --git a/pyccel/codegen/printing/fcode.py b/pyccel/codegen/printing/fcode.py
index 7d794dab..fc6d8414 100644
--- a/pyccel/codegen/printing/fcode.py
+++ b/pyccel/codegen/printing/fcode.py
@@ -3126,19 +3126,12 @@ class FCodePrinter(CodePrinter):
         return code
 
     def _print_PyccelFloorDiv(self, expr):
-
-        code     = self._print(expr.args[0])
-        adtype   = expr.args[0].dtype.primitive_type
-        is_float = isinstance(expr.dtype.primitive_type, PrimitiveFloatingPointType)
-        for b in expr.args[1:]:
-            bdtype    = b.dtype.primitive_type
-            if all(isinstance(dtype, PrimitiveIntegerType) for dtype in (adtype, bdtype)):
-                b = NumpyFloat(b)
-            c = self._print(b)
-            adtype = bdtype
-            code = 'FLOOR({}/{},{})'.format(code, c, self.print_kind(expr))
-            if is_float:
-                code = 'real({}, {})'.format(code, self.print_kind(expr))
+        new_args = [self._apply_cast(expr.dtype, arg) for arg in expr.args]
+        args = [self._print(arg) for arg in new_args]
+        if all(isinstance(arg.dtype.primitive_type, (PrimitiveBooleanType, PrimitiveIntegerType)) for arg in expr.args):
+            self.add_import(Import('pyc_math_f90', Module('pyc_math_f90',(),())))
+            return f'pyc_floor_div({args[0]}, {args[1]})'
+        code = f'real(FLOOR({args[0]} / {args[1]}, {self.print_kind(expr)}), {self.print_kind(expr)})'
         return code
 
     def _print_PyccelRShift(self, expr):
diff --git a/pyccel/stdlib/math/pyc_math_c.c b/pyccel/stdlib/math/pyc_math_c.c
index 9b589e84..9094facf 100644
--- a/pyccel/stdlib/math/pyc_math_c.c
+++ b/pyccel/stdlib/math/pyc_math_c.c
@@ -6,9 +6,9 @@
 #include "pyc_math_c.h"
 
 /*---------------------------------------------------------------------------*/
-int64_t                     pyc_factorial(int64_t n)
+int64_t pyc_factorial(int64_t n)
 {
-    int64_t    res = 1;
+    int64_t res = 1;
 
     /* ValueError: factorial() not defined for negative values */
     if (n < 0)
@@ -18,11 +18,12 @@ int64_t                     pyc_factorial(int64_t n)
     return (res);
 }
 /*---------------------------------------------------------------------------*/
-int64_t                     pyc_gcd (int64_t a, int64_t b)
+int64_t pyc_gcd(int64_t a, int64_t b)
 {
-    a = a > 0 ? a : - a;
-    b = b > 0 ? b : - b;
-    while (b) {
+    a = a > 0 ? a : -a;
+    b = b > 0 ? b : -b;
+    while (b)
+    {
         a %= b;
         /* swap a and b*/
         a = a ^ b;
@@ -32,10 +33,10 @@ int64_t                     pyc_gcd (int64_t a, int64_t b)
     return a;
 }
 /*---------------------------------------------------------------------------*/
-int64_t                     pyc_lcm (int64_t a, int64_t b)
+int64_t pyc_lcm(int64_t a, int64_t b)
 {
-    a = a > 0 ? a : - a;
-    b = b > 0 ? b : - b;
+    a = a > 0 ? a : -a;
+    b = b > 0 ? b : -b;
     return a / pyc_gcd(a, b) * b;
 }
 /*---------------------------------------------------------------------------*/
@@ -43,13 +44,13 @@ int64_t                     pyc_lcm (int64_t a, int64_t b)
 /* numpy.sign for float, double and integers */
 long long int isign(long long int x)
 {
-    return (x>0) - (x<0);
+    return (x > 0) - (x < 0);
 }
 
 /* numpy.sign for float, double and integers */
-double  fsign(double x)
+double fsign(double x)
 {
-    return (double)((x>0) - (x<0));
+    return (double)((x > 0) - (x < 0));
 }
 
 /* numpy.sign for complex for NumPy v1 */
@@ -77,8 +78,9 @@ double fpyc_bankers_round(double arg, int64_t ndigits)
 
     double diff = arg - rnd;
 
-    if (ndigits <= 0 && (diff == 0.5 || diff == -0.5)) {
-        rnd = ((int64_t)(arg*0.5 + nearest_int_fix))*2.0;
+    if (ndigits <= 0 && (diff == 0.5 || diff == -0.5))
+    {
+        rnd = ((int64_t)(arg * 0.5 + nearest_int_fix)) * 2.0;
     }
 
     return rnd / factor;
@@ -86,18 +88,25 @@ double fpyc_bankers_round(double arg, int64_t ndigits)
 
 int64_t ipyc_bankers_round(int64_t arg, int64_t ndigits)
 {
-    if (ndigits >= 0) {
+    if (ndigits >= 0)
+    {
         return arg;
-    } else {
+    }
+    else
+    {
         int64_t mul_fact = 1;
-        for (int i = 0; i< -ndigits; ++i) mul_fact *= 10;
+        for (int i = 0; i < -ndigits; ++i)
+            mul_fact *= 10;
 
-        int64_t pivot_point = copysign(5*mul_fact/10, arg);
+        int64_t pivot_point = copysign(5 * mul_fact / 10, arg);
         int64_t remainder = arg % mul_fact;
-        if ( remainder == pivot_point ) {
+        if (remainder == pivot_point)
+        {
             int64_t val = (mul_fact - remainder) / mul_fact;
             return (val + (val & 1)) * mul_fact;
-        } else {
+        }
+        else
+        {
             return ((arg + pivot_point) / mul_fact) * mul_fact;
         }
     }
@@ -109,3 +118,7 @@ extern inline int64_t      pyc_modulo(int64_t a, int64_t b);
 extern inline double        pyc_fmodulo(double a, double b);
 extern inline double complex complex_min(double complex a, double complex b);
 extern inline double complex complex_max(double complex a, double complex b);
+extern inline int8_t py_floor_div_int8_t(int8_t x, int8_t y);
+extern inline int16_t py_floor_div_int16_t(int16_t x, int16_t y);
+extern inline int32_t py_floor_div_int32_t(int32_t x, int32_t y);
+extern inline int64_t py_floor_div_int64_t(int64_t x, int64_t y);
diff --git a/pyccel/stdlib/math/pyc_math_c.h b/pyccel/stdlib/math/pyc_math_c.h
index 1224402d..62e7db90 100644
--- a/pyccel/stdlib/math/pyc_math_c.h
+++ b/pyccel/stdlib/math/pyc_math_c.h
@@ -52,6 +52,16 @@ 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);
 
+#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)))); \
+    }
+
+PY_FLOOR_DIV_TYPE(int8_t)
+PY_FLOOR_DIV_TYPE(int16_t)
+PY_FLOOR_DIV_TYPE(int32_t)
+PY_FLOOR_DIV_TYPE(int64_t)
+
 inline double complex complex_min(double complex a, double complex b) {
     bool lt = creal(a) == creal(b) ? cimag(a) < cimag(b) : creal(a) < creal(b);
     return lt ? a : b;
diff --git a/pyccel/stdlib/math/pyc_math_f90.f90 b/pyccel/stdlib/math/pyc_math_f90.f90
index 7ea9496e..e75e373d 100644
--- a/pyccel/stdlib/math/pyc_math_f90.f90
+++ b/pyccel/stdlib/math/pyc_math_f90.f90
@@ -11,9 +11,10 @@ module pyc_math_f90
          f64 => C_DOUBLE, &
          c64 => C_DOUBLE_COMPLEX, &
          c32 => C_FLOAT_COMPLEX
-
 implicit none
 
+private
+
 real(f64), parameter, private :: pi = 4.0_f64 * DATAN(1.0_f64)
 
 interface pyc_gcd
@@ -56,6 +57,13 @@ interface pyc_bankers_round
     module procedure pyc_bankers_round_int
 end interface pyc_bankers_round
 
+interface pyc_floor_div
+    module procedure pyc_floor_div_i8
+    module procedure pyc_floor_div_i16
+    module procedure pyc_floor_div_i32
+    module procedure pyc_floor_div_i64
+end interface pyc_floor_div
+
 public :: pyc_gcd, &
           pyc_factorial, &
           pyc_lcm, &
@@ -65,8 +73,8 @@ public :: pyc_gcd, &
           amin, &
           csgn, &
           csign, &
-          pyc_bankers_round
-
+          pyc_bankers_round, &
+          pyc_floor_div
 contains
 
 ! Implementation of math factorial function
@@ -421,4 +429,29 @@ pure function pyc_bankers_round_int(arg, ndigits) result(rnd)
 
 end function pyc_bankers_round_int
 
+elemental pure integer(kind=1) function pyc_floor_div_i8(x, y) result(res)
+  implicit none
+  integer(kind=1), intent(in) :: x, y
+  res = x / y - merge(1, 0, mod(x, y) /= 0 .and. ((x < 0) .neqv. (y < 0)))
+end function pyc_floor_div_i8
+
+elemental pure integer(kind=2) function pyc_floor_div_i16(x, y) result(res)
+  implicit none
+  integer(kind=2), intent(in) :: x, y
+  res = x / y - merge(1, 0, mod(x, y) /= 0 .and. ((x < 0) .neqv. (y < 0)))
+end function pyc_floor_div_i16
+
+elemental pure integer(kind=4) function pyc_floor_div_i32(x, y) result(res)
+  implicit none
+  integer(kind=4), intent(in) :: x, y
+  res = x / y - merge(1, 0, mod(x, y) /= 0 .and. ((x < 0) .neqv. (y < 0)))
+end function pyc_floor_div_i32
+
+elemental pure integer(kind=8) function pyc_floor_div_i64(x, y) result(res)
+  implicit none
+  integer(kind=8), intent(in) :: x, y
+  res = x / y - merge(1, 0, mod(x, y) /= 0 .and. ((x < 0) .neqv. (y < 0)))
+end function pyc_floor_div_i64
+
+
 end module pyc_math_f90


Test patch:
diff --git a/tests/epyccel/test_epyccel_division.py b/tests/epyccel/test_epyccel_division.py
index 11aec7a2..10cda6ee 100644
--- a/tests/epyccel/test_epyccel_division.py
+++ b/tests/epyccel/test_epyccel_division.py
@@ -184,6 +184,19 @@ def test_call_fdiv_i_i_i(language):
     assert (f(-x, -y, z) == fdiv_i_i_i(-x, -y, z))
     assert isinstance(f(x, y, z), type(fdiv_i_i_i(x, y, z)))
 
+def test_call_fdiv_b_b(language):
+    def fdiv_b_b(x : 'bool', y : 'bool'):
+        return x // y
+
+    fflags = "-Werror -Wconversion"
+
+    f = epyccel(fdiv_b_b, language=language, fflags=fflags)
+    assert (f(True, True) == fdiv_b_b(True, True))
+    assert (f(False, True) == fdiv_b_b(False, True))
+    assert isinstance(f(True, True), type(fdiv_b_b(True, True)))
+
+
+
 def test_call_fdiv_i_r(language):
     def fdiv_i_r(x : int, y : 'float'):
         return x // y