/*
  Copyright (C) 2013 Andrew Magill

  This software is provided 'as-is', without any express or implied
  warranty.  In no event will the authors be held liable for any damages
  arising from the use of this software.

  Permission is granted to anyone to use this software for any purpose,
  including commercial applications, and to alter it and redistribute it
  freely, subject to the following restrictions:

  1. The origin of this software must not be misrepresented; you must not
     claim that you wrote the original software. If you use this software
     in a product, an acknowledgment in the product documentation would be
     appreciated but is not required.
  2. Altered source versions must be plainly marked as such, and must not be
     misrepresented as being the original software.
  3. This notice may not be removed or altered from any source distribution.

*/

part of vector_math_64;

/*
 * This is based on the implementation of Simplex Noise by Stefan Gustavson
 * found at: http://webstaff.itn.liu.se/~stegu/simplexnoise/SimplexNoise.java
 */

class SimplexNoise {
  static final List<List<double>> _grad3 = <List<double>>[
    <double>[1.0, 1.0, 0.0],
    <double>[-1.0, 1.0, 0.0],
    <double>[1.0, -1.0, 0.0],
    <double>[-1.0, -1.0, 0.0],
    <double>[1.0, 0.0, 1.0],
    <double>[-1.0, 0.0, 1.0],
    <double>[1.0, 0.0, -1.0],
    <double>[-1.0, 0.0, -1.0],
    <double>[0.0, 1.0, 1.0],
    <double>[0.0, -1.0, 1.0],
    <double>[0.0, 1.0, -1.0],
    <double>[0.0, -1.0, -1.0]
  ];

  static final List<List<double>> _grad4 = <List<double>>[
    <double>[0.0, 1.0, 1.0, 1.0],
    <double>[0.0, 1.0, 1.0, -1.0],
    <double>[0.0, 1.0, -1.0, 1.0],
    <double>[0.0, 1.0, -1.0, -1.0],
    <double>[0.0, -1.0, 1.0, 1.0],
    <double>[0.0, -1.0, 1.0, -1.0],
    <double>[0.0, -1.0, -1.0, 1.0],
    <double>[0.0, -1.0, -1.0, -1.0],
    <double>[1.0, 0.0, 1.0, 1.0],
    <double>[1.0, 0.0, 1.0, -1.0],
    <double>[1.0, 0.0, -1.0, 1.0],
    <double>[1.0, 0.0, -1.0, -1.0],
    <double>[-1.0, 0.0, 1.0, 1.0],
    <double>[-1.0, 0.0, 1.0, -1.0],
    <double>[-1.0, 0.0, -1.0, 1.0],
    <double>[-1.0, 0.0, -1.0, -1.0],
    <double>[1.0, 1.0, 0.0, 1.0],
    <double>[1.0, 1.0, 0.0, -1.0],
    <double>[1.0, -1.0, 0.0, 1.0],
    <double>[1.0, -1.0, 0.0, -1.0],
    <double>[-1.0, 1.0, 0.0, 1.0],
    <double>[-1.0, 1.0, 0.0, -1.0],
    <double>[-1.0, -1.0, 0.0, 1.0],
    <double>[-1.0, -1.0, 0.0, -1.0],
    <double>[1.0, 1.0, 1.0, 0.0],
    <double>[1.0, 1.0, -1.0, 0.0],
    <double>[1.0, -1.0, 1.0, 0.0],
    <double>[1.0, -1.0, -1.0, 0.0],
    <double>[-1.0, 1.0, 1.0, 0.0],
    <double>[-1.0, 1.0, -1.0, 0.0],
    <double>[-1.0, -1.0, 1.0, 0.0],
    <double>[-1.0, -1.0, -1.0, 0.0]
  ];

  // To remove the need for index wrapping, double the permutation table length
  List<int> _perm;
  List<int> _permMod12;

  // Skewing and unskewing factors for 2, 3, and 4 dimensions
  static final double _F2 = 0.5 * (math.sqrt(3.0) - 1.0);
  static final double _G2 = (3.0 - math.sqrt(3.0)) / 6.0;
  static const double _F3 = 1.0 / 3.0;
  static const double _G3 = 1.0 / 6.0;
  static final double _F4 = (math.sqrt(5.0) - 1.0) / 4.0;
  static final double _G4 = (5.0 - math.sqrt(5.0)) / 20.0;

