p2-colomon

use v6;
use Test;

# borrowed from List::Utils
sub sliding-window-wrapped(@a, \$n) {
my \$a-list = @a.iterator.list;
my @values;
gather {
while defined(my \$a = \$a-list.shift) {
@values.push(\$a);
@values.shift if +@values > \$n;
take @values if +@values == \$n;
}

for ^(\$n-1) {
@values.push(@a[\$_]);
@values.shift if +@values > \$n;
take @values if +@values == \$n;
}
}
}

# solves
# ax+by = c
# dx+ey = f
sub solve-two-variables(\$a, \$b, \$c, \$d, \$e, \$f)
{
my (\$x, \$y);
if \$a.abs < 1e-10 {
return Any if \$b.abs < 1e-10 || \$d.abs < 1e-10;
\$y = \$c / \$b;
\$x = (\$f - \$e * \$y) / \$d;
return \$x, \$y;
}
if \$b.abs < 1e-10 {
return Any if \$a.abs < 1e-10 || \$e.abs < 1e-10;
\$x = \$c / \$a;
\$y = (\$f - \$d * \$x) / \$e;
return \$x, \$y;
}
if \$d.abs < 1e-10 {
return Any if \$e.abs < 1e-10 || \$a.abs < 1e-10;
\$y = \$f / \$e;
\$x = (\$c - \$b * \$y) / \$a;
return \$x, \$y;
}
if \$e.abs < 1e-10 {
return Any if \$d.abs < 1e-10 || \$b.abs < 1e-10;
\$x = \$f / \$d;
\$y = (\$c - \$a * \$x) / \$b;
return \$x, \$y;
}

my \$denom = \$b - \$a * \$e / \$d;
return Any if \$denom.abs < 1e-10;

\$y = (\$c - \$a * \$f / \$d) / \$denom;
\$x = (\$c - \$b * \$y) / \$a;
\$x, \$y;
}

# x * d1 - y * d2 = b2 - b1
sub lines-intersect(Complex \$base1, Complex \$direction1, Complex \$base2, Complex \$direction2) {
solve-two-variables(\$direction1.re, -\$direction2.re, \$base2.re - \$base1.re,
\$direction1.im, -\$direction2.im, \$base2.im - \$base1.im);
}

sub segments-intersect(Complex \$start1, Complex \$end1, Complex \$start2, Complex \$end2) {
my (\$s, \$t) = lines-intersect(\$start1, \$end1 - \$start1, \$start2, \$end2 - \$start2);
\$s.defined ?? \$s ~~ 0 ^..^ 1 && \$t ~~ 0 ^..^ 1 !! Bool;
}

sub point-distance-to-line(Complex \$point, Complex \$start, Complex \$end) {
((\$end.re - \$start.re) * (\$start.im - \$point.im)
- (\$start.re - \$point.re) * (\$end.im - \$start.im)).abs
/ (\$end - \$start).abs;
}

sub intersections-line-and-loop(\$p1, \$p2, @loop) {
# if the line we are intersecting comes too close to one of the loop's
# vertices, this algorithm cannot be trusted, so we just return Any
for @loop -> \$point {
return Any if point-distance-to-line(\$point, \$p1, \$p2) < 1e-15;
}

my \$intersections = 0;
for sliding-window-wrapped(@loop, 2) -> \$start, \$end {
# segments-intersect should only return undefined if the segments are parallel
\$intersections += segments-intersect(\$p1, \$p2, \$start, \$end) // 0;
}
\$intersections;
}

