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[0] ~~ Real && @a[1] ~~ 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[0] * \$x + @coefficients[1] * \$y, @coefficients[2], "ax + by = c"; is @coefficients[3] * \$x + @coefficients[4] * \$y, @coefficients[5], "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[0] == @coefficients[1] * @coefficients[4] / @coefficients[3], "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; }