Add note on filtering

Flake lock file updates:

• Updated input 'nixpkgs':
    'github:nixos/nixpkgs/cd749f58ba83f7155b7062dd49d08e5e47e44d50' (2023-04-17)
  → 'github:nixos/nixpkgs/897876e4c484f1e8f92009fd11b7d988a121a4e7' (2023-05-06)

COPYING

@@ -1,14 +0,0 @@
-Copyright (C) 2023 Simon Bruder
-
-This program is free software: you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation, either version 3 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License
-along with this program.  If not, see <https://www.gnu.org/licenses/>.

README.md

@@ -8,6 +8,19 @@ Each exercise includes a separate `readme.txt`
 that describes things specific to how I implemented the task.
 However, those documents are in German, which is a requirement.
 
+**Note**:
+Due to copyright restrictions,
+I am not allowed to publish many of my solutions.
+The repository `ecg-prog-filtered` only includes files solely written by me.
+It is filtered with the [`filter.sh`](./filter.sh) script.
+Please **don’t** rely on the history of this repository,
+as rewriting is part of how it can exist.
+See [Licence](#licence) for more information on the licencing.
+
+If you want access to all of my solutions,
+please write me an email (or contact me in another way),
+I’ll see what I can do.
+
 ## Usage
 
 ### Build
@@ -39,5 +52,16 @@ ln -sf build/compile_commands.json ..
 
 ## Licence
 
-[GPL-3.0-or-later](https://spdx.org/licenses/GPL-3.0-or-later.html).
-See [COPYING](./COPYING) and [gpl-3.0.txt](./gpl-3.0.txt) for details.
+Depending on the file,
+[GPL-3.0-or-later](https://spdx.org/licenses/GPL-3.0-or-later.html)
+or [LGPL-3.0-or-later](https://spdx.org/licenses/LGPL-3.0-or-later.html)
+or “Copyright (C) CGV TU Dresden - All Rights Reserved”.
+
+The published version only includes files under free software licences.
+Sadly, that makes them very incomplete,
+because most implementations done by me are in files from the template,
+which makes my additions a derivative work.
+
+All files not marked differently are released under GPL-3.0-or-later.
+
+See [gpl-3.0.txt](./gpl-3.0.txt) and [lgpl-3.0.txt](./lgpl-3.0.txt) for details.

filter.sh

@@ -0,0 +1,9 @@
+#!/usr/bin/env bash
+# SPDX-License-Identifier: GPL-3.0-or-later
+set -euo pipefail
+export FILTER_BRANCH_SQUELCH_WARNING=1
+git branch -D filtered
+git switch -c filtered
+git filter-branch -f --tree-filter "rm -f $(rg --files-without-match 'SPDX-License-Identifier: (L)?GPL-3.0-or-later' -g '/u??/**/*' | tr '\n' ' ')" HEAD
+git push -f -u filtered filtered
+git switch -

flake.lock

@@ -20,16 +20,16 @@
     },
     "nixpkgs": {
       "locked": {
-        "lastModified": 1681759395,
-        "narHash": "sha256-7aaRtLxLAy8qFVIA26ulB+Q5nDVzuQ71qi0s0wMjAws=",
+        "lastModified": 1683408522,
+        "narHash": "sha256-9kcPh6Uxo17a3kK3XCHhcWiV1Yu1kYj22RHiymUhMkU=",
         "owner": "nixos",
         "repo": "nixpkgs",
-        "rev": "cd749f58ba83f7155b7062dd49d08e5e47e44d50",
+        "rev": "897876e4c484f1e8f92009fd11b7d988a121a4e7",
         "type": "github"
       },
       "original": {
         "owner": "nixos",
-        "ref": "nixos-22.11",
+        "ref": "nixos-unstable",
         "repo": "nixpkgs",
         "type": "github"
       }

flake.nix

@@ -1,6 +1,6 @@
 {
   inputs = {
-    nixpkgs.url = "github:nixos/nixpkgs/nixos-22.11";
+    nixpkgs.url = "github:nixos/nixpkgs/nixos-unstable";
     flake-utils.url = "github:numtide/flake-utils";
   };
 
@@ -21,10 +21,22 @@
             doCheck = true;
           })
           { };
+
+        u02 = pkgs.callPackage
+          ({ stdenv, catch2_3, cmake, freeglut, libGL, libGLU }: stdenv.mkDerivation {
+            name = "ecg-u02";
+
+            src = ./u02;
+
+            nativeBuildInputs = [ catch2_3 cmake freeglut libGL libGLU ];
+
+            doCheck = true;
+          })
+          { };
       };
 
       devShells.default = pkgs.mkShell {
-        nativeBuildInputs = with pkgs; [ catch2_3 cmake ];
+        nativeBuildInputs = with pkgs; [ catch2_3 cmake freeglut libGL libGLU ];
       };
     });
 }

lgpl-3.0.txt