  double _dot2(List<double> g, double x, double y) => g[0] * x + g[1] * y;

  double _dot3(List<double> g, double x, double y, double z) =>
      g[0] * x + g[1] * y + g[2] * z;

  double _dot4(List<double> g, double x, double y, double z, double w) =>
      g[0] * x + g[1] * y + g[2] * z + g[3] * w;

  SimplexNoise([math.Random r]) {
    r ??= math.Random();
    final List<int> p =
        List<int>.generate(256, (_) => r.nextInt(256), growable: false);
    _perm = List<int>.generate(p.length * 2, (int i) => p[i % p.length],
        growable: false);
    _permMod12 = List<int>.generate(_perm.length, (int i) => _perm[i] % 12,
        growable: false);
  }

  double noise2D(double xin, double yin) {
    double n0, n1, n2; // Noise contributions from the three corners
    // Skew the input space to determine which simplex cell we're in
    final double s = (xin + yin) * _F2; // Hairy factor for 2D
    final int i = (xin + s).floor();
    final int j = (yin + s).floor();
    final double t = (i + j) * _G2;
    final double X0 = i - t; // Unskew the cell origin back to (x,y) space
    final double Y0 = j - t;
    final double x0 = xin - X0; // The x,y distances from the cell origin
    final double y0 = yin - Y0;
    // For the 2D case, the simplex shape is an equilateral triangle.
    // Determine which simplex we are in.
    int i1, j1; // Offsets for second (middle) corner of simplex in (i,j) coords
    if (x0 > y0) {
      i1 = 1;
      j1 = 0;
    } // lower triangle, XY order: (0,0)->(1,0)->(1,1)
    else {
      i1 = 0;
      j1 = 1;
    } // upper triangle, YX order: (0,0)->(0,1)->(1,1)
    // A step of (1,0) in (i,j) means a step of (1-c,-c) in (x,y), and
    // a step of (0,1) in (i,j) means a step of (-c,1-c) in (x,y), where
    // c = (3-sqrt(3))/6
    final double x1 =
        x0 - i1 + _G2; // Offsets for middle corner in (x,y) unskewed coords
    final double y1 = y0 - j1 + _G2;
    final double x2 = x0 -
        1.0 +
        2.0 * _G2; // Offsets for last corner in (x,y) unskewed coords
    final double y2 = y0 - 1.0 + 2.0 * _G2;
    // Work out the hashed gradient indices of the three simplex corners
    final int ii = i & 255;
    final int jj = j & 255;
    final int gi0 = _permMod12[ii + _perm[jj]];
    final int gi1 = _permMod12[ii + i1 + _perm[jj + j1]];
    final int gi2 = _permMod12[ii + 1 + _perm[jj + 1]];
    // Calculate the contribution from the three corners
    double t0 = 0.5 - x0 * x0 - y0 * y0;
    if (t0 < 0)
      n0 = 0.0;
    else {
      t0 *= t0;
      n0 = t0 *
          t0 *
          _dot2(_grad3[gi0], x0, y0); // (x,y) of grad3 used for 2D gradient
    }
    double t1 = 0.5 - x1 * x1 - y1 * y1;
    if (t1 < 0)
      n1 = 0.0;
    else {
      t1 *= t1;
      n1 = t1 * t1 * _dot2(_grad3[gi1], x1, y1);
    }
    double t2 = 0.5 - x2 * x2 - y2 * y2;
    if (t2 < 0)
      n2 = 0.0;
    else {
      t2 *= t2;
      n2 = t2 * t2 * _dot2(_grad3[gi2], x2, y2);
    }
    // Add contributions from each corner to get the final noise value.
    // The result is scaled to return values in the interval [-1,1].
    return 70.0 * (n0 + n1 + n2);
  }

