#ifndef MATH_FUNCTIONS_H
#define MATH_FUNCTIONS_H

#include "SPlisHSPlasH/Common.h"

// ------------------------------------------------------------------------------------
namespace SPH
{
	class MathFunctions
	{
	public:
		/** Implementation of the paper: \n
		 * Matthias Müller, Jan Bender, Nuttapong Chentanez and Miles Macklin,
		 * "A Robust Method to Extract the Rotational Part of Deformations", 
		 * ACM SIGGRAPH Motion in Games, 2016
		 */
		static void extractRotation(const Matrix3r &A, Quaternionr &q, const unsigned int maxIter);

		static void pseudoInverse(const Matrix3r &a, Matrix3r &res);
		static void svdWithInversionHandling(const Matrix3r &A, Vector3r &sigma, Matrix3r &U, Matrix3r &VT);
		static void eigenDecomposition(const Matrix3r &A, Matrix3r &eigenVecs, Vector3r &eigenVals);
		static void jacobiRotate(Matrix3r &A, Matrix3r &R, int p, int q);

		/** Returns two orthogonal vectors to vec which are also orthogonal to each other.
		*/
		static void getOrthogonalVectors(const Vector3r &vec, Vector3r &x, Vector3r &y);

		/** computes the APD of 8 deformation gradients. (Alg. 3 from the paper: Kugelstadt et al. "Fast Corotated FEM using Operator Splitting", CGF 2018)
		*/
		static void APD_Newton(const Matrix3r& F, Quaternionr& q);

		/** Computes rotation matrix R from deformation gradient F using iARAP approach.
		 * (Lin et al. "Isotropic ARAP energy using Cauchy-Green invariants", 2022)
		 * Returns R directly as a 3x3 matrix computed from Cauchy-Green invariants.
		 */
		static void iARAP(const Matrix3r& F, Matrix3r& R);

		/** Computes singular values from deformation gradient F using iARAP quartic.
		 * Extracts singular values directly from quartic roots.
		 * sigma is returned in ascending order (sigma[0] = smallest).
		 */
		static void ARAP_eigenvalues(const Matrix3r& F, Vector3r& sigma);

		/** Computes full iARAP decomposition: R, sigma, and V.
		 * Uses quartic for sigma, iARAP for R, and analytical eigenvectors for V.
		 * sigma in ascending order. U = R * V.
		 */
		static void ARAP_decomposition(const Matrix3r& F, Matrix3r& R, Vector3r& sigma, Matrix3r& V);
	};
}

#endif