using System;

namespace SquishyTerry.Audio;

public sealed partial class FaceSynth
{
	// musl (fdlibm) single precision powf, the implementation the original binary was built with.
	// Porting it keeps the patch's [pow] / mtof results bit identical to the original.

	static readonly float[] PowBp = { 1.0f, 1.5f };
	static readonly float[] PowDpH = { 0.0f, 5.84960938e-01f };
	static readonly float[] PowDpL = { 0.0f, 1.56322085e-06f };

	static float Bits( int i ) => BitConverter.Int32BitsToSingle( i );
	static int Word( float f ) => BitConverter.SingleToInt32Bits( f );

	/// <summary>musl powf.</summary>
	static float PowF( float x, float y )
	{
		const float two24 = 16777216.0f, huge = 1.0e30f, tiny = 1.0e-30f;
		const float L1 = 6.0000002384e-01f, L2 = 4.2857143283e-01f, L3 = 3.3333334327e-01f, L4 = 2.7272811532e-01f,
			L5 = 2.3066075146e-01f, L6 = 2.0697501302e-01f;
		const float P1 = 1.6666667163e-01f, P2 = -2.7777778450e-03f, P3 = 6.6137559770e-05f, P4 = -1.6533901999e-06f,
			P5 = 4.1381369442e-08f;
		const float lg2 = 6.9314718246e-01f, lg2_h = 6.93145752e-01f, lg2_l = 1.42860654e-06f, ovt = 4.2995665694e-08f;
		const float cp = 9.6179670095e-01f, cp_h = 9.6191406250e-01f, cp_l = -1.1736857402e-04f;
		const float ivln2 = 1.4426950216e+00f, ivln2_h = 1.4426879883e+00f, ivln2_l = 7.0526075433e-06f;

		float z, ax, z_h, z_l, p_h, p_l, y1, t1, t2, r, s, sn, t, u, v, w;
		int i, j, k, yisint, n, hx, hy, ix, iy, isw;

		hx = Word( x );
		hy = Word( y );
		ix = hx & 0x7fffffff;
		iy = hy & 0x7fffffff;
		if ( hx == 0x3f800000 || iy == 0 ) return 1.0f;
		if ( ix > 0x7f800000 || iy > 0x7f800000 ) return x + y;

		yisint = 0;
		if ( hx < 0 )
		{
			if ( iy >= 0x4b800000 ) yisint = 2;
			else if ( iy >= 0x3f800000 )
			{
				k = (iy >> 23) - 0x7f;
				j = iy >> (23 - k);
				if ( (j << (23 - k)) == iy ) yisint = 2 - (j & 1);
			}
		}

		if ( iy == 0x7f800000 )
		{
			if ( ix == 0x3f800000 ) return 1.0f;
			if ( ix > 0x3f800000 ) return hy >= 0 ? y : 0.0f;
			return hy >= 0 ? 0.0f : -y;
		}
		if ( iy == 0x3f800000 ) return hy >= 0 ? x : 1.0f / x;
		if ( hy == 0x40000000 ) return x * x;
		if ( hy == 0x3f000000 && hx >= 0 ) return MathF.Sqrt( x );

		ax = MathF.Abs( x );
		if ( ix == 0x7f800000 || ix == 0 || ix == 0x3f800000 )
		{
			z = ax;
			if ( hy < 0 ) z = 1.0f / z;
			if ( hx < 0 )
			{
				if ( ((ix - 0x3f800000) | yisint) == 0 ) z = (z - z) / (z - z);
				else if ( yisint == 1 ) z = -z;
			}
			return z;
		}

		sn = 1.0f;
		if ( hx < 0 )
		{
			if ( yisint == 0 ) return (x - x) / (x - x);
			if ( yisint == 1 ) sn = -1.0f;
		}

		if ( iy > 0x4d000000 )
		{
			if ( ix < 0x3f7ffff8 ) return hy < 0 ? sn * huge * huge : sn * tiny * tiny;
			if ( ix > 0x3f800007 ) return hy > 0 ? sn * huge * huge : sn * tiny * tiny;
			t = ax - 1;
			w = (t * t) * (0.5f - t * (0.333333333333f - t * 0.25f));
			u = ivln2_h * t;
			v = t * ivln2_l - w * ivln2;
			t1 = u + v;
			t1 = Bits( Word( t1 ) & unchecked((int)0xfffff000) );
			t2 = v - (t1 - u);
		}
		else
		{
			float s2, s_h, s_l, t_h, t_l;
			n = 0;
			if ( ix < 0x00800000 )
			{
				ax *= two24;
				n -= 24;
				ix = Word( ax );
			}
			n += (ix >> 23) - 0x7f;
			j = ix & 0x007fffff;
			ix = j | 0x3f800000;
			if ( j <= 0x1cc471 ) k = 0;
			else if ( j < 0x5db3d7 ) k = 1;
			else
			{
				k = 0;
				n += 1;
				ix -= 0x00800000;
			}
			ax = Bits( ix );

			u = ax - PowBp[k];
			v = 1.0f / (ax + PowBp[k]);
			s = u * v;
			s_h = Bits( Word( s ) & unchecked((int)0xfffff000) );
			isw = ((ix >> 1) & unchecked((int)0xfffff000)) | 0x20000000;
			t_h = Bits( isw + 0x00400000 + (k << 21) );
			t_l = ax - (t_h - PowBp[k]);
			s_l = v * ((u - s_h * t_h) - s_h * t_l);
			s2 = s * s;
			r = s2 * s2 * (L1 + s2 * (L2 + s2 * (L3 + s2 * (L4 + s2 * (L5 + s2 * L6)))));
			r += s_l * (s_h + s);
			s2 = s_h * s_h;
			t_h = 3.0f + s2 + r;
			t_h = Bits( Word( t_h ) & unchecked((int)0xfffff000) );
			t_l = r - ((t_h - 3.0f) - s2);
			u = s_h * t_h;
			v = s_l * t_h + t_l * s;
			p_h = u + v;
			p_h = Bits( Word( p_h ) & unchecked((int)0xfffff000) );
			p_l = v - (p_h - u);
			z_h = cp_h * p_h;
			z_l = cp_l * p_h + p_l * cp + PowDpL[k];
			t = n;
			t1 = ((z_h + z_l) + PowDpH[k]) + t;
			t1 = Bits( Word( t1 ) & unchecked((int)0xfffff000) );
			t2 = z_l - (((t1 - t) - PowDpH[k]) - z_h);
		}

		y1 = Bits( Word( y ) & unchecked((int)0xfffff000) );
		p_l = (y - y1) * t1 + y * t2;
		p_h = y1 * t1;
		z = p_l + p_h;
		j = Word( z );
		if ( j > 0x43000000 ) return sn * huge * huge;
		else if ( j == 0x43000000 )
		{
			if ( p_l + ovt > z - p_h ) return sn * huge * huge;
		}
		else if ( (j & 0x7fffffff) > 0x43160000 ) return sn * tiny * tiny;
		else if ( j == unchecked((int)0xc3160000) )
		{
			if ( p_l <= z - p_h ) return sn * tiny * tiny;
		}

		i = j & 0x7fffffff;
		k = (i >> 23) - 0x7f;
		n = 0;
		if ( i > 0x3f000000 )
		{
			n = j + (0x00800000 >> (k + 1));
			k = ((n & 0x7fffffff) >> 23) - 0x7f;
			t = Bits( n & ~(0x007fffff >> k) );
			n = ((n & 0x007fffff) | 0x00800000) >> (23 - k);
			if ( j < 0 ) n = -n;
			p_h -= t;
		}
		t = p_l + p_h;
		t = Bits( Word( t ) & unchecked((int)0xffff8000) );
		u = t * lg2_h;
		v = (p_l - (t - p_h)) * lg2 + t * lg2_l;
		z = u + v;
		w = v - (z - u);
		t = z * z;
		t1 = z - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
		r = (z * t1) / (t1 - 2.0f) - (w + z * w);
		z = 1.0f - (r - z);
		j = Word( z );
		j += n << 23;
		if ( (j >> 23) <= 0 ) z = ScalbnF( z, n );
		else z = Bits( j );
		return sn * z;
	}

	static float ScalbnF( float x, int n )
	{
		float y = x;
		if ( n > 127 )
		{
			y *= 1.7014118e38f;
			n -= 127;
			if ( n > 127 )
			{
				y *= 1.7014118e38f;
				n -= 127;
				if ( n > 127 ) n = 127;
			}
		}
		else if ( n < -126 )
		{
			y *= 1.1754944e-38f;
			n += 126;
			if ( n < -126 )
			{
				y *= 1.1754944e-38f;
				n += 126;
				if ( n < -126 ) n = -126;
			}
		}
		return y * Bits( (0x7f + n) << 23 );
	}
}