@@ -0,0 +1,165 @@
+                   GNU LESSER GENERAL PUBLIC LICENSE
+                       Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc. <https://fsf.org/>
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+
+  This version of the GNU Lesser General Public License incorporates
+the terms and conditions of version 3 of the GNU General Public
+License, supplemented by the additional permissions listed below.
+
+  0. Additional Definitions.
+
+  As used herein, "this License" refers to version 3 of the GNU Lesser
+General Public License, and the "GNU GPL" refers to version 3 of the GNU
+General Public License.
+
+  "The Library" refers to a covered work governed by this License,
+other than an Application or a Combined Work as defined below.
+
+  An "Application" is any work that makes use of an interface provided
+by the Library, but which is not otherwise based on the Library.
+Defining a subclass of a class defined by the Library is deemed a mode
+of using an interface provided by the Library.
+
+  A "Combined Work" is a work produced by combining or linking an
+Application with the Library.  The particular version of the Library
+with which the Combined Work was made is also called the "Linked
+Version".
+
+  The "Minimal Corresponding Source" for a Combined Work means the
+Corresponding Source for the Combined Work, excluding any source code
+for portions of the Combined Work that, considered in isolation, are
+based on the Application, and not on the Linked Version.
+
+  The "Corresponding Application Code" for a Combined Work means the
+object code and/or source code for the Application, including any data
+and utility programs needed for reproducing the Combined Work from the
+Application, but excluding the System Libraries of the Combined Work.
+
+  1. Exception to Section 3 of the GNU GPL.
+
+  You may convey a covered work under sections 3 and 4 of this License
+without being bound by section 3 of the GNU GPL.
+
+  2. Conveying Modified Versions.
+
+  If you modify a copy of the Library, and, in your modifications, a
+facility refers to a function or data to be supplied by an Application
+that uses the facility (other than as an argument passed when the
+facility is invoked), then you may convey a copy of the modified
+version:
+
+   a) under this License, provided that you make a good faith effort to
+   ensure that, in the event an Application does not supply the
+   function or data, the facility still operates, and performs
+   whatever part of its purpose remains meaningful, or
+
+   b) under the GNU GPL, with none of the additional permissions of
+   this License applicable to that copy.
+
+  3. Object Code Incorporating Material from Library Header Files.
+
+  The object code form of an Application may incorporate material from
+a header file that is part of the Library.  You may convey such object
+code under terms of your choice, provided that, if the incorporated
+material is not limited to numerical parameters, data structure
+layouts and accessors, or small macros, inline functions and templates
+(ten or fewer lines in length), you do both of the following:
+
+   a) Give prominent notice with each copy of the object code that the
+   Library is used in it and that the Library and its use are
+   covered by this License.
+
+   b) Accompany the object code with a copy of the GNU GPL and this license
+   document.
+
+  4. Combined Works.
+
+  You may convey a Combined Work under terms of your choice that,
+taken together, effectively do not restrict modification of the
+portions of the Library contained in the Combined Work and reverse
+engineering for debugging such modifications, if you also do each of
+the following:
+
+   a) Give prominent notice with each copy of the Combined Work that
+   the Library is used in it and that the Library and its use are
+   covered by this License.
+
+   b) Accompany the Combined Work with a copy of the GNU GPL and this license
+   document.
+
+   c) For a Combined Work that displays copyright notices during
+   execution, include the copyright notice for the Library among
+   these notices, as well as a reference directing the user to the
+   copies of the GNU GPL and this license document.
+
+   d) Do one of the following:
+
+       0) Convey the Minimal Corresponding Source under the terms of this
+       License, and the Corresponding Application Code in a form
+       suitable for, and under terms that permit, the user to
+       recombine or relink the Application with a modified version of
+       the Linked Version to produce a modified Combined Work, in the
+       manner specified by section 6 of the GNU GPL for conveying
+       Corresponding Source.
+
+       1) Use a suitable shared library mechanism for linking with the
+       Library.  A suitable mechanism is one that (a) uses at run time
+       a copy of the Library already present on the user's computer
+       system, and (b) will operate properly with a modified version
+       of the Library that is interface-compatible with the Linked
+       Version.
+
+   e) Provide Installation Information, but only if you would otherwise
+   be required to provide such information under section 6 of the
+   GNU GPL, and only to the extent that such information is
+   necessary to install and execute a modified version of the
+   Combined Work produced by recombining or relinking the
+   Application with a modified version of the Linked Version. (If
+   you use option 4d0, the Installation Information must accompany
+   the Minimal Corresponding Source and Corresponding Application
+   Code. If you use option 4d1, you must provide the Installation
+   Information in the manner specified by section 6 of the GNU GPL
+   for conveying Corresponding Source.)
+
+  5. Combined Libraries.
+
+  You may place library facilities that are a work based on the
+Library side by side in a single library together with other library
+facilities that are not Applications and are not covered by this
+License, and convey such a combined library under terms of your
+choice, if you do both of the following:
+
+   a) Accompany the combined library with a copy of the same work based
+   on the Library, uncombined with any other library facilities,
+   conveyed under the terms of this License.
+
+   b) Give prominent notice with the combined library that part of it
+   is a work based on the Library, and explaining where to find the
+   accompanying uncombined form of the same work.
+
+  6. Revised Versions of the GNU Lesser General Public License.
+
+  The Free Software Foundation may publish revised and/or new versions
+of the GNU Lesser General Public License from time to time. Such new
+versions will be similar in spirit to the present version, but may
+differ in detail to address new problems or concerns.
+
+  Each version is given a distinguishing version number. If the
+Library as you received it specifies that a certain numbered version
+of the GNU Lesser General Public License "or any later version"
+applies to it, you have the option of following the terms and
+conditions either of that published version or of any later version
+published by the Free Software Foundation. If the Library as you
+received it does not specify a version number of the GNU Lesser
+General Public License, you may choose any version of the GNU Lesser
+General Public License ever published by the Free Software Foundation.
+
+  If the Library as you received it specifies that a proxy can decide
+whether future versions of the GNU Lesser General Public License shall
+apply, that proxy's public statement of acceptance of any version is
+permanent authorization for you to choose that version for the
+Library.

u02/include/sweep_line_tool.h

@@ -0,0 +1,24 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+#pragma once
+
+#include "tool_base.h"
+
+class sweep_line_tool : public tool_base {
+public:
+  sweep_line_tool(canvas_buffer &canvas);
+
+  // Draw example triangle
+  void draw();
+  // Compatibility for main application (only handles draw methods with one or two points)
+  void draw(int _x, int _y);
+  // Draw triangle provided by three given points
+  void draw(int x0, int y0, int x1, int y1, int x2, int y2);
+
+  void set_text(std::stringstream &stream);
+
+private:
+  // Draw every pixel on the specified y coordinate,
+  // in the interval given by the boundaries b1 and b2.
+  // The boundaries do not need to be sorted.
+  void draw_interval(int b1, int b2, int y);
+};

u02/include/util.h