  // 3D simplex noise
  double noise3D(double xin, double yin, double zin) {
    double n0, n1, n2, n3; // Noise contributions from the four corners
    // Skew the input space to determine which simplex cell we're in
    final double s =
        (xin + yin + zin) * _F3; // Very nice and simple skew factor for 3D
    final int i = (xin + s).floor();
    final int j = (yin + s).floor();
    final int k = (zin + s).floor();
    final double t = (i + j + k) * _G3;
    final double X0 = i - t; // Unskew the cell origin back to (x,y,z) space
    final double Y0 = j - t;
    final double Z0 = k - t;
    final double x0 = xin - X0; // The x,y,z distances from the cell origin
    final double y0 = yin - Y0;
    final double z0 = zin - Z0;
    // For the 3D case, the simplex shape is a slightly irregular tetrahedron.
    // Determine which simplex we are in.
    int i1, j1, k1; // Offsets for second corner of simplex in (i,j,k) coords
    int i2, j2, k2; // Offsets for third corner of simplex in (i,j,k) coords
    if (x0 >= y0) {
      if (y0 >= z0) {
        i1 = 1;
        j1 = 0;
        k1 = 0;
        i2 = 1;
        j2 = 1;
        k2 = 0;
      } // X Y Z order
      else if (x0 >= z0) {
        i1 = 1;
        j1 = 0;
        k1 = 0;
        i2 = 1;
        j2 = 0;
        k2 = 1;
      } // X Z Y order
      else {
        i1 = 0;
        j1 = 0;
        k1 = 1;
        i2 = 1;
        j2 = 0;
        k2 = 1;
      } // Z X Y order
    } else {
      // x0<y0
      if (y0 < z0) {
        i1 = 0;
        j1 = 0;
        k1 = 1;
        i2 = 0;
        j2 = 1;
        k2 = 1;
      } // Z Y X order
      else if (x0 < z0) {
        i1 = 0;
        j1 = 1;
        k1 = 0;
        i2 = 0;
        j2 = 1;
        k2 = 1;
      } // Y Z X order
      else {
        i1 = 0;
        j1 = 1;
        k1 = 0;
        i2 = 1;
        j2 = 1;
        k2 = 0;
      } // Y X Z order
    }
    // A step of (1,0,0) in (i,j,k) means a step of (1-c,-c,-c) in (x,y,z),
    // a step of (0,1,0) in (i,j,k) means a step of (-c,1-c,-c) in (x,y,z), and
    // a step of (0,0,1) in (i,j,k) means a step of (-c,-c,1-c) in (x,y,z), where
    // c = 1/6.
    final double x1 =
        x0 - i1 + _G3; // Offsets for second corner in (x,y,z) coords
    final double y1 = y0 - j1 + _G3;
    final double z1 = z0 - k1 + _G3;
    final double x2 =
        x0 - i2 + 2.0 * _G3; // Offsets for third corner in (x,y,z) coords
    final double y2 = y0 - j2 + 2.0 * _G3;
    final double z2 = z0 - k2 + 2.0 * _G3;
    final double x3 =
        x0 - 1.0 + 3.0 * _G3; // Offsets for last corner in (x,y,z) coords
    final double y3 = y0 - 1.0 + 3.0 * _G3;
    final double z3 = z0 - 1.0 + 3.0 * _G3;
    // Work out the hashed gradient indices of the four simplex corners
    final int ii = i & 255;
    final int jj = j & 255;
    final int kk = k & 255;
    final int gi0 = _permMod12[ii + _perm[jj + _perm[kk]]];
    final int gi1 = _permMod12[ii + i1 + _perm[jj + j1 + _perm[kk + k1]]];
    final int gi2 = _permMod12[ii + i2 + _perm[jj + j2 + _perm[kk + k2]]];
    final int gi3 = _permMod12[ii + 1 + _perm[jj + 1 + _perm[kk + 1]]];
    // Calculate the contribution from the four corners
    double t0 = 0.6 - x0 * x0 - y0 * y0 - z0 * z0;
    if (t0 < 0)
      n0 = 0.0;
    else {
      t0 *= t0;
      n0 = t0 * t0 * _dot3(_grad3[gi0], x0, y0, z0);
    }
    double t1 = 0.6 - x1 * x1 - y1 * y1 - z1 * z1;
    if (t1 < 0)
      n1 = 0.0;
    else {
      t1 *= t1;
      n1 = t1 * t1 * _dot3(_grad3[gi1], x1, y1, z1);
    }
    double t2 = 0.6 - x2 * x2 - y2 * y2 - z2 * z2;
    if (t2 < 0)
      n2 = 0.0;
    else {
      t2 *= t2;
      n2 = t2 * t2 * _dot3(_grad3[gi2], x2, y2, z2);
    }
    double t3 = 0.6 - x3 * x3 - y3 * y3 - z3 * z3;
    if (t3 < 0)
      n3 = 0.0;
    else {
      t3 *= t3;
      n3 = t3 * t3 * _dot3(_grad3[gi3], x3, y3, z3);
    }
    // Add contributions from each corner to get the final noise value.
    // The result is scaled to stay just inside [-1,1]
    return 32.0 * (n0 + n1 + n2 + n3);
  }