sub is-point-inside-loop(\$point, @loop) {
# If the point is really close to one of the loop points,
# then we declare it is not in the loop.
for @loop -> \$p {
return Bool::False if (\$point - \$p).abs < 1e-15;
}

# find a circle enclosing the loop
my \$u-range = @loop.map(*.re).minmax;
my \$v-range = @loop.map(*.im).minmax;
my \$safe-distance = \$u-range.max - \$u-range.min + \$v-range.max - \$v-range.min + 1;
my \$center = Complex.new(\$u-range.max + \$u-range.min / 2, \$v-range.max + \$v-range.min / 2);

# now intersect the line from the given point to random points on that circle until
# we find one that does not intersect one of the loop's vertices
my \$intersections;
loop {
my \$outside = \$center + \$safe-distance.unpolar((0 .. 2 * pi).rand);
\$intersections = intersections-line-and-loop(\$point, \$outside, @loop);
last if \$intersections.defined;
}

# if the number of intersections is odd, we have a point inside the loop
\$intersections % 2 == 1;
}

sub validate-point-string(\$point-string) {
my @a = \$point-string.eval;
+@a == 2 && @a ~~ Real && @a ~~ Real;
}

sub validate-point-strings(@point-strings) {
my @invalid-point-strings = @point-strings.grep({ !validate-point-string(\$_) });
+@point-strings > 2 && !@invalid-point-strings;
}

multi MAIN() {
my @point-strings = \$*IN.get.comb(/\S+/);
unless validate-point-strings(@point-strings) {
say "invalid input";
return;
}
my @loop = @point-strings.map({ Complex.new(|\$_.eval); });

my \$point-string = \$*IN.get;
unless validate-point-string(\$point-string) {
say "invalid input";
return;
}
my \$point = Complex.new(|\$point-string.eval);

say is-point-inside-loop(\$point, @loop) ?? "yes" !! "no";
}