@@ -0,0 +1,51 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+#pragma once
+#include <cstdint>
+#include <tuple>
+
+typedef uint8_t Transformation;
+
+const Transformation TRANSFORM_MIRROR_X = 1 << 0;
+const Transformation TRANSFORM_MIRROR_Y = 1 << 1;
+const Transformation TRANSFORM_ROTATE_CW = 1 << 2;
+const Transformation TRANSFORM_ROTATE_CCW = 1 << 3;
+
+// Applies the provided transformation to the point, mutating it in place.
+// Rotation is done before mirroring.
+void transform_mut(Transformation transformation, int &x, int &y);
+// Applies the provided transformation to the point,
+// returning the transformed point.
+// Rotation is done before mirroring.
+std::pair<int, int> transform(Transformation transformation, int x, int y);
+
+// Applies the inverse transformation to the point, mutating it in place.
+// Composition of this with the transformation is the identity function.
+void transform_inv_mut(Transformation transformation, int &x, int &y);
+// Applies the inverse transformation to the point,
+// returning the transformed point.
+// Composition of this with the transformation is the identity function.
+std::pair<int, int> transform_inv(Transformation transformation, int x, int y);
+
+// Returns the transformation required
+// to make the given endpoints of a line conform
+// to the standard case for rasterization.
+Transformation transformation_to_standard_case(int x0, int y0, int x1, int y1);
+
+// Returns the barycentric coordinates of the point given by (xp, yp)
+// in the triangle given by the three points (x0, y0), (x1, y1), (x2, y2).
+// May return non-real values when the points are on a straight line.
+std::tuple<float, float, float> barycentric_coordinates(int x0, int y0, int x1,
+                                                        int y1, int x2, int y2,
+                                                        int xp, int yp);
+
+// Checks if the point given by (xp, yp) is inside the triangle
+// given by the three points (x0, y0), (x1, y1), (x2, y2).
+bool point_in_triangle(int x0, int y0, int x1, int y1, int x2, int y2, int xp,
+                       int yp);
+
+// Sorts the points of a triangle to be in ascending order (y0 < y1 < y2).
+void sort_triangle_points(int &x0, int &y0, int &x1, int &y1, int &x2, int &y2);
+
+// This calculates the slope of a line from (x0, y0) to (x1, y1).
+// It handles special cases to ensure the return value will be real.
+float slope(int x0, int y0, int x1, int y1);

u02/src/sweep_line_tool.cpp

@@ -0,0 +1,100 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+#include "sweep_line_tool.h"
+#include "dda_line_tool.h"
+#include "util.h"
+#include <cmath>
+#include <iostream>
+#include <vector>
+
+// Calculate the inverse of the DDA function.
+int dda_inv(int x0, int y0, float m, int y) {
+  // This uses the regular function that is the basis of DDA
+  //
+  //     y_i = y_0 + m·(x_i - x_0)
+  //
+  // but rearranges it to be the inverse function:
+  //
+  //   y_i = y_0 + m·(x_i - x_0)
+  // ⇔ y_i - y_0 + x_0·m = x_i·m
+  // ⇔ x_i = (y_i - y_0)/m + x_0
+  //
+  // This returns a valid x coordinate on the line
+  // starting from (x0, y0) with the slope m.
+
+  // Handle special case of flat line
+  if (m == 0)
+    return x0;
+  else
+    return round((y - y0) / m + x0);
+}
+
+sweep_line_tool::sweep_line_tool(canvas_buffer &canvas) : tool_base(canvas) {
+  shape = TS_NONE;
+  is_draggable = false;
+}
+
+void sweep_line_tool::draw_interval(int b1, int b2, int y) {
+  for (int x = std::min(b1, b2); x <= std::max(b1, b2); x++) {
+    canvas.set_pixel(x, y);
+  }
+}
+
+void sweep_line_tool::draw() { draw(10, 10, 90, 30, 30, 90); }
+
+void sweep_line_tool::draw(int _x, int _y) { draw(); }
+
+void sweep_line_tool::draw(int x0, int y0, int x1, int y1, int x2, int y2) {
+  // Terminology:
+  //
+  //       (x0, y0)
+  //          +
+  //          | \
+  //          |   \ m_1
+  //          |first\
+  //          | pass  \
+  // m_shared |---------+ (x1, y1)
+  //          |second /
+  //          |pass /
+  //          |   / m_2
+  //          | /
+  //          +
+  //       (x2, y2)
+
+  // Sort triangle points (in place) so that y0 < y1 < y2
+  sort_triangle_points(x0, y0, x1, y1, x2, y2);
+
+  // Slope of the side limiting the first pass (only)
+  float m_1 = slope(x0, y0, x1, y1);
+  // Slope of the side limiting the second pass (only)
+  float m_2 = slope(x1, y1, x2, y2);
+  // Slope of the side limiting both passes
+  float m_shared = slope(x0, y0, x2, y2);
+
+  // First pass
+  if (y0 == y1) {
+    // If the first two points are on the same height, only draw one line.
+    // This is only needed for the first interval,
+    // because in the case that y1 == y2,
+    // the problematic line would have already been handled in the first pass.
+    draw_interval(x0, x1, y0);
+  } else {
+    for (int y = y0; y <= y1; y++) {
+      int b1 = dda_inv(x0, y0, m_1, y);
+      int b2 = dda_inv(x0, y0, m_shared, y);
+      draw_interval(b1, b2, y);
+    }
+  }
+
+  // Second pass
+  // it can start iterating at y1 + 1,
+  // because y1 is already included in the first pass.
+  for (int y = y1 + 1; y <= y2; y++) {
+    int b1 = dda_inv(x1, y1, m_2, y);
+    int b2 = dda_inv(x0, y0, m_shared, y);
+    draw_interval(b1, b2, y);
+  }
+}
+
+void sweep_line_tool::set_text(std::stringstream &stream) {
+  stream << "Tool: Sweep-Line";
+}

u02/src/tests.cpp

