Suggested corrections to mod function

PRECISION.md

@@ -439,36 +439,36 @@ The `intervalMissing.cpp` file implement placeholders for primitives that have n
 
 Mod is the modulo operator. 
 
-It is implemented using C++ function `std::fmod`, which has the property that $fmod(x, y) = x - y \cdot \lfloor \frac{x}{y} \rfloor$. 
-$\lfloor \frac{x}{y}\rfloor$ will be called the quotient of $x$ by $y$.
+It is implemented using C++ function `std::fmod`, which has the property that $\mathrm{fmod}(x, y) = x - y \cdot \lfloor \frac{x}{y} \rfloor$. 
+The value $\lfloor \frac{x}{y}\rfloor$ will be called the quotient of $x$ by $y$.
 
 Let's denote the input intervals $X = [\underline{x}; \overline{x}]$ for the dividend and $Y =[\underline{y}; \overline{y}]$ for the divisor.
-The computation of the result interval is split between two functions, `positiveFMod`, dealing with the case where intervals for both arguments are positive, and `Mod`, which reduces the general case to the first function by splitting the arguments intervals between positive and negative values. 
+The case where both intervals are positive is a particular case, implemented with function `positiveFMod`, and the general case can be reduced to the positive case by splitting the arguments intervals between positive and negative values, which is done in function `Mod`.
 
 In the positive case, the result interval is automatically included in $[0; \overline{y}[$.
 This is a direct consequence of the definition of $\lfloor~\rfloor$:
 $\lfloor \frac{x}{y} \rfloor \le \frac{x}{y} < \lfloor \frac{x}{y} \rfloor + 1$,
-so $y\cdot \lfloor \frac{x}{y} \rfloor \le x < y \lfloor \frac{x}{y} \rfloor + y$
+so $y\cdot \lfloor \frac{x}{y} \rfloor \le x < y \cdot \lfloor \frac{x}{y} \rfloor + y$
 and $0 \le x - y\cdot \lfloor \frac{x}{y} \rfloor  < y \le \overline{y}$.
 
-The $fmod$ function is a piecewise linear function, with discontinuities along the lines of equation $x = k\cdot y$, $k \in \mathbb{Z}$. 
+The $\mathrm{fmod}$ function is a piecewise linear function, with discontinuities along the lines of equation $x = k\cdot y$, $k \in \mathbb{Z}$. 
 These discontinuities correspond to the places where $x$ is a multiple of $y$ (and the remainder is thus null).
-- On the $x \ge ky$ side of the line, the limit of $fmod(x, y)$ is $0$ as we approach the line and the quotient $\lfloor\frac{x}{y}\rfloor = k$.
-- On the $x < ky$ side, the limit of $fmod(x, y)$ is $y$ and $\lfloor \frac{x}{y} \rfloor = k-1$.
+- On the $x \ge ky$ side of the line, the limit of $\mathrm{fmod}(x, y)$ is $0$ as we approach the line and the quotient $\lfloor\frac{x}{y}\rfloor = k$.
+- On the $x < ky$ side, the limit of $\mathrm{fmod}(x, y)$ is $y$ and $\lfloor \frac{x}{y} \rfloor = k-1$.
 
 We have to distinguish two cases, depending on whether the interval product $X \times Y$ spans a discontinuity or not.
 Let us pose $n = \lfloor \frac{\underline{x}}{\overline{y}}\rfloor$, which is the smallest value that can be taken by the integral quotient between $x$ and $y$.
 
 
 If it doesn't span a discontinuity, the quotient $\lfloor \frac{x}{y}\rfloor$ is constant with value $n$.
-The $fmod$ function is thus linear over the domain: $fmod(x, y) = x - n\cdot y$, and has maximum value $\overline{x} - n \cdot\underline{y}$ and minimum $\underline{x} - n \cdot \overline{y}$.
+The $\mathrm{fmod}$ function is thus linear over the domain: $\mathrm{fmod}(x, y) = x - n\cdot y$, and has maximum value $\overline{x} - n \cdot\underline{y}$ and minimum $\underline{x} - n \cdot \overline{y}$.
 This is the case represented by the last `return` statement in the implementation.
 
 If the domain does span discontinuities, the minimum value is $0$ and is attained on the discontinuities.
-Regarding the maximum, the function $fmod$ presents local suprema along these discontinuities of equation $x = k\cdot y, ~k \ge n+1$.
+Regarding the maximum, the function $\mathrm{fmod}$ presents local suprema along these discontinuities of equation $x = k\cdot y, ~k \ge n+1$.
 We exclude the discontinuity $x = n \cdot y$ because $n = \lfloor \frac{\underline{x}}{\overline{y}} \rfloor$ (by definition of $n$) means that $\underline{x} \ge n \overline{y}$ (by definition of $\lfloor ~\rfloor$), from which we can deduce that $\forall x \in X, y \in Y, x > n \cdot y$, except maybe in the case $x = \underline{x}, y = \overline{y}$, where the discontinuity line intersects the rectangle domain in the single point $(\underline{x}, \overline{y})$.
 
-We are looking for the largest value of $fmod(x, y)$. 
+We are looking for the largest value of $\mathrm{fmod}(x, y)$. 
 We have seen that the limit of this value is $y$ along lines of equation $x = k\cdot y, ~k\ge n+1$.
 We are thus looking for the largest $y$ that verify $\exists x \in X, k \ge n+1, y = \frac{1}{k} x$.
 
@@ -483,6 +483,11 @@ Thus, the upper bound of the output interval is:
 
 ### Precision
 
+Let us assume that $x$ is represented with precision $l_x$, from which the ulp $u_x$, and $y$ has precision $l_y$ and ulp $u_y$.
+If we assume that $u_x < u_y$, we can write $x = p \cdot u_x$ and $y = q \cdot u_x$, where $p, q \in \mathbb{Z}$.
+Then, $\mathrm{fmod}(x, y) = x - \lfloor\frac{x}{y}\rfloor \cdot y = (p - \lfloor \frac{x}{y} \rfloor \cdot q)\cdot u_x$, thus $\mathrm{fmod}(x, y)$ can be represented with precision $l_x$.
+A similar reasoning takes place when $u_y < u_x$.
+
 Output precision is the finest precision between those of the two arguments: this will be the one to determine which discrete grid the result aligns with.
 
 ### Remark on the sampling tests

interval/intervalMod.cpp

@@ -76,7 +76,7 @@ std::pair<interval, interval> splitnz(interval x)
     if (x.hi() < 0) {
         return {x, empty()};
     }
-    return {interval{x.lo(), nexttoward(0.0, -1.0)}, interval{nexttoward(0.0, 1.0), x.hi()}};
+    return {interval{x.lo(), nexttoward(0.0, -1.0), x.lsb()}, interval{nexttoward(0.0, 1.0), x.hi(), x.lsb()}};
 }
 
 //------------------------------------------------------------------------------------------
@@ -97,9 +97,11 @@ interval positiveFMod(const interval& x, const interval& y)
         return empty();
     }
     int n = int(x.lo() / y.hi());