multi MAIN("test") {
plan *;

{
# ax+by = c
# dx+ey = f

is ~solve-two-variables(0, 1, 2, 1, 0, 3), ~(3, 2), "y = 2, x = 3";
is ~solve-two-variables(1, 0, 3, 0, 1, 2), ~(3, 2), "x = 3, y = 2";
is ~solve-two-variables(0, 1, 2, 1, 1, 3), ~(1, 2), "y = 2, x + y = 3";
is ~solve-two-variables(1, 1, 3, 0, 1, 2), ~(1, 2), "x + y = 3, y = 2";

nok solve-two-variables(0, 0, 3, 0, 1, 2).defined, "No solution for 0x + 0y = 3, y = 2";
nok solve-two-variables(0, 12, 3, 0, 1, 2).defined, "No solution for 12y = 3, y = 2";
nok solve-two-variables(1, 0, 3, 23, 0, 2).defined, "No solution for x = 3, 23x = 2";
nok solve-two-variables(2, 2, 3, 1, 1, 2).defined, "No solution for 2x + 2y = 3, x + y = 2";

for ^5 {
my @coefficients = (^6).map({ (-10..10).roll / (1..10).roll });

my (\$x, \$y) = solve-two-variables(|@coefficients);
if \$x.defined {
# is solution correct?
is @coefficients * \$x + @coefficients * \$y, @coefficients, "ax + by = c";
is @coefficients * \$x + @coefficients * \$y, @coefficients, "dx + ey = f";

is ~solve-two-variables(|@coefficients[3,4,5,0,1,2]), ~(\$x, \$y),
"Swapping abc and def doesn't change answer";
is ~solve-two-variables(|@coefficients[1,0,2,4,3,5]), ~(\$y, \$x),
"Swapping a for b and d for e switches x & y";
is ~solve-two-variables(|@coefficients[4,3,5,1,0,2]), ~(\$y, \$x),
"Swapping bac and edf switches x & y";
} else {
nok \$y.defined, "x and y are both undefined";
ok +@coefficients.grep(* == 0)
|| @coefficients == @coefficients * @coefficients / @coefficients,
"Found a reason to be undefined";
}
}
}

{
for 1 + 0i, 1i, 1 + 1i -> \$dir {
nok lines-intersect(Complex.new(0, 0), \$dir, Complex.new(-1, 1), \$dir).defined,
"Parallel lines return undefined";
nok lines-intersect(Complex.new(0, 0), \$dir, Complex.new(0, 0), \$dir).defined,
"Idential lines return undefined";

my (\$s, \$t) = lines-intersect(Complex.new(0, 0), \$dir, Complex.new(0, 1/2), 1 - 1i);
is_approx Complex.new(0, 0) + \$s * \$dir, Complex.new(0, 1/2) + \$t * (1 - 1i),
"Lines intersect as reported";
ok segments-intersect(Complex.new(0, 0), \$dir, Complex.new(-1, 1.5), Complex.new(1, -0.5)),
"Segments do intersect";
nok segments-intersect(Complex.new(0, 0), \$dir / 10, Complex.new(-1, 1.5), Complex.new(1, -0.5)),
"Segments do not intersect";
}

nok segments-intersect(Complex.new(0, 0), Complex.new(1, 0), Complex.new(1/2, 0), Complex.new(1/2, 1)),
"Segments that touch exactly do not intersect";
nok segments-intersect(Complex.new(0, 0), Complex.new(1, 0), Complex.new(1/2, 1), Complex.new(1/2, 0)),
"Segments that touch exactly do not intersect";
nok segments-intersect(Complex.new(1/2, 0), Complex.new(1/2, 1), Complex.new(0, 0), Complex.new(1, 0)),
"Segments that touch exactly do not intersect";
nok segments-intersect(Complex.new(1/2, 1), Complex.new(1/2, 0), Complex.new(0, 0), Complex.new(1, 0)),
"Segments that touch exactly do not intersect";
}

{
is_approx point-distance-to-line(Complex.new(0, 0),
Complex.new(-1, 0), Complex.new(1, 0)),
0, "Point on line is distance 0";
is_approx point-distance-to-line(Complex.new(0, 1),
Complex.new(-1, 0), Complex.new(1, 0)),
1, "Point (0, 1) is distance 1";
is_approx point-distance-to-line(Complex.new(0, 5/2),
Complex.new(-1, 0), Complex.new(1, 0)),
5/2, "Point (0, 5/2) is distance 5/2";
is_approx point-distance-to-line(Complex.new(0, -pi),
Complex.new(-1, 0), Complex.new(1, 0)),
pi, "Point (0, -pi) is distance pi";

my \$a = Complex.new(2, 3);
my \$b = Complex.new(-13, 2);
my \$v = Complex.new(\$b.im - \$a.im, \$a.re - \$b.re);
\$v /= \$v.abs;

is_approx point-distance-to-line(\$a, \$a, \$b), 0, "Distance to first point is 0";
is_approx point-distance-to-line(\$b, \$a, \$b), 0, "Distance to second point is 0";
is_approx point-distance-to-line((\$a + \$b) / 2, \$a, \$b), 0, "Distance to mid point is 0";
is_approx point-distance-to-line(\$a + \$v, \$a, \$b), 1, "Distance to first point plus perp is 1";
is_approx point-distance-to-line(\$a + 3/2 * \$v, \$a, \$b), 3/2, "Distance to first point plus perp is 3/2";
is_approx point-distance-to-line(\$a - pi * \$v, \$a, \$b), pi, "Distance to first point plus perp is pi";
is_approx point-distance-to-line(\$b + \$v, \$a, \$b), 1, "Distance to second point plus perp is 1";
is_approx point-distance-to-line(\$b + 3/2 * \$v, \$a, \$b), 3/2, "Distance to second point plus perp is 3/2";
is_approx point-distance-to-line(\$b - pi * \$v, \$a, \$b), pi, "Distance to second point plus perp is pi";
}

{
my @loop = Complex.new(1.0,1.0),
Complex.new(1.0,-1.0),
Complex.new(-1.0,-1.0),
Complex.new(-1.0,1.0);
is intersections-line-and-loop(Complex.new(0, 0), Complex.new(2,0), @loop), 1,
"Got one intersection";
is intersections-line-and-loop(Complex.new(-2, 0), Complex.new(2,0), @loop), 2,
"Got two intersections";
nok intersections-line-and-loop(Complex.new(0, 0), Complex.new(2,2), @loop).defined,
"Got undefined, because line intersects a vertex";
}

for 0, 10 + 10i, 100 - 100i -> \$shift {
my @loop = Complex.new(1.0,1.0) + \$shift,
Complex.new(1.0,-1.0) + \$shift,
Complex.new(-1.0,-1.0) + \$shift,
Complex.new(-1.0,1.0) + \$shift;
ok is-point-inside-loop(Complex.new(0, 0) + \$shift, @loop), "(0, 0) is in loop";
ok is-point-inside-loop(Complex.new(0.0, 0.0) + \$shift, @loop), "(0.0, 0.0) is in loop";
ok is-point-inside-loop(Complex.new(0.Num, 0.Num) + \$shift, @loop), "(0.Num, 0.Num) is in loop";
ok is-point-inside-loop(Complex.new(0.9, 0.9) + \$shift, @loop), "(0.9, 0.9) is in loop";
ok is-point-inside-loop(Complex.new(-0.9, -0.9) + \$shift, @loop), "(-0.9, -0.9) is in loop";

# points on the vertices of the loop
nok is-point-inside-loop(Complex.new(1.0, 1.0) + \$shift, @loop), "(1.0, 1.0) is not in loop";
nok is-point-inside-loop(Complex.new(-1.0, 1.0) + \$shift, @loop), "(-1.0, 1.0) is not in loop";
nok is-point-inside-loop(Complex.new(1.0, -1.0) + \$shift, @loop), "(1.0, -1.0) is not in loop";
nok is-point-inside-loop(Complex.new(-1.0, -1.0) + \$shift, @loop), "(-1.0, -1.0) is not in loop";

# clearly outside the loop
nok is-point-inside-loop(Complex.new(1.1, 1.1) + \$shift, @loop), "(1.1, 1.1) is not in loop";
nok is-point-inside-loop(Complex.new(-1.1, 1.1) + \$shift, @loop), "(-1.1, 1.1) is not in loop";
nok is-point-inside-loop(Complex.new(1.1, -1.1) + \$shift, @loop), "(1.1, -1.1) is not in loop";
nok is-point-inside-loop(Complex.new(-1.1, -1.1) + \$shift, @loop), "(-1.1, -1.1) is not in loop";
nok is-point-inside-loop(Complex.new(1.1, 0.0) + \$shift, @loop), "(1.1, 0.0) is not in loop";
nok is-point-inside-loop(Complex.new(-1.1, 0.0) + \$shift, @loop), "(-1.1, 0.0) is not in loop";
nok is-point-inside-loop(Complex.new(0.0, -1.1) + \$shift, @loop), "(0.0, -1.1) is not in loop";
nok is-point-inside-loop(Complex.new(0.0, -1.1) + \$shift, @loop), "(0.0, -1.1) is not in loop";
}

{
my @loop = Complex.new(1.0, 1.0),
Complex.new(1.0, 1/4),
Complex.new(3/4, 1/4),
Complex.new(3/4, 3/4),
Complex.new(-3/4, 3/4),
Complex.new(-3/4, -3/4),
Complex.new(3/4, -3/4),
Complex.new(3/4, -1/4),
Complex.new(1, -1/4),
Complex.new(1.0, -1.0),
Complex.new(-1.0, -1.0),
Complex.new(-1.0, 1.0);

# actually print out the inside of this loop
# for -17/16, * + 1/8 ... * > 1 -> \$y {
#     for -17/16, * + 1/8 ... * > 1 -> \$x {
#         if is-point-inside-loop(\$x + \$y\i, @loop) {
#             print "*";
#         } else {
#             print " ";
#         }
#     }
#     say " \$y";
# }

nok is-point-inside-loop(Complex.new(0, 0), @loop), "(0, 0) is not in loop";

nok is-point-inside-loop(Complex.new(1/2, 0), @loop), "(1/2, 0) is not in loop";
nok is-point-inside-loop(Complex.new(-1/2, 0), @loop), "(-1/2, 0) is not in loop";
nok is-point-inside-loop(Complex.new(0, 1/2), @loop), "(0, 1/2) is not in loop";
nok is-point-inside-loop(Complex.new(0, -1/2), @loop), "(0, -1/2) is not in loop";

nok is-point-inside-loop(Complex.new(7/8, 0), @loop), "(7/8, 0) is not in loop";
ok is-point-inside-loop(Complex.new(-7/8, 0), @loop), "(-7/8, 0) is in loop";
ok is-point-inside-loop(Complex.new(0, 7/8), @loop), "(0, 7/8) is in loop";
ok is-point-inside-loop(Complex.new(0, -7/8), @loop), "(0, -7/8) is in loop";

ok is-point-inside-loop(Complex.new(0.9, 0.9), @loop), "(0.9, 0.9) is in loop";
ok is-point-inside-loop(Complex.new(-0.9, -0.9), @loop), "(-0.9, -0.9) is in loop";

for @loop -> \$vertex {
nok is-point-inside-loop(\$vertex, @loop), "Vertex \$vertex is not in loop";
}

# clearly outside the loop
nok is-point-inside-loop(Complex.new(1.1, 1.1), @loop), "(1.1, 1.1) is not in loop";
nok is-point-inside-loop(Complex.new(-1.1, 1.1), @loop), "(-1.1, 1.1) is not in loop";
nok is-point-inside-loop(Complex.new(1.1, -1.1), @loop), "(1.1, -1.1) is not in loop";
nok is-point-inside-loop(Complex.new(-1.1, -1.1), @loop), "(-1.1, -1.1) is not in loop";
nok is-point-inside-loop(Complex.new(1.1, 0.0), @loop), "(1.1, 0.0) is not in loop";
nok is-point-inside-loop(Complex.new(-1.1, 0.0), @loop), "(-1.1, 0.0) is not in loop";
nok is-point-inside-loop(Complex.new(0.0, -1.1), @loop), "(0.0, -1.1) is not in loop";
nok is-point-inside-loop(Complex.new(0.0, -1.1), @loop), "(0.0, -1.1) is not in loop";
}

done;
}

