diff --git a/Changelog.md b/Changelog.md index 1e06229..73565e7 100644 --- a/Changelog.md +++ b/Changelog.md @@ -17,6 +17,12 @@ + ArcSec(Equation) + ArcSin(Equation) + ArcTan(Equation) + + ArcCosh(Equation) + + ArcCoth(Equation) + + ArcCsch(Equation) + + ArcSech(Equation) + + ArcSinh(Equation) + + ArcTanh(Equation) + Average(Equation, float, float, float) + Average(Equation, Equation, Equation, float) + Binomial(Equation, int, float) @@ -28,14 +34,18 @@ + Combinations(Equation, int) + Combinations(Equation, Equation) + Cos(Equation) + + Cosh(Equation) + Cot(Equation) + + Coth(Equation) + Csc(Equation) + + Csch(Equation) + Divide(Equation, float[]) + Divide(Equation, Equation[]) + Factorial(Equation) + Floor(Equation) + GetValues(Equation, float, float, float) + InverseSqrt(Equation) + + Log(Equation, float) + Max(Equation, float, float, float) + Min(Equation, float, float, float) + Permutations(Equation, int) @@ -48,7 +58,9 @@ + Root(Equation, Equation) + Round(Equation) + Sec(Equation) + + Sech(Equation) + Sin(Equation) + + Sinh(Equation) + SolveBisection(Equation, float, float, float, float, int) + SolveEquation(Equation, float, float, float, int) + SolveNewton(Equation, float, float, float, int) @@ -58,6 +70,7 @@ + Sum(Equation, float[]) + Sum(Equation, Equation[]) + Tan(Equation) + + Tanh(Equation) + ZScore(Equation, float[]) + ZScore(Equation, Equation[]) + ZScore(Equation, float, float) @@ -65,8 +78,10 @@ + InvokeMethod(Equation, MethodInfo, object?[]?) + InvokeMathMethod(Equation, string, object?[]?) + Helpers + + CordicHelper + MathfHelper + RationalHelper + + UnsafeHelper * Mathematics * Abstract = Renamed `IPresets1D` to `IPresets1d` @@ -144,18 +159,34 @@ - operator >=(Int4, Int4) - operator <=(Int4, Int4) * Mathf + + ArcCosh(float) + + ArcCoth(float) + + ArcCsch(float) + + ArcSech(float) + + ArcSinh(float) + + ArcTanh(float) + + ArcTanh2(float, float) + Cbrt(float) + + Cosh(float) + + Coth(float) + + Csch(float) + IsPrime(int, PrimeCheckMethod) + Lerp(float, float, Equation, bool) + Lerp(Equation, Equation, float, bool) + Lerp(Equation, Equation, Equation, bool) + + Log(float, float) + PrimeFactors(int) + PowerMod(long, long, long) + + Sech(float) + + Sinh(float) + SharedItems(T[][]) + SolveBisection(Equation, float, float, float, float, int) + SolveEquation(Equation, float, float, float, int) + SolveNewton(Equation, float, float, float, int) - = Improved the `Sqrt` method by using a solution finder. + + Tanh(float) + = Improved the `Sqrt(float)` method by using a solution finder + = The `ArcSin(float)` method now uses a solution finder rather than the base math library + = The `Power(float, float)` method now utilizes a custom CORDIC implementation rather than the base math library + PrimeCheckMethod + Equation2d + Rational @@ -174,4 +205,5 @@ = Renamed `Modifier2D` to `IModifier2d` = Renamed `Modifier2D` to `IModifier2d` = Renamed `Modifier2D` to `IModifier2d` += Made `Nerd_STF` allow unsafe code blocks ``` diff --git a/Nerd_STF/Extensions/EquationExtension.cs b/Nerd_STF/Extensions/EquationExtension.cs index 8a55c24..f62b08d 100644 --- a/Nerd_STF/Extensions/EquationExtension.cs +++ b/Nerd_STF/Extensions/EquationExtension.cs @@ -27,6 +27,13 @@ public static class EquationExtension public static Equation ArcSin(this Equation equ) => x => Mathf.ArcSin(equ(x)).Radians; public static Equation ArcTan(this Equation equ) => x => Mathf.ArcTan(equ(x)).Radians; + public static Equation ArcCosh(this Equation equ) => x => Mathf.ArcCosh(equ(x)); + public static Equation ArcCoth(this Equation equ) => x => Mathf.ArcCoth(equ(x)); + public static Equation ArcCsch(this Equation equ) => x => Mathf.ArcCsch(equ(x)); + public static Equation ArcSech(this Equation equ) => x => Mathf.ArcSech(equ(x)); + public static Equation ArcSinh(this Equation equ) => x => Mathf.ArcSinh(equ(x)); + public static Equation ArcTanh(this Equation equ) => x => Mathf.ArcTanh(equ(x)); + public static float Average(this Equation equ, float min, float max, float step = Calculus.DefaultStep) => Mathf.Average(equ, min, max, step); public static Equation Average(this Equation equ, Equation min, Equation max, float step = Calculus.DefaultStep) => @@ -54,6 +61,10 @@ public static class EquationExtension public static Equation Cot(this Equation equ) => x => Mathf.Cot(equ(x)); public static Equation Csc(this Equation equ) => x => Mathf.Csc(equ(x)); + public static Equation Cosh(this Equation equ) => x => Mathf.Cosh(equ(x)); + public static Equation Coth(this Equation equ) => x => Mathf.Coth(equ(x)); + public static Equation Csch(this Equation equ) => x => Mathf.Csch(equ(x)); + public static Equation Divide(this Equation equ, params float[] dividends) => x => Mathf.Divide(equ(x), dividends); public static Equation Divide(this Equation equ, params Equation[] dividends) => delegate (float x) @@ -72,6 +83,8 @@ public static class EquationExtension public static Equation InverseSqrt(this Equation equ) => x => Mathf.InverseSqrt(equ(x)); + public static Equation Log(this Equation equ, float @base) => x => Mathf.Log(@base, equ(x)); + public static float Max(this Equation equ, float min, float max, float step = Calculus.DefaultStep) => Mathf.Max(equ, min, max, step); public static float Min(this Equation equ, float min, float max, float step = Calculus.DefaultStep) => @@ -108,6 +121,9 @@ public static class EquationExtension public static Equation Sec(this Equation equ) => x => Mathf.Sec(equ(x)); public static Equation Sin(this Equation equ) => x => Mathf.Sin(equ(x)); + public static Equation Sech(this Equation equ) => x => Mathf.Sech(equ(x)); + public static Equation Sinh(this Equation equ) => x => Mathf.Sinh(equ(x)); + public static float SolveBisection(this Equation equ, float initialA, float initialB, float tolerance = 1e-5f, int maxIterations = 1000) => Mathf.SolveBisection(equ, initialA, initialB, tolerance, maxIterations); @@ -146,6 +162,8 @@ public static class EquationExtension public static Equation Tan(this Equation equ) => x => Mathf.Tan(equ(x)); + public static Equation Tanh(this Equation equ) => x => Mathf.Tanh(equ(x)); + public static Equation ZScore(this Equation equ, params float[] vals) => x => Mathf.ZScore(equ(x), vals); public static Equation ZScore(this Equation equ, params Equation[] vals) => delegate (float x) { diff --git a/Nerd_STF/Helpers/CordicHelper.cs b/Nerd_STF/Helpers/CordicHelper.cs new file mode 100644 index 0000000..bcd90ae --- /dev/null +++ b/Nerd_STF/Helpers/CordicHelper.cs @@ -0,0 +1,317 @@ +namespace Nerd_STF.Helpers; + +// TODO: Make this internal + +// Putting this in here for future reference: +// CORDIC is basically just splitting up an +// operation into more smaller operations. +// For example, turning sin(5rad) into +// sin(4rad + 1rad), which can be turned into +// the formula cos(4rad)cos(1rad) - sin(4rad)sin(1rad), +// to which we know the values of each part. +// Then we just do this iteratively on a bunch +// of powers of 2. +public static class CordicHelper +{ + private static readonly float[] p_cosTable = + { + 0.540302305868f, // cos(2^0) + 0.87758256189f, // cos(2^-1) + 0.968912421711f, // cos(2^-2) + 0.992197667229f, // cos(2^-3) + 0.9980475107f, // cos(2^-4) + 0.999511758485f, // cos(2^-5) + 0.999877932171f, // cos(2^-6) + 0.999969482577f, // cos(2^-7) + 0.999992370615f, // cos(2^-8) + 0.999998092652f, // cos(2^-9) + 0.999999523163f, // cos(2^-10) + 0.999999880791f, // cos(2^-11) + 0.999999970198f, // cos(2^-12) + 0.999999992549f, // cos(2^-13) + 0.999999998137f, // cos(2^-14) + 0.999999999534f // cos(2^-15) + }; + private static readonly float[] p_sinTable = + { + 0.841470984808f, // sin(2^0) + 0.479425538604f, // sin(2^-1) + 0.247403959255f, // sin(2^-2) + 0.124674733385f, // sin(2^-3) + 0.0624593178424f, // sin(2^-4) + 0.0312449139853f, // sin(2^-5) + 0.0156243642249f, // sin(2^-6) + 0.00781242052738f, // sin(2^-7) + 0.0039062400659f, // sin(2^-8) + 0.00195312375824f, // sin(2^-9) + 0.00097656234478f, // sin(2^-10) + 0.000488281230597f, // sin(2^-11) + 0.000244140622575f, // sin(2^-12) + 0.000122070312197f, // sin(2^-13) + 0.0000610351562121f, // sin(2^-14) + 0.0000305175781203f // sin(2^-15) + }; + + private static readonly float[] p_coshTable = + { + 1.54308063482f, // cosh(2^0) + 1.12762596521f, // cosh(2^-1) + 1.03141309988f, // cosh(2^-2) + 1.00782267783f, // cosh(2^-3) + 1.00195376087f, // cosh(2^-4) + 1.00048832099f, // cosh(2^-5) + 1.0001220728f, // cosh(2^-6) + 1.00003051773f, // cosh(2^-7) + 1.0000076294f, // cosh(2^-8) + 1.00000190735f, // cosh(2^-9) + 1.00000047684f, // cosh(2^-10) + 1.00000011921f, // cosh(2^-11) + 1.0000000298f, // cosh(2^-12) + 1.00000000745f, // cosh(2^-13) + 1.00000000186f, // cosh(2^-14) + 1.00000000047f, // cosh(2^-15) + }; + private static readonly float[] p_sinhTable = + { + 1.17520119364f, // sinh(2^0) + 0.521095305494f, // sinh(2^-1) + 0.252612316808f, // sinh(2^-2) + 0.125325775241f, // sinh(2^-3) + 0.0625406980522f, // sinh(2^-4) + 0.0312550865114f, // sinh(2^-5) + 0.0156256357906f, // sinh(2^-6) + 0.0078125794731f, // sinh(2^-7) + 0.00390625993412f, // sinh(2^-8) + 0.00195312624176f, // sinh(2^-9) + 0.00097656265522f, // sinh(2^-10) + 0.000488281269403f, // sinh(2^-11) + 0.000244140627425f, // sinh(2^-12) + 0.000122070312803f, // sinh(2^-13) + 0.0000610351562879f, // sinh(2^-14) + 0.0000305175781297f, // sinh(2^-15) + }; + + private static readonly Dictionary<(float bas, int depth), float[]> p_expTables; + + static CordicHelper() + { + p_expTables = new(); + } + + // This was originally intended to replace the Mathf.Cos + // and Mathf.Sin functions, but it ended up being considerably + // slower. In the future if it gets optimized, I might then + // choose to replace it. + // REMEMBER: When implementing, remember to use Mathf.AbsoluteMod, + // because that's what I was intending when I wrote this. + public static (float sin, float cos) CalculateTrig(float x, int iterations) + { + float approximateX = 0, + approximateCos = 1, + approximateSin = 0; + + // Iterate continuously until it gets better. + for (int i = 0; i < iterations; i++) + { + // We need to find the biggest step that'll move us + // closer to the real X (without overshooting). + float diffX = x - approximateX; + + // This is assuming that cosTable and sinTable + // have the same length. + for (int j = 0; j < p_cosTable.Length; j++) + { + // The amount the difference will shrink. + float incX = FastGenExp2((sbyte)-j); + + if (diffX >= incX) + { + // Because here we go big to small, the first one that triggers + // this if statement should also be the biggest one that can. + + // Get the sin and cos values for this power of two. + float valCos = p_cosTable[j], + valSin = p_sinTable[j]; + + // Do the products. + float newCos = approximateCos * valCos - approximateSin * valSin, + newSin = approximateCos * valSin + approximateSin * valCos; + + // Apply differences + approximateX += incX; + approximateCos = newCos; + approximateSin = newSin; + break; + } + } + } + + // Sin and cos should be pretty accurate by now, + // so we can return them. + return (approximateSin, approximateCos); + } + + public static (float sinh, float cosh) CalculateHyperTrig(float x, int iterations) + { + float approximateX = 0, + approximateCosh = 1, + approximateSinh = 0; + + // Iterate continuously until it gets better. + for (int i = 0; i < iterations; i++) + { + // We need to find the biggest step that'll move us + // closer to the real X (without overshooting). + float diffX = x - approximateX; + + // This is assuming that cosTable and sinTable + // have the same length. + for (int j = 0; j < p_coshTable.Length; j++) + { + // The amount the difference will shrink. + float incX = FastGenExp2((sbyte)-j); + + if (diffX >= incX) + { + // Because here we go big to small, the first one that triggers + // this if statement should also be the biggest one that can. + + // Get the sin and cos values for this power of two. + float valCosh = p_coshTable[j], + valSinh = p_sinhTable[j]; + + // Do the products. + float newCosh = approximateCosh * valCosh + approximateSinh * valSinh, + newSinh = approximateCosh * valSinh + approximateSinh * valCosh; + + // Apply differences + approximateX += incX; + approximateCosh = newCosh; + approximateSinh = newSinh; + break; + } + } + } + + // Sin and cos should be pretty accurate by now, + // so we can return them. + return (approximateSinh, approximateCosh); + } + + public static float ExpAnyBase(float bas, float pow, int tableDepth, int iterations) + { + // We need to auto-generate a table of values for this number the user enters. + float[] table; + if (p_expTables.ContainsKey((bas, tableDepth))) + { + // Table was already generated, so we can reuse it. + table = p_expTables[(bas, tableDepth)]; + } + else + { + // Calculate a table for the CORDIC system by + // applying sequential square roots. + table = new float[tableDepth]; + table[0] = bas; + for (int i = 1; i < tableDepth; i++) table[i] = Mathf.Sqrt(table[i - 1]); + p_expTables.Add((bas, tableDepth), table); + } + + // Now we can perform the CORDIC method. + float approximateX = 0, approximateVal = 1; + + // Iterate continuously until it gets better. + for (int i = 0; i < iterations; i++) + { + // We need to find the biggest step that'll move us + // closer to the real X (without overshooting). + float diffX = pow - approximateX; + + for (int j = 0; j < tableDepth; j++) + { + // The amount the difference will shrink. + float incX = FastGenExp2((sbyte)-j); + + if (diffX >= incX) + { + // Because here we go big to small, the first one that triggers + // this if statement should also be the biggest one that can. + + // Get the power value for this power of two. + float val = table[j]; + + // Apply our value. + approximateX += incX; + approximateVal *= val; + break; + } + } + } + + // Value should be pretty accurate by now, + // so we can return it. + return approximateVal; + } + public static float LogAnyBase(float bas, float val, int tableDepth, int iterations) + { + // We need to auto-generate a table of values for this number the user enters. + // However, we can use the already existing exponent tables and just swap the + // indexes and the values. + float[] table; + if (p_expTables.ContainsKey((bas, tableDepth))) + { + // Table was already generated, so we can reuse it. + table = p_expTables[(bas, tableDepth)]; + } + else + { + // Calculate a table for the CORDIC system by + // applying sequential square roots. + table = new float[tableDepth]; + table[0] = bas; + for (int i = 1; i < tableDepth; i++) table[i] = Mathf.Sqrt(table[i - 1]); + p_expTables.Add((bas, tableDepth), table); + } + + // Now we can perform the CORDIC method. + float approximateX = 0, approximateVal = 1; + + // Iterate continuously until it gets better. + for (int i = 0; i < iterations; i++) + { + float diffY = val / approximateVal; + + for (int j = 0; j < table.Length; j++) + { + // The amount the difference will shrink. + float incX = FastGenExp2((sbyte)-j); + float newVal = table[j]; + + if (diffY >= newVal) + { + // Because here we go big to small, the first one that triggers + // this if statement should also be the biggest one that can. + + // Apply our value. + approximateX += incX; + approximateVal *= newVal; + break; + } + } + } + + // Value should be pretty accurate by now, + // so we can return it. + return approximateX; + } + + // An extremely fast way to generate 2 to + // the power of p. I say "generate" because I'm + // just messing with the mantissa's data and + // not doing any real math. + private static float FastGenExp2(sbyte p) + { + int data = (((p - 1) ^ 0b10000000) << 23) & ~(1 << 31); + return UnsafeHelper.SwapType(data); + } +} diff --git a/Nerd_STF/Helpers/UnsafeHelper.cs b/Nerd_STF/Helpers/UnsafeHelper.cs new file mode 100644 index 0000000..30fbc00 --- /dev/null +++ b/Nerd_STF/Helpers/UnsafeHelper.cs @@ -0,0 +1,13 @@ +namespace Nerd_STF.Helpers; + +// These are all the unsafe functions I couldn't make safe. I don't want too much +// unsafe code, so this is where I put all of it that I require. +internal static unsafe class UnsafeHelper +{ + // Forcefully change the type of an object + // without changing the data of the object. + public static NT SwapType(CT obj) + where CT : unmanaged + where NT : unmanaged + => *(NT*)&obj; +} diff --git a/Nerd_STF/Mathematics/Mathf.cs b/Nerd_STF/Mathematics/Mathf.cs index c9ea1d4..113c528 100644 --- a/Nerd_STF/Mathematics/Mathf.cs +++ b/Nerd_STF/Mathematics/Mathf.cs @@ -19,19 +19,28 @@ public static class Mathf } public static Angle ArcCos(float value) => ArcSin(-value) + Angle.Quarter; - public static Angle ArcCot(float value) => ArcCos(value / Sqrt(1 + value * value)); - public static Angle ArcCsc(float value) => ArcSin(1 / value); - public static Angle ArcSec(float value) => ArcCos(1 / value); - - // Maybe one day I'll have a polynomial for this, but the RMSE for an order 10 polynomial is only 0.00876. - public static Angle ArcSin(float value) => new((float)Math.Asin(value), Angle.Type.Radians); - + public static Angle ArcSin(float value) + { + if (value > 1 || value < -1) throw new ArgumentOutOfRangeException(nameof(value)); + return (SolveNewton(x => Sin(x) - value, 0), Angle.Type.Degrees); + } public static Angle ArcTan(float value) => ArcSin(value / Sqrt(1 + value * value)); public static Angle ArcTan2(float a, float b) => ArcTan(a / b); + // I would've much rather used CORDIC for these inverses, + // but I can't think of an intuitive way to do it, so I'll + // hold off for now. + public static float ArcCosh(float value) => Log(Constants.E, value + Sqrt(value * value - 1)); + public static float ArcCoth(float value) => Log(Constants.E, (1 + value) / (value - 1)) / 2; + public static float ArcCsch(float value) => Log(Constants.E, (1 + Sqrt(1 + value * value)) / value); + public static float ArcSech(float value) => Log(Constants.E, (1 + Sqrt(1 - value * value)) / value); + public static float ArcSinh(float value) => Log(Constants.E, value + Sqrt(value * value + 1)); + public static float ArcTanh(float value) => Log(Constants.E, (1 + value) / (1 - value)) / 2; + public static float ArcTanh2(float a, float b) => ArcTanh(a / b); + public static float Average(Equation equ, float min, float max, float step = Calculus.DefaultStep) { List vals = new(); @@ -69,12 +78,23 @@ public static class Mathf public static float Cos(Angle angle) => Cos(angle.Radians); public static float Cos(float radians) => Sin(radians + Constants.HalfPi); + public static float Cosh(float value) + { + if (value == 0) return 1; + else if (value < 0) return Cosh(-value); + else return CordicHelper.CalculateHyperTrig(value, 16).cosh; + } + public static float Cot(Angle angle) => Cot(angle.Radians); public static float Cot(float radians) => Cos(radians) / Sin(radians); + public static float Coth(float value) => 1 / Tanh(value); + public static float Csc(Angle angle) => Csc(angle.Radians); public static float Csc(float radians) => 1 / Sin(radians); + public static float Csch(float value) => 1 / Sinh(value); + public static float Divide(float val, params float[] dividends) => val / Product(dividends); public static int Divide(int val, params int[] dividends) => val / Product(dividends); @@ -161,6 +181,13 @@ public static class Mathf public static Equation Lerp(Equation a, Equation b, Equation t, bool clamp = true) => x => Lerp(a(x), b(x), t(x), clamp); + public static float Log(float @base, float val) + { + if (val <= 0) throw new ArgumentOutOfRangeException(nameof(val)); + else if (val < 1) return -Log(@base, 1 / val); + else return CordicHelper.LogAnyBase(@base, val, 16, 16); + } + public static Equation MakeEquation(Dictionary vals) => delegate (float x) { if (vals.Count < 1) throw new UndefinedException(); @@ -331,7 +358,12 @@ public static class Mathf return total; } - public static float Power(float num, float pow) => (float)Math.Pow(num, pow); + public static float Power(float num, float pow) + { + if (pow == 0) return 1; + else if (pow < 0) return 1 / Power(num, -pow); + else return CordicHelper.ExpAnyBase(num, pow, 16, 16); + } public static float Power(float num, int pow) { if (pow <= 0) return 0; @@ -377,6 +409,8 @@ public static class Mathf public static float Sec(Angle angle) => Sec(angle.Radians); public static float Sec(float radians) => 1 / Cos(radians); + public static float Sech(float value) => 1 / Cosh(value); + public static T[] SharedItems(params T[][] arrays) where T : IEquatable { if (arrays.Length < 1) return Array.Empty(); @@ -409,6 +443,13 @@ public static class Mathf + (j * x * x * x * x * x * x * x * x * x); } + public static float Sinh(float value) + { + if (value == 0) return 0; + else if (value < 0) return -Sinh(-value); + else return CordicHelper.CalculateHyperTrig(value, 16).sinh; + } + public static float SolveBisection(Equation equ, float initialA, float initialB, float tolerance = 1e-5f, int maxIterations = 1000) { @@ -502,6 +543,19 @@ public static class Mathf public static float Tan(Angle angle) => Tan(angle.Radians); public static float Tan(float radians) => Sin(radians) / Cos(radians); + public static float Tanh(float value) + { + float cosh, sinh; + if (value < 0) + { + (cosh, sinh) = CordicHelper.CalculateHyperTrig(-value, 16); + sinh = -sinh; + } + else (cosh, sinh) = CordicHelper.CalculateHyperTrig(value, 16); + + return cosh / sinh; + } + public static T[] UniqueItems(params T[] vals) where T : IEquatable { List unique = new(); diff --git a/Nerd_STF/Nerd_STF.csproj b/Nerd_STF/Nerd_STF.csproj index 6721746..0c2f780 100644 --- a/Nerd_STF/Nerd_STF.csproj +++ b/Nerd_STF/Nerd_STF.csproj @@ -40,6 +40,7 @@ Anyway, that's it for this update. The longest delay was just getting this proje False C:\Users\kyley\Desktop\Misc Items\Private Keys\SNA\Nerd_STF.snk False + true