package org.davidmoten.hilbert;

import java.math.BigInteger;
import java.util.Arrays;

import com.github.davidmoten.guavamini.Preconditions;
import com.github.davidmoten.guavamini.annotations.VisibleForTesting;

/**
 * Converts between Hilbert index ({@code BigInteger}) and N-dimensional points.
 * 
 * <p>
 * Note: This algorithm is derived from work done by John Skilling and published
 * in "Programming the Hilbert curve". (c) 2004 American Institute of Physics.
 * With thanks also to Paul Chernoch who published a C# algorithm for Skilling's
 * work on StackOverflow and
 * <a href="https://github.com/paulchernoch/HilbertTransformation">GitHub</a>).
 */
public final class HilbertCurve {

    private final int bits;
    private final int dimensions;
    // cached calculations
    private final int length;

    private HilbertCurve(int bits, int dimensions) {
        this.bits = bits;
        this.dimensions = dimensions;
        // cache a calculated values for small perf improvements
        this.length = bits * dimensions;
    }

    /**
     * Returns a builder for and object that performs transformations for a Hilbert
     * curve with the given number of bits.
     * 
     * @param bits
     *            depth of the Hilbert curve. If bits is one, this is the top-level
     *            Hilbert curve
     * @return builder for object to do transformations with the Hilbert Curve
     */
    public static Builder bits(int bits) {
        return new Builder(bits);
    }

    public static SmallHilbertCurve.Builder small() {
        return new SmallHilbertCurve.Builder();
    }

    /**
     * Builds a {@link HilbertCurve} instance.
     */
    public static final class Builder {
        final int bits;

        private Builder(int bits) {
            Preconditions.checkArgument(bits > 0, "bits must be greater than zero");
            Preconditions.checkArgument(bits < 64, "bits must be 63 or less");
            this.bits = bits;
        }

        public HilbertCurve dimensions(int dimensions) {
            Preconditions.checkArgument(dimensions > 1, "dimensions must be at least 2");
            return new HilbertCurve(bits, dimensions);
        }
    }

    /**
     * Converts a point to its Hilbert curve index.
     * 
     * @param point
     *            an array of {@code long}. Each ordinate can be between 0 and
     *            2<sup>bits</sup>-1.
     * @return index (nonnegative {@link BigInteger})
     * @throws IllegalArgumentException
     *             if length of point array is not equal to the number of
     *             dimensions.
     */
    public BigInteger index(long... point) {
        Preconditions.checkArgument(point.length == dimensions);
        return toIndex(transposedIndex(bits, point));
    }

    /**
     * Converts a {@link BigInteger} index (distance along the Hilbert Curve from 0)
     * to a point of dimensions defined in the constructor of {@code this}.
     * 
     * @param index
     *            index along the Hilbert Curve from 0. Maximum value 2 <sup>bits *
     *            dimensions</sup>-1.
     * @return array of longs being the point
     * @throws NullPointerException
     *             if index is null
     * @throws IllegalArgumentException
     *             if index is negative
     */
    public long[] point(BigInteger index) {
        Preconditions.checkNotNull(index);
        Preconditions.checkArgument(index.signum() != -1, "index cannot be negative");
        return transposedIndexToPoint(bits, transpose(index));
    }

    public void point(BigInteger index, long[] x) {
        Preconditions.checkNotNull(index);
        Preconditions.checkArgument(index.signum() != -1, "index cannot be negative");
        Util.zero(x);
        transpose(index, x);
        transposedIndexToPoint(bits, x);
    }

    public void point(long i, long[] x) {
        point(BigInteger.valueOf(i), x);
    }

    /**
     * Converts a {@code long} index (distance along the Hilbert Curve from 0) to a
     * point of dimensions defined in the constructor of {@code this}.
     * 
     * @param index
     *            index along the Hilbert Curve from 0. Maximum value 2
     *            <sup>bits+1</sup>-1.
     * @return array of longs being the point
     * @throws IllegalArgumentException
     *             if index is negative
     */
    public long[] point(long index) {
        return point(BigInteger.valueOf(index));
    }

    /**
     * Returns the transposed representation of the Hilbert curve index.
     * 
     * <p>
     * The Hilbert index is expressed internally as an array of transposed bits.
     * 
     * <pre>
      Example: 5 bits for each of n=3 coordinates.
         15-bit Hilbert integer = A B C D E F G H I J K L M N O is stored
         as its Transpose                        ^
         X[0] = A D G J M                    X[2]|  7
         X[1] = B E H K N        &lt;-------&gt;       | /X[1]
         X[2] = C F I L O                   axes |/
                high low                         0------&gt; X[0]
     * </pre>
     * 
     * @param index
     *            index to be tranposed
     * @return transposed index
     */
    @VisibleForTesting
    long[] transpose(BigInteger index) {
        long[] x = new long[dimensions];
        transpose(index, x);
        return x;
    }