+    std::cout << "n = " << n << std::endl;
     int precision = std::min(x.lsb(), y.lsb());
 
-    if (n == 0) {
+    // n == 0 obeys the same rules as the general case
+    /* if (n == 0) {
         // prop: x.lo() < y.hi()
         if (y.hi() > x.hi()) {
             if (y.lo() > x.hi()) {
@@ -111,7 +113,7 @@ interval positiveFMod(const interval& x, const interval& y)
         }
         // prop: x.lo() < y.hi() <= x.hi()
         return interval{0.0, nexttoward(y.hi(), 0), precision};
-    }
+    }*/
 
     // prop: n > 0 && y.hi() <= x.lo()
     double hi = x.hi() / (n + 1);
@@ -123,16 +125,8 @@ interval positiveFMod(const interval& x, const interval& y)
         return interval{0.0, nexttoward(hi, 0), precision};
     }
     // prop : y.lo() > hi
-    // we compute the intersections of y with the segment
-    // (hi,hi)->(zhi,0) and (lo,lo)->(zlo,0)
-    double lo  = x.lo() / (n + 1);
-    double zhi = x.hi() / n;
-    double zlo = x.lo() / n;
-
-    double v1 = hi - (y.lo() - hi) / (zhi - hi) * hi;
-    double v0 = lo - (y.hi() - lo) / (zlo - lo) * lo;
-
-    return interval{v0, v1, precision};
+    // in that case, the quotient between x and y is constant and equal
+    return interval{x.lo() - n*y.hi(), x.hi() - n*y.lo(), precision};
 }
 
 // fmod of two signed intervals
