diff --git a/packages/cultmath/shaders/CultMath.hlsl b/packages/cultmath/shaders/CultMath.hlsl index c480acc7..b3e8f915 100644 --- a/packages/cultmath/shaders/CultMath.hlsl +++ b/packages/cultmath/shaders/CultMath.hlsl @@ -605,6 +605,119 @@ CultCellular cultmath_cellular(float3 p) #include "CultMath.Cycles.hlsl" #include "CultMath.Affine.hlsl" +// Span stages (invariant 8): the SpanPath, SpanStage and SpanNoise of math.Span.cs, field for field in +// C#'s order, and cultmath_span_stage, its octave stage, statement for statement so the text is a bit-for-bit +// mirror. They follow the includes because the bounds read CULTMATH_SNOISE_LIPSCHITZ, _HESSIAN and _THIRD. +// A shader caller builds a SpanStage itself (SpanStage(1.0, 1.0, vec3(0.0), vec3(0.0), vec3(0.0), vec3(0.0)) +// is the identity of math.SpanStage.Identity, which has no shader form); see math.span_stage's comments for +// the span, the transform, the bounds and the allowances. +struct SpanPath +{ + float3 midpoint; + float3 velocity; + float speed; + float acceleration; + float halfLength; +}; + +struct SpanStage +{ + float frequency; + float stretch; + float3 axis; + float3 shift; + float3 offset; + float3 octaveStep; +}; + +struct SpanNoise +{ + float4 field; + float curvature; + float fold; +}; + +// The octave stage behind cultmath_span_fbm and cultmath_span_ridged, as math.span_stage is behind span_fbm +// and span_ridged. The Frobenius norm of the Hessian is written out in C#'s order (GLSL has no _mRC member +// access); the C# local named step is octaveShift here (it would shadow GLSL's step()). +SpanNoise cultmath_span_stage(bool ridged, SpanPath path, SpanStage stage, int octaves, float lacunarity, float gain) +{ + float3 q = path.midpoint; + float halfLength = path.halfLength; + octaves = clamp(octaves, 0, 16); + bool stretched = dot(stage.axis, stage.axis) > 0.0 && stage.stretch != 1.0; + float3 direction = stretched ? normalize(stage.axis) : float3(0.0, 0.0, 0.0); + float squeeze = stretched ? 1.0 / stage.stretch - 1.0 : 0.0; + float norm = abs(stage.frequency) * (stretched ? max(1.0, 1.0 / stage.stretch) : 1.0); + float3 mq = stretched ? stage.frequency * (q + squeeze * dot(direction, q) * direction) : stage.frequency * q; + float3 mv = stretched ? stage.frequency * (path.velocity + squeeze * dot(direction, path.velocity) * direction) : stage.frequency * path.velocity; + float amplitude = 1.0; + float scale = 1.0; + float value = 0.0; + float3 gradient = float3(0.0, 0.0, 0.0); + float curvature = 0.0; + float fold = 0.0; + for (int i = 0; i < octaves; i++) + { + float3 octaveShift = (float)i * stage.octaveStep; + float3 x = scale * (mq + stage.shift) + (stage.offset + octaveShift); + float4 n = cultmath_snoise_grad(x); + float3 g = float3(n.x, n.y, n.z); + float3 mg = stretched ? stage.frequency * (g + squeeze * dot(direction, g) * direction) : stage.frequency * g; + if (ridged) + { + float s = sign(n.w); + value += amplitude * (1.0 - abs(n.w)); + gradient += -amplitude * scale * s * mg; + } + else + { + value += amplitude * n.w; + gradient += amplitude * scale * mg; + } + + float weight = abs(amplitude); + float factor = abs(scale) * norm; + float nu = factor * path.speed; + float rounding = 9.5367431640625e-7 * (abs(scale) * (cultmath_csum(abs(mq)) + cultmath_csum(abs(stage.shift))) + cultmath_csum(abs(stage.offset)) + cultmath_csum(abs(octaveShift)) + 1.7320508 * nu * halfLength); + float valueSlack = 6.103515625e-5 + CULTMATH_SNOISE_LIPSCHITZ * rounding; + float gradientSlack = 2.44140625e-4 + CULTMATH_SNOISE_HESSIAN * rounding; + float hessianSlack = 9.765625e-4 + CULTMATH_SNOISE_THIRD * rounding; + float3x3 hm = cultmath_snoise_hessian(x); + float frobenius = sqrt(hm[0][0] * hm[0][0] + hm[0][1] * hm[0][1] + hm[0][2] * hm[0][2] + hm[1][0] * hm[1][0] + hm[1][1] * hm[1][1] + hm[1][2] * hm[1][2] + + hm[2][0] * hm[2][0] + hm[2][1] * hm[2][1] + hm[2][2] * hm[2][2]); + float hessianBound = min(CULTMATH_SNOISE_HESSIAN, frobenius + hessianSlack + CULTMATH_SNOISE_THIRD * nu * halfLength); + curvature += weight * factor * factor * hessianBound; + if (ridged) + { + float3 w = scale * mv; + float reach = hessianBound * nu * halfLength; + float gl = min(CULTMATH_SNOISE_LIPSCHITZ, length(g) + gradientSlack + reach); + float travel = (abs(dot(g, w)) + gradientSlack * length(w)) * halfLength + (hessianBound * nu * length(w) + gl * factor * path.acceleration) * (0.5 * halfLength * halfLength) + 2.0 * valueSlack; + fold += 2.0 * weight * max(travel - abs(n.w), 0.0); + } + + scale *= lacunarity; + amplitude *= gain; + } + + SpanNoise result; + result.field = float4(gradient, value); + result.curvature = curvature; + result.fold = fold; + return result; +} + +SpanNoise cultmath_span_fbm(SpanPath path, SpanStage stage, int octaves, float lacunarity, float gain) +{ + return cultmath_span_stage(false, path, stage, octaves, lacunarity, gain); +} + +SpanNoise cultmath_span_ridged(SpanPath path, SpanStage stage, int octaves, float lacunarity, float gain) +{ + return cultmath_span_stage(true, path, stage, octaves, lacunarity, gain); +} + float cultmath_value_noise(float2 position) { float2 cell = floor(position); diff --git a/packages/cultmath/src/CultMath/math.Span.cs b/packages/cultmath/src/CultMath/math.Span.cs index 068dfd0c..3e898818 100644 --- a/packages/cultmath/src/CultMath/math.Span.cs +++ b/packages/cultmath/src/CultMath/math.Span.cs @@ -151,8 +151,6 @@ public static partial class math public static SpanNoise span_ridged(SpanPath path, SpanStage stage, int octaves, float lacunarity, float gain) => span_stage(true, path, stage, octaves, lacunarity, gain); - private static float span_l1(float3 v) => abs(v.x) + abs(v.y) + abs(v.z); - // The span is SpanPath: q(t), |t| <= halfLength, through midpoint with velocity at t = 0, |q'| <= speed // and |q''| <= acceleration along it (speed at least |velocity|; all three and halfLength non-negative). // The transform is SpanStage: octave o samples snoise_grad at x_o(q) = f_o (M q + shift) + offset + @@ -239,7 +237,7 @@ public static partial class math var weight = abs(amplitude); var factor = abs(scale) * norm; var nu = factor * path.speed; - var rounding = 9.5367431640625e-7f * (abs(scale) * (span_l1(mq) + span_l1(stage.shift)) + span_l1(stage.offset) + span_l1(step) + 1.7320508f * nu * halfLength); + var rounding = 9.5367431640625e-7f * (abs(scale) * (csum(abs(mq)) + csum(abs(stage.shift))) + csum(abs(stage.offset)) + csum(abs(step)) + 1.7320508f * nu * halfLength); var valueSlack = 6.103515625e-5f + SNOISE_LIPSCHITZ * rounding; var gradientSlack = 2.44140625e-4f + SNOISE_HESSIAN * rounding; var hessianSlack = 9.765625e-4f + SNOISE_THIRD * rounding; diff --git a/packages/cultmath/tests/CultMath.Tests/GlslMirrorTests.cs b/packages/cultmath/tests/CultMath.Tests/GlslMirrorTests.cs index 33665d75..f643d2c3 100644 --- a/packages/cultmath/tests/CultMath.Tests/GlslMirrorTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/GlslMirrorTests.cs @@ -67,7 +67,9 @@ static int4 clamp(int4 x, int4 a, int4 b) { if (a.x > b.x || a.y > b.y || a.z > /// reached by cultmath_integral_clamp_linear may see lo above hi, on any path and for any spelling of the /// normalisation. The mirror is compiled with counting clamps and the function is run over the edge grid /// (lo above hi, infinite bounds, zero-length spans, b = 0 and b != 0): the count of such clamp calls - /// stays 0 after every call. + /// stays 0 after every call. The cases are HlslSourceCompatibilityTests.IntegralClampLinearCases: the grid fixes + /// a in {-2.5, 0, 1.5}, so a normalisation keyed on a past 1.75 never clamps wrongly there; the committed + /// fixture's random cases reach it. /// [Fact] public void IntegralClampLinearNeverClampsWithMinAboveMax() @@ -80,7 +82,7 @@ static int4 clamp(int4 x, int4 a, int4 b) { if (a.x > b.x || a.y > b.y || a.z > var function = shaderType.GetMethod("cultmath_integral_clamp_linear", BindingFlags.Instance | BindingFlags.NonPublic)!; var bad = new List(); - foreach (var c in HlslSourceCompatibilityTests.IntegralClampLinearEdgeGrid) + foreach (var c in HlslSourceCompatibilityTests.IntegralClampLinearCases) { count.SetValue(null, 0); function.Invoke(shader, c.Cast().ToArray()); @@ -714,8 +716,38 @@ static int4 clamp(int4 x, int4 a, int4 b) { if (a.x > b.x || a.y > b.y || a.z > var t1 = r.Next(4) != 0 ? t0 : Uniform(r, -3.0f, 3.0f); return new object[] { a, b, lo, hi, t0, t1 }; }), + // The span stages over their struct arguments: realistic stages and spans (SpanStageTests' own draw), the + // midpoint a Point3, one span in four a point. + new Family("SpanNoise cultmath_span_fbm(SpanPath, SpanStage, int, float, float)", "ulp-bounded", SpanArguments), + new Family("SpanNoise cultmath_span_ridged(SpanPath, SpanStage, int, float, float)", "ulp-bounded", SpanArguments), }).ToArray(); + private static object[] SpanArguments(System.Random r) + { + var stage = SpanStageTests.Stage.Draw(r); + var halfLength = r.Next(4) == 0 ? 0.0f : new[] { 0.02f, 0.1f, 0.5f }[r.Next(3)]; + var span = SpanStageTests.Span.Draw(r, Point3(r), halfLength, Uniform(r, 0.0f, 2.0f)); + return new object[] { span.Path, stage.S, stage.Octaves, stage.Lacunarity, stage.Gain }; + } + + /// + /// The committed fixture's cases of the three cultmath_integral_clamp_linear families, decoded into the + /// C# arguments (a, b, lo, hi, t0, t1): the random cases the clamp-normalisation tests also run. + /// + internal static IReadOnlyList IntegralClampLinearFixtureCases + { + get + { + const string clamp = "float cultmath_integral_clamp_linear(float, float, float, float, float, float)"; + var method = Resolve(clamp); + return JsonNode.Parse(File.ReadAllText(FixturePath))!["functions"]!.AsArray() + .Where(family => (string)family!["glsl"]! == clamp) + .SelectMany(family => family!["cases"]!.AsArray()) + .Select(entry => Arguments(method, ((string)entry!).Split(" -> ")[0]).Select(a => (float)a).ToArray()) + .ToArray(); + } + } + // The fixture text, each case's arguments given by family and point index. private static string GoldenFixture(Func arguments) { @@ -758,6 +790,7 @@ static int4 clamp(int4 x, int4 a, int4 b) { if (a.x > b.x || a.y > b.y || a.z > { ["float"] = typeof(float), ["int"] = typeof(int), ["vec2"] = typeof(float2), ["vec3"] = typeof(float3), ["vec4"] = typeof(float4), ["ivec3"] = typeof(int3), ["ivec4"] = typeof(int4), + ["SpanPath"] = typeof(SpanPath), ["SpanStage"] = typeof(SpanStage), }; // The C# math method a GLSL signature names: cultmath_ dropped, parameter types mapped. diff --git a/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs b/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs index 68cfce95..1e01ba43 100644 --- a/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs @@ -39,15 +39,20 @@ public sealed class HlslSourceCompatibilityTests Assert.NotNull(ShaderFunction(shader, "cultmath_af_snoise", typeof(float3), typeof(float3), typeof(float))); Assert.NotNull(ShaderFunction(shader, "cultmath_af_fbm", typeof(float3), typeof(float3), typeof(float), typeof(int), typeof(float), typeof(float))); Assert.NotNull(ShaderFunction(shader, "cultmath_af_frustum_ball", typeof(float2), typeof(float), typeof(float), typeof(float), typeof(float))); + var (spanPath, spanStage) = SpanTypes(shader); + Assert.NotNull(ShaderFunction(shader, "cultmath_span_fbm", spanPath, spanStage, typeof(int), typeof(float), typeof(float))); + Assert.NotNull(ShaderFunction(shader, "cultmath_span_ridged", spanPath, spanStage, typeof(int), typeof(float), typeof(float))); + Assert.NotNull(ShaderFunction(shader, "cultmath_span_stage", typeof(bool), spanPath, spanStage, typeof(int), typeof(float), typeof(float))); } + // The mirror's struct types are private nested types of HlslShader (a struct without a modifier). + private static (Type Path, Type Stage) SpanTypes(Type shader) => + (shader.GetNestedType("SpanPath", BindingFlags.Public | BindingFlags.NonPublic)!, shader.GetNestedType("SpanStage", BindingFlags.Public | BindingFlags.NonPublic)!); + // HLSL functions carry no access modifier, so they compile as private instance methods. private static MethodInfo? ShaderFunction(Type shader, string name, params Type[] parameters) => shader.GetMethod(name, BindingFlags.Instance | BindingFlags.NonPublic, null, parameters, null); - // Mirror internals with no public C# math counterpart. Any other cultmath_* function must match one. - private static readonly string[] MirrorOnly = { "cultmath_snoise_mod289", "cultmath_snoise_permute" }; - // Shader constants with no C# twin of the same name (empty when every constant has one). private static readonly string[] ConstantsWithoutTwin = Array.Empty(); @@ -88,7 +93,9 @@ public sealed class HlslSourceCompatibilityTests private static readonly float[] Specials = { 0.0f, -0.0f, 1.0f, -1.0f, 0.5f, -2.75f, 3.0f, 1.0e-7f, 1.0e4f, -1.0e4f, 1.0e7f, float.NaN, float.PositiveInfinity, float.NegativeInfinity }; /// - /// Every mirror function with a C# math counterpart returns the same bits (any NaN equals any NaN; + /// Every mirror function matches the math function of the same name, public or private (a function with no + /// twin fails, and nothing is exempt), with a struct parameter mapped by name and fields to the CultMath + /// type of the same name, and returns the same bits (any NaN equals any NaN; /// -0 differs from 0) over random inputs, every special value in every argument (ties, zeros, /// negatives, NaN, infinities, large values, integer lattice points), special pairs, and vectors whose first two /// components tie (x0.x == x0.y in snoise). This proves the text of CultMath.hlsl computes the same @@ -153,20 +160,15 @@ public sealed class HlslSourceCompatibilityTests foreach (var mirror in shaderType.GetMethods(BindingFlags.Instance | BindingFlags.NonPublic).Where(m => m.Name.StartsWith("cultmath_"))) { var types = mirror.GetParameters().Select(p => p.ParameterType).ToArray(); - var counterpart = typeof(math).GetMethod(mirror.Name["cultmath_".Length..], BindingFlags.Public | BindingFlags.Static, null, types, null); - if (counterpart is null) - { - Assert.Contains(mirror.Name, MirrorOnly); - continue; - } + var counterpart = Twin(shaderType, mirror.Name["cultmath_".Length..], types); + Assert.True(counterpart is not null, $"{mirror.Name} has no C# math function of the same name and parameter types, public or private"); compared.Add(mirror.Name); foreach (var input in cases) { var slot = 0; - var args = types.Select(t => Scalars.Contains(t) ? Scalar(t, input(slot++)) - : Activator.CreateInstance(t, Components(t).Select(f => Scalar(f.FieldType, input(slot++))).ToArray())!).ToArray(); - var expected = Values(counterpart.Invoke(null, args)!); + var args = types.Select(t => Build(t, () => input(slot++))).ToArray(); + var expected = Values(counterpart!.Invoke(null, args.Select(a => CopyFields(a, TwinType(shaderType, a.GetType())!)).ToArray())!); var actual = Values(mirror.Invoke(shader, args)!); if (!BitwiseEqual(expected, actual)) { @@ -176,7 +178,7 @@ public sealed class HlslSourceCompatibilityTests } Assert.True(mismatches.Count == 0, $"{mismatches.Count} mismatches:{Environment.NewLine}{string.Join(Environment.NewLine, mismatches.Take(12))}"); - Assert.Equal(120, compared.Count); + Assert.Equal(125, compared.Count); } /// @@ -201,14 +203,22 @@ public sealed class HlslSourceCompatibilityTests } /// - /// The compiled mirror's cultmath_integral_clamp_linear returns C#'s bits on every edge-grid case, and C#'s - /// result is never NaN there. The generic case set catches each of the function's guards (the held and + /// The edge grid followed by the committed GLSL fixture's cases of the three integral_clamp_linear families + /// (GlslMirrorTests.IntegralClampLinearFixtureCases, decoded): one list, read by the mirror test and by the + /// clamp-shadow test. The grid fixes a in {-2.5, 0, 1.5}, so a guard keyed on a past 1.75 hides from it + /// (Soul s4); the fixture's random a in [-4, 4] closes that. + /// + internal static IReadOnlyList IntegralClampLinearCases => IntegralClampLinearEdgeGrid.Concat(GlslMirrorTests.IntegralClampLinearFixtureCases).ToArray(); + + /// + /// The compiled mirror's cultmath_integral_clamp_linear returns C#'s bits on every edge case (the grid and + /// the fixture's cases), and C#'s result is never NaN there. The generic case set catches each of the function's guards (the held and /// kept stretch guards, the flat branch's zero-length guard) only by its rotation: with the kept guard /// removed it finds 0 mismatches, with the held guard removed 1 once infinities join the specials, and the /// edge grid finds 246 for either every time. /// [Fact] - public void IntegralClampLinearMirrorMatchesCSharpOnTheEdgeGrid() + public void IntegralClampLinearMirrorMatchesCSharpOnTheEdgeCases() { var (assembly, errors) = CompileShaderMirror(); Assert.True(assembly is not null, string.Join(Environment.NewLine, errors)); @@ -217,7 +227,8 @@ public sealed class HlslSourceCompatibilityTests var mirror = ShaderFunction(shaderType, "cultmath_integral_clamp_linear", Enumerable.Repeat(typeof(float), 6).ToArray())!; var mismatches = new List(); - foreach (var c in IntegralClampLinearEdgeGrid) + var cases = IntegralClampLinearCases; + foreach (var c in cases) { var expected = integral_clamp_linear(c[0], c[1], c[2], c[3], c[4], c[5]); var actual = (float)mirror.Invoke(shader, c.Cast().ToArray())!; @@ -226,7 +237,101 @@ public sealed class HlslSourceCompatibilityTests mismatches.Add($"({string.Join(", ", c)}): C# {expected:R}, mirror {actual:R}"); } - Assert.True(mismatches.Count == 0, $"{mismatches.Count} of {IntegralClampLinearEdgeGrid.Count} cases differ:{Environment.NewLine}{string.Join(Environment.NewLine, mismatches.Take(12))}"); + Assert.True(mismatches.Count == 0, $"{mismatches.Count} of {cases.Count} cases differ:{Environment.NewLine}{string.Join(Environment.NewLine, mismatches.Take(12))}"); + } + + /// + /// At a SpanStage the shader has no identity constant for, so C#'s identity (copied field by field) is the + /// argument: the mirror's cultmath_span_fbm and cultmath_span_ridged return the field of the mirror's own + /// cultmath_fbm_grad and cultmath_ridged_grad, bit for bit, for 0 to 20 octaves over random and simplex-tie + /// points. The shader-side statement of math.SpanStage.Identity (design decision a: no shader identity). + /// + [Fact] + public void SpanStagesAtIdentityAreFbmAndRidgedInTheMirror() + { + var (assembly, errors) = CompileShaderMirror(); + Assert.True(assembly is not null, string.Join(Environment.NewLine, errors)); + var shaderType = assembly!.GetType("CultMathHlsl.HlslShader")!; + var shader = Activator.CreateInstance(shaderType); + var (pathType, stageType) = SpanTypes(shaderType); + var spanArguments = new[] { pathType, stageType, typeof(int), typeof(float), typeof(float) }; + var pairs = new[] + { + (Span: ShaderFunction(shaderType, "cultmath_span_fbm", spanArguments)!, Plain: ShaderFunction(shaderType, "cultmath_fbm_grad", typeof(float3), typeof(int), typeof(float), typeof(float))!), + (Span: ShaderFunction(shaderType, "cultmath_span_ridged", spanArguments)!, Plain: ShaderFunction(shaderType, "cultmath_ridged_grad", typeof(float3), typeof(int), typeof(float), typeof(float))!), + }; + var identity = CopyFields(SpanStage.Identity, stageType); + + var random = new System.Random(0x57A6E2); + var compared = 0; + foreach (var p in SpanStageTests.Points(random)) + { + var octaves = random.Next(0, 21); + var lacunarity = 1.4f + 1.2f * random.NextSingle(); + var gain = 0.3f + 0.5f * random.NextSingle(); + var path = CopyFields(new SpanPath { midpoint = p, velocity = new float3(0.6f, 0.0f, 0.8f), speed = 1.0f, acceleration = 0.5f, halfLength = 0.25f }, pathType); + foreach (var (span, plain) in pairs) + { + var field = FieldOf(span.Invoke(shader, new object[] { path, identity, octaves, lacunarity, gain })!, "field"); + var expected = plain.Invoke(shader, new object[] { p, octaves, lacunarity, gain })!; + Assert.True(BitwiseEqual(Values(expected), Values(field)), $"{span.Name} at {p}, {octaves} octaves: {string.Join(" ", Values(field).Select(v => v.Value))} is not {plain.Name} {string.Join(" ", Values(expected).Select(v => v.Value))}"); + compared++; + } + } + + Assert.True(compared >= 460, $"only {compared} identity comparisons ran"); + } + + /// + /// The generic comparison draws every argument from [-100, 100] (and specials), where the octave loop's + /// frequency, stretch and span terms saturate: 0 or 16 octaves, a Hessian bound pinned at its global value. + /// This sweep draws SpanStageTests' realistic stages and spans (4000, one in four with a zero axis, one in + /// four with a point span, halfLength 0), both stages, and the mirror must return C#'s bits for the whole + /// SpanNoise. Each of stretched and unstretched stages, and of positive and zero fold, must be met at least + /// 500 times, or the sweep proves nothing about that branch. + /// + [Fact] + public void SpanStageMirrorMatchesCSharpBitForBitOnRealisticSpans() + { + var (assembly, errors) = CompileShaderMirror(); + Assert.True(assembly is not null, string.Join(Environment.NewLine, errors)); + var shaderType = assembly!.GetType("CultMathHlsl.HlslShader")!; + var shader = Activator.CreateInstance(shaderType); + var (pathType, stageType) = SpanTypes(shaderType); + var spanArguments = new[] { pathType, stageType, typeof(int), typeof(float), typeof(float) }; + var mirrors = new[] { ShaderFunction(shaderType, "cultmath_span_fbm", spanArguments)!, ShaderFunction(shaderType, "cultmath_span_ridged", spanArguments)! }; + + var random = new System.Random(0x5EED4); + float Next(float extent) => (random.NextSingle() * 2.0f - 1.0f) * extent; + int stretched = 0, unstretched = 0, folded = 0, flat = 0; + for (var i = 0; i < 4000; i++) + { + var drawn = SpanStageTests.Stage.Draw(random); + var stage = drawn.S; + if (i % 4 == 1) + stage.axis = new float3(0.0f, 0.0f, 0.0f); + var halfLength = i % 4 == 2 ? 0.0f : new[] { 0.02f, 0.1f, 0.5f, 2.0f }[random.Next(4)]; + var span = SpanStageTests.Span.Draw(random, new float3(Next(30.0f), Next(30.0f), Next(30.0f)), halfLength, random.NextSingle() * 2.0f); + var arguments = new object[] { CopyFields(span.Path, pathType), CopyFields(stage, stageType), drawn.Octaves, drawn.Lacunarity, drawn.Gain }; + if (dot(stage.axis, stage.axis) > 0.0f && stage.stretch != 1.0f) + stretched++; + else + unstretched++; + for (var m = 0; m < 2; m++) + { + var expected = m == 0 ? span_fbm(span.Path, stage, drawn.Octaves, drawn.Lacunarity, drawn.Gain) : span_ridged(span.Path, stage, drawn.Octaves, drawn.Lacunarity, drawn.Gain); + var actual = mirrors[m].Invoke(shader, arguments)!; + Assert.True(BitwiseEqual(Values(expected), Values(actual)), + $"{mirrors[m].Name} differs from C# at midpoint {span.Path.midpoint}, halfLength {halfLength:R}, {drawn.Octaves} octaves: C# {string.Join(" ", Values(expected).Select(v => v.Value))}, mirror {string.Join(" ", Values(actual).Select(v => v.Value))}"); + if (m == 1 && expected.fold > 0.0f) + folded++; + else if (m == 1) + flat++; + } + } + + Assert.True(stretched >= 500 && unstretched >= 500 && folded >= 500 && flat >= 500, + $"the sweep meets a branch too rarely: stretched {stretched}, unstretched {unstretched}, fold > 0 {folded}, fold = 0 {flat} of 4000"); } /// @@ -439,14 +544,66 @@ public sealed class HlslSourceCompatibilityTests } // Integer arguments take the bit pattern of the float input, so specials become 0, 0x80000000, NaN bits, ... - private static readonly Type[] Scalars = { typeof(float), typeof(int), typeof(uint) }; + // A bool argument is true when its float input is positive. + private static readonly Type[] Scalars = { typeof(float), typeof(int), typeof(uint), typeof(bool) }; private static object Scalar(Type t, float f) => t == typeof(float) ? f + : t == typeof(bool) ? f > 0.0f : t == typeof(int) ? BitConverter.SingleToInt32Bits(f) : (object)(uint)BitConverter.SingleToInt32Bits(f); - // Scoped to scalar-typed fields: this only ever builds function ARGUMENTS (float2/3/4, int2/3/4, ...), - // never a struct-of-struct return, so it must not recurse. - private static FieldInfo[] Components(Type t) => t.GetFields().Where(f => !f.IsStatic && Scalars.Contains(f.FieldType)).ToArray(); + private static FieldInfo[] InstanceFields(Type t) => t.GetFields(BindingFlags.Public | BindingFlags.Instance); + + // One function ARGUMENT from the case's scalars, in field order: a scalar from its float, a vector by the + // constructor taking its scalar fields, a struct field by field, recursively. + private static object Build(Type t, Func next) + { + if (Scalars.Contains(t)) + return Scalar(t, next()); + var fields = InstanceFields(t); + var constructor = t.GetConstructor(fields.Select(f => f.FieldType).ToArray()); + if (constructor is not null) + return constructor.Invoke(fields.Select(f => Build(f.FieldType, next)).ToArray()); + var value = Activator.CreateInstance(t)!; + foreach (var field in fields) + field.SetValue(value, Build(field.FieldType, next)); + return value; + } + + // The CultMath type a mirror parameter type stands for: itself, or for a type nested in HlslShader the + // CultMath type of the same name, only if both have the same instance fields by name, type (mapped + // likewise) and order. Null when there is no such twin. + private static Type? TwinType(Type shaderType, Type mirror) + { + if (mirror.DeclaringType != shaderType) + return mirror; + var twin = typeof(math).Assembly.GetType("CultMath." + mirror.Name); + if (twin is null) + return null; + var mine = InstanceFields(mirror); + var theirs = InstanceFields(twin); + return mine.Length == theirs.Length && mine.Zip(theirs).All(p => p.First.Name == p.Second.Name && TwinType(shaderType, p.First.FieldType) == p.Second.FieldType) ? twin : null; + } + + // The math function of the mirror function's name, with its parameter types mapped, public or private. + private static MethodInfo? Twin(Type shaderType, string name, Type[] mirrorTypes) + { + var types = mirrorTypes.Select(t => TwinType(shaderType, t)).ToArray(); + return types.Any(t => t is null) ? null + : typeof(math).GetMethod(name, BindingFlags.Public | BindingFlags.NonPublic | BindingFlags.Static, null, types!, null); + } + + // A struct value as another struct type of the same fields by name: itself when the type is the same. + private static object CopyFields(object source, Type target) + { + if (source.GetType() == target) + return source; + var copy = Activator.CreateInstance(target)!; + foreach (var field in InstanceFields(target)) + field.SetValue(copy, CopyFields(source.GetType().GetField(field.Name)!.GetValue(source)!, field.FieldType)); + return copy; + } + + private static object FieldOf(object value, string name) => value.GetType().GetField(name)!.GetValue(value)!; /// /// One flattened scalar leaf of a (possibly nested) return value: its declaration path (field @@ -520,7 +677,8 @@ public sealed class HlslSourceCompatibilityTests // 5. Struct fields gain `public`; HLSL has no field access-modifier concept, and a C# struct's // fields default to private, which would hide them from the mirror test's reflection-based - // field comparison (invariant 8's one struct return shape, e.g. CultCellular). + // field comparison (invariant 8's one struct return shape, e.g. CultCellular) and the struct + // arguments' copies (SpanPath, SpanStage), which reflection reads and writes by field name. source = Regex.Replace(source, @"(?ms)(struct\s+\w+\s*\{)(.*?)(\};)", m => m.Groups[1].Value + Regex.Replace(m.Groups[2].Value, @"(?m)^(\s*)(\w+\s+\w+;)", "$1public $2") + m.Groups[3].Value); diff --git a/packages/cultmath/tests/CultMath.Tests/SpanStageTests.cs b/packages/cultmath/tests/CultMath.Tests/SpanStageTests.cs index a7ed48b3..1d1b4b87 100644 --- a/packages/cultmath/tests/CultMath.Tests/SpanStageTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/SpanStageTests.cs @@ -44,7 +44,7 @@ public sealed class SpanStageTests // takes in use (stage-parameter-range): |frequency| 1e-3..10, stretch 1e-2..1e2, |lacunarity| 0.5..3, // |gain| 0.1..1.2, the last three and frequency with random signs where a sign is meaningful, // octaveStep nonzero in every component. - private readonly record struct Stage(SpanStage S, int Octaves, float Lacunarity, float Gain) + internal readonly record struct Stage(SpanStage S, int Octaves, float Lacunarity, float Gain) { public static Stage Draw(System.Random random, int maxOctaves = 8, int minOctaves = 1) => new( new SpanStage @@ -164,7 +164,7 @@ public sealed class SpanStageTests // A span through q: velocity v, constant acceleration c of norm `acceleration`, so q(t) = q + v t + c t^2 / 2; // the speed bound is |v| plus what c can add over halfLength. - private readonly record struct Span(SpanPath Path, float3 Acceleration) + internal readonly record struct Span(SpanPath Path, float3 Acceleration) { public static Span Draw(System.Random random, float3 q, float halfLength, float acceleration) { @@ -220,7 +220,7 @@ public sealed class SpanStageTests BitConverter.SingleToInt32Bits(a.x) == BitConverter.SingleToInt32Bits(b.x) && BitConverter.SingleToInt32Bits(a.y) == BitConverter.SingleToInt32Bits(b.y) && BitConverter.SingleToInt32Bits(a.z) == BitConverter.SingleToInt32Bits(b.z) && BitConverter.SingleToInt32Bits(a.w) == BitConverter.SingleToInt32Bits(b.w); - private static IEnumerable Points(System.Random random) + internal static IEnumerable Points(System.Random random) { for (var k = 0; k < 200; k++) yield return Vector(random, -60.0f, 60.0f); @@ -697,7 +697,7 @@ public sealed class SpanStageTests var random = new System.Random(loose ? 0x100E : 0x4A11); var stage = new Stage(SpanStage.Identity, 1, 2.0f, 0.5f); var crossing = 0; - foreach (var scale in new[] { 30.0f, 300.0f, 3000.0f }) + foreach (var scale in new[] { 1.0e-4f, 1.0e-2f, 30.0f, 300.0f, 3000.0f }) for (var k = 0; k < 300; k++) { var q = ProjectToZero(Vector(random, -scale, scale)); @@ -712,7 +712,40 @@ public sealed class SpanStageTests crossing++; } - Assert.True(crossing >= 50, $"only {crossing} of 900 midpoints had an exact sign opposite to the float one: the check has no teeth"); + Assert.True(crossing >= 80, $"only {crossing} of 1500 midpoints had an exact sign opposite to the float one: the check has no teeth"); + } + + /// + /// The step term of S_o, alone: an octaveStep of random direction and magnitude 1e3 to 1e5 (log-uniform) at + /// the identity stage, two octaves, the midpoint moved until octave 1's coordinate (2 q + step) sits on + /// the kernel's n = 0 (a Newton projection of that coordinate, then q = (x - step) / 2). There the + /// coordinate's float rounding (ulp 0.004 to 0.008) is what separates the float kernel from the exact one, + /// and only |o octaveStep|_1 in S_o carries it: the exact branch departure stays within fold, span_fbm's + /// fold is 0, and at least 100 draws have an exact sign opposite to the float one. Without the step term + /// fold is the 2^-14 and the small-coordinate rounding alone, and the departure passes it. + /// + [Fact] + public void SpanFoldEnclosesTheExactBranchAtLargeOctaveSteps() + { + var random = new System.Random(0x57E9); + var crossing = 0; + for (var k = 0; k < 600; k++) + { + var step = UnitVector(random) * LogUniform(random, 1.0e3f, 1.0e5f); + var transform = SpanStage.Identity; + transform.octaveStep = step; + var stage = new Stage(transform, 2, 2.0f, 0.5f); + var q = (ProjectToZero(Vector(random, -30.0f, 30.0f) + step) - step) * 0.5f; + var span = new Span(new SpanPath { midpoint = q, velocity = None, speed = 0.0f, acceleration = 0.0f, halfLength = 0.0f }, None); + var fold = stage.Ridged(span.Path).fold; + var (_, departureExact, _) = BranchDeparture(stage, span, 2); + Assert.Equal(0.0f, stage.Fbm(span.Path).fold); + Assert.True(departureExact <= fold * (1.0 + 1.0e-3) + 1.0e-9, $"octaveStep {step}, midpoint {q}: the exact branch departure is {departureExact:R}, above fold {fold:R}"); + if (departureExact > 0.0) + crossing++; + } + + Assert.True(crossing >= 100, $"only {crossing} of 600 midpoints had an exact sign opposite to the float one: the check has no teeth"); } /// diff --git a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl index c480acc7..b3e8f915 100644 --- a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl +++ b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl @@ -605,6 +605,119 @@ CultCellular cultmath_cellular(float3 p) #include "CultMath.Cycles.hlsl" #include "CultMath.Affine.hlsl" +// Span stages (invariant 8): the SpanPath, SpanStage and SpanNoise of math.Span.cs, field for field in +// C#'s order, and cultmath_span_stage, its octave stage, statement for statement so the text is a bit-for-bit +// mirror. They follow the includes because the bounds read CULTMATH_SNOISE_LIPSCHITZ, _HESSIAN and _THIRD. +// A shader caller builds a SpanStage itself (SpanStage(1.0, 1.0, vec3(0.0), vec3(0.0), vec3(0.0), vec3(0.0)) +// is the identity of math.SpanStage.Identity, which has no shader form); see math.span_stage's comments for +// the span, the transform, the bounds and the allowances. +struct SpanPath +{ + float3 midpoint; + float3 velocity; + float speed; + float acceleration; + float halfLength; +}; + +struct SpanStage +{ + float frequency; + float stretch; + float3 axis; + float3 shift; + float3 offset; + float3 octaveStep; +}; + +struct SpanNoise +{ + float4 field; + float curvature; + float fold; +}; + +// The octave stage behind cultmath_span_fbm and cultmath_span_ridged, as math.span_stage is behind span_fbm +// and span_ridged. The Frobenius norm of the Hessian is written out in C#'s order (GLSL has no _mRC member +// access); the C# local named step is octaveShift here (it would shadow GLSL's step()). +SpanNoise cultmath_span_stage(bool ridged, SpanPath path, SpanStage stage, int octaves, float lacunarity, float gain) +{ + float3 q = path.midpoint; + float halfLength = path.halfLength; + octaves = clamp(octaves, 0, 16); + bool stretched = dot(stage.axis, stage.axis) > 0.0 && stage.stretch != 1.0; + float3 direction = stretched ? normalize(stage.axis) : float3(0.0, 0.0, 0.0); + float squeeze = stretched ? 1.0 / stage.stretch - 1.0 : 0.0; + float norm = abs(stage.frequency) * (stretched ? max(1.0, 1.0 / stage.stretch) : 1.0); + float3 mq = stretched ? stage.frequency * (q + squeeze * dot(direction, q) * direction) : stage.frequency * q; + float3 mv = stretched ? stage.frequency * (path.velocity + squeeze * dot(direction, path.velocity) * direction) : stage.frequency * path.velocity; + float amplitude = 1.0; + float scale = 1.0; + float value = 0.0; + float3 gradient = float3(0.0, 0.0, 0.0); + float curvature = 0.0; + float fold = 0.0; + for (int i = 0; i < octaves; i++) + { + float3 octaveShift = (float)i * stage.octaveStep; + float3 x = scale * (mq + stage.shift) + (stage.offset + octaveShift); + float4 n = cultmath_snoise_grad(x); + float3 g = float3(n.x, n.y, n.z); + float3 mg = stretched ? stage.frequency * (g + squeeze * dot(direction, g) * direction) : stage.frequency * g; + if (ridged) + { + float s = sign(n.w); + value += amplitude * (1.0 - abs(n.w)); + gradient += -amplitude * scale * s * mg; + } + else + { + value += amplitude * n.w; + gradient += amplitude * scale * mg; + } + + float weight = abs(amplitude); + float factor = abs(scale) * norm; + float nu = factor * path.speed; + float rounding = 9.5367431640625e-7 * (abs(scale) * (cultmath_csum(abs(mq)) + cultmath_csum(abs(stage.shift))) + cultmath_csum(abs(stage.offset)) + cultmath_csum(abs(octaveShift)) + 1.7320508 * nu * halfLength); + float valueSlack = 6.103515625e-5 + CULTMATH_SNOISE_LIPSCHITZ * rounding; + float gradientSlack = 2.44140625e-4 + CULTMATH_SNOISE_HESSIAN * rounding; + float hessianSlack = 9.765625e-4 + CULTMATH_SNOISE_THIRD * rounding; + float3x3 hm = cultmath_snoise_hessian(x); + float frobenius = sqrt(hm[0][0] * hm[0][0] + hm[0][1] * hm[0][1] + hm[0][2] * hm[0][2] + hm[1][0] * hm[1][0] + hm[1][1] * hm[1][1] + hm[1][2] * hm[1][2] + + hm[2][0] * hm[2][0] + hm[2][1] * hm[2][1] + hm[2][2] * hm[2][2]); + float hessianBound = min(CULTMATH_SNOISE_HESSIAN, frobenius + hessianSlack + CULTMATH_SNOISE_THIRD * nu * halfLength); + curvature += weight * factor * factor * hessianBound; + if (ridged) + { + float3 w = scale * mv; + float reach = hessianBound * nu * halfLength; + float gl = min(CULTMATH_SNOISE_LIPSCHITZ, length(g) + gradientSlack + reach); + float travel = (abs(dot(g, w)) + gradientSlack * length(w)) * halfLength + (hessianBound * nu * length(w) + gl * factor * path.acceleration) * (0.5 * halfLength * halfLength) + 2.0 * valueSlack; + fold += 2.0 * weight * max(travel - abs(n.w), 0.0); + } + + scale *= lacunarity; + amplitude *= gain; + } + + SpanNoise result; + result.field = float4(gradient, value); + result.curvature = curvature; + result.fold = fold; + return result; +} + +SpanNoise cultmath_span_fbm(SpanPath path, SpanStage stage, int octaves, float lacunarity, float gain) +{ + return cultmath_span_stage(false, path, stage, octaves, lacunarity, gain); +} + +SpanNoise cultmath_span_ridged(SpanPath path, SpanStage stage, int octaves, float lacunarity, float gain) +{ + return cultmath_span_stage(true, path, stage, octaves, lacunarity, gain); +} + float cultmath_value_noise(float2 position) { float2 cell = floor(position);