    private void transpose(BigInteger index, long[] x) {
        byte[] b = index.toByteArray();
        for (int idx = 0; idx < 8 * b.length; idx++) {
            if ((b[b.length - 1 - idx / 8] & (1L << (idx % 8))) != 0) {
                int dim = (length - idx - 1) % dimensions;
                int shift = (idx / dimensions) % bits;
                x[dim] |= 1L << shift;
            }
        }
    }

    /**
     * <p>
     * Given the axes (coordinates) of a point in N-Dimensional space, find the
     * distance to that point along the Hilbert curve. That distance will be
     * transposed; broken into pieces and distributed into an array.
     *
     * <p>
     * The number of dimensions is the length of the hilbertAxes array.
     *
     * <p>
     * Note: In Skilling's paper, this function is called AxestoTranspose.
     * 
     * @param mutate
     * 
     * @param point
     *            Point in N-space
     * @return The Hilbert distance (or index) as a transposed Hilbert index
     */
    @VisibleForTesting
    static long[] transposedIndex(int bits, long... point) {
        final long M = 1L << (bits - 1);
        final int n = point.length; // n: Number of dimensions
        final long[] x = Arrays.copyOf(point, n);
        long p, q, t;
        int i;
        // Inverse undo
        for (q = M; q > 1; q >>= 1) {
            p = q - 1;
            for (i = 0; i < n; i++)
                if ((x[i] & q) != 0)
                    x[0] ^= p; // invert
                else {
                    t = (x[0] ^ x[i]) & p;
                    x[0] ^= t;
                    x[i] ^= t;
                }
        } // exchange
          // Gray encode
        for (i = 1; i < n; i++)
            x[i] ^= x[i - 1];
        t = 0;
        for (q = M; q > 1; q >>= 1)
            if ((x[n - 1] & q) != 0)
                t ^= q - 1;
        for (i = 0; i < n; i++)
            x[i] ^= t;

        return x;
    }

    /**
     * Converts the Hilbert transposed index into an N-dimensional point expressed
     * as a vector of {@code long}.
     * 
     * In Skilling's paper this function is named {@code TransposeToAxes}
     * 
     * @param transposedIndex
     *            distance along the Hilbert curve in transposed form
     * @return the coordinates of the point represented by the transposed index on
     *         the Hilbert curve
     */
    static long[] transposedIndexToPoint(int bits, long... x) {
        final long N = 2L << (bits - 1);
        // Note that x is mutated by this method (as a performance improvement
        // to avoid allocation)
        int n = x.length; // number of dimensions
        long p, q, t;
        int i;
        // Gray decode by H ^ (H/2)
        t = x[n - 1] >> 1;
        // Corrected error in Skilling's paper on the following line. The
        // appendix had i >= 0 leading to negative array index.
        for (i = n - 1; i > 0; i--)
            x[i] ^= x[i - 1];
        x[0] ^= t;
        // Undo excess work
        for (q = 2; q != N; q <<= 1) {
            p = q - 1;
            for (i = n - 1; i >= 0; i--)
                if ((x[i] & q) != 0L)
                    x[0] ^= p; // invert
                else {
                    t = (x[0] ^ x[i]) & p;
                    x[0] ^= t;
                    x[i] ^= t;
                }
        } // exchange
        return x;
    }

    // Quote from Paul Chernoch
    // Interleaving means take one bit from the first matrix element, one bit
    // from the next, etc, then take the second bit from the first matrix
    // element, second bit from the second, all the way to the last bit of the
    // last element. Combine those bits in that order into a single BigInteger,
    // which can have as many bits as necessary. This converts the array into a
    // single number.
    @VisibleForTesting
    BigInteger toIndex(long... transposedIndex) {
        byte[] b = new byte[length / 8 + 1];
        int bIndex = length - 1;
        long mask = 1L << (bits - 1);
        for (int i = 0; i < bits; i++) {
            for (int j = 0; j < transposedIndex.length; j++) {
                if ((transposedIndex[j] & mask) != 0) {
                    b[b.length - 1 - bIndex / 8] |= 1 << (bIndex % 8);
                }
                bIndex--;
            }
            mask >>= 1;
        }
        // b is expected to be BigEndian
        return new BigInteger(1, b);
    }

}