@@ -148,27 +142,34 @@ interval interval_algebra::Mod(const interval& x, const interval& y)
     auto xpyp = positiveFMod(xp, yp);
 
     // Make sure these 4 values are in the resulting interval
-    auto bb = singleton(fmod(x.hi(), y.hi())) + singleton(fmod(x.lo(), y.hi())) + singleton(fmod(x.hi(), y.lo())) +
-              singleton(fmod(x.lo(), y.lo()));
+    auto bb = singleton(fmod(x.hi(), y.hi())) 
+            + singleton(fmod(x.lo(), y.hi())) 
+            + singleton(fmod(x.hi(), y.lo())) 
+            + singleton(fmod(x.lo(), y.lo()));
+
+    bb = interval{bb.lo(), bb.hi(), std::min(x.lsb(), y.lsb())};
 
     // combine all the intervals
     return bb + xnyn + xnyp + xpyn + xpyp;
 }
 
 void interval_algebra::testMod()
 {
-    // check("test algebra Mod", Mod(interval(-100, 100), 1.0), interval(nextafter(-1.0, 0.0), nextafter(1.0, 0.0)));
-    // check("test algebra Mod", Mod(interval(0, 100), 2), interval(0, nextafter(2.0, 0)));
-    // check("test algebra Mod", Mod(interval(0, 100), -1.0), interval(0, nextafter(1.0, 0)));
-    /* check("test algebra Mod", Mod(interval(5, 7), interval(4, 4.5)), interval(0.5, 3));
+    check("test algebra Mod", Mod(interval(-100, 100), 1.0), interval(nextafter(-1.0, 0.0), nextafter(1.0, 0.0)));
+    check("test algebra Mod", Mod(interval(0, 100), 2), interval(0, nextafter(2.0, 0)));
+    check("test algebra Mod", Mod(interval(0, 100), -1.0), interval(0, nextafter(1.0, 0)));
+    check("test algebra Mod", Mod(interval(5, 7), interval(4, 4.5)), interval(0.5, 3));
     check("test algebra Mod", Mod(interval(5, 7), interval(8, 10)), interval(5, 7));
     check("test algebra Mod", Mod(interval(-7, 7), interval(8, 10)), interval(-7, 7));
-    check("test algebra Mod", Mod(interval(0, 100), interval(7, 7)), interval(0, nextafter(7.0, 0.0))); */
+    check("test algebra Mod", Mod(interval(0, 100), interval(7, 7)), interval(0, nextafter(7.0, 0.0)));
 
     analyzeBinaryMethod(10, 10000, "mod", interval(0, 10, -5), interval(0, 10, -5), fmod, &interval_algebra::Mod);
     analyzeBinaryMethod(10, 10000, "mod", interval(0, 10, 1), interval(0, 10, 0), fmod, &interval_algebra::Mod);
     analyzeBinaryMethod(10, 10000, "mod", interval(0, 10, 0), interval(0, 10, 0), fmod, &interval_algebra::Mod);
     analyzeBinaryMethod(10, 10000, "mod", interval(0, 10, 0), interval(0, 10, -5), fmod, &interval_algebra::Mod);
+
+    analyzeBinaryMethod(10, 100000000, "mod", interval(3, 4, -3), interval(1.2, 1.4, -3), fmod, &interval_algebra::Mod);
+    // analyzeBinaryMethod(10, 10000, "mod", interval(-10, 10, -5), interval(-10, 10, -5), fmod, &interval_algebra::Mod);
 }
 
 }  // namespace itv