// 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 <iostream>

#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 "star_tool.h"
#include "sweep_line_tool.h"
#include "util.h"

using Catch::Matchers::WithinAbs;

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:
  default: // required to make static compiler warnings happy
    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);
  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:
  default: // required to make static compiler warnings happy
    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)) {
    // catch2’s SKIP macro does not exactly do what I want here,
    // so this just returns from the function and prints a message.
    std::cerr
        << "All coordinates have extreme value, skipping (avoid rounding error)"
        << std::endl;
    return;
  }

  const float 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++) {
      float 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)));

  // If all points are on a straight line, the property does not hold.
  // Checking this is equivalent to checking if the area of the triangle is 0.
  if ((x0 - x1) * (y0 - y2) - (y0 - y1) * (x0 - x2) == 0) {
    // see above
    std::cerr << "Points do not form reasonable triangle, skipping"
              << std::endl;
    return;
  }

  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);

  REQUIRE_THAT(b1 + b2 + b3, WithinAbs(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 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 deviate.
  // This ist not accurate (false negatives 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));
}

TEST_CASE("Star (prop: all points inside circle)") {
  // See test “Bresenham circle (prop: √(x²+y²)-r<ε)” for what this does.
  const int size = 100;
  const int max_c = std::floor(size / (3 - std::sqrt(2)));
  const int min_c = std::ceil(((4 - std::sqrt(2)) / 7) * size);

  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)) {
    // see above
    std::cerr
        << "All coordinates have extreme value, skipping (avoid rounding error)"
        << std::endl;
    return;
  }

  const float r = std::sqrt(std::pow((x1 - x0), 2) + std::pow((y1 - y0), 2));

  // correction factor for very small stars
  const float correction = r < 2 ? 1 : 0;

  canvas_buffer *canvas = new canvas_buffer(size, size);
  star_tool *tool = new star_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 (canvas->get_pixel(x, y) &&
          (std::sqrt(std::pow(x0 - x, 2) + std::pow(y0 - y, 2)) >
           r * r + correction)) {
        pass = false;
      }
    }
  }
  REQUIRE(pass);
}

TEST_CASE("Star (prop: average of all points is centre)") {
  const int x0 = GENERATE(take(5, random(-100, 100)));
  const int y0 = GENERATE(take(5, random(-100, 100)));
  const int n = GENERATE(take(5, random(3, 1000)));
  const float r2 = GENERATE(take(5, random(0, 100)));
  const float r_factor = GENERATE(take(5, random(0, 1)));
  const float r1 = r_factor * r2;
  const float angle = GENERATE(take(5, random(-16.0 * atan(1), 16 * atan(1))));

  std::vector<std::pair<float, float>> points = star(n, r1, r2, x0, y0, angle);

  float x = 0, y = 0;
  for (std::pair<float, float> point : points) {
    x += point.first;
    y += point.second;
  }

  REQUIRE_THAT(x / points.size(), WithinAbs(x0, 0.5));
  REQUIRE_THAT(y / points.size(), WithinAbs(y0, 0.5));
}

// The Haskell implementation of regular_polygon_mod has many more tests,
// but they are not implemented here for brevity.