euler165.cpp

#include <cstdint>
#include <iostream>
#include <set>
#include <vector>

#include <boost/rational.hpp>

using ext_coord = std::int64_t;
using int_coord = boost::rational<ext_coord>;

struct ext_point { ext_coord x, y; };
struct int_point { int_coord x, y; };
struct segment   { ext_point p, q; };

ext_point operator + (const ext_point & p, const ext_point & q)
{
    return ext_point { p.x + q.x, p.y + q.y };
}

ext_point operator - (const ext_point & p, const ext_point & q)
{
    return ext_point { p.x - q.x, p.y - q.y };
}

ext_point operator * (const ext_point & p, ext_coord k)
{
    return ext_point { p.x * k, p.y * k };
}

int_point operator / (const ext_point & p, ext_coord k)
{
    return int_point { int_coord(p.x, k), int_coord(p.y, k) };
}

bool operator < (const int_point & p, const int_point & q)
{
    if (p.x < q.x) return true;
    if (p.x > q.x) return false;
    return p.y < q.y;
}

bool unaligned(ext_coord b, ext_coord c)
{
    if (c < 0) return b >= 0 || b <= c;
    if (c > 0) return b <= 0 || b >= c;
    return true;
}

void intersect(const segment & s, const segment & t, std::set<int_point> & ps)
{
    ext_point dp1 = s.p - s.q;  // a,d
    ext_point dp2 = t.p - t.q;  // b,e
    ext_point del = s.p - t.p;  // c,f
    
    // ae - bd
    ext_coord den = dp1.x * dp2.y - dp1.y * dp2.x;
    if (den == 0)
        return;
    
    // ce - bf
    ext_coord num1 = dp2.y * del.x - dp2.x * del.y;
    if (unaligned(num1, den))
        return;
    
    // cd - af
    ext_coord num2 = dp1.y * del.x - dp1.x * del.y;
    if (unaligned(num2, den))
        return;
    
    ps.insert((s.p * (den-num1) + s.q * num1) / den);
}

ext_coord generate(ext_coord & state)
{
    state *= state;
    state %= 50515093;
    return state % 500;
}

int main()
{
    std::vector<segment> ss;
    std::set<int_point> ps;
    ext_coord state = 290797;
    ss.reserve(5000);
    
    for (int i = 0; i < 5000; ++i)
    {
        segment t { { generate(state), generate(state) },
                    { generate(state), generate(state) } };
        for (auto & s : ss)
            intersect(t, s, ps);
        ss.emplace_back(t);
    }
    
    std::cout << ps.size() << std::endl;
    
    return 0;
}

euler165.hs

import Data.Int
import Data.List
import Data.Maybe
import Data.Ratio
import qualified Data.Set as S

type ExtCoord = Int64
type ExtPoint = [ExtCoord]      -- exactly two coordinates: x and y
type Segment  = [ExtPoint]      -- exactly two points: start and end

type IntCoord = Ratio ExtCoord
type IntPoint = [IntCoord]      -- exactly two coordinates: x and y

coords :: [ExtCoord]
coords = unfoldr next 290797
  where
    next s = Just (s' `mod` 500, s')
      where
        s' = s * s `mod` 50515093

segments :: [ExtCoord] -> [Segment]
segments (a:b:c:d:ts) = [[a,b], [c,d]] : segments ts
segments _            = []

unaligned :: ExtCoord -> ExtCoord -> Bool
unaligned b c
  | 0 < b && b < c = False
  | 0 > b && b > c = False
  | otherwise      = True

add, sub :: Segment -> ExtPoint
add [p1,p2] = zipWith (+) p1 p2
sub [p1,p2] = zipWith (-) p1 p2

mult :: ExtPoint -> Segment -> Segment
mult = zipWith $ map . (*)

