//! Space-filling curves and other related functionality.

#![allow(dead_code)]

pub mod hilbert {
    const N: u32 = 1 << 16;

    /// Convert (x,y) to hilbert curve index.
    ///
    /// x: The x coordinate.  Must be no greater than 2^16-1.
    /// y: The y coordinate.  Must be no greater than 2^16-1.
    ///
    /// Returns the hilbert curve index corresponding to the (x,y) coordinates given.
    pub fn encode(x: u32, y: u32) -> u32 {
        assert!(x < N);
        assert!(y < N);

        let (mut x, mut y) = (x, y);
        let mut d = 0;
        let mut s = N >> 1;
        while s > 0 {
            let rx = if (x & s) > 0 { 1 } else { 0 };
            let ry = if (y & s) > 0 { 1 } else { 0 };
            d += s * s * ((3 * rx) ^ ry);
            (x, y) = hilbert_rotate(s, rx, ry, x, y);

            s >>= 1
        }

        d
    }

    /// Convert hilbert curve index to (x,y).
    ///
    /// d: The hilbert curve index.
    ///
    /// Returns the (x, y) coords at the given index.
    pub fn decode(d: u32) -> (u32, u32) {
        let (mut x, mut y) = (0, 0);
        let mut s = 1;
        let mut t = d;
        while s < N {
            let rx = 1 & (t >> 1);
            let ry = 1 & (t ^ rx);
            (x, y) = hilbert_rotate(s, rx, ry, x, y);
            x += s * rx;
            y += s * ry;
            t >>= 2;

            s <<= 1;
        }

        (x, y)
    }

    //------------
    // Utilities.

    fn hilbert_rotate(n: u32, rx: u32, ry: u32, x: u32, y: u32) -> (u32, u32) {
        if ry == 0 {
            if rx == 1 {
                ((n - 1).wrapping_sub(y), (n - 1).wrapping_sub(x))
            } else {
                (y, x)
            }
        } else {
            (x, y)
        }
    }
}

pub mod morton {
    const N: u32 = 1 << 16;

    /// Convert (x,y) to morton curve index.
    ///
    /// x: The x coordinate.  Should be no greater than 2^16-1.
    /// y: The y coordinate.  Should be no greater than 2^16-1.
    ///
    /// Returns the morton curve index corresponding to the (x,y) coordinates given.
    pub fn encode(x: u32, y: u32) -> u32 {
        debug_assert!(x < N);
        debug_assert!(y < N);

        part_1_by_1(x) | (part_1_by_1(y) << 1)
    }

    /// Convert morton curve index to (x,y).
    ///
    /// i: The morton curve index.
    ///
    /// Returns the (x, y) coords at the given index.
    pub fn decode(i: u32) -> (u32, u32) {
        (compact_1_by_1(i), compact_1_by_1(i >> 1))
    }

    //------------
    // Utilities.

    #[inline(always)]
    fn part_1_by_1(mut x: u32) -> u32 {
        x &= 0x0000ffff;
        x = (x ^ (x << 8)) & 0x00ff00ff;
        x = (x ^ (x << 4)) & 0x0f0f0f0f;
        x = (x ^ (x << 2)) & 0x33333333;
        x = (x ^ (x << 1)) & 0x55555555;
        x
    }

    #[inline(always)]
    fn compact_1_by_1(mut x: u32) -> u32 {
        x &= 0x55555555;
        x = (x ^ (x >> 1)) & 0x33333333;
        x = (x ^ (x >> 2)) & 0x0f0f0f0f;
        x = (x ^ (x >> 4)) & 0x00ff00ff;
        x = (x ^ (x >> 8)) & 0x0000ffff;
        x
    }
}

//-------------------------------------------------------------

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn hilbert_reversible() {
        let i = 0x4c8587a2;
        let (x, y) = hilbert::decode(i);
        let i2 = hilbert::encode(x, y);

        assert_eq!(i, i2);
    }

    #[test]
    fn morton_reversible() {
        let i = 0x4c8587a2;
        let (x, y) = morton::decode(i);
        let i2 = morton::encode(x, y);

        assert_eq!(i, i2);
    }
}