  // 4D simplex noise, better simplex rank ordering method 2012-03-09
  double noise4D(double x, double y, double z, double w) {
    double n0, n1, n2, n3, n4; // Noise contributions from the five corners
    // Skew the (x,y,z,w) space to determine which cell of 24 simplices we're in
    final double s = (x + y + z + w) * _F4; // Factor for 4D skewing
    final int i = (x + s).floor();
    final int j = (y + s).floor();
    final int k = (z + s).floor();
    final int l = (w + s).floor();
    final double t = (i + j + k + l) * _G4; // Factor for 4D unskewing
    final double X0 = i - t; // Unskew the cell origin back to (x,y,z,w) space
    final double Y0 = j - t;
    final double Z0 = k - t;
    final double W0 = l - t;
    final double x0 = x - X0; // The x,y,z,w distances from the cell origin
    final double y0 = y - Y0;
    final double z0 = z - Z0;
    final double w0 = w - W0;
    // For the 4D case, the simplex is a 4D shape I won't even try to describe.
    // To find out which of the 24 possible simplices we're in, we need to
    // determine the magnitude ordering of x0, y0, z0 and w0.
    // Six pair-wise comparisons are performed between each possible pair
    // of the four coordinates, and the results are used to rank the numbers.
    int rankx = 0;
    int ranky = 0;
    int rankz = 0;
    int rankw = 0;
    if (x0 > y0)
      rankx++;
    else
      ranky++;
    if (x0 > z0)
      rankx++;
    else
      rankz++;
    if (x0 > w0)
      rankx++;
    else
      rankw++;
    if (y0 > z0)
      ranky++;
    else
      rankz++;
    if (y0 > w0)
      ranky++;
    else
      rankw++;
    if (z0 > w0)
      rankz++;
    else
      rankw++;
    int i1, j1, k1, l1; // The integer offsets for the second simplex corner
    int i2, j2, k2, l2; // The integer offsets for the third simplex corner
    int i3, j3, k3, l3; // The integer offsets for the fourth simplex corner
    // simplex[c] is a 4-vector with the numbers 0, 1, 2 and 3 in some order.
    // Many values of c will never occur, since e.g. x>y>z>w makes x<z, y<w and x<w
    // impossible. Only the 24 indices which have non-zero entries make any sense.
    // We use a thresholding to set the coordinates in turn from the largest magnitude.
    // Rank 3 denotes the largest coordinate.
    i1 = rankx >= 3 ? 1 : 0;
    j1 = ranky >= 3 ? 1 : 0;
    k1 = rankz >= 3 ? 1 : 0;
    l1 = rankw >= 3 ? 1 : 0;
    // Rank 2 denotes the second largest coordinate.
    i2 = rankx >= 2 ? 1 : 0;
    j2 = ranky >= 2 ? 1 : 0;
    k2 = rankz >= 2 ? 1 : 0;
    l2 = rankw >= 2 ? 1 : 0;
    // Rank 1 denotes the second smallest coordinate.
    i3 = rankx >= 1 ? 1 : 0;
    j3 = ranky >= 1 ? 1 : 0;
    k3 = rankz >= 1 ? 1 : 0;
    l3 = rankw >= 1 ? 1 : 0;
    // The fifth corner has all coordinate offsets = 1, so no need to compute that.
    final double x1 =
        x0 - i1 + _G4; // Offsets for second corner in (x,y,z,w) coords
    final double y1 = y0 - j1 + _G4;
    final double z1 = z0 - k1 + _G4;
    final double w1 = w0 - l1 + _G4;
    final double x2 =
        x0 - i2 + 2.0 * _G4; // Offsets for third corner in (x,y,z,w) coords
    final double y2 = y0 - j2 + 2.0 * _G4;
    final double z2 = z0 - k2 + 2.0 * _G4;
    final double w2 = w0 - l2 + 2.0 * _G4;
    final double x3 =
        x0 - i3 + 3.0 * _G4; // Offsets for fourth corner in (x,y,z,w) coords
    final double y3 = y0 - j3 + 3.0 * _G4;
    final double z3 = z0 - k3 + 3.0 * _G4;
    final double w3 = w0 - l3 + 3.0 * _G4;
    final double x4 =
        x0 - 1.0 + 4.0 * _G4; // Offsets for last corner in (x,y,z,w) coords
    final double y4 = y0 - 1.0 + 4.0 * _G4;
    final double z4 = z0 - 1.0 + 4.0 * _G4;
    final double w4 = w0 - 1.0 + 4.0 * _G4;
    // Work out the hashed gradient indices of the five simplex corners
    final int ii = i & 255;
    final int jj = j & 255;
    final int kk = k & 255;
    final int ll = l & 255;
    final int gi0 = _perm[ii + _perm[jj + _perm[kk + _perm[ll]]]] % 32;
    final int gi1 =
        _perm[ii + i1 + _perm[jj + j1 + _perm[kk + k1 + _perm[ll + l1]]]] % 32;
    final int gi2 =
        _perm[ii + i2 + _perm[jj + j2 + _perm[kk + k2 + _perm[ll + l2]]]] % 32;
    final int gi3 =
        _perm[ii + i3 + _perm[jj + j3 + _perm[kk + k3 + _perm[ll + l3]]]] % 32;
    final int gi4 =
        _perm[ii + 1 + _perm[jj + 1 + _perm[kk + 1 + _perm[ll + 1]]]] % 32;
    // Calculate the contribution from the five corners
    double t0 = 0.6 - x0 * x0 - y0 * y0 - z0 * z0 - w0 * w0;
    if (t0 < 0)
      n0 = 0.0;
    else {
      t0 *= t0;
      n0 = t0 * t0 * _dot4(_grad4[gi0], x0, y0, z0, w0);
    }
    double t1 = 0.6 - x1 * x1 - y1 * y1 - z1 * z1 - w1 * w1;
    if (t1 < 0)
      n1 = 0.0;
    else {
      t1 *= t1;
      n1 = t1 * t1 * _dot4(_grad4[gi1], x1, y1, z1, w1);
    }
    double t2 = 0.6 - x2 * x2 - y2 * y2 - z2 * z2 - w2 * w2;
    if (t2 < 0)
      n2 = 0.0;
    else {
      t2 *= t2;
      n2 = t2 * t2 * _dot4(_grad4[gi2], x2, y2, z2, w2);
    }
    double t3 = 0.6 - x3 * x3 - y3 * y3 - z3 * z3 - w3 * w3;
    if (t3 < 0)
      n3 = 0.0;
    else {
      t3 *= t3;
      n3 = t3 * t3 * _dot4(_grad4[gi3], x3, y3, z3, w3);
    }
    double t4 = 0.6 - x4 * x4 - y4 * y4 - z4 * z4 - w4 * w4;
    if (t4 < 0)
      n4 = 0.0;
    else {
      t4 *= t4;
      n4 = t4 * t4 * _dot4(_grad4[gi4], x4, y4, z4, w4);
    }
    // Sum up and scale the result to cover the range [-1,1]
    return 27.0 * (n0 + n1 + n2 + n3 + n4);
  }
}