intersection :: Segment -> Segment -> Maybe IntPoint
intersection s1@(p1:_) s2@(p2:_) =
  | unaligned num1 den = Nothing        -- Outside open segment s1
  | unaligned num2 den = Nothing        -- Outside open segment s2
  | otherwise          = Just $ (% den) `map` add nums
  where
    [a,d] = sub s1                      -- Solve for s,t:
    [b,e] = sub s2                      -- 
    [c,f] = sub [p1,p2]                 --   as - bt = c
    num1  = c*e - b*f                   --   ds - et = f
    num2  = c*d - a*f                   --   0 < s,t < 1
    den   = a*e - b*d                   -- 
    nums  = mult [den-num1, num1] s1    -- Solution must be unique

main :: IO ()
main = print . S.size . S.fromList . next . segments $ take 20000 coords
  where
    next []     = []                    -- 
    next (x:xs) = current ++ next xs    -- 
      where
        current = catMaybes $ intersection x `map` xs

euler165.rs

extern crate collections;
extern crate num;

use collections::treemap::TreeSet;
use num::rational::Ratio;

type ExtCoord = int;
type IntCoord = Ratio<ExtCoord>;

struct ExtPoint {
    x: ExtCoord,
    y: ExtCoord
}

#[deriving(Eq, Ord, TotalEq, TotalOrd)]
struct IntPoint {
    x: IntCoord,
    y: IntCoord
}

struct Segment {
    p: ExtPoint,
    q: ExtPoint
}

struct Generator {
    state: ExtCoord
}

impl Add<ExtPoint, ExtPoint> for ExtPoint {
    fn add(&self, rhs: &ExtPoint) -> ExtPoint {
        ExtPoint {
            x: self.x + rhs.x,
            y: self.y + rhs.y
        }
    }
}

impl Sub<ExtPoint, ExtPoint> for ExtPoint {
    fn sub(&self, rhs: &ExtPoint) -> ExtPoint {
        ExtPoint {
            x: self.x - rhs.x,
            y: self.y - rhs.y
        }
    }
}

impl Mul<ExtCoord, ExtPoint> for ExtPoint {
    fn mul(&self, &rhs: &ExtCoord) -> ExtPoint {
        ExtPoint {
            x: self.x * rhs,
            y: self.y * rhs
        }
    }
}

impl Div<ExtCoord, IntPoint> for ExtPoint {
    fn div(&self, &rhs: &ExtCoord) -> IntPoint {
        IntPoint {
            x: Ratio::new(self.x, rhs),
            y: Ratio::new(self.y, rhs)
        }
    }
}

impl Generator {
    fn new() -> Generator {
        Generator { state: 290797 }
    }
    
    fn next(&mut self) -> ExtCoord {
        self.state *= self.state;
        self.state %= 50515093;
        self.state % 500
    }
}

impl ExtPoint {
    fn new(gen: &mut Generator) -> ExtPoint {
        ExtPoint {
            x: gen.next(),
            y: gen.next()
        }
    }
}

fn unaligned(b: ExtCoord, c: ExtCoord) -> bool {
    if c < 0 { b >= 0 || b <= c } else
    if c > 0 { b <= 0 || b >= c } else { true }
}

impl Segment {
    fn new(gen: &mut Generator) -> Segment {
        Segment {
            p: ExtPoint::new(gen),
            q: ExtPoint::new(gen)
        }
    }
    
    fn intersect(&self, that: &Segment) -> Option<IntPoint> {
        let dp1 = self.p - self.q;
        let dp2 = that.p - that.q;
        let del = self.p - that.p;
        
        let den = dp1.x * dp2.y - dp1.y * dp2.x;
        if den == 0 { return None }
        
        let num1 = dp2.y * del.x - dp2.x * del.y;
        if unaligned(num1, den) { return None }
        
        let num2 = dp1.y * del.x - dp1.x * del.y;
        if unaligned(num2, den) { return None }
        
        Some((self.p * (den - num1) + self.q * num1) / den)
    }
}


fn main() {
    let mut gen = Generator::new();
    let mut ss  = Vec::with_capacity(5000);
    let mut ps  = TreeSet::new();
    
    for _ in range(0, 5000) {
        let t = Segment::new(&mut gen);
        for s in ss.iter() {
            for p in t.intersect(s).move_iter() {
                ps.insert(p);
            }
        }
        ss.push(t)
    }
    
    println!("{}", ps.len())
}