@@ -0,0 +1,524 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+#include <algorithm>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/generators/catch_generators.hpp>
+#include <catch2/generators/catch_generators_adapters.hpp>
+#include <catch2/generators/catch_generators_random.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+
+#include "bresenham_circle_tool.h"
+#include "bresenham_line_tool.h"
+#include "canvas_buffer.h"
+#include "dda_line_tool.h"
+#include "non_recursive_fill_tool.h"
+#include "rectangle_tool.h"
+#include "recursive_fill_tool.h"
+#include "sweep_line_tool.h"
+#include "util.h"
+
+using Catch::Matchers::WithinRel;
+
+TEST_CASE("Transform Mirror") {
+  // elementary operations
+  REQUIRE(transform(TRANSFORM_MIRROR_X, 10, 20) == std::make_pair(10, -20));
+  REQUIRE(transform(TRANSFORM_MIRROR_Y, 10, 20) == std::make_pair(-10, 20));
+
+  // composite operations
+  REQUIRE(transform(TRANSFORM_MIRROR_X | TRANSFORM_MIRROR_Y, 10, 20) ==
+          std::make_pair(-10, -20));
+}
+
+TEST_CASE("Transform Rotate") {
+  REQUIRE(transform(TRANSFORM_ROTATE_CW, 10, 20) ==
+          std::make_pair(-20, 10)); // 4th quadrant
+  REQUIRE(transform(TRANSFORM_ROTATE_CW, -20, 10) ==
+          std::make_pair(-10, -20)); // 3rd quadrant
+  REQUIRE(transform(TRANSFORM_ROTATE_CW, -10, -20) ==
+          std::make_pair(20, -10)); // 2nd quadrant
+  REQUIRE(transform(TRANSFORM_ROTATE_CW, 20, -10) ==
+          std::make_pair(10, 20)); // 1st quadrant
+
+  REQUIRE(transform(TRANSFORM_ROTATE_CCW, 20, -10) ==
+          std::make_pair(-10, -20)); // 1nd quadrant
+  REQUIRE(transform(TRANSFORM_ROTATE_CCW, -10, -20) ==
+          std::make_pair(-20, 10)); // 2rd quadrant
+  REQUIRE(transform(TRANSFORM_ROTATE_CCW, -20, 10) ==
+          std::make_pair(10, 20)); // 3th quadrant
+  REQUIRE(transform(TRANSFORM_ROTATE_CCW, 10, 20) ==
+          std::make_pair(20, -10)); // 4st quadrant
+}
+
+TEST_CASE("Transform Rotate + Mirror") {
+  REQUIRE(transform(TRANSFORM_ROTATE_CW | TRANSFORM_MIRROR_X, 10, 20) ==
+          std::make_pair(-20, -10));
+  REQUIRE(transform(TRANSFORM_ROTATE_CW | TRANSFORM_MIRROR_Y, 10, 20) ==
+          std::make_pair(20, 10));
+  REQUIRE(
+      transform(TRANSFORM_ROTATE_CW | TRANSFORM_MIRROR_X | TRANSFORM_MIRROR_Y,
+                10, 20) == std::make_pair(20, -10));
+
+  REQUIRE(transform(TRANSFORM_ROTATE_CCW | TRANSFORM_MIRROR_X, 10, 20) ==
+          std::make_pair(20, 10));
+  REQUIRE(transform(TRANSFORM_ROTATE_CCW | TRANSFORM_MIRROR_Y, 10, 20) ==
+          std::make_pair(-20, -10));
+  REQUIRE(
+      transform(TRANSFORM_ROTATE_CCW | TRANSFORM_MIRROR_X | TRANSFORM_MIRROR_Y,
+                10, 20) == std::make_pair(-20, 10));
+}
+
+TEST_CASE("Transform = Inverse Transform ○ Transform") {
+  const int x = GENERATE(take(10, random(-100, 100)));
+  const int y = GENERATE(take(10, random(-100, 100)));
+  // this iterates over all possible transformations,
+  // even bogus ones (like rotating cw and ccw)
+  for (Transformation transformation = 0; transformation < 0b10000;
+       transformation++) {
+    int xt, yt;
+    std::tie(xt, yt) = transform(transformation, x, y);
+    int xti, yti;
+    std::tie(xti, yti) = transform_inv(transformation, xt, yt);
+    REQUIRE(x == xti);
+    REQUIRE(y == yti);
+  }
+}
+
+TEST_CASE("Transformation to standard case") {
+  REQUIRE(transformation_to_standard_case(5, 20, 20, 10) == 0);
+  REQUIRE(transformation_to_standard_case(5, 5, 20, 15) == TRANSFORM_MIRROR_X);
+  REQUIRE(transformation_to_standard_case(20, 15, 5, 5) == TRANSFORM_MIRROR_Y);
+  REQUIRE(transformation_to_standard_case(20, 10, 5, 20) ==
+          (TRANSFORM_MIRROR_X | TRANSFORM_MIRROR_Y));
+
+  REQUIRE(transformation_to_standard_case(5, 20, 15, 5) ==
+          (TRANSFORM_ROTATE_CW | TRANSFORM_MIRROR_X));
+  REQUIRE(transformation_to_standard_case(5, 5, 15, 20) ==
+          TRANSFORM_ROTATE_CCW);
+  REQUIRE(transformation_to_standard_case(15, 5, 5, 20) ==
+          (TRANSFORM_ROTATE_CCW | TRANSFORM_MIRROR_X));
+  REQUIRE(transformation_to_standard_case(15, 20, 5, 5) == TRANSFORM_ROTATE_CW);
+}
+
+TEST_CASE("Transformation to standard case (prop)") {
+  int x0 = GENERATE(take(10, random(-100, 100)));
+  int y0 = GENERATE(take(10, random(-100, 100)));
+  int x1 = GENERATE(take(10, random(-100, 100)));
+  int y1 = GENERATE(take(10, random(-100, 100)));
+
+  const Transformation transformation =
+      transformation_to_standard_case(x0, y0, x1, y1);
+
+  transform_mut(transformation, x0, y0);
+  transform_mut(transformation, x1, y1);
+
+  REQUIRE(x0 <= x1);
+  REQUIRE(y0 >= y1);
+}
+
+TEST_CASE("Bresenham/DDA line tool (prop: for every row/column, only one pixel "
+          "is set)") {
+  const int size = 100;
+  canvas_buffer *canvas = new canvas_buffer(size, size);
+
+  bresenham_line_tool *tool_bresenham = new bresenham_line_tool(*canvas);
+  dda_line_tool *tool_dda = new dda_line_tool(*canvas);
+
+  tool_base *tool;
+  const int tool_idx = GENERATE(0, 1);
+  switch (tool_idx) {
+  case 0:
+    tool = tool_bresenham;
+    break;
+  case 1:
+    tool = tool_dda;
+    break;
+  }
+
+  const int x0 = GENERATE(take(10, random(0, size - 1)));
+  const int y0 = GENERATE(take(10, random(0, size - 1)));
+  const int x1 = GENERATE(take(10, random(0, size - 1)));
+  const int y1 = GENERATE(take(10, random(0, size - 1)));
+
+  tool->draw(x0, y0, x1, y1);
+
+  const int x_min = std::min(x0, x1);
+  const int x_max = std::max(x0, x1);
+  const int y_min = std::min(y0, y1);
+  const int y_max = std::max(y0, y1);
+
+  // Depending on what the direction of the line (rounded to the next 90°) is,
+  // either every row or column has only one pixel set.
+  bool vertical = false;
+  int draw_direction_min;
+  int draw_direction_max;
+  int unique_direction_min;
+  int unique_direction_max;
+  if (abs(y1 - y0) > abs(x1 - x0)) {
+    vertical = true;
+    draw_direction_min = y_min;
+    draw_direction_max = y_max;
+    unique_direction_min = x_min;
+    unique_direction_max = x_max;
+  } else {
+    draw_direction_min = x_min;
+    draw_direction_max = x_max;
+    unique_direction_min = y_min;
+    unique_direction_max = y_max;
+  }
+
+  bool all_sums_are_one = true;
+  int sum;
+  for (int dd = draw_direction_min; dd <= draw_direction_max; dd++) {
+    sum = 0;
+    for (int ud = unique_direction_min; ud <= unique_direction_max; ud++) {
+      int x, y;
+      if (vertical) {
+        x = ud;
+        y = dd;
+      } else {
+        x = dd;
+        y = ud;
+      }
+      if (canvas->get_pixel(x, y))
+        sum++;
+    }
+    if (sum != 1)
+      all_sums_are_one = false;
+  }
+
+  REQUIRE(all_sums_are_one);
+}
+
+TEST_CASE("Fill (recursive and non recursive) test shape") {
+  canvas_buffer *canvas = new canvas_buffer(100, 100);
+  bresenham_line_tool *tool_line = new bresenham_line_tool(*canvas);
+  recursive_fill_tool *tool_fill_recursive = new recursive_fill_tool(*canvas);
+  non_recursive_fill_tool *tool_fill_non_recursive =
+      new non_recursive_fill_tool(*canvas);
+
+  tool_base *tool_fill;
+  const int tool_fill_idx = GENERATE(0, 1);
+  switch (tool_fill_idx) {
+  case 0:
+    tool_fill = tool_fill_recursive;
+    break;
+  case 1:
+    tool_fill = tool_fill_non_recursive;
+    break;
+  }
+
+  canvas->draw_test_shape();
+
+  REQUIRE_FALSE(canvas->get_pixel(50, 49));
+  REQUIRE_FALSE(canvas->get_pixel(50, 25));
+  REQUIRE_FALSE(canvas->get_pixel(50, 75));
+  tool_fill->draw(50, 25);
+  REQUIRE(canvas->get_pixel(50, 49));
+  REQUIRE(canvas->get_pixel(50, 25));
+  REQUIRE(canvas->get_pixel(50, 75));
+
+  REQUIRE_FALSE(canvas->get_pixel(75, 40));
+  REQUIRE_FALSE(canvas->get_pixel(75, 60));
+  tool_fill->draw(75, 50);
+  REQUIRE(canvas->get_pixel(75, 40));
+  REQUIRE(canvas->get_pixel(75, 60));
+
+  REQUIRE_FALSE(canvas->get_pixel(0, 0));
+  REQUIRE_FALSE(canvas->get_pixel(99, 99));
+  tool_fill->draw(25, 50);
+  REQUIRE(canvas->get_pixel(0, 0));
+  REQUIRE(canvas->get_pixel(99, 99));
+}
+
+TEST_CASE("Fill recursive == Fill non recursive (prop, 5 random lines)") {
+  const int size = 100;
+
+  canvas_buffer *canvas_recursive = new canvas_buffer(size, size);
+  canvas_buffer *canvas_non_recursive = new canvas_buffer(size, size);
+
+  bresenham_line_tool *tool_line_recursive =
+      new bresenham_line_tool(*canvas_recursive);
+  bresenham_line_tool *tool_line_non_recursive =
+      new bresenham_line_tool(*canvas_non_recursive);
+
+  recursive_fill_tool *tool_fill_recursive =
+      new recursive_fill_tool(*canvas_recursive);
+  non_recursive_fill_tool *tool_fill_non_recursive =
+      new non_recursive_fill_tool(*canvas_non_recursive);
+
+  for (int i = 0; i < 5; i++) {
+    const int x0 = GENERATE(take(1, random(0, size - 1)));
+    const int y0 = GENERATE(take(1, random(0, size - 1)));
+    const int x1 = GENERATE(take(1, random(0, size - 1)));
+    const int y1 = GENERATE(take(1, random(0, size - 1)));
+    tool_line_recursive->draw(x0, y0, x1, y1);
+    tool_line_non_recursive->draw(x0, y0, x1, y1);
+  }
+
+  const int x = GENERATE(take(3, random(0, size - 1)));
+  const int y = GENERATE(take(3, random(0, size - 1)));
+  tool_fill_recursive->draw(x, y);
+  tool_fill_non_recursive->draw(x, y);
+
+  bool equal = true;
+  for (int x = 0; x < size; x++) {
+    for (int y = 0; y < size; y++) {
+      if (canvas_recursive->get_pixel(x, y) !=
+          canvas_non_recursive->get_pixel(x, y))
+        equal = false;
+    }
+  }
+  REQUIRE(equal);
+}
+
+TEST_CASE("Rectangle (prop)") {
+  const int size = 100;
+
+  const int x0 = GENERATE(take(5, random(0, size - 1)));
+  const int y0 = GENERATE(take(5, random(0, size - 1)));
+  const int x1 = GENERATE(take(5, random(0, size - 1)));
+  const int y1 = GENERATE(take(5, random(0, size - 1)));
+
+  const int x_min = std::min(x0, x1);
+  const int x_max = std::max(x0, x1);
+  const int y_min = std::min(y0, y1);
+  const int y_max = std::max(y0, y1);
+
+  canvas_buffer *canvas = new canvas_buffer(size, size);
+  rectangle_tool *tool = new rectangle_tool(*canvas);
+
+  tool->draw(x0, y0, x1, y1);
+
+  bool pass = true;
+  for (int x = 0; x < size; x++) {
+    for (int y = 0; y < size; y++) {
+      if (((x == x0 || x == x1) && (y >= y_min && y <= y_max)) ||
+          ((y == y0 || y == y1) && (x >= x_min && x <= x_max))) {
+        if (!canvas->get_pixel(x, y))
+          pass = false;
+      } else {
+        if (canvas->get_pixel(x, y))
+          pass = false;
+      }
+    }
+  }
+  REQUIRE(pass);
+}
+
+TEST_CASE("Bresenham circle (prop: √(x²+y²)-r<ε)") {
+  // Let s be the size of the canvas (s,s).
+  // Let m be the smallest coordinate (x and y) for random points
+  // and M be the largest coordinate (x and y) for random points.
+  //
+  // The largest radius that can fit on canvas (for arbitrary centres) is m.
+  // The largest radius that is possible to create is √(2)(M-m).
+  //
+  // ⇒ √(2)(M-m)   ≥ m
+  // ⇔ √(2)M-√(2)m ≥ m
+  // ⇔ √(2)M       ≥ (1+√(2))m
+  // ⇔ m           ≤ (√(2)/(1+√(2)))m
+  // ⇔ m           ≤ (2-√(2))M    (1)
+  //
+  // Additionally, to have a centered point field,
+  // s-M=m must hold
+  //
+  //   s-M = m
+  // ⇒ s-M ≤ (2-√(2))M
+  // ⇔ s   ≤ (3-√(2))M
+  // ⇔ M   ≥ s/(3-√(2))    (2)
+  //
+  // With this, it now is possible to express M and m in terms of s:
+  //
+  //   m ≤ (2-√(2))M    (1)
+  // ⇔ m ≤ ((2-√(2))/(3-√(2)))s
+  // ⇔ m ≤ ((4-√(2))/7)s (3)
+  //
+  // When the points are rounded to the nearest integer,
+  // M must be rounded down and m rounded down.
+  //
+  // An interactive version of this can be found here:
+  // https://www.desmos.com/calculator/kn19qhue20
+  const int size = 100;                                         // s
+  const int max_c = std::floor(size / (3 - std::sqrt(2)));      // M (2)
+  const int min_c = std::ceil(((4 - std::sqrt(2)) / 7) * size); // m (3)
+
+  const int x0 = GENERATE_COPY(take(10, random(min_c, max_c)));
+  const int y0 = GENERATE_COPY(take(10, random(min_c, max_c)));
+  const int x1 = GENERATE_COPY(take(10, random(min_c, max_c)));
+  const int y1 = GENERATE_COPY(take(10, random(min_c, max_c)));
+
+  if ((x0 == min_c || x0 == max_c) && (x1 == min_c || x1 == max_c) &&
+      (y0 == min_c || y0 == max_c) && (y1 == min_c || y1 == max_c)) {
+    SKIP("All coordinates have extreme value, skipping (avoid rounding error)");
+  }
+
+  const int r =
+      round(std::sqrt(std::pow((x1 - x0), 2) + std::pow((y1 - y0), 2)));
+
+  canvas_buffer *canvas = new canvas_buffer(size, size);
+  bresenham_circle_tool *tool = new bresenham_circle_tool(*canvas);
+
+  tool->draw(x0, y0, x1, y1);
+
+  bool pass = true;
+  for (int x = 0; x < size; x++) {
+    for (int y = 0; y < size; y++) {
+      double distance =
+          std::abs(std::sqrt(std::pow(x0 - x, 2) + std::pow(y0 - y, 2)) - r);
+      // Because of rounding errors, an exact test (for all pixels) is not
+      // feasible.
+      // Therefore, it is only tested if set pixels have a distance <= 1.
+      if (canvas->get_pixel(x, y) && distance > 1) {
+        pass = false;
+      }
+    }
+  }
+  REQUIRE(pass);
+}
+
+TEST_CASE("Barycentric coordinates: Edge cases") {
+  int x0 = 0, y0 = 0, x1 = 0, y1 = 10, x2 = 10, y2 = 0;
+  float b1, b2, b3;
+
+  // point on vertex
+  std::tie(b1, b2, b3) = barycentric_coordinates(x0, y0, x1, y1, x2, y2, 0, 0);
+  REQUIRE(b1 == 1);
+  REQUIRE(b2 == 0);
+  REQUIRE(b3 == 0);
+
+  std::tie(b1, b2, b3) = barycentric_coordinates(x0, y0, x1, y1, x2, y2, 0, 10);
+  REQUIRE(b1 == 0);
+  REQUIRE(b2 == 1);
+  REQUIRE(b3 == 0);
+
+  std::tie(b1, b2, b3) = barycentric_coordinates(x0, y0, x1, y1, x2, y2, 10, 0);
+  REQUIRE(b1 == 0);
+  REQUIRE(b2 == 0);
+  REQUIRE(b3 == 1);
+
+  // point on edge
+  std::tie(b1, b2, b3) = barycentric_coordinates(x0, y0, x1, y1, x2, y2, 0, 5);
+  REQUIRE(b1 == 0.5);
+  REQUIRE(b2 == 0.5);
+  REQUIRE(b3 == 0);
+
+  std::tie(b1, b2, b3) = barycentric_coordinates(x0, y0, x1, y1, x2, y2, 5, 0);
+  REQUIRE(b1 == 0.5);
+  REQUIRE(b2 == 0);
+  REQUIRE(b3 == 0.5);
+
+  std::tie(b1, b2, b3) = barycentric_coordinates(x0, y0, x1, y1, x2, y2, 5, 5);
+  REQUIRE(b1 == 0);
+  REQUIRE(b2 == 0.5);
+  REQUIRE(b3 == 0.5);
+
+  // All points on straight line
+  std::tie(b1, b2, b3) = barycentric_coordinates(0, y0, 0, y1, 0, y2, 0, 0);
+  REQUIRE(std::isnan(b1));
+  REQUIRE(std::isnan(b2));
+  REQUIRE(std::isnan(b3));
+
+  std::tie(b1, b2, b3) = barycentric_coordinates(x0, 0, x1, 0, x2, 0, 0, 0);
+  REQUIRE(std::isnan(b1));
+  REQUIRE(std::isnan(b2));
+  REQUIRE(std::isnan(b3));
+}
+
+TEST_CASE("Barycentric coordinates (prop: Σ = 1)") {
+  int x0 = GENERATE(take(2, random(-100, 100)));
+  int y0 = GENERATE(take(2, random(-100, 100)));
+  int x1 = GENERATE(take(2, random(-100, 100)));
+  int y1 = GENERATE(take(2, random(-100, 100)));
+  int x2 = GENERATE(take(2, random(-100, 100)));
+  int y2 = GENERATE(take(2, random(-100, 100)));
+
+  int x = GENERATE(take(5, random(-100, 100)));
+  int y = GENERATE(take(5, random(-100, 100)));
+  float b1, b2, b3;
+
+  std::tie(b1, b2, b3) = barycentric_coordinates(x0, y0, x1, y1, x2, y2, x, y);
+
+  // If all points are on a straight line, the property does not hold
+  if (!(x0 == x1 && x1 == x2) && !(y0 == y1 && y1 == y2)) {
+    REQUIRE_THAT(b1 + b2 + b3, WithinRel(1.0, 0.01));
+  }
+}
+
+TEST_CASE("Sort triangle points") {
+  int x0 = 10, y0 = 60, x1 = 50, y1 = 90, x2 = 40, y2 = 30;
+  sort_triangle_points(x0, y0, x1, y1, x2, y2);
+  REQUIRE(x0 == 40);
+  REQUIRE(y0 == 30);
+  REQUIRE(x1 == 10);
+  REQUIRE(y1 == 60);
+  REQUIRE(x2 == 50);
+  REQUIRE(y2 == 90);
+}
+
+TEST_CASE("Sort triangle points (prop: y0 < y1 < y2)") {
+  int x0 = GENERATE(take(3, random(-100, 100)));
+  int y0 = GENERATE(take(3, random(-100, 100)));
+  int x1 = GENERATE(take(3, random(-100, 100)));
+  int y1 = GENERATE(take(3, random(-100, 100)));
+  int x2 = GENERATE(take(3, random(-100, 100)));
+  int y2 = GENERATE(take(3, random(-100, 100)));
+
+  sort_triangle_points(x0, y0, x1, y1, x2, y2);
+
+  REQUIRE(y0 <= y1);
+  REQUIRE(y1 <= y2);
+}
+
+TEST_CASE("Slope") {
+  REQUIRE(slope(5, 10, 20, 10) == 0.0);
+  REQUIRE(slope(0, 0, 10, 10) == 1.0);
+  REQUIRE(slope(0, 0, 10, -10) == -1.0);
+  REQUIRE(slope(0, 0, 10, 5) == 0.5);
+  REQUIRE(slope(0, 0, 10, -5) == -0.5);
+  REQUIRE(slope(0, 10, 10, 40) == 3.0);
+  REQUIRE(slope(0, 10, 10, -40) == -5.0);
+
+  // Special case: Infinite slope, must be normalized
+  REQUIRE(slope(10, 10, 10, 40) == 0.0);
+}
+
+TEST_CASE("Sweep line (prop: Barycentric coordinates)") {
+  const int size = 100;
+  int x0 = GENERATE(take(3, random(0, size - 1)));
+  int y0 = GENERATE(take(3, random(0, size - 1)));
+  int x1 = GENERATE(take(3, random(0, size - 1)));
+  int y1 = GENERATE(take(3, random(0, size - 1)));
+  int x2 = GENERATE(take(3, random(0, size - 1)));
+  int y2 = GENERATE(take(3, random(0, size - 1)));
+
+  canvas_buffer *canvas = new canvas_buffer(size, size);
+  sweep_line_tool *tool = new sweep_line_tool(*canvas);
+
+  tool->draw(x0, y0, x1, y1, x2, y2);
+
+  int deviating = 0;
+  bool pass = true;
+  for (int x = 0; x < size; x++) {
+    for (int y = 0; y < size; y++) {
+      if (point_in_triangle(x0, y0, x1, y1, x2, y2, x, y)) {
+        if (!canvas->get_pixel(x, y)) {
+          // Barycentric coordinates say, point is not in triangle,
+          // but point is not set.
+          // This must not happen → fail test.
+          pass = false;
+        }
+      } else if (canvas->get_pixel(x, y)) {
+        // Barycentric coordinates say, point is not in triangle,
+        // but point is set.
+        // The point is most likely on edge → mark it as deviating.
+        deviating++;
+      }
+    }
+  }
+  REQUIRE(pass);
+  // Crude heuristic:
+  // No more than differences of all edge point coordinates can fail.
+  // This ist not accurate (false positives possible) on small/spiky triangles,
+  // but overall it gives an okayish result.
+  REQUIRE(deviating < abs(y1 - y0) + abs(y2 - y1) + abs(y0 - y2) +
+                          abs(x1 - x0) + abs(x2 - x1) + abs(x0 - x2));
+}