The code is a pleasant read, especially is-point-inside-loop starting on line 95. Many of the subroutines are well-named and straightforward.

Consistency

Line 122 contains an eval of user input, as part of validating said input. Oops.

segments-intersect is a case of a method possibly doing more than it claims. It has three return values: True (the segments intersect), False (the segments don't intersect), and Bool (the segments don't intersect because they are parallel). The last kind of return value isn't used for anything in the program.

Clarity of intent

The sub solve-two-variables has 13 occurrences of the same numeric literal. That's 12 infractions of DRY. :) In fact, the literal occurs in the same syntactic surroundings each time, so maybe a whole subroutine is missing.

There's another, smaller numeric literal occuring twice in two different subs. Is it "the same" measurement? Hard to tell, which is why it's mentioned here.

No harm in using Complex for representing 2D coordinate tuples (as opposed to, say, using the type for its arithmetic properties) -- but it wouldn't hurt to document the "mapping" somewhere along with the code. (.unpolar is used once, but not in any way that uses the features of complex numbers.)

Algorithmic efficiency

Efficiency isn't really a big issue in this one. I see nothing slower than linear.

Idiomatic use of Perl 6

.comb(/\S+/) on line 132 can also be spelled .words.

(0 .. 2 * pi).rand on line 112 doesn't cause the distribution intended. Ranges don't have a dedicated .rand method, so the range numifies to 7 before .randing:

\$ perl6 -e 'say (0..2 * pi).rand'
6.55095617280338

Brevity

This is not the shortest solution imaginable. It makes up for it somewhat by being methodical. The test-to-algorithm ratio is admirable.