#ifndef MATH_H
#define MATH_H


#include <random>
#define restrict  __restrict


namespace Math
{
  float dotProduct(const float *v1, const float *v2, int stride1 = 1, int stride2 = 1, int count = 4);
  void crossProduct(float *restrict out, const float *restrict v1, const float *restrict v2);
  void transformVector(float *restrict vout, const float *restrict m, const float *restrict vin);
  void normalizeVector(float *vInOut, int count = 4);

  // Functions to make a new matrix
  void makeIdentityMatrix4x4(float* out);
  void makePerspectiveMatrix4x4(float* out, float a_fov, float a_aspect, float a_near, float a_far);
  void makeOrthographicMatrix4x4(float* out, float , float , float, float, float a_near, float a_far);

  // Apply transfomations to a matrix (pass in an identity matrix if need to make a matrix)
  void translateMatrix4x4(float *restrict inOut, const float *restrict vec3);
  void scaleMatrix4x4(float *restrict inOut, const float *restrict vec3);
  void rotateMatrix4x4(float *restrict inOut, float radians, float *restrict vec3); // side-effect: vec3 will be normalized

  void translateMatrix4x4(float* inOut, float x, float y, float z);
  void scaleMatrix4x4(float* inOut, float x, float y, float z);
  void rotateMatrix4x4(float* inOut, float rad, float x, float y, float z);

  // Helpers
  void multiplyMatrix4x4(float *restrict out, const float *restrict m1, const float *restrict m2);
  bool matrix4x4ToNormalMatrix3x3(float *restrict out, const float *restrict m);
  bool inverseMatrix4x4(float *restrict out, const float *restrict m);
  void transposeMatrix4x4(float *restrict out, const float *restrict m);

  template <typename T>
  inline T lerp(T t0, T t1, float ratio) {
    return t0 + ratio * (t1 - t0);
  }

  template <typename T>
  inline T clamp(T t, T t1, T t2) {
    T max = t > t1 ? t : t1;
    return max > t2 ? t2 : max;
  }

/*
  template <typename T>
  inline T floor2(T t1, T t2) {
    // 0.5 -> 0.0      -0.5 -> -1.0
    return T(int(((t1 - ((t2+1)>>1))/t2)));
  }
*/

  // different result to x%y, result here should always be positive
  template <typename T>
  T modulus(T x, T y)
  {
    // return x - y * ((x-y+T(1))/y);
    return x - y * int(floor(float(x)/y)); // inefficient at the moment, probably can do much better if test and use something like the floor2 function above
  }

  // different results to fmod, this function matches what glsl shaders mod function does
  // this is important for predicting where a model will be rendered when allowing the vertex shader
  // to wrap positions
  template <typename T>
  float modulus(float x, float y)
  {
    return x - y * floor(x/y);
  }

  class RandomNumberGenerator
  {
    public:
      RandomNumberGenerator() :  gen(rd()) {}
      template <typename T>
      T generate(T min = 0, T max = std::numeric_limits<T>::max()) {
        std::uniform_int_distribution<T> dis(min, max);
        return dis(gen);
      }
    private:
      std::random_device rd;
      std::mt19937 gen;
  };

}


#endif // MATH_H