u02/src/util.cpp

@@ -0,0 +1,207 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+#include "util.h"
+#include <cmath>
+
+void transform_mut(Transformation transformation, int &x, int &y) {
+  if (transformation & TRANSFORM_ROTATE_CW) {
+    std::swap(x, y);
+    x = -x;
+  }
+  if (transformation & TRANSFORM_ROTATE_CCW) {
+    std::swap(x, y);
+    y = -y;
+  }
+  if (transformation & TRANSFORM_MIRROR_X) {
+    y = -y;
+  }
+  if (transformation & TRANSFORM_MIRROR_Y) {
+    x = -x;
+  }
+}
+
+std::pair<int, int> transform(Transformation transformation, int x, int y) {
+  transform_mut(transformation, x, y);
+  return std::make_pair(x, y);
+}
+
+void transform_inv_mut(Transformation transformation, int &x, int &y) {
+  if (transformation & TRANSFORM_MIRROR_Y) {
+    x = -x;
+  }
+  if (transformation & TRANSFORM_MIRROR_X) {
+    y = -y;
+  }
+  if (transformation & TRANSFORM_ROTATE_CCW) {
+    // does clockwise rotation
+    std::swap(x, y);
+    x = -x;
+  }
+  if (transformation & TRANSFORM_ROTATE_CW) {
+    // does counterclockwise rotation
+    std::swap(x, y);
+    y = -y;
+  }
+}
+
+std::pair<int, int> transform_inv(Transformation transformation, int x, int y) {
+  transform_inv_mut(transformation, x, y);
+  return std::make_pair(x, y);
+}
+
+/*
+ * After it took me many hours to get this right,
+ * I at least want to document how I got to it:
+ *
+ * N.B. In the following,
+ * Modulo is *not* the remainder of the euclidean division,
+ * but instead the remainder of truncated division
+ * (i.e., negative quotients produce negative results).
+ *
+ * There are two main cases:
+ * The simple one is where angle % 90° ≤ 45°.
+ * To transform this into the standard case,
+ * only mirrors are needed.
+ * The more complicated is when angle % 90° ≥ 45°.
+ * To transform this into the standard case,
+ * a rotation has to be done, followed by a mirror in some cases.
+ *
+ * The following matrices show what must be done when.
+ * The ASCII art arrows show the line as it should be drawn,
+ * the the column/row headings show how they can be identified in code,
+ * the capital letters in the field show what needs to be done
+ * to reach the standard case
+ * (X/Y: mirror X/Y; CW/CCW: rotate CW/CCW).
+ * Because there is no nice way to draw arrows with an angle < 45°,
+ * they are just differentiated by the heading.
+ *
+ * Let (x_0, y_0) be the starting point and (x_1, y_1) the end point.
+ * Let Δx = x_1 - x_0, Δy = y_1 - y_0.
+ * Let m = Δy/Δx, α = atan(m).
+ *
+ * α ≤ 45°:
+ *
+ *      Δx>0   Δx<0
+ *
+ *         A | A
+ * Δy<0   /  |  \
+ *       /   |   \
+ *      /  - | Y  \
+ *     ------+------
+ *      \  X | XY /
+ *       \   |   /
+ * Δy>0   \  |  /
+ *         V | V
+ *
+ * α ≥ 45°:
+ *
+ *      Δx>0   Δx<0                       Δx>0   Δx<0
+ *
+ *         A | A                          \    |    A
+ * Δy<0   /  |  \                    Δy<0  \   |   /
+ *       /   |   \                          \  |  /
+ *      / CW | CW \                       X  V | /
+ *    -------+------- → after rotation → ------+-----
+ *     \ CCW | CCW /                         A | \
+ *      \    |    /                         /  |  \
+ * Δy>0  \   |   /                   Δy>0  /   |   \
+ *        V  |  V                         /    | X  V
+ */
+Transformation transformation_to_standard_case(int x0, int y0, int x1, int y1) {
+  Transformation transformation = 0;
+
+  int delta_y = y1 - y0;
+  int delta_x = x1 - x0;
+  // checks if angle ∈ (-90°, 90°) is ≥ 45°
+  // this is a simplified version of atan(Δy/Δx) > π/4:
+  //   atan(Δy/Δx) > π/4 | tan(…)
+  // ⇔ Δx/Δy       > 1 | Δy
+  // ⇔ Δx          > Δy
+  if (abs(delta_y) > abs(delta_x)) {
+    // if-else is needed, because of special case Δy = 0
+    if (delta_y < 0) {
+      transformation |= TRANSFORM_ROTATE_CW;
+    } else if (delta_y > 0) {
+      transformation |= TRANSFORM_ROTATE_CCW;
+    }
+    // the sign of Δx and Δy (pre-rotation!) differ,
+    // an additional mirror is needed
+    if (delta_x * delta_y < 0) {
+      transformation |= TRANSFORM_MIRROR_X;
+    }
+  } else {
+    if (delta_x < 0) {
+      transformation |= TRANSFORM_MIRROR_Y;
+    }
+    if (delta_y > 0) {
+      transformation |= TRANSFORM_MIRROR_X;
+    }
+  }
+
+  return transformation;
+}
+
+std::tuple<float, float, float> barycentric_coordinates(int x0, int y0, int x1,
+                                                        int y1, int x2, int y2,
+                                                        int xp, int yp) {
+  // Source:
+  // https://en.wikipedia.org/wiki/Barycentric_coordinate_system#Vertex_approach
+  float b1 = x1 * y2 - x2 * y1 + xp * (y1 - y2) + yp * (x2 - x1);
+  float b2 = x2 * y0 - x0 * y2 + xp * (y2 - y0) + yp * (x0 - x2);
+  float b3 = x0 * y1 - x1 * y0 + xp * (y0 - y1) + yp * (x1 - x0);
+
+  // reciprocal computed directly for performance
+  float area_factor =
+      1 / static_cast<float>(x0 * (y1 - y2) + x1 * (y2 - y0) + x2 * (y0 - y1));
+
+  b1 *= area_factor;
+  b2 *= area_factor;
+  b3 *= area_factor;
+
+  return {b1, b2, b3};
+}
+
+bool point_in_triangle(int x0, int y0, int x1, int y1, int x2, int y2, int xp,
+                       int yp) {
+  float b1, b2, b3;
+  std::tie(b1, b2, b3) =
+      barycentric_coordinates(x0, y0, x1, y1, x2, y2, xp, yp);
+  return b1 >= 0.0 && b1 <= 1.0 && b2 >= 0.0 && b2 <= 1.0 && b3 >= 0.0 &&
+         b3 <= 1.0;
+}
+
+void sort_triangle_points(int &x0, int &y0, int &x1, int &y1, int &x2,
+                          int &y2) {
+  // Bubble sort is not really ideal in general.
+  // It could be changed to use a more efficient algorithm,
+  // but for only 3 values, it should suffice.
+  // Moreover, implementing sorting on an array/vector of tuples
+  // is probably more overhead.
+  if (y0 > y1) {
+    std::swap(x0, x1);
+    std::swap(y0, y1);
+  }
+  if (y0 > y2) {
+    std::swap(x0, x2);
+    std::swap(y0, y2);
+  }
+  if (y1 > y2) {
+    std::swap(x1, x2);
+    std::swap(y1, y2);
+  }
+}
+
+float slope(int x0, int y0, int x1, int y1) {
+  float m = ((float)(y1 - y0)) / ((float)(x1 - x0));
+  if (std::isinf(m) || std::isnan(m)) {
+    // This is a special case for two things:
+    //
+    // IEEE 754 specifies ∞ × 0 / 0 × ∞ to be an invalid operation,
+    // and therefore return NaN.
+    // That makes the computation of Δy fail when x0 == x1.
+    //
+    // In the case that additionally y0 == y1,
+    // the expression is 0/0, also defined in IEEE 754 as invalid.
+    m = 0;
+  }
+  return m;
+}