From 789549ed2fb28133230a45625d6e2893d85be687 Mon Sep 17 00:00:00 2001 From: Yingjie Wang Date: Fri, 9 Oct 2026 00:14:54 -0400 Subject: [PATCH] Fix: Validate asymptotic entries with precise roots and fallback Use scaled long-double quadratic arithmetic without explicit FMA. Validate entry candidates against backend geometry and localize uncertain entries along the original exterior trajectory. Preserve conservative miss semantics and propagate concrete entry failures. Add production-sample and numerical regression coverage. --- Makefile | 17 +- nr_spacetime_movie_renderer_design.md | 29 + src/asymptotic.c | 676 +++++++++++++++++++--- src/asymptotic.h | 10 + src/asymptotic_entry.c | 326 +++++++++++ src/asymptotic_entry.h | 136 +++++ src/asymptotic_schwarzschild.c | 62 +- src/asymptotic_schwarzschild.h | 13 + src/geodesic.c | 8 +- src/geodesic.h | 1 + src/ray.c | 7 +- tests/test_asymptotic.c | 261 ++++++++- tests/test_asymptotic_entry.c | 781 ++++++++++++++++++++++++++ tests/test_asymptotic_quadratic.c | 65 +++ tests/test_asymptotic_schwarzschild.c | 120 ++++ tests/test_geodesic.c | 2 + 16 files changed, 2421 insertions(+), 93 deletions(-) create mode 100644 src/asymptotic_entry.c create mode 100644 src/asymptotic_entry.h create mode 100644 tests/test_asymptotic_entry.c create mode 100644 tests/test_asymptotic_quadratic.c diff --git a/Makefile b/Makefile index 87ba0bf..9cb7970 100644 --- a/Makefile +++ b/Makefile @@ -108,6 +108,8 @@ TEST_OUT_DIR := $(OBJECT_DIR)/$(HDR_BUILD_TAG) TEST_TARGET := $(TEST_OUT_DIR)/test_geodesic ADAPTIVE_GEODESIC_TEST_TARGET := $(TEST_OUT_DIR)/test_geodesic_adaptive ASYMPTOTIC_TEST_TARGET := $(TEST_OUT_DIR)/test_asymptotic +ASYMPTOTIC_ENTRY_TEST_TARGET := $(TEST_OUT_DIR)/test_asymptotic_entry +ASYMPTOTIC_QUADRATIC_TEST_TARGET := $(TEST_OUT_DIR)/test_asymptotic_quadratic ASYMPTOTIC_SCHWARZSCHILD_TEST_TARGET := $(TEST_OUT_DIR)/test_asymptotic_schwarzschild TERMINATION_ORACLE_TEST_TARGET := $(TEST_OUT_DIR)/test_termination_oracle FRAME_TEST_TARGET := $(TEST_OUT_DIR)/test_frame @@ -217,9 +219,20 @@ $(ADAPTIVE_GEODESIC_TEST_TARGET): tests/test_geodesic_adaptive.c $(COMMON_SOURCE $(ASYMPTOTIC_TEST_TARGET): tests/test_asymptotic.c $(CORE_MINKOWSKI_SOURCES) $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@ +# The test includes asymptotic.c to cover its private floating-point kernel. +$(ASYMPTOTIC_QUADRATIC_TEST_TARGET): tests/test_asymptotic_quadratic.c $(CORE_MINKOWSKI_SOURCES) $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) + $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $(filter-out src/asymptotic.c,$^) $(LDLIBS) -o $@ + $(ASYMPTOTIC_SCHWARZSCHILD_TEST_TARGET): tests/test_asymptotic_schwarzschild.c $(COMMON_SOURCES) src/spacetime_schwarzschild.c $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -DSPACETIME_SCHWARZSCHILD -Isrc $^ $(LDLIBS) -o $@ +# Independent core regression for the backend-free numerical entry localizer. +# It links only the new module and the shared spacetime dispatch wrapper: no +# analytic backend, no geodesic integrator and no asymptotic.c are required, +# so it stays exercisable independently of the route integration. +$(ASYMPTOTIC_ENTRY_TEST_TARGET): tests/test_asymptotic_entry.c src/asymptotic_entry.c src/spacetime_common.c src/asymptotic_entry.h src/asymptotic.h src/geodesic.h src/spacetime.h src/observer.h | $(TEST_OUT_DIR) + $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc tests/test_asymptotic_entry.c src/asymptotic_entry.c src/spacetime_common.c $(LDLIBS) -o $@ + $(FRAME_TEST_TARGET): tests/test_frame.c $(CORE_MINKOWSKI_SOURCES) $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@ @@ -280,12 +293,14 @@ FAST_PSF_FFTW_TEST_DEP := FAST_PSF_FFTW_TEST_RUN := endif -test: $(CAMERA_TEST_TARGETS) $(TEST_TARGET) $(ADAPTIVE_GEODESIC_TEST_TARGET) $(ASYMPTOTIC_TEST_TARGET) $(ASYMPTOTIC_SCHWARZSCHILD_TEST_TARGET) $(TERMINATION_ORACLE_TEST_TARGET) $(FRAME_TEST_TARGET) $(SCHWARZSCHILD_TEST_TARGET) $(ALCUBIERRE_TEST_TARGET) $(OBSERVER_TRACK_TEST_TARGET) $(CATALOG_PREFETCH_TEST_TARGET) $(FAST_PSF_FFTW_TEST_DEP) $(TONE_MAP_TEST_TARGET) $(MOVIE_OUTPUT_TEST_TARGET) $(SENSOR_BLOOM_TEST_TARGET) +test: $(CAMERA_TEST_TARGETS) $(TEST_TARGET) $(ADAPTIVE_GEODESIC_TEST_TARGET) $(ASYMPTOTIC_TEST_TARGET) $(ASYMPTOTIC_ENTRY_TEST_TARGET) $(ASYMPTOTIC_QUADRATIC_TEST_TARGET) $(ASYMPTOTIC_SCHWARZSCHILD_TEST_TARGET) $(TERMINATION_ORACLE_TEST_TARGET) $(FRAME_TEST_TARGET) $(SCHWARZSCHILD_TEST_TARGET) $(ALCUBIERRE_TEST_TARGET) $(OBSERVER_TRACK_TEST_TARGET) $(CATALOG_PREFETCH_TEST_TARGET) $(FAST_PSF_FFTW_TEST_DEP) $(TONE_MAP_TEST_TARGET) $(MOVIE_OUTPUT_TEST_TARGET) $(SENSOR_BLOOM_TEST_TARGET) $(TEST_OUT_DIR)/test_observer_minkowski $(TEST_OUT_DIR)/test_observer_schwarzschild $(TEST_TARGET) $(ADAPTIVE_GEODESIC_TEST_TARGET) $(ASYMPTOTIC_TEST_TARGET) + $(ASYMPTOTIC_ENTRY_TEST_TARGET) + $(ASYMPTOTIC_QUADRATIC_TEST_TARGET) $(ASYMPTOTIC_SCHWARZSCHILD_TEST_TARGET) $(TERMINATION_ORACLE_TEST_TARGET) $(FRAME_TEST_TARGET) diff --git a/nr_spacetime_movie_renderer_design.md b/nr_spacetime_movie_renderer_design.md index fe8e868..f900ad3 100644 --- a/nr_spacetime_movie_renderer_design.md +++ b/nr_spacetime_movie_renderer_design.md @@ -1007,6 +1007,35 @@ residual、Chebyshev 表或解析主项；运行期不得建表。 - Minkowski entry quadratic 用稳定根公式（$q=-\tfrac12(b+\mathrm{copysign} (\sqrt\Delta,b))$，取最小正根）；$c=0$（相机在边界）时按 $b=\mathrm dF/ \mathrm ds$ 分类。worldtube sample 必须有限且 $R>0$，否则 `INVALID`。 + Alcubierre 等平直外区共用此解法；相对位置/速度、系数、判别式和根采用 + `long double` 中间量，最后才转换为 backend 的 double 状态。某些平台的 + `long double` 不提供额外有效位，因此精度提升不能替代入口验证。 + 送入根公式前，以系数最大绝对值的二进制指数共同缩放 $a,b,c$（最大值进入 + $[1/2,1)$），避免判别式乘积溢出；不改变传播参数及根，也不覆盖系数构造前 + 已发生的误差。非有限系数或缩放后非零系数落入 subnormal 范围，转入不确定 + 路径而非当作退化线性式。判别式 $b^2-4ac$ 和入口斜率采用普通 `long double` + 运算，不使用显式 FMA；精度/成本对照见 `benchmarks/quadratic_precision/`。 + 本机软件 `fmal` 成本明显，测试未显示足以抵偿成本的精度收益；不确定带、 + 几何验证及后备定位仍保留，不能把此样本结论当作全参数可靠性保证。 +- **入口快路径与后备定位分离**：解析相交只提供候选入口；采用前须在实际 backend + 坐标、实际入口时间重新检查 worldtube 残差。通过原有几何舍入容差的候选沿快路径 + 激活；未通过的候选沿同一外区 geodesic，从原相机事件重新求值并用确定在外/ + 严格在内的区间定位首次入口。公共 localizer 不依赖具体 metric backend， + 各已支持的外区提供轨迹求值与首次入口 bracket；不能把离散采样无交点当作 miss。 + miss 判定必须排除真实入口；擦边或求根病态造成的数值不确定性不能作为 + `ESCAPED` 的依据，允许保守地生成候选入口并进行后备验证。 + 浮点判别式等于零不构成精确擦边证明；重建最小点的单次正残差也不构成 + miss 证明（大时间/坐标的舍入会移动该最小点）。精确擦边可由独立几何证据 + 排除入口，例如固定半径段内某个不变坐标的精确距离已不小于球半径； + 无法获得这种证据且无法构造可信 bracket 时仍返回 `ENTRY_UNCONFIRMED`。 + 后备路径返回 bracket 收敛后的可表示内侧轨迹状态，不做径向位置投影，也不放宽 + 原 worldtube 检查容差。未能确认首次入口须报告 `INCOMPLETE/ENTRY_UNCONFIRMED`， + callback、历史及几何错误仍传播各自具体 reason，不能改写为 capture 或 escape。 + 定位在 pre-route 中完成，不在 slab sweep 中倒回相机时间；不重置相机 `L0`， + 内区 accepted-step/lookback 预算仍从最终激活事件起算。成功的后备入口由 + `entry_fallback_evaluations` 记录 localizer 的轨迹求值次数（不含构造 bracket + 的探测），不混入内区 RHS 计数；临时状态由 + 调用者独占，不引入共享可变缓存。此标记不改变 lens-map 二进制布局。 - **构造期验证优先**：`spacetime_create_*()` 成功即承诺该 source 已可安全光追。 每个 constructor 在安装 ops/context 后调用公共 `spacetime_source_finalize()`： 检查 ops/context 完整、`end_id` 唯一且非 `NONE`、exterior kind 受支持、 diff --git a/src/asymptotic.c b/src/asymptotic.c index 8e27f2e..76fa55d 100644 --- a/src/asymptotic.c +++ b/src/asymptotic.c @@ -1,5 +1,6 @@ #include "asymptotic.h" +#include "asymptotic_entry.h" #include "asymptotic_schwarzschild.h" #include @@ -188,73 +189,190 @@ int asymptotic_backend_from_canonical(const SpacetimeSource *source, return 0; } -/* Solve |d + q s|^2 = (R0 - rr s)^2 for the smallest s >= 0 with outside -> - * inside crossing. Returns 1 on entry (sets s), 0 on miss, -1 on error. */ -static int solve_entry_quadratic(const double d[3], const double q[3], - double R0, double rr, double *s_out) { - const double qq = dot3(q, q); - const double a = qq - rr * rr; - const double b = 2.0 * (dot3(d, q) + R0 * rr); - const double c = dot3(d, d) - R0 * R0; - if (c < 0.0) { +/* Long-double coefficients of the relative-distance quadratic + * F(sigma) = |d + q sigma|^2 - (R0 - rr sigma)^2 = a sigma^2 + b sigma + c, + * where d = x_cur - c_frame and q = w_frame + v_frame. Keeping them in long + * double preserves the cancellation-prone grazing entries; the caller hands + * the coefficients to both the root solve and the conservative fallback. */ +typedef struct { + long double a, b, c; +} EntryQuadratic; + +static EntryQuadratic entry_quadratic_coeffs( + const double x_cur[3], const double c_frame[3], const double w_frame[3], + const double v_frame[3], double R0, double rr) { + long double d[3], q[3]; + for (int i = 0; i < 3; ++i) { + d[i] = (long double)x_cur[i] - (long double)c_frame[i]; + q[i] = (long double)w_frame[i] + (long double)v_frame[i]; + } + long double qq = 0.0L, dot_dq = 0.0L, dot_dd = 0.0L; + for (int i = 0; i < 3; ++i) { + qq += q[i] * q[i]; + dot_dq += d[i] * q[i]; + dot_dd += d[i] * d[i]; + } + const long double R0_ld = (long double)R0; + const long double rr_ld = (long double)rr; + EntryQuadratic k; + k.a = qq - rr_ld * rr_ld; + k.b = 2.0L * (dot_dq + R0_ld * rr_ld); + k.c = dot_dd - R0_ld * R0_ld; + return k; +} + +typedef enum { + ENTRY_SOLVE_MISS = 0, /* proven no outside->inside crossing */ + ENTRY_SOLVE_ENTRY = 1, /* a first-crossing candidate parameter was found */ + ENTRY_SOLVE_UNCERTAIN = 2 /* discriminant/degeneracy at the resolution floor */ +} EntrySolveResult; + +/* For normalized coefficients, compute b*b - 4*a*c without overflowing. + * Near-cancellation is handled conservatively by the uncertainty band. */ +static long double entry_discriminant(const EntryQuadratic *k, + long double *scale) { + const long double four_a = 4.0L * k->a; + const long double ac = four_a * k->c; + *scale = k->b * k->b + fabsl(ac); + return k->b * k->b - ac; +} + +/* Solve the quadratic for the smallest sigma >= 0 with an outside -> inside + * crossing. The stable roots are evaluated in long double. A discriminant + * that is negative but within its own rounding bound, an exact double root, or + * a root that does not move inward is not a proof: it is reported as + * ENTRY_SOLVE_UNCERTAIN so the caller attempts a strict-inside bracket instead + * of fabricating a miss. Exact algebra (a == 0) stays the only linear branch; + * a tiny-but-nonzero `a` always goes through the discriminant, so a future + * entry at a huge parameter is never discarded by an approximate threshold. */ +static EntrySolveResult entry_solve(const EntryQuadratic *k, double *s_out) { + if (!isfinite(k->a) || !isfinite(k->b) || !isfinite(k->c)) + return ENTRY_SOLVE_UNCERTAIN; + /* A common power-of-two scale preserves roots and coefficient signs without + * adding division rounding. Normalize before squaring/products, retaining + * the original coefficients in the caller for geometric fallback. Do not + * silently turn an underflowed coefficient into a linear/boundary case. */ + EntryQuadratic scaled = *k; + const long double magnitude = fmaxl(fabsl(k->a), + fmaxl(fabsl(k->b), fabsl(k->c))); + if (magnitude > 0.0L) { + int exponent; + (void)frexpl(magnitude, &exponent); + scaled.a = scalbnl(k->a, -exponent); + scaled.b = scalbnl(k->b, -exponent); + scaled.c = scalbnl(k->c, -exponent); + if ((k->a != 0.0L && fabsl(scaled.a) < LDBL_MIN) || + (k->b != 0.0L && fabsl(scaled.b) < LDBL_MIN) || + (k->c != 0.0L && fabsl(scaled.c) < LDBL_MIN)) + return ENTRY_SOLVE_UNCERTAIN; + } + k = &scaled; + if (k->c < 0.0L) { /* Strictly inside; the lifecycle normally handles this as INSIDE. */ *s_out = 0.0; - return 1; + return ENTRY_SOLVE_ENTRY; } - if (c == 0.0) { + if (k->c == 0.0L) { /* On the boundary: classify by dF/ds = b. Past-inward enters at once; * outward/tangent rays may still re-enter later when the sphere shrinks * (a < 0), so do not declare a permanent miss on the zero root. */ - if (b < 0.0) { + if (k->b < 0.0L) { *s_out = 0.0; - return 1; + return ENTRY_SOLVE_ENTRY; } - if (b == 0.0) { - if (a < 0.0) { + if (k->b == 0.0L) { + if (k->a < 0.0L) { *s_out = 0.0; - return 1; + return ENTRY_SOLVE_ENTRY; } - return 0; + return ENTRY_SOLVE_MISS; } - if (a < 0.0) { - *s_out = -b / a; - return 1; + if (k->a < 0.0L) { + *s_out = (double)(-k->b / k->a); + return ENTRY_SOLVE_ENTRY; } - return 0; + return ENTRY_SOLVE_MISS; } - /* Compare the quadratic coefficient against the velocity-squared scale it - * is built from; mixing in R0^2 would let a large radius misclassify a - * genuinely quadratic entry as linear. */ - const double scale = qq + rr * rr; - if (fabs(a) <= 32.0 * DBL_EPSILON * scale) { - if (!isfinite(b) || b >= 0.0) - return 0; - const double s = -c / b; - if (s <= 0.0) - return 0; - *s_out = s; - return 1; + /* c > 0: the ray starts outside. */ + if (k->a == 0.0L) { + /* Exact linear branch only. b >= 0 never crosses for sigma > 0. */ + if (!(k->b < 0.0L)) + return ENTRY_SOLVE_MISS; + const long double s = -k->c / k->b; + if (!(s > 0.0L)) + return ENTRY_SOLVE_MISS; + *s_out = (double)s; + return ENTRY_SOLVE_ENTRY; } - const double disc = b * b - 4.0 * a * c; - if (!isfinite(disc) || disc <= 0.0) - return 0; - /* Numerically stable quadratic roots: q avoids cancellation in the root + /* Ordinary long-double products; near-zero differences need fallback. */ + long double disc_scale; + const long double disc = entry_discriminant(k, &disc_scale); + /* Conservative discriminant resolution: the long-double evaluation error + * plus the rounding the double inputs already carry through the frame + * rotation/translation into d and q. The double term dominates and keeps a + * near-tangent discriminant from being read as a proven miss. */ + const long double disc_err = + 64.0L * ((long double)DBL_EPSILON + (long double)LDBL_EPSILON) * + disc_scale; + if (!isfinite(disc)) + return ENTRY_SOLVE_UNCERTAIN; + if (disc < -disc_err) + return ENTRY_SOLVE_MISS; + if (fabsl(disc) <= disc_err) + return ENTRY_SOLVE_UNCERTAIN; + /* Numerically stable quadratic roots: qq2 avoids cancellation in the root * with the same sign as b, which is exactly the small entry root when the * camera sits just outside a large sphere. */ - const double root = sqrt(disc); - const double qq2 = -0.5 * (b + copysign(root, b)); - const double r1 = qq2 / a; - const double r2 = c / qq2; + const long double root = sqrtl(disc); + const long double qq2 = -0.5L * (k->b + copysignl(root, k->b)); + if (qq2 == 0.0L) + return ENTRY_SOLVE_UNCERTAIN; + const long double r1 = qq2 / k->a; + const long double r2 = k->c / qq2; /* The first outside->inside crossing is the smallest positive root. */ - double s = INFINITY; - if (r1 > 0.0) + long double s = INFINITY; + if (r1 > 0.0L) s = r1; - if (r2 > 0.0 && r2 < s) + if (r2 > 0.0L && r2 < s) s = r2; if (!(s < INFINITY)) - return 0; - *s_out = s; - return 1; + return (r1 > 0.0L || r2 > 0.0L) ? ENTRY_SOLVE_UNCERTAIN : ENTRY_SOLVE_MISS; + /* First crossing must move inward (dF/dsigma < 0). A nonnegative slope + * means the stable-root selection picked the exit root or the roots merged; + * that is uncertain, not a proof of a miss. */ + const long double slope = 2.0L * k->a * s + k->b; + if (!(slope < 0.0L)) + return ENTRY_SOLVE_UNCERTAIN; + *s_out = (double)s; + return ENTRY_SOLVE_ENTRY; +} + +static void minkowski_route_entry(const SpacetimeAsymptoticEnd *end, double t0, + const double x_frame[3], + const double w_frame[3], double s_entry, + SpacetimeEndId end_id, + AsymptoticRoute *route); + +/* Opaque evaluator context: repropagate constant-velocity motion from the + * ORIGINAL camera state to a total past parameter, never from a nearby root. */ +typedef struct { + const SpacetimeAsymptoticEnd *end; + double t0; + const double *x_frame; + const double *w_frame; + double log_alpha_p0; + SpacetimeEndId end_id; +} MinkowskiEntryContext; + +static AsymptoticStatus minkowski_entry_evaluate(void *opaque, double parameter, + AsymptoticRoute *state) { + const MinkowskiEntryContext *ctx = opaque; + *state = (AsymptoticRoute){0}; + minkowski_route_entry(ctx->end, ctx->t0, ctx->x_frame, ctx->w_frame, + parameter, ctx->end_id, state); + state->log_alpha_p0 = ctx->log_alpha_p0; + state->log_alpha_p0_camera = ctx->log_alpha_p0; + return ASYMPTOTIC_OK; } static void minkowski_route_escaped(const SpacetimeAsymptoticEnd *end, @@ -287,13 +405,248 @@ static void minkowski_route_entry(const SpacetimeAsymptoticEnd *end, double t0, route->Pi[i] = -w_backend[i]; } +/* Recover a first-entry bracket inside the CURRENT constant-motion segment + * when the closed-form candidate fails geometric validation or the + * discriminant is uncertain. The outside endpoint is the segment start (the + * previous segments produced no entry), and the inside endpoint is either the + * convex minimum or a modest, geometrically grown step past the candidate + * root. Both are confirmed with the actual worldtube callback. On success + * `route` carries the localized ENTRY and a nonzero fallback evaluation count. + * Any failure leaves `route->failure_reason` set and never reports an escape; + * one reconstructed probe cannot certify a miss. */ +static AsymptoticStatus minkowski_fallback_entry( + const SpacetimeSource *source, const SpacetimeAsymptoticEnd *end, double t0, + const double x_frame[3], const double w_frame[3], double log_alpha_p0, + SpacetimeEndId end_id, double s_base, double s_segment, + const EntryQuadratic *k, EntrySolveResult solve, double candidate_sigma, + AsymptoticRoute *route) { + MinkowskiEntryContext ctx = {.end = end, + .t0 = t0, + .x_frame = x_frame, + .w_frame = w_frame, + .log_alpha_p0 = log_alpha_p0, + .end_id = end_id}; + + /* Outside endpoint: the segment start must still be outside. */ + AsymptoticRoute outside_state; + AsymptoticStatus status = + minkowski_entry_evaluate(&ctx, s_base, &outside_state); + if (status != ASYMPTOTIC_OK) { + route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED; + route->end_id = end_id; + return status; + } + RayReason why = RAY_REASON_NONE; + double f_out = 0.0, tol_out = 0.0; + status = asymptotic_entry_geometry(source, end_id, outside_state.activate_t, + outside_state.x, &f_out, &tol_out, &why); + if (status != ASYMPTOTIC_OK) { + route->failure_reason = why; + route->end_id = end_id; + return status; + } + if (!(f_out >= 0.0)) { + /* Already inside at the segment start: an earlier segment missed the + * crossing. Do not fabricate a bracket from it. */ + route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED; + route->end_id = end_id; + return ASYMPTOTIC_INVALID; + } + + int have_inside = 0; + double inside_sigma = 0.0; + + if (k->a > 0.0L) { + /* Convex: only the quadratic minimum can be strictly inside, and F is + * monotonically decreasing from the segment start to that minimum, so the + * bracket still straddles the first crossing. */ + const long double sigma_min_ld = -k->b / (2.0L * k->a); + if (!(sigma_min_ld > 0.0L)) { + route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED; + route->end_id = end_id; + return ASYMPTOTIC_INVALID; + } + double probe = (double)sigma_min_ld; + if (probe > s_segment) + probe = s_segment; + if (!(probe > 0.0)) { + route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED; + route->end_id = end_id; + return ASYMPTOTIC_INVALID; + } + AsymptoticRoute probe_state; + status = minkowski_entry_evaluate(&ctx, s_base + probe, &probe_state); + if (status != ASYMPTOTIC_OK) { + route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED; + route->end_id = end_id; + return status; + } + double f_probe = 0.0, tol_probe = 0.0; + status = asymptotic_entry_geometry(source, end_id, probe_state.activate_t, + probe_state.x, &f_probe, &tol_probe, + &why); + if (status != ASYMPTOTIC_OK) { + route->failure_reason = why; + route->end_id = end_id; + return status; + } + if (f_probe < 0.0) { + have_inside = 1; + inside_sigma = probe; + } else { + /* One reconstructed probe is not a miss proof: coordinate-time rounding + * can shift the minimum and hide an inside point at a nearby parameter. */ + route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED; + route->end_id = end_id; + return ASYMPTOTIC_INVALID; + } + } else { + /* Monotone (a == 0 linear) or concave after the first crossing: step + * modestly past the candidate root and grow geometrically, staying inside + * the declared segment and the positive-radius domain. */ + if (solve != ENTRY_SOLVE_ENTRY || !(candidate_sigma >= 0.0)) { + route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED; + route->end_id = end_id; + return ASYMPTOTIC_INVALID; + } + double sigma = candidate_sigma; + if (sigma > s_segment) + sigma = s_segment; + AsymptoticRoute probe_state; + status = minkowski_entry_evaluate(&ctx, s_base + sigma, &probe_state); + if (status != ASYMPTOTIC_OK) { + route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED; + route->end_id = end_id; + return status; + } + double f_sigma = 0.0, tol_sigma = 0.0; + status = asymptotic_entry_geometry(source, end_id, probe_state.activate_t, + probe_state.x, &f_sigma, &tol_sigma, + &why); + if (status != ASYMPTOTIC_OK) { + route->failure_reason = why; + route->end_id = end_id; + return status; + } + if (f_sigma < 0.0) { + have_inside = 1; + inside_sigma = sigma; + } else { + const long double slope = + 2.0L * k->a * (long double)sigma + k->b; /* < 0 for an entry */ + double delta = 0.0; + if (slope < 0.0L) + delta = 2.0 * fabs(f_sigma) / fabs((double)slope); + const double ulp_term = + 16.0 * DBL_EPSILON * fmax(1.0, fabs(s_base + sigma)); + if (!(delta > ulp_term)) + delta = ulp_term; + for (int attempt = 0; attempt < 64 && !have_inside; ++attempt) { + const double probe = sigma + delta; + if (!(probe > sigma) || probe > s_segment) + break; + AsymptoticRoute grown_state; + status = minkowski_entry_evaluate(&ctx, s_base + probe, &grown_state); + if (status != ASYMPTOTIC_OK) { + route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED; + route->end_id = end_id; + return status; + } + double f_grown = 0.0, tol_grown = 0.0; + status = asymptotic_entry_geometry(source, end_id, + grown_state.activate_t, + grown_state.x, &f_grown, &tol_grown, + &why); + if (status != ASYMPTOTIC_OK) { + route->failure_reason = why; + route->end_id = end_id; + return status; + } + if (f_grown < 0.0) { + have_inside = 1; + inside_sigma = probe; + } else { + delta *= 2.0; + } + } + if (!have_inside) { + route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED; + route->end_id = end_id; + return ASYMPTOTIC_INVALID; + } + } + } + + /* Bracket confirmed: hand it to the common, exterior-independent driver. */ + AsymptoticRoute out; + unsigned int evaluations = 0; + RayReason localize_reason = RAY_REASON_NONE; + status = asymptotic_entry_localize( + source, end_id, minkowski_entry_evaluate, &ctx, s_base, + s_base + inside_sigma, &out, &evaluations, &localize_reason); + if (status == ASYMPTOTIC_OK) { + *route = out; + route->failure_reason = RAY_REASON_NONE; + route->entry_fallback_evaluations = evaluations; + return ASYMPTOTIC_OK; + } + if (status == ASYMPTOTIC_TIME_RANGE_EXHAUSTED || + status == ASYMPTOTIC_UNSUPPORTED) { + route->failure_reason = localize_reason; + route->end_id = end_id; + return status; + } + route->failure_reason = + (localize_reason == RAY_REASON_ESCAPE_LOCALIZATION_FAILED || + localize_reason == RAY_REASON_NONE || + localize_reason == RAY_REASON_PROTOCOL_ERROR) + ? RAY_REASON_ENTRY_UNCONFIRMED + : localize_reason; + route->end_id = end_id; + return ASYMPTOTIC_INVALID; +} + +/* A fixed coordinate outside the sphere is an independent algebraic miss + * certificate, including exact tangency. Restrict this cheap certificate to + * identity axes and zero origin so frame reconstruction cannot change the + * original camera component. Sterbenz's lemma certifies the subtraction when + * both nonzero operands have the same sign and are within a factor of two. */ +static int minkowski_coordinate_miss( + const SpacetimeAsymptoticEnd *end, const double x_cur[3], + const double w_frame[3], const SpacetimeEscapeWorldtubeSample *sample) { + if (sample->radius_rate != 0.0) + return 0; + for (int i = 0; i < 3; ++i) + if (end->frame_origin[i] != 0.0) + return 0; + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + if (end->frame_axes[i][j] != (i == j ? 1.0 : 0.0)) + return 0; + for (int i = 0; i < 3; ++i) { + if (w_frame[i] != 0.0 || sample->velocity[i] != 0.0) + continue; + const double x = x_cur[i] + end->frame_origin[i]; + const double c = sample->center[i]; + const int exact = x == 0.0 || c == 0.0 || + (signbit(x) == signbit(c) && fabs(x) * 0.5 <= fabs(c) && + fabs(c) * 0.5 <= fabs(x)); + const double d = x - c; + if (isfinite(x) && isfinite(d) && exact && fabs(d) >= sample->radius) + return 1; + } + return 0; +} + static AsymptoticStatus minkowski_preroute( const SpacetimeSource *source, const SpacetimeAsymptoticEnd *end, double t0, - const double x_frame[3], const double w_frame[3], SpacetimeEndId end_id, - AsymptoticRoute *route) { + const double x_frame[3], const double w_frame[3], double log_alpha_p0, + SpacetimeEndId end_id, AsymptoticRoute *route) { /* Walk constant-velocity motion segments. A quadratic root is only valid * inside the current segment and while the radius stays positive; otherwise * advance to the next segment boundary and re-sample. */ + route->end_id = end_id; + route->failure_reason = RAY_REASON_NONE; double s = 0.0; for (int segment = 0; segment < 1000000; ++segment) { const double t = t0 - s; @@ -313,32 +666,68 @@ static AsymptoticStatus minkowski_preroute( * path below. */ return ASYMPTOTIC_UNSUPPORTED; } - double c_frame[3], v_frame[3], d[3], q[3]; + double c_frame[3], v_frame[3]; backend_position_to_frame(end, sample.center, c_frame); backend_vector_to_frame(end, sample.velocity, v_frame); - for (int i = 0; i < 3; ++i) { - d[i] = x_cur[i] - c_frame[i]; - q[i] = w_frame[i] + v_frame[i]; - } const double boundary = spacetime_escape_worldtube_next_segment( source, end->end_id, t); const double s_segment = isfinite(boundary) ? (t - boundary) : INFINITY; if (!(s_segment >= 0.0)) return ASYMPTOTIC_INVALID; - double sigma; - const int hit = solve_entry_quadratic(d, q, sample.radius, - sample.radius_rate, &sigma); + const EntryQuadratic k = entry_quadratic_coeffs( + x_cur, c_frame, w_frame, v_frame, sample.radius, sample.radius_rate); + double sigma = 0.0; + EntrySolveResult solve = entry_solve(&k, &sigma); + if (solve == ENTRY_SOLVE_UNCERTAIN && + minkowski_coordinate_miss(end, x_cur, w_frame, &sample)) + solve = ENTRY_SOLVE_MISS; /* The backend constructor guarantees R > 0 throughout every segment, so * a root inside the segment is a real entry. A root past the segment * boundary is not adopted here; the next segment is sampled instead. * The cheap R > 0 test at the root guards against a backend that * bypasses its constructor. */ - if (hit && sigma >= 0.0 && sigma <= s_segment) { + if (solve == ENTRY_SOLVE_ENTRY && sigma >= 0.0 && sigma <= s_segment) { if (sample.radius - sample.radius_rate * sigma <= 0.0) return ASYMPTOTIC_INVALID; + AsymptoticRoute candidate = {.kind = ASYMPTOTIC_ROUTE_INVALID}; minkowski_route_entry(end, t0, x_frame, w_frame, s + sigma, end_id, - route); - return ASYMPTOTIC_OK; + &candidate); + candidate.log_alpha_p0 = log_alpha_p0; + candidate.log_alpha_p0_camera = log_alpha_p0; + candidate.failure_reason = RAY_REASON_NONE; + int valid = 0; + RayReason why = RAY_REASON_NONE; + status = asymptotic_entry_validate(source, end_id, &candidate, &valid, + &why); + if (status == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) { + route->failure_reason = why; + return status; + } + if (status != ASYMPTOTIC_OK) { + route->failure_reason = why; + route->end_id = end_id; + return status; + } + if (valid) { + /* Fast path: the reconstructed entry is already on the boundary + * within the geometric ULP band. */ + *route = candidate; + return ASYMPTOTIC_OK; + } + /* The closed-form candidate lands off the reconstructed boundary: + * fall back to the common numerical localizer inside this segment. */ + status = minkowski_fallback_entry(source, end, t0, x_frame, w_frame, + log_alpha_p0, end_id, s, s_segment, &k, + solve, sigma, route); + return status; + } else if (solve == ENTRY_SOLVE_UNCERTAIN) { + /* Near-tangent/degenerate discriminant: not a proof of a miss. Attempt + * a strict-inside bracket; a positive reconstructed probe cannot prove + * that the continuous trajectory misses. */ + status = minkowski_fallback_entry(source, end, t0, x_frame, w_frame, + log_alpha_p0, end_id, s, s_segment, &k, + solve, 0.0, route); + return status; } if (!isfinite(s_segment)) { /* Open final segment with no entry: a genuine miss. */ @@ -355,6 +744,127 @@ static AsymptoticStatus minkowski_preroute( return ASYMPTOTIC_INVALID; } +/* Evaluate the exact Schwarzschild inward orbit from the original camera at + * the radius parameter p = -rho. Recomputed per call (rotate + integral), not + * projected from an earlier state, and using no backend metric. */ +typedef struct { + const SpacetimeAsymptoticEnd *end; + const SchwarzschildCanonical *camera; + double log_alpha_p0_camera; + SpacetimeEndId end_id; +} SchwarzschildEntryContext; + +static AsymptoticStatus schwarzschild_entry_evaluate(void *opaque, + double parameter, + AsymptoticRoute *state) { + const SchwarzschildEntryContext *ctx = opaque; + const double rho = -parameter; + double x[3], Pi[3], log_alpha_p0 = 0.0, activate_t = 0.0; + if (asymptotic_schwarzschild_inward_state_at_radius( + ctx->end, ctx->camera, rho, x, Pi, &log_alpha_p0, &activate_t)) + return ASYMPTOTIC_INVALID; + *state = (AsymptoticRoute){0}; + state->kind = ASYMPTOTIC_ROUTE_ENTRY; + state->end_id = ctx->end_id; + state->activate_t = activate_t; + for (int i = 0; i < 3; ++i) { + state->x[i] = x[i]; + state->Pi[i] = Pi[i]; + } + state->log_alpha_p0 = log_alpha_p0; + state->log_alpha_p0_camera = ctx->log_alpha_p0_camera; + state->failure_reason = RAY_REASON_NONE; + state->entry_fallback_evaluations = 0; + return ASYMPTOTIC_OK; +} + +/* Bracket the first inward radius crossing when the closed-form entry state + * lands off the reconstructed worldtube boundary. The parameter p = -rho + * increases inward from the original camera radius; the inside end is nudged + * just below the worldtube radius, and expanded inward only while staying + * above the turning radius and the rho > 2 state domain. Failure is an + * explicit unconfirmed entry, never an escape. */ +static AsymptoticStatus schwarzschild_fallback_entry( + const SpacetimeSource *source, const SpacetimeAsymptoticEnd *end, + const SchwarzschildCanonical *camera, double worldtube_radius, + double log_alpha_p0_camera, AsymptoticRoute *route) { + SchwarzschildEntryContext ctx = {.end = end, + .camera = camera, + .log_alpha_p0_camera = log_alpha_p0_camera, + .end_id = end->end_id}; + double floor = 2.0 + 1e-12 * fmax(1.0, worldtube_radius); + const double rho_turn = asymptotic_schwarzschild_turning_rho(camera->beta); + if (isfinite(rho_turn) && rho_turn > floor) + floor = rho_turn; + if (!(floor < worldtube_radius)) { + route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED; + route->end_id = end->end_id; + return ASYMPTOTIC_INVALID; + } + double rho_inside = nextafter(worldtube_radius, -INFINITY); + if (!(rho_inside > floor)) + rho_inside = 0.5 * (worldtube_radius + floor); + double decrement = 0.0; + int found = 0; + for (int attempt = 0; attempt < 64; ++attempt) { + if (!(rho_inside > floor)) + break; + double x[3], Pi[3], log_alpha_p0 = 0.0, activate_t = 0.0; + if (asymptotic_schwarzschild_inward_state_at_radius( + end, camera, rho_inside, x, Pi, &log_alpha_p0, &activate_t)) + break; + double F = 0.0, tol = 0.0; + RayReason why = RAY_REASON_NONE; + const AsymptoticStatus st = asymptotic_entry_geometry( + source, end->end_id, activate_t, x, &F, &tol, &why); + if (st != ASYMPTOTIC_OK) { + route->failure_reason = why; + route->end_id = end->end_id; + return st; + } + if (F < 0.0) { + found = 1; + break; + } + if (decrement == 0.0) + decrement = (worldtube_radius - rho_inside) * 2.0; + else + decrement *= 2.0; + rho_inside = worldtube_radius - decrement; + } + if (!found) { + route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED; + route->end_id = end->end_id; + return ASYMPTOTIC_INVALID; + } + AsymptoticRoute out; + unsigned int evaluations = 0; + RayReason localize_reason = RAY_REASON_NONE; + const AsymptoticStatus status = asymptotic_entry_localize( + source, end->end_id, schwarzschild_entry_evaluate, &ctx, -camera->rho, + -rho_inside, &out, &evaluations, &localize_reason); + if (status == ASYMPTOTIC_OK) { + *route = out; + route->failure_reason = RAY_REASON_NONE; + route->entry_fallback_evaluations = evaluations; + return ASYMPTOTIC_OK; + } + if (status == ASYMPTOTIC_TIME_RANGE_EXHAUSTED || + status == ASYMPTOTIC_UNSUPPORTED) { + route->failure_reason = localize_reason; + route->end_id = end->end_id; + return status; + } + route->failure_reason = + (localize_reason == RAY_REASON_ESCAPE_LOCALIZATION_FAILED || + localize_reason == RAY_REASON_NONE || + localize_reason == RAY_REASON_PROTOCOL_ERROR) + ? RAY_REASON_ENTRY_UNCONFIRMED + : localize_reason; + route->end_id = end->end_id; + return ASYMPTOTIC_INVALID; +} + static AsymptoticStatus schwarzschild_route( const SpacetimeSource *source, const SpacetimeAsymptoticEnd *end, const MetricData *metric, const GeodesicRayState *state, @@ -397,14 +907,38 @@ static AsymptoticStatus schwarzschild_route( } route->end_id = end->end_id; if (kind == SCH_ROUTE_ENTRY) { - route->kind = ASYMPTOTIC_ROUTE_ENTRY; - route->activate_t = activate_t; + AsymptoticRoute candidate = {.kind = ASYMPTOTIC_ROUTE_ENTRY, + .end_id = end->end_id, + .activate_t = activate_t, + .log_alpha_p0 = log_alpha_p0, + .log_alpha_p0_camera = state->log_alpha_p0, + .failure_reason = RAY_REASON_NONE}; for (int i = 0; i < 3; ++i) { - route->x[i] = x[i]; - route->Pi[i] = Pi[i]; + candidate.x[i] = x[i]; + candidate.Pi[i] = Pi[i]; } - route->log_alpha_p0 = log_alpha_p0; - return ASYMPTOTIC_OK; + int valid = 0; + RayReason why = RAY_REASON_NONE; + const AsymptoticStatus validate_status = asymptotic_entry_validate( + source, end->end_id, &candidate, &valid, &why); + if (validate_status == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) { + route->failure_reason = why; + route->kind = ASYMPTOTIC_ROUTE_TIME_RANGE_EXHAUSTED; + return validate_status; + } + if (validate_status != ASYMPTOTIC_OK) { + route->failure_reason = why; + return validate_status; + } + if (valid) { + *route = candidate; + return ASYMPTOTIC_OK; + } + /* The closed-form entry is off the reconstructed boundary: localize it in + * the radius parameter against the same exact inward transfer. */ + return schwarzschild_fallback_entry(source, end, &camera, + sample.radius / end->mass, + state->log_alpha_p0, route); } route->kind = ASYMPTOTIC_ROUTE_ESCAPED; for (int i = 0; i < 3; ++i) @@ -521,14 +1055,20 @@ AsymptoticStatus asymptotic_route_camera(const SpacetimeSource *source, AsymptoticRoute candidate = {.kind = ASYMPTOTIC_ROUTE_INVALID}; const AsymptoticStatus status = minkowski_preroute( source, &end, state.coordinate_time, canonical.x, canonical.w, - end.end_id, &candidate); + state.log_alpha_p0, end.end_id, &candidate); if (status == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) { route->kind = ASYMPTOTIC_ROUTE_TIME_RANGE_EXHAUSTED; route->end_id = end.end_id; + route->failure_reason = candidate.failure_reason; return ASYMPTOTIC_TIME_RANGE_EXHAUSTED; } - if (status != ASYMPTOTIC_OK) + if (status != ASYMPTOTIC_OK) { + /* Propagate the specific validation/fallback failure the segment walk + * refused with, not a generic preroute error. */ + route->failure_reason = candidate.failure_reason; + route->end_id = candidate.end_id; return status; + } if (candidate.kind == ASYMPTOTIC_ROUTE_ENTRY) { const double s = state.coordinate_time - candidate.activate_t; if (!have_entry || s < best_s) { diff --git a/src/asymptotic.h b/src/asymptotic.h index 8ab6e4f..fe1d43b 100644 --- a/src/asymptotic.h +++ b/src/asymptotic.h @@ -48,6 +48,16 @@ typedef struct { /* Terminal infinity endpoint for ESCAPED. */ double n_infinity[3]; double frequency_ratio; + /* Diagnostic reason for an ASYMPTOTIC_INVALID return: set specifically by the + * validation/fallback failure that refused the route, so the lifecycle can + * report the concrete cause instead of a generic preroute failure. NONE on + * success. */ + RayReason failure_reason; + /* Nonzero when this route was produced by the generic bracketed first-entry + * localizer rather than a closed-form/fast entry solve. It records the + * number of evaluator calls the localizer spent; zero on the fast path. This + * is private in-memory provenance only and is not serialized. */ + unsigned int entry_fallback_evaluations; } AsymptoticRoute; /* Pre-route one camera ray against every declared end's worldtube. */ diff --git a/src/asymptotic_entry.c b/src/asymptotic_entry.c new file mode 100644 index 0000000..8600b64 --- /dev/null +++ b/src/asymptotic_entry.c @@ -0,0 +1,326 @@ +#include "asymptotic_entry.h" + +#include +#include +#include + +/* See asymptotic_entry.h for the contract. This module deliberately keeps no + * global mutable state: every cache/scratch value lives on the stack of the + * calling trace, so it stays thread-safe under the coarse-grained OpenMP ray + * parallelism of the renderer. */ + +static void set_failure(RayReason *failure, RayReason reason) { + if (failure != NULL) + *failure = reason; +} + +/* Left-to-right double accumulation, matching the geodesic event layer's + * worldtube F. */ +static double dot3(const double a[3], const double b[3]) { + return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; +} + +static int entry_state_finite(const AsymptoticRoute *state) { + if (!isfinite(state->activate_t) || !isfinite(state->log_alpha_p0) || + !isfinite(state->log_alpha_p0_camera)) + return 0; + for (int i = 0; i < 3; ++i) + if (!isfinite(state->x[i]) || !isfinite(state->Pi[i])) + return 0; + return 1; +} + +AsymptoticStatus asymptotic_entry_geometry(const SpacetimeSource *source, + SpacetimeEndId end_id, double t, + const double x[3], double *F, + double *tol, RayReason *failure) { + set_failure(failure, RAY_REASON_PROTOCOL_ERROR); + if (source == NULL || x == NULL || F == NULL || tol == NULL) { + set_failure(failure, RAY_REASON_INVALID_ARGUMENT); + return ASYMPTOTIC_INVALID; + } + if (!isfinite(t)) { + set_failure(failure, RAY_REASON_WORLDTUBE_GEOMETRY_INVALID); + return ASYMPTOTIC_INVALID; + } + SpacetimeEscapeWorldtubeSample sample; + if (spacetime_escape_worldtube_sample(source, end_id, t, &sample)) { + set_failure(failure, RAY_REASON_WORLDTUBE_SAMPLE_FAILED); + return ASYMPTOTIC_INVALID; + } + if (!sample.valid) { + /* The backend cannot describe the worldtube at this time; this is history + * exhaustion, never a miss. */ + set_failure(failure, RAY_REASON_TIME_RANGE_EXHAUSTED); + return ASYMPTOTIC_TIME_RANGE_EXHAUSTED; + } + if (!(sample.radius > 0.0) || !isfinite(sample.radius) || + !isfinite(sample.radius_rate)) { + set_failure(failure, RAY_REASON_WORLDTUBE_GEOMETRY_INVALID); + return ASYMPTOTIC_INVALID; + } + for (int i = 0; i < 3; ++i) { + if (!isfinite(sample.center[i]) || !isfinite(sample.velocity[i])) { + set_failure(failure, RAY_REASON_WORLDTUBE_GEOMETRY_INVALID); + return ASYMPTOTIC_INVALID; + } + } + double d[3]; + for (int i = 0; i < 3; ++i) { + d[i] = x[i] - sample.center[i]; + if (!isfinite(d[i])) { + set_failure(failure, RAY_REASON_WORLDTUBE_GEOMETRY_INVALID); + return ASYMPTOTIC_INVALID; + } + } + /* Exact same grouping as the geodesic event layer: F uses + * dot3(d,d) - radius^2, while the geometric ULP band accumulates d2 with an + * explicit left-to-right loop. Keeping both identical means a candidate that + * passes this validator at the tolerance threshold is grouped exactly like + * the geodesic's own f_before. */ + const double value = dot3(d, d) - sample.radius * sample.radius; + const double r2 = sample.radius * sample.radius; + double d2 = 0.0; + for (int i = 0; i < 3; ++i) + d2 += d[i] * d[i]; + const double geometry_tol = 128.0 * DBL_EPSILON * fmax(r2, d2); + /* Overflow to inf and NaN propagation both land here. */ + if (!isfinite(value) || !isfinite(geometry_tol)) { + set_failure(failure, RAY_REASON_WORLDTUBE_GEOMETRY_INVALID); + return ASYMPTOTIC_INVALID; + } + *F = value; + *tol = geometry_tol; + set_failure(failure, RAY_REASON_NONE); + return ASYMPTOTIC_OK; +} + +AsymptoticStatus asymptotic_entry_validate(const SpacetimeSource *source, + SpacetimeEndId end_id, + const AsymptoticRoute *candidate, + int *valid, RayReason *failure) { + set_failure(failure, RAY_REASON_PROTOCOL_ERROR); + if (source == NULL || candidate == NULL || valid == NULL) { + set_failure(failure, RAY_REASON_INVALID_ARGUMENT); + return ASYMPTOTIC_INVALID; + } + *valid = 0; + if (candidate->kind != ASYMPTOTIC_ROUTE_ENTRY) { + set_failure(failure, RAY_REASON_PROTOCOL_ERROR); + return ASYMPTOTIC_INVALID; + } + if (!entry_state_finite(candidate)) { + set_failure(failure, RAY_REASON_PROTOCOL_ERROR); + return ASYMPTOTIC_INVALID; + } + double value = 0.0, tol = 0.0; + const AsymptoticStatus status = asymptotic_entry_geometry( + source, end_id, candidate->activate_t, candidate->x, &value, &tol, + failure); + if (status != ASYMPTOTIC_OK) + return status; + /* Boundary-near means within the geometric ULP band on either side. Outside + * and arbitrary-deep-inside are both non-candidates, not protocol errors. */ + *valid = fabs(value) <= tol; + set_failure(failure, RAY_REASON_NONE); + return ASYMPTOTIC_OK; +} + +/* Map an evaluator's own failure to a diagnostic reason. The worldtube + * callback/history/geometry reasons come from asymptotic_entry_geometry; this + * only covers the case where the evaluator itself refuses to produce a state. */ +static AsymptoticStatus entry_evaluator_failure(AsymptoticStatus status, + RayReason *failure) { + switch (status) { + case ASYMPTOTIC_TIME_RANGE_EXHAUSTED: + set_failure(failure, RAY_REASON_TIME_RANGE_EXHAUSTED); + break; + case ASYMPTOTIC_UNSUPPORTED: + set_failure(failure, RAY_REASON_UNSUPPORTED); + break; + default: + /* The evaluator refused to produce a state, so the first entry is not + * confirmed. */ + set_failure(failure, RAY_REASON_ENTRY_UNCONFIRMED); + break; + } + return status; +} + +AsymptoticStatus asymptotic_entry_localize( + const SpacetimeSource *source, SpacetimeEndId end_id, + AsymptoticEntryEvaluator evaluate, void *context, + double outside_parameter, double inside_parameter, AsymptoticRoute *out, + unsigned int *evaluations, RayReason *failure) { + set_failure(failure, RAY_REASON_ENTRY_UNCONFIRMED); + if (evaluations != NULL) + *evaluations = 0; + if (source == NULL || evaluate == NULL || out == NULL) { + set_failure(failure, RAY_REASON_INVALID_ARGUMENT); + return ASYMPTOTIC_INVALID; + } + if (!isfinite(outside_parameter) || !isfinite(inside_parameter) || + !(outside_parameter < inside_parameter)) { + set_failure(failure, RAY_REASON_INVALID_ARGUMENT); + return ASYMPTOTIC_INVALID; + } + + unsigned int count = 0; + AsymptoticRoute lo_state, hi_state; + double f_lo = 0.0, f_hi = 0.0, tol_lo = 0.0, tol_hi = 0.0; + + AsymptoticStatus status = evaluate(context, outside_parameter, &lo_state); + ++count; + if (status != ASYMPTOTIC_OK) { + if (evaluations != NULL) + *evaluations = count; + return entry_evaluator_failure(status, failure); + } + if (!entry_state_finite(&lo_state)) { + if (evaluations != NULL) + *evaluations = count; + set_failure(failure, RAY_REASON_ENTRY_UNCONFIRMED); + return ASYMPTOTIC_INVALID; + } + status = asymptotic_entry_geometry(source, end_id, lo_state.activate_t, + lo_state.x, &f_lo, &tol_lo, failure); + if (status != ASYMPTOTIC_OK) { + if (evaluations != NULL) + *evaluations = count; + return status; + } + + status = evaluate(context, inside_parameter, &hi_state); + ++count; + if (status != ASYMPTOTIC_OK) { + if (evaluations != NULL) + *evaluations = count; + return entry_evaluator_failure(status, failure); + } + if (!entry_state_finite(&hi_state)) { + if (evaluations != NULL) + *evaluations = count; + set_failure(failure, RAY_REASON_ENTRY_UNCONFIRMED); + return ASYMPTOTIC_INVALID; + } + status = asymptotic_entry_geometry(source, end_id, hi_state.activate_t, + hi_state.x, &f_hi, &tol_hi, failure); + if (status != ASYMPTOTIC_OK) { + if (evaluations != NULL) + *evaluations = count; + return status; + } + + if (!(f_lo >= 0.0) || !(f_hi < 0.0)) { + /* The caller owns the first-entry bracket; a bracket that does not straddle + * the boundary is an unconfirmed entry, never an escape. */ + set_failure(failure, RAY_REASON_ENTRY_UNCONFIRMED); + if (evaluations != NULL) + *evaluations = count; + return ASYMPTOTIC_INVALID; + } + + /* Parameter increases backward, so the inside end must not be later in + * coordinate time than the outside end. */ + if (!(hi_state.activate_t <= lo_state.activate_t)) { + set_failure(failure, RAY_REASON_ENTRY_UNCONFIRMED); + if (evaluations != NULL) + *evaluations = count; + return ASYMPTOTIC_INVALID; + } + + if (f_lo == 0.0) { + /* Exact boundary at the outside end with a strictly inside other end: this + * is the legitimate inward first-entry state already on the worldtube, so + * no bisection is needed. */ + *out = lo_state; + out->kind = ASYMPTOTIC_ROUTE_ENTRY; + out->end_id = end_id; + if (evaluations != NULL) + *evaluations = count; + set_failure(failure, RAY_REASON_NONE); + return ASYMPTOTIC_OK; + } + + double lo = outside_parameter; + double hi = inside_parameter; + unsigned int iteration = 0; + for (; iteration < ASYMPTOTIC_ENTRY_BISECTION_LIMIT; ++iteration) { + const double span = hi - lo; + const double mid = isfinite(span) ? lo + span * 0.5 + : lo * 0.5 + hi * 0.5; + if (!(mid > lo && mid < hi)) + break; /* Parameter midpoint cannot be represented; bracket is adjacent. */ + AsymptoticRoute mid_state; + status = evaluate(context, mid, &mid_state); + ++count; + if (status != ASYMPTOTIC_OK) { + if (evaluations != NULL) + *evaluations = count; + return entry_evaluator_failure(status, failure); + } + if (!entry_state_finite(&mid_state) || + !(hi_state.activate_t <= mid_state.activate_t && + mid_state.activate_t <= lo_state.activate_t)) { + if (evaluations != NULL) + *evaluations = count; + set_failure(failure, RAY_REASON_ENTRY_UNCONFIRMED); + return ASYMPTOTIC_INVALID; + } + double f_mid = 0.0, tol_mid = 0.0; + status = asymptotic_entry_geometry(source, end_id, mid_state.activate_t, + mid_state.x, &f_mid, &tol_mid, failure); + if (status != ASYMPTOTIC_OK) { + if (evaluations != NULL) + *evaluations = count; + return status; + } + if (f_mid >= 0.0) { + lo = mid; + lo_state = mid_state; + } else { + hi = mid; + hi_state = mid_state; + f_hi = f_mid; + tol_hi = tol_mid; + } + } + + if (iteration >= ASYMPTOTIC_ENTRY_BISECTION_LIMIT) { + /* The parameter interval never contracted to adjacent doubles within the + * implementation guard; do not fabricate an entry. */ + set_failure(failure, RAY_REASON_ENTRY_UNCONFIRMED); + if (evaluations != NULL) + *evaluations = count; + return ASYMPTOTIC_INVALID; + } + + /* Final consistency: keep the strict-inside adjacent endpoint. Its negative + * residual need not fit the fast-path band at coarse coordinate resolution; + * the bracket, rather than a radial displacement, establishes the entry. */ + if (!(f_hi < 0.0) || + !(hi_state.activate_t <= lo_state.activate_t) || + !isfinite(hi_state.activate_t) || !isfinite(hi_state.log_alpha_p0) || + !isfinite(hi_state.log_alpha_p0_camera)) { + set_failure(failure, RAY_REASON_ENTRY_UNCONFIRMED); + if (evaluations != NULL) + *evaluations = count; + return ASYMPTOTIC_INVALID; + } + for (int i = 0; i < 3; ++i) { + if (!isfinite(hi_state.x[i]) || !isfinite(hi_state.Pi[i])) { + set_failure(failure, RAY_REASON_ENTRY_UNCONFIRMED); + if (evaluations != NULL) + *evaluations = count; + return ASYMPTOTIC_INVALID; + } + } + + *out = hi_state; + out->kind = ASYMPTOTIC_ROUTE_ENTRY; + out->end_id = end_id; + if (evaluations != NULL) + *evaluations = count; + set_failure(failure, RAY_REASON_NONE); + return ASYMPTOTIC_OK; +} diff --git a/src/asymptotic_entry.h b/src/asymptotic_entry.h new file mode 100644 index 0000000..b675e3e --- /dev/null +++ b/src/asymptotic_entry.h @@ -0,0 +1,136 @@ +#ifndef ASYMPTOTIC_ENTRY_H +#define ASYMPTOTIC_ENTRY_H + +#include "asymptotic.h" + +/* Backend-independent numerical entry localizer for the common asymptotic + * exterior. + * + * The camera pre-route must decide, without touching any backend metric outside + * a worldtube, whether a past-directed camera ray crosses an escape worldtube + * from outside to inside and where the FIRST such entry lies. A supported + * exterior model may provide a closed-form (quadratic/analytic) entry; when it + * cannot, the caller supplies a path-parameter bracket that is already known to + * straddle the first entry and this module refines it numerically. + * + * The module never evaluates a metric, never loads a slab, never sweeps the + * movie in reverse and never snaps a position onto a radial shell. It only + * consumes the worldtube `escape_worldtube_sample` callback through + * `asymptotic_entry_geometry`, and an opaque evaluator callback that + * repropagates the exact supported exterior geodesic from its camera state to a + * path parameter. The driver is therefore independent of the exterior model + * and does not solve the entry equation itself. + * + * Parameter and time convention (18A.5/18A.6): the evaluator parameter + * increases along the renderer's backward propagation; the returned state's + * `activate_t` is the coordinate time at that parameter, so it is + * nonincreasing as the parameter increases. `outside_parameter` is the + * smaller-parameter end (worldtube outside, or exactly on the boundary) and + * `inside_parameter` is the larger-parameter end that is strictly inside. + * + * The caller is responsible for providing a correct first-entry bracket. This + * driver does not search arbitrary samples for a crossing; a bracket that does + * not straddle the boundary is reported as an unconfirmed entry + * (RAY_REASON_ENTRY_UNCONFIRMED), never as an escape. */ +typedef AsymptoticStatus (*AsymptoticEntryEvaluator)(void *context, + double parameter, + AsymptoticRoute *state); + +/* Convenience bundle for callers that want to keep the callback and its opaque + * context together. Not required by any entry point. */ +typedef struct { + AsymptoticEntryEvaluator evaluate; + void *context; +} AsymptoticEntryPropagator; + +/* Hard implementation guard on the number of bisection refinements. This + * bounds the double-parameter bisection; it is not a physical parameter. Any + * bracket near a finite nonzero entry contracts in roughly 60 halvings; + * very wide exponent ranges may instead exhaust the guard explicitly. */ +#define ASYMPTOTIC_ENTRY_BISECTION_LIMIT 256u + +/* Evaluate the worldtube function + * + * F(t, x) = |x - center(t)|^2 - radius(t)^2 + * + * at one state through the worldtube sample callback alone (no metric + * evaluation). `*F` is accumulated with the same left-to-right double + * arithmetic as the geodesic event layer, and `*tol` is the matching geometric + * ULP band 128 * DBL_EPSILON * max(radius^2, |x-center|^2). + * + * Status / failure reason: + * ASYMPTOTIC_OK -> *failure = RAY_REASON_NONE + * ASYMPTOTIC_TIME_RANGE_EXHAUSTED -> RAY_REASON_TIME_RANGE_EXHAUSTED + * ASYMPTOTIC_INVALID (callback failed)-> RAY_REASON_WORLDTUBE_SAMPLE_FAILED + * ASYMPTOTIC_INVALID (bad geometry) -> RAY_REASON_WORLDTUBE_GEOMETRY_INVALID + * + * A non-positive/non-finite radius, non-finite radius_rate, centers or + * velocities, and any non-finite (overflow/NaN) F or tolerance are geometry + * failures. `*failure` may be NULL. */ +AsymptoticStatus asymptotic_entry_geometry(const SpacetimeSource *source, + SpacetimeEndId end_id, double t, + const double x[3], double *F, + double *tol, RayReason *failure); + +/* Validate one caller-produced entry candidate. The candidate must be finite + * (activate_t, x, Pi, log_alpha_p0, log_alpha_p0_camera), carry kind + * ASYMPTOTIC_ROUTE_ENTRY, and sit + * on the worldtube boundary within the geometric ULP band: + * + * |F| <= tol -> *valid = 1 + * + * A candidate outside the worldtube (F > tol) or arbitrarily deep inside + * (F < -tol) is NOT an entry candidate: it yields *valid = 0 but still returns + * ASYMPTOTIC_OK, because the caller owns the exterior solve and a non-candidate + * is not a protocol error. A wrong kind or non-finite state is a protocol + * error. Worldtube callback/history/geometry failures propagate with the same + * reasons as asymptotic_entry_geometry. + * + * This function does not fabricate a velocity or re-derive the entry: the + * caller already guarantees that its exterior solve produced an inward entry. */ +AsymptoticStatus asymptotic_entry_validate(const SpacetimeSource *source, + SpacetimeEndId end_id, + const AsymptoticRoute *candidate, + int *valid, RayReason *failure); + +/* Numerically locate the first outside -> inside entry inside a known bracket. + * + * `evaluate(context, parameter, state)` repropagates the exact exterior geodesic + * from its camera state to `parameter`; it must fill `state->activate_t` (the + * coordinate time at that parameter), `state->x`, `state->Pi`, + * `state->log_alpha_p0`, `state->log_alpha_p0_camera` and `state->end_id`. + * Every call is counted into `*evaluations` (may be NULL). + * + * The two bracket endpoints are evaluated first. The driver requires + * F(outside_parameter) >= 0 and F(inside_parameter) < 0; a bracket with an + * exact boundary at the outside end (F == 0) and a strictly inside other end is + * returned directly as the legitimate inward entry. Otherwise it bisects the + * parameter. Each midpoint reconstructs a fresh state through the evaluator + * (no propagation from a prior midpoint, so no accumulated rounding or + * projection error), keeps the low end outside (F >= 0) and the high end + * strictly inside (F < 0), and never turns an F >= 0 midpoint into an escape. + * + * Bisection stops only when the parameter midpoint can no longer be represented + * strictly between the two ends (adjacent doubles), not on an F tolerance, so a + * curved trajectory is not stopped early by a coarser coordinate-time + * resolution. On success the strictly-inside endpoint adjacent to the entry in + * path parameter is returned: its Pi/L/L_camera are copied from the evaluator + * state unchanged, with kind forced to ENTRY and end_id set to `end_id`. The + * result is a representable bracketing of the entry (the true entry lies + * between the final outside and inside endpoints), not an absolute positional + * error claim; the returned inside state may legitimately have F < -tol. + * + * `ASYMPTOTIC_ENTRY_BISECTION_LIMIT` is an implementation guard on parameter + * representability, not a physical parameter. If the bracket never contracts, + * if the endpoints do not straddle the boundary, or if the endpoint coordinate + * times are not ordered (inside time <= outside time), the result is + * ASYMPTOTIC_INVALID with RAY_REASON_ENTRY_UNCONFIRMED. Callback, history and + * geometry failures from the evaluator/geometry propagate their exact status + * and reason; nothing is silently reported as a successful entry or escape. */ +AsymptoticStatus asymptotic_entry_localize( + const SpacetimeSource *source, SpacetimeEndId end_id, + AsymptoticEntryEvaluator evaluate, void *context, + double outside_parameter, double inside_parameter, AsymptoticRoute *out, + unsigned int *evaluations, RayReason *failure); + +#endif diff --git a/src/asymptotic_schwarzschild.c b/src/asymptotic_schwarzschild.c index 7f6a439..40af79e 100644 --- a/src/asymptotic_schwarzschild.c +++ b/src/asymptotic_schwarzschild.c @@ -408,6 +408,47 @@ int asymptotic_schwarzschild_state_from_canonical( return 0; } +int asymptotic_schwarzschild_inward_state_at_radius( + const SpacetimeAsymptoticEnd *end, const SchwarzschildCanonical *camera, + double rho, double x[3], double Pi[3], double *log_alpha_p0, + double *activate_t) { + if (end == NULL || camera == NULL || x == NULL || Pi == NULL) + return -1; + if (!(rho > 2.0) || !(rho <= camera->rho)) + return -1; + /* The inward branch only exists while the orbit has not turned before the + * requested radius. Allow a few ULP at the grazing limit so rounding in the + * turning root does not reject a legitimate boundary radius; a genuine + * inside-the-turning radius still fails through the Q >= 0 check below. */ + const double rho_turn = asymptotic_schwarzschild_turning_rho(camera->beta); + if (isfinite(rho_turn) && + rho < rho_turn - 16.0 * DBL_EPSILON * fmax(1.0, rho_turn)) + return -1; + const double dphi = asymptotic_schwarzschild_phi(rho, camera->beta) - + asymptotic_schwarzschild_phi(camera->rho, camera->beta); + if (!isfinite(dphi)) + return -1; + double rhat_rho[3]; + if (camera->beta > 0.0) + rotate_axis(camera->rhat, camera->Lhat, -dphi, rhat_rho); + else + for (int i = 0; i < 3; ++i) + rhat_rho[i] = camera->rhat[i]; + SchwarzschildCanonical state = *camera; + state.rho = rho; + for (int i = 0; i < 3; ++i) + state.rhat[i] = rhat_rho[i]; + state.radial_sign = -1; + if (asymptotic_schwarzschild_state_from_canonical(end, &state, x, Pi, + log_alpha_p0)) + return -1; + if (activate_t != NULL) { + const double T = sch_time_transfer(camera->rho, rho, camera->beta); + *activate_t = camera->t - end->mass * T; + } + return 0; +} + int asymptotic_schwarzschild_finish(const SpacetimeAsymptoticEnd *end, const SchwarzschildCanonical *canonical, double n_infinity[3], @@ -476,25 +517,10 @@ int asymptotic_schwarzschild_preroute( } if (camera->beta < beta_R) { - const double dphi = asymptotic_schwarzschild_phi(R, camera->beta) - - asymptotic_schwarzschild_phi(camera->rho, camera->beta); - double rhat_entry[3]; - if (camera->beta > 0.0) - rotate_axis(camera->rhat, camera->Lhat, -dphi, rhat_entry); - else - for (int i = 0; i < 3; ++i) - rhat_entry[i] = camera->rhat[i]; - SchwarzschildCanonical entry = *camera; - entry.rho = R; - for (int i = 0; i < 3; ++i) - entry.rhat[i] = rhat_entry[i]; - entry.radial_sign = -1; - if (asymptotic_schwarzschild_state_from_canonical(end, &entry, x, Pi, - log_alpha_p0)) + /* Shared exact inward transfer, also used by the common entry fallback. */ + if (asymptotic_schwarzschild_inward_state_at_radius( + end, camera, R, x, Pi, log_alpha_p0, activate_t)) return -1; - const double T = - sch_time_transfer(camera->rho, R, camera->beta); - *activate_t = camera->t - end->mass * T; *kind = SCH_ROUTE_ENTRY; return 0; } diff --git a/src/asymptotic_schwarzschild.h b/src/asymptotic_schwarzschild.h index 66232a0..af328b7 100644 --- a/src/asymptotic_schwarzschild.h +++ b/src/asymptotic_schwarzschild.h @@ -42,6 +42,19 @@ int asymptotic_schwarzschild_state_from_canonical( const SpacetimeAsymptoticEnd *end, const SchwarzschildCanonical *canonical, double x[3], double Pi[3], double *log_alpha_p0); +/* Repropagate the exact inward orbit from the ORIGINAL camera canonical state + * to `rho >= turning_rho` (and > 2) by recomputing the swept angle and the + * coordinate-time integral -- never by projecting a nearby state. Fills the + * backend entry state, its local L = ln(alpha p^0), and the activation + * coordinate time. Returns -1 when `rho` lies outside the reachable inward + * domain. This is the radius-parameter evaluator used by the common first + * entry localizer; `R` may dip slightly below the worldtube radius because the + * common driver only needs a strictly-inside bracket. */ +int asymptotic_schwarzschild_inward_state_at_radius( + const SpacetimeAsymptoticEnd *end, const SchwarzschildCanonical *camera, + double rho, double x[3], double Pi[3], double *log_alpha_p0, + double *activate_t); + /* Infinity endpoint for an outward crossing at the canonical radius. */ int asymptotic_schwarzschild_finish(const SpacetimeAsymptoticEnd *end, const SchwarzschildCanonical *canonical, diff --git a/src/geodesic.c b/src/geodesic.c index edb8eb0..c79df84 100644 --- a/src/geodesic.c +++ b/src/geodesic.c @@ -83,6 +83,8 @@ const char *ray_reason_name(RayReason reason) { return "THRESHOLD_EVENT_UNCONFIRMED"; case RAY_REASON_SLAB_LOAD_FAILED: return "SLAB_LOAD_FAILED"; + case RAY_REASON_ENTRY_UNCONFIRMED: + return "ENTRY_UNCONFIRMED"; case RAY_REASON_COUNT: break; } @@ -136,6 +138,7 @@ RayReason ray_reason_category(RayReason reason) { case RAY_REASON_SUBINTEGRATION_TARGET_MISSED: case RAY_REASON_ESCAPE_EVENT_UNCONFIRMED: case RAY_REASON_THRESHOLD_EVENT_UNCONFIRMED: + case RAY_REASON_ENTRY_UNCONFIRMED: return RAY_REASON_INTEGRATION_ERROR; case RAY_REASON_SLAB_LOAD_FAILED: return RAY_REASON_IO_ERROR; @@ -2169,7 +2172,10 @@ RayEndpoint geodesic_trace_past(const SpacetimeSource *source, } if (route_status != ASYMPTOTIC_OK) { out.outcome = RAY_OUTCOME_INCOMPLETE; - out.reason = RAY_REASON_CAMERA_PREROUTE_FAILED; + out.reason = ray_reason_valid(route.failure_reason) && + route.failure_reason != RAY_REASON_NONE + ? route.failure_reason : RAY_REASON_CAMERA_PREROUTE_FAILED; + out.end_id = route.end_id; return out; } if (route.kind == ASYMPTOTIC_ROUTE_ESCAPED) { diff --git a/src/geodesic.h b/src/geodesic.h index 0267b10..937710c 100644 --- a/src/geodesic.h +++ b/src/geodesic.h @@ -64,6 +64,7 @@ typedef enum { RAY_REASON_THRESHOLD_EVENT_UNCONFIRMED, /* threshold event retry exhausted */ /* I/O-derived. */ RAY_REASON_SLAB_LOAD_FAILED, /* spacetime_load_slab failed */ + RAY_REASON_ENTRY_UNCONFIRMED, /* no trustworthy first-entry bracket */ RAY_REASON_COUNT /* sentinel: valid ids are < COUNT */ } RayReason; diff --git a/src/ray.c b/src/ray.c index c005972..e92f749 100644 --- a/src/ray.c +++ b/src/ray.c @@ -149,8 +149,11 @@ void ray_pool_preroute(RayPool *p, const SpacetimeSource *source) { if (status == ASYMPTOTIC_UNSUPPORTED || status == ASYMPTOTIC_INVALID) { p->endpoint[i].outcome = RAY_OUTCOME_INCOMPLETE; p->endpoint[i].reason = status == ASYMPTOTIC_UNSUPPORTED - ? RAY_REASON_UNSUPPORTED - : RAY_REASON_CAMERA_PREROUTE_FAILED; + ? RAY_REASON_UNSUPPORTED + : (ray_reason_valid(route.failure_reason) && + route.failure_reason != RAY_REASON_NONE + ? route.failure_reason + : RAY_REASON_CAMERA_PREROUTE_FAILED); p->endpoint[i].end_id = route.end_id; p->status[i] = RAY_POOL_FAILED; continue; diff --git a/tests/test_asymptotic.c b/tests/test_asymptotic.c index 95cdd06..0d64b3d 100644 --- a/tests/test_asymptotic.c +++ b/tests/test_asymptotic.c @@ -3,6 +3,7 @@ #include "ray.h" #include "spacetime.h" +#include #include #include @@ -50,6 +51,7 @@ typedef struct { int sample_nonpositive_radius; int fail_on_sample_call; /* 1-based callback invocation to fail. */ int sample_call_count; + double frame_origin[3]; } SyntheticContext; static SpacetimePointStatus synthetic_eval(const SpacetimeSource *source, @@ -95,7 +97,8 @@ static int synthetic_end(const SpacetimeSource *source, size_t index, .end_id = 0, .exterior_kind = kind, .mass = mass, - .frame_origin = {0.0, 0.0, 0.0}, + .frame_origin = {context->frame_origin[0], context->frame_origin[1], + context->frame_origin[2]}, .frame_axes = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}}; return 0; } @@ -134,12 +137,17 @@ static int synthetic_worldtube(const SpacetimeSource *source, *out = (SpacetimeEscapeWorldtubeSample){.valid = 0}; return 0; } + /* The sample contract is the worldtube value at this coordinate time, so the + * reported radius must carry its own time dependence: R(t) = R0 + rr t, with + * dR/dt = radius_rate. A constant `radius` with a nonzero rate would make + * the closed-form segment model and the callback geometry disagree. */ + const double radius_t = context->radius + context->radius_rate * t; if (context->has_segment && t < context->segment_t) { /* Second segment: center moves toward +x as t decreases. */ *out = (SpacetimeEscapeWorldtubeSample){ .center = {context->segment_t - t, 0.0, 0.0}, .velocity = {-1.0, 0.0, 0.0}, - .radius = context->radius, + .radius = radius_t, .radius_rate = context->radius_rate, .velocity_constant = context->constant, .valid = 1}; @@ -148,7 +156,7 @@ static int synthetic_worldtube(const SpacetimeSource *source, *out = (SpacetimeEscapeWorldtubeSample){ .center = {context->vx * t + 0.5 * context->accel * t * t, 0.0, 0.0}, .velocity = {context->vx + context->accel * t, 0.0, 0.0}, - .radius = context->radius, + .radius = radius_t, .radius_rate = context->radius_rate, .velocity_constant = context->constant, .valid = 1}; @@ -710,6 +718,166 @@ static void test_moving_sphere(void) { "co-moving ray misses"); } +/* Fixed-observer tetrad used by the production Alcubierre observer-track rows + * 63/64, with spatial axes (e1, e2, e3) = (y-hat, z-hat, x-hat). The literal + * values are embedded here so this regression does not depend on the + * untracked observer CSV. */ +static ObserverState track_observer(double coordinate_time) { + ObserverState o = {0}; + o.coordinate_time = coordinate_time; + o.coordinate_position[0] = 0.0; + o.coordinate_position[1] = -24.0; + o.coordinate_position[2] = 0.0; + o.tetrad[0][0] = 1.0; + o.tetrad[1][2] = 1.0; + o.tetrad[2][3] = 1.0; + o.tetrad[3][1] = 1.0; + return o; +} + +/* Two exact production RayPool pre-route samples (frame 63 sample 12315 and + * frame 64 sample 3994). They are grazing (disc/b^2 ~ 1e-4), so the plain + * double root solve left the reconstructed entry state at F ~ 1.0-1.2 x + * geom_tol, which the event layer rejected as OUTSIDE_WORLDTUBE. The moving + * sphere fixture reproduces the Alcubierre worldtube (center = 2 t, radius 5); + * the observer time/position/tetrad and the camera direction are the exact raw + * production values. The route must land inside the geometric tolerance with + * the entry direction (Pi) unchanged. */ +static void test_grazing_production_entries(void) { + SyntheticContext context = {.vx = 2.0, + .accel = 0.0, + .radius = 5.0, + .radius_rate = 0.0, + .valid_t_min = -1.0e30, + .constant = 1}; + SpacetimeSource source = {.ops = &synthetic_ops, .context = &context}; + struct { + double time; + double direction[3]; + double pi[3]; + } cases[2] = { + {0x1.fa8f5c28f5c29p+3, + {0x1.c7378f8e872d1p-1, 0x1.15bad4e30e8ddp-8, -0x1.d4afba4704cap-2}, + {0x1.d4afba4704cap-2, -0x1.c7378f8e872d1p-1, -0x1.15bad4e30e8ddp-8}}, + {0x1.fb17e4b17e4b1p+3, + {0x1.bd3bb364ac492p-1, 0x1.102d2a1c6ac74p-7, -0x1.f98ae1a782104p-2}, + {0x1.f98ae1a782104p-2, -0x1.bd3bb364ac492p-1, -0x1.102d2a1c6ac74p-7}}, + }; + for (int c = 0; c < 2; ++c) { + ObserverState observer = track_observer(cases[c].time); + AsymptoticRoute route; + CHECK(asymptotic_route_camera(&source, &observer, cases[c].direction, + &route) == ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ENTRY, + "grazing production entry found"); + CHECK(route.entry_fallback_evaluations == 0, + "grazing production entry stays on the fast path"); + for (int i = 0; i < 3; ++i) + CHECK(route.Pi[i] == cases[c].pi[i], + "grazing entry direction is unchanged"); + SpacetimeEscapeWorldtubeSample sample; + CHECK(spacetime_escape_worldtube_sample(&source, route.end_id, + route.activate_t, &sample) == 0 && + sample.valid && sample.radius > 0.0, + "grazing entry sample valid"); + double d2 = 0.0; + for (int k = 0; k < 3; ++k) { + const double dk = route.x[k] - sample.center[k]; + d2 += dk * dk; + } + const double r2 = sample.radius * sample.radius; + double value; + CHECK(asymptotic_worldtube_value(&source, route.end_id, route.activate_t, + route.x, &value) == 0, + "grazing entry worldtube value"); + const double geom_tol = 128.0 * DBL_EPSILON * fmax(r2, d2); + /* The event layer rejects the entry (OUTSIDE_WORLDTUBE) exactly when + * F > geom_tol; require the residual to sit inside the tolerance band + * rather than accepting an arbitrary sign. */ + CHECK(value <= geom_tol && value >= -geom_tol, + "grazing entry F within geometric tolerance"); + } +} + +/* Linear (a == 0) entry: a growing sphere whose radius rate cancels the + * relative closing speed, so qq == rr^2 and the quadratic degenerates. The + * stable solver must still take the smallest positive root. The fixture's + * sample() reports the consistent radius R(t) = 10 - t, i.e. R(s) = 10 + s + * along the past parameter s = -t, so the contact point is on the true + * ruled-surface boundary. */ +static void test_linear_a_zero_entry(void) { + SyntheticContext context = {.vx = 0.0, + .accel = 0.0, + .radius = 10.0, + .radius_rate = -1.0, /* R(t) = 10 - t */ + .valid_t_min = -1.0e30, + .constant = 1}; + SpacetimeSource source = {.ops = &synthetic_ops, .context = &context}; + const ObserverState camera = flat_observer(100.0, 0.0, 0.0); + AsymptoticRoute route; + CHECK(asymptotic_route_camera(&source, &camera, (double[]){-1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ENTRY, + "linear a~0 entry"); + CHECK(fabs(route.activate_t + 45.0) < 1e-12, "linear entry time"); + CHECK(fabs(route.x[0] - 55.0) < 1e-12 && fabs(route.x[1]) < 1e-12 && + fabs(route.x[2]) < 1e-12, + "linear entry position"); +} + +/* Large-coordinate-time cancellation: the long-double closed-form root is + * accurate in the frame, but the entry state reconstructed in double at a huge + * t0 loses the sub-ULP part of the event and lands far outside the geometric + * ULP band (F ~ 0.3 >> tol). The common fallback must repropagate from the + * original camera, localize the first entry numerically, and return a + * strict-inside endpoint (F < 0) while preserving the camera direction Pi and + * reference L exactly. The fixture is a legitimate constant-velocity + * worldtube, not a nonlinearity injection. */ +static void test_fallback_reconstruction_cancellation(void) { + SyntheticContext context = {.vx = 0x1.999999999999ap-4, /* 0.1 */ + .accel = 0.0, + .radius = 10.0, + .radius_rate = 0.0, + .valid_t_min = -1.0e30, + .constant = 1}; + SpacetimeSource source = {.ops = &synthetic_ops, .context = &context}; + const double t0 = 1.0e15; + /* Exactly 0.1 * 1e15 + 100.123456789, pinned as a hex literal. The offset + * is not aligned to the t0 ULP, so activate_t = t0 - s rounds and the + * reconstructed boundary residual exceeds the tolerance band. */ + const double camera_x = 0x1.6bcc41e901908p+46; + ObserverState observer = flat_observer(camera_x, 0.0, 0.0); + observer.coordinate_time = t0; + const double direction[3] = {-1.0, 0.0, 0.0}; + AsymptoticRoute route; + CHECK(asymptotic_route_camera(&source, &observer, direction, &route) == + ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ENTRY, + "cancellation fallback still produces an entry"); + CHECK(route.entry_fallback_evaluations > 0, + "cancellation entry used the common fallback"); + CHECK(route.failure_reason == RAY_REASON_NONE, + "successful fallback has no failure reason"); + double value; + CHECK(asymptotic_worldtube_value(&source, route.end_id, route.activate_t, + route.x, &value) == 0 && + value <= 0.0, + "fallback entry state is inside the worldtube"); + + MetricData metric; + GeodesicRayState camera_state; + CHECK(spacetime_eval(&source, t0, observer.coordinate_position, &metric) == + 0 && + geodesic_initialize_past_ray_metric(&metric, &observer, direction, + &camera_state) == 0, + "cancellation camera state"); + for (int i = 0; i < 3; ++i) + CHECK(route.Pi[i] == camera_state.Pi[i], + "fallback preserves the entry direction exactly"); + CHECK(route.log_alpha_p0_camera == camera_state.log_alpha_p0, + "fallback preserves the camera reference L exactly"); +} + static void test_accelerated_worldtube_unsupported(void) { /* A genuinely accelerating (non-constant velocity) worldtube has no strict * relative-motion interval bound, so the route is explicitly unsupported. @@ -789,11 +957,98 @@ static void test_ray_pool_lifecycle(void) { spacetime_destroy(&source); } +static void test_zero_discriminant_is_not_miss(void) { + SpacetimeSource source; + CHECK(spacetime_create_minkowski(&source, 1.0) == 0, + "zero-discriminant source"); + const ObserverState observer = flat_observer(1e10, 0.5, 0.0); + AsymptoticRoute route; + /* Forming c = 1e20 + 0.25 - 1 loses the transverse contribution even in + * 80-bit arithmetic; b*b - 4*a*c then rounds to zero despite a real entry. */ + CHECK(asymptotic_route_camera(&source, &observer, + (double[]){-1.0, 0.0, 0.0}, &route) == + ASYMPTOTIC_OK && route.kind == ASYMPTOTIC_ROUTE_ENTRY && + route.entry_fallback_evaluations > 0, + "rounded zero discriminant uses entry fallback, not escape"); + double value = 0.0; + CHECK(asymptotic_worldtube_value(&source, 0, route.activate_t, route.x, + &value) == ASYMPTOTIC_OK && value < 0.0, + "zero-discriminant fallback produces actual inside state"); + const ObserverState tangent = flat_observer(1e10, 1.0, 0.0); + CHECK(asymptotic_route_camera(&source, &tangent, + (double[]){-1.0, 0.0, 0.0}, &route) == + ASYMPTOTIC_OK && route.kind == ASYMPTOTIC_ROUTE_ESCAPED, + "fixed transverse coordinate independently certifies exact tangency"); + spacetime_destroy(&source); +} + +static void test_positive_reconstructed_minimum_is_not_miss(void) { + SyntheticContext ctx = {.vx = 0.1, .radius = 10.0, .constant = 1}; + SpacetimeSource source = {.ops = &synthetic_ops, .context = &ctx}; + ObserverState observer = flat_observer(0x1.6bcc41e901904p+46, + 0x1.3ffffde7210bfp+3, 0.0); + observer.coordinate_time = 1e15; + const double parameter = 0x1.bca8814065f1ep+6; + const double witness[3] = {observer.coordinate_position[0] - parameter, + observer.coordinate_position[1], 0.0}; + double value = 0.0; + CHECK(asymptotic_worldtube_value(&source, 0, + observer.coordinate_time - parameter, + witness, &value) == ASYMPTOTIC_OK && value < 0, + "strict-inside witness exists despite positive reconstructed minimum"); + AsymptoticRoute route; + CHECK(asymptotic_route_camera(&source, &observer, + (double[]){-1.0, 0.0, 0.0}, &route) == + ASYMPTOTIC_INVALID && + route.failure_reason == RAY_REASON_ENTRY_UNCONFIRMED, + "positive probe without a miss certificate fails explicitly"); + RayPool pool; + CHECK(ray_pool_init(&pool, 1) == 0, "ambiguous entry pool"); + CHECK(ray_pool_append(&pool, &observer, (double[]){-1.0, 0.0, 0.0}, + 0, 0) == 0, "ambiguous entry ray"); + ray_pool_preroute(&pool, &source); + CHECK(pool.status[0] == RAY_POOL_FAILED && + pool.endpoint[0].outcome == RAY_OUTCOME_INCOMPLETE && + pool.endpoint[0].reason == RAY_REASON_ENTRY_UNCONFIRMED, + "pool propagates unconfirmed entry instead of fabricating escape"); + ray_pool_destroy(&pool); +} + +static void test_translated_frame_cannot_certify_miss(void) { + SyntheticContext ctx = {.radius = 1.0, .constant = 1, + .valid_t_min = -1e100, + .frame_origin = {0.0, 1e10, 0.0}}; + SpacetimeSource source = {.ops = &synthetic_ops, .context = &ctx}; + const ObserverState observer = flat_observer(1e10, 1.0 - 0x1p-22, 0.0); + double value = 0.0; + const double witness[3] = {0.0, observer.coordinate_position[1], 0.0}; + CHECK(asymptotic_worldtube_value(&source, 0, -1e10, witness, &value) == + ASYMPTOTIC_OK && value < 0.0, + "original untranslated trajectory has an inside witness"); + AsymptoticRoute route; + CHECK(asymptotic_route_camera(&source, &observer, + (double[]){-1.0, 0.0, 0.0}, &route) == + ASYMPTOTIC_INVALID && + route.failure_reason == RAY_REASON_ENTRY_UNCONFIRMED, + "lossy translated frame must not certify a miss"); + ctx.frame_origin[1] = 0.0; + CHECK(asymptotic_route_camera(&source, &observer, + (double[]){-1.0, 0.0, 0.0}, &route) == + ASYMPTOTIC_OK && route.kind == ASYMPTOTIC_ROUTE_ENTRY, + "same unshifted trajectory confirms an entry"); +} + int main(void) { test_fixed_sphere(); test_large_radius_quadratic(); test_round_trip(); test_moving_sphere(); + test_grazing_production_entries(); + test_linear_a_zero_entry(); + test_fallback_reconstruction_cancellation(); + test_zero_discriminant_is_not_miss(); + test_positive_reconstructed_minimum_is_not_miss(); + test_translated_frame_cannot_certify_miss(); test_accelerated_worldtube_unsupported(); test_piecewise_segment_entry(); test_boundary_semantics_minkowski(); diff --git a/tests/test_asymptotic_entry.c b/tests/test_asymptotic_entry.c new file mode 100644 index 0000000..69b1421 --- /dev/null +++ b/tests/test_asymptotic_entry.c @@ -0,0 +1,781 @@ +#include "asymptotic_entry.h" +#include "spacetime.h" + +#include +#include +#include + +/* Backend-independent core regression for the numerical entry localizer. It + * deliberately links no analytic backend and no geodesic integrator: the fake + * SpacetimeSource exposes only `escape_worldtube_sample` (plus a deliberately + * trapped `eval`), and the evaluator is an analytic path-parameter callback. + * + * Nothing here depends on an untracked production track, CSV or binary. */ + +static int failures = 0; + +#define CHECK(condition, message) \ + do { \ + if (!(condition)) { \ + fprintf(stderr, "FAIL %s:%d: %s\n", __FILE__, __LINE__, message); \ + ++failures; \ + } \ + } while (0) + +#ifndef TEST_PI +#define TEST_PI 3.14159265358979323846 +#endif + +/* ------------------------------------------------------------------ */ +/* Fake worldtube source */ +/* ------------------------------------------------------------------ */ + +typedef struct { + double center0[3]; + double center_vel[3]; /* dc/dt */ + double radius0; + double radius_rate; /* dR/dt */ + double valid_t_min; /* sample is valid for t >= valid_t_min */ + int callback_fails; /* always return -1 */ + int fail_at_call; /* 1-based sample-call index to fail, 0 disabled */ + int nan_radius; + int zero_radius; + double hole_center; /* isolated invalid time window */ + double hole_halfwidth; /* 0 disables the window */ + int call_count; + int eval_calls; /* trap: how often the metric eval callback ran */ +} EntryWorldtube; + +static SpacetimePointStatus entry_eval_trap(const SpacetimeSource *source, + double t, const double x[3], + MetricData *metric) { + EntryWorldtube *wt = source->context; + ++wt->eval_calls; + (void)t; + (void)x; + (void)metric; + /* This source is deliberately outside the metric domain. The localizer must + * never reach here because it does no metric evaluation. */ + return SPACETIME_POINT_OUT_OF_DOMAIN; +} + +static int entry_worldtube_cb(const SpacetimeSource *source, + SpacetimeEndId end_id, double t, + SpacetimeEscapeWorldtubeSample *out) { + EntryWorldtube *wt = source->context; + if (end_id != 0) + return -1; + ++wt->call_count; + if (wt->fail_at_call > 0 && wt->call_count == wt->fail_at_call) + return -1; + if (wt->callback_fails) + return -1; + if (!isfinite(t) || t < wt->valid_t_min) { + *out = (SpacetimeEscapeWorldtubeSample){.valid = 0}; + return 0; + } + if (wt->hole_halfwidth > 0.0 && + fabs(t - wt->hole_center) <= wt->hole_halfwidth) { + *out = (SpacetimeEscapeWorldtubeSample){.valid = 0}; + return 0; + } + if (wt->nan_radius) { + *out = (SpacetimeEscapeWorldtubeSample){.radius = NAN, .valid = 1}; + return 0; + } + if (wt->zero_radius) { + *out = (SpacetimeEscapeWorldtubeSample){.radius = 0.0, .valid = 1}; + return 0; + } + *out = (SpacetimeEscapeWorldtubeSample){ + .center = {wt->center0[0] + wt->center_vel[0] * t, + wt->center0[1] + wt->center_vel[1] * t, + wt->center0[2] + wt->center_vel[2] * t}, + .velocity = {wt->center_vel[0], wt->center_vel[1], wt->center_vel[2]}, + .radius = wt->radius0 + wt->radius_rate * t, + .radius_rate = wt->radius_rate, + .velocity_constant = 1, + .valid = 1}; + return 0; +} + +static const SpacetimeOps entry_ops = { + .eval = entry_eval_trap, + .escape_worldtube_sample = entry_worldtube_cb, +}; + +static SpacetimeSource entry_source(EntryWorldtube *wt) { + return (SpacetimeSource){.ops = &entry_ops, .context = wt}; +} + +/* ------------------------------------------------------------------ */ +/* Analytic path-parameter evaluator */ +/* ------------------------------------------------------------------ */ + +typedef struct { + double camera_t; + double camera_x[3]; + double w[3]; /* unit past direction (straight mode) */ + int arc_mode; + double arc_center[3]; + double arc_radius; + double arc_theta0; + double L0; + double L0camera; + int evaluator_fails_at; + AsymptoticStatus fail_status; + int evaluator_call_count; + int nonfinite_at; + int reversed_time_at; +} EntryEvaluator; + +static void entry_trajectory(const EntryEvaluator *c, double parameter, + double x[3], double w[3], double *t) { + if (c->arc_mode) { + /* Circular analytic arc: parameter is arc length. Not a physical + * geodesic, but a generic curved callback that exercises the driver beyond + * straight lines. */ + const double theta = c->arc_theta0 + parameter / c->arc_radius; + x[0] = c->arc_center[0] + c->arc_radius * cos(theta); + x[1] = c->arc_center[1] + c->arc_radius * sin(theta); + x[2] = c->arc_center[2]; + w[0] = -sin(theta); + w[1] = cos(theta); + w[2] = 0.0; + } else { + for (int i = 0; i < 3; ++i) { + x[i] = c->camera_x[i] + parameter * c->w[i]; + w[i] = c->w[i]; + } + } + *t = c->camera_t - parameter; +} + +static AsymptoticStatus entry_evaluator_cb(void *context, double parameter, + AsymptoticRoute *state) { + EntryEvaluator *c = context; + ++c->evaluator_call_count; + if (c->evaluator_fails_at > 0 && + c->evaluator_call_count == c->evaluator_fails_at) + return c->fail_status; + double x[3], w[3], t; + entry_trajectory(c, parameter, x, w, &t); + *state = (AsymptoticRoute){0}; + state->kind = ASYMPTOTIC_ROUTE_ENTRY; + state->end_id = 0; + state->activate_t = t; + for (int i = 0; i < 3; ++i) { + state->x[i] = x[i]; + state->Pi[i] = -w[i]; + } + state->log_alpha_p0 = c->L0; + state->log_alpha_p0_camera = c->L0camera; + if (c->evaluator_call_count == c->nonfinite_at) + state->log_alpha_p0_camera = NAN; + if (c->evaluator_call_count == c->reversed_time_at) + state->activate_t = c->camera_t + 1.0; + return ASYMPTOTIC_OK; +} + +/* Independent test-side oracle: the same worldtube F the driver sees, but + * computed directly from the analytic trajectory. Used only to find the true + * first entry for comparison. */ +typedef struct { + const EntryEvaluator *ev; + const EntryWorldtube *wt; +} EntryOracle; + +static double entry_oracle_F(void *context, double parameter) { + const EntryOracle *o = context; + double x[3], w[3], t; + entry_trajectory(o->ev, parameter, x, w, &t); + double d[3]; + for (int i = 0; i < 3; ++i) + d[i] = x[i] - (o->wt->center0[i] + o->wt->center_vel[i] * t); + const double d2 = d[0] * d[0] + d[1] * d[1] + d[2] * d[2]; + const double radius = o->wt->radius0 + o->wt->radius_rate * t; + return d2 - radius * radius; +} + +static double entry_oracle_root(const EntryEvaluator *ev, + const EntryWorldtube *wt, double lo, + double hi) { + EntryOracle o = {.ev = ev, .wt = wt}; + if (!(entry_oracle_F(&o, lo) >= 0.0 && entry_oracle_F(&o, hi) < 0.0)) + return NAN; + for (int i = 0; i < 200; ++i) { + const double mid = 0.5 * (lo + hi); + if (!(mid > lo && mid < hi)) + break; + if (entry_oracle_F(&o, mid) >= 0.0) + lo = mid; + else + hi = mid; + } + return 0.5 * (lo + hi); +} + +static double path_parameter(const EntryEvaluator *ev, + const AsymptoticRoute *state) { + return ev->camera_t - state->activate_t; +} + +/* ------------------------------------------------------------------ */ +/* Tests */ +/* ------------------------------------------------------------------ */ + +static void test_geometry_contract(void) { + EntryWorldtube wt = {.radius0 = 10.0, .valid_t_min = -1.0e300}; + SpacetimeSource source = entry_source(&wt); + double F = NAN, tol = NAN; + RayReason reason = RAY_REASON_COUNT; + + CHECK(asymptotic_entry_geometry(&source, 0, 0.0, (double[]){10.0, 0.0, 0.0}, + &F, &tol, &reason) == ASYMPTOTIC_OK, + "boundary geometry status"); + CHECK(F == 0.0, "boundary F is exactly zero"); + CHECK(tol > 0.0 && isfinite(tol), "boundary tolerance finite positive"); + CHECK(reason == RAY_REASON_NONE, "boundary reason none"); + + CHECK(asymptotic_entry_geometry(&source, 0, 0.0, (double[]){20.0, 0.0, 0.0}, + &F, &tol, &reason) == ASYMPTOTIC_OK && + F == 300.0, + "outside F is positive 300"); + CHECK(asymptotic_entry_geometry(&source, 0, 0.0, (double[]){5.0, 0.0, 0.0}, + &F, &tol, &reason) == ASYMPTOTIC_OK && + F == -75.0, + "inside F is negative 75"); + + /* History exhaustion beats a miss. */ + EntryWorldtube hole = {.radius0 = 10.0, .valid_t_min = 0.0}; + source = entry_source(&hole); + CHECK(asymptotic_entry_geometry(&source, 0, -1.0, (double[]){10.0, 0.0, 0.0}, + &F, &tol, &reason) == + ASYMPTOTIC_TIME_RANGE_EXHAUSTED && + reason == RAY_REASON_TIME_RANGE_EXHAUSTED, + "valid=0 is history exhaustion"); + + /* Callback failure is distinct from an invalid geometry. */ + EntryWorldtube fail = {.radius0 = 10.0, .callback_fails = 1}; + source = entry_source(&fail); + CHECK(asymptotic_entry_geometry(&source, 0, 0.0, (double[]){10.0, 0.0, 0.0}, + &F, &tol, &reason) == ASYMPTOTIC_INVALID && + reason == RAY_REASON_WORLDTUBE_SAMPLE_FAILED, + "callback failure reason"); + + EntryWorldtube nanr = {.radius0 = 10.0, .nan_radius = 1}; + source = entry_source(&nanr); + CHECK(asymptotic_entry_geometry(&source, 0, 0.0, (double[]){10.0, 0.0, 0.0}, + &F, &tol, &reason) == ASYMPTOTIC_INVALID && + reason == RAY_REASON_WORLDTUBE_GEOMETRY_INVALID, + "NaN radius is invalid geometry"); + + EntryWorldtube zeror = {.radius0 = 10.0, .zero_radius = 1}; + source = entry_source(&zeror); + CHECK(asymptotic_entry_geometry(&source, 0, 0.0, (double[]){10.0, 0.0, 0.0}, + &F, &tol, &reason) == ASYMPTOTIC_INVALID && + reason == RAY_REASON_WORLDTUBE_GEOMETRY_INVALID, + "non-positive radius is invalid geometry"); + + source = entry_source(&wt); + CHECK(asymptotic_entry_geometry(&source, 0, NAN, (double[]){10.0, 0.0, 0.0}, + &F, &tol, &reason) == ASYMPTOTIC_INVALID && + reason == RAY_REASON_WORLDTUBE_GEOMETRY_INVALID, + "NaN time is invalid geometry"); + CHECK(asymptotic_entry_geometry(NULL, 0, 0.0, (double[]){10.0, 0.0, 0.0}, + &F, &tol, &reason) == ASYMPTOTIC_INVALID && + reason == RAY_REASON_INVALID_ARGUMENT, + "NULL source rejected"); + CHECK(asymptotic_entry_geometry(&source, 0, 0.0, NULL, &F, &tol, + &reason) == ASYMPTOTIC_INVALID, + "NULL position rejected"); +} + +static void test_validate_contract(void) { + EntryWorldtube wt = {.radius0 = 10.0, .valid_t_min = -1.0e300}; + SpacetimeSource source = entry_source(&wt); + int valid = -1; + RayReason reason = RAY_REASON_COUNT; + AsymptoticRoute candidate = {.kind = ASYMPTOTIC_ROUTE_ENTRY, + .end_id = 0, + .activate_t = 0.0, + .x = {10.0, 0.0, 0.0}, + .Pi = {-1.0, 0.0, 0.0}, + .log_alpha_p0 = 0.5, + .log_alpha_p0_camera = 0.25}; + + CHECK(asymptotic_entry_validate(&source, 0, &candidate, &valid, &reason) == + ASYMPTOTIC_OK && + valid == 1, + "boundary candidate is valid"); + + candidate.x[0] = 11.0; /* F = 21 > tol */ + CHECK(asymptotic_entry_validate(&source, 0, &candidate, &valid, &reason) == + ASYMPTOTIC_OK && + valid == 0 && reason == RAY_REASON_NONE, + "outside candidate is valid=0 with OK status"); + + candidate.x[0] = 5.0; /* F = -75, far inside */ + CHECK(asymptotic_entry_validate(&source, 0, &candidate, &valid, &reason) == + ASYMPTOTIC_OK && + valid == 0, + "deep-inside candidate is valid=0 with OK status"); + + candidate.x[0] = 10.0; + candidate.kind = ASYMPTOTIC_ROUTE_ESCAPED; + CHECK(asymptotic_entry_validate(&source, 0, &candidate, &valid, &reason) == + ASYMPTOTIC_INVALID && + valid == 0 && reason == RAY_REASON_PROTOCOL_ERROR, + "wrong candidate kind is a protocol error"); + + candidate.kind = ASYMPTOTIC_ROUTE_ENTRY; + candidate.x[0] = NAN; + CHECK(asymptotic_entry_validate(&source, 0, &candidate, &valid, &reason) == + ASYMPTOTIC_INVALID && + reason == RAY_REASON_PROTOCOL_ERROR, + "non-finite candidate is a protocol error"); + + /* History exhaustion propagates through validation. */ + candidate.x[0] = 10.0; + candidate.activate_t = -1.0; + EntryWorldtube hole = {.radius0 = 10.0, .valid_t_min = 0.0}; + source = entry_source(&hole); + CHECK(asymptotic_entry_validate(&source, 0, &candidate, &valid, &reason) == + ASYMPTOTIC_TIME_RANGE_EXHAUSTED && + reason == RAY_REASON_TIME_RANGE_EXHAUSTED, + "validation propagates history exhaustion"); +} + +static void test_fixed_sphere_localize(void) { + EntryWorldtube wt = {.radius0 = 10.0, .valid_t_min = -1.0e300}; + SpacetimeSource source = entry_source(&wt); + EntryEvaluator ev = {.camera_t = 0.0, + .camera_x = {50.0, 0.0, 0.0}, + .w = {-1.0, 0.0, 0.0}, + .L0 = 0.75, + .L0camera = 0.5}; + AsymptoticRoute out; + unsigned int evaluations = 0; + RayReason reason = RAY_REASON_COUNT; + + CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 30.0, + 45.0, &out, &evaluations, + &reason) == ASYMPTOTIC_OK, + "fixed-sphere localize succeeds"); + CHECK(out.kind == ASYMPTOTIC_ROUTE_ENTRY && out.end_id == 0, + "localized kind and end"); + const double s = path_parameter(&ev, &out); + const double s_true = entry_oracle_root(&ev, &wt, 30.0, 45.0); + CHECK(isfinite(s_true), "oracle found the same bracket"); + CHECK(s >= s_true && s - s_true <= 1e-9, + "localized just past first entry"); + CHECK(fabs(s - 40.0) <= 1e-9, "fixed-sphere entry at s=40"); + CHECK(out.Pi[0] == 1.0 && out.Pi[1] == 0.0 && out.Pi[2] == 0.0, + "direction preserved exactly"); + CHECK(out.log_alpha_p0 == 0.75 && out.log_alpha_p0_camera == 0.5, + "L and camera L preserved exactly"); + CHECK(evaluations >= 2 && evaluations <= 260, "evaluation count bounded"); + CHECK(ev.evaluator_call_count == (int)evaluations, + "evaluator calls counted once each"); + CHECK(wt.eval_calls == 0, "no metric evaluation outside the worldtube"); +} + +static void test_too_early_hint(void) { + /* Outside endpoint is the camera (a deliberately too-early, corrupted + * bracket); the localizer still returns the true first entry, not the + * inside hint and not the camera. */ + EntryWorldtube wt = {.radius0 = 10.0, .valid_t_min = -1.0e300}; + SpacetimeSource source = entry_source(&wt); + EntryEvaluator ev = {.camera_t = 0.0, + .camera_x = {50.0, 0.0, 0.0}, + .w = {-1.0, 0.0, 0.0}, + .L0 = 0.1, + .L0camera = 0.2}; + AsymptoticRoute out; + unsigned int evaluations = 0; + RayReason reason = RAY_REASON_COUNT; + + CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 0.0, + 45.0, &out, &evaluations, + &reason) == ASYMPTOTIC_OK, + "too-early hint still localizes"); + const double s = path_parameter(&ev, &out); + CHECK(fabs(s - 40.0) <= 1e-9, "returns actual first entry, not the hint"); + CHECK(s > 1.0 && s < 45.0, "not the camera and not the inside hint"); + CHECK(evaluations <= 260, "hint evaluation budget"); +} + +static void test_moving_sphere_localize(void) { + EntryWorldtube wt = {.radius0 = 10.0, + .center_vel = {0.5, 0.0, 0.0}, + .valid_t_min = -1.0e300}; + SpacetimeSource source = entry_source(&wt); + EntryEvaluator ev = {.camera_t = 0.0, + .camera_x = {100.0, 0.0, 0.0}, + .w = {-1.0, 0.0, 0.0}, + .L0 = 0.3, + .L0camera = 0.4}; + AsymptoticRoute out; + unsigned int evaluations = 0; + RayReason reason = RAY_REASON_COUNT; + + CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 150.0, + 200.0, &out, &evaluations, + &reason) == ASYMPTOTIC_OK, + "translated+ moving sphere localize"); + const double s = path_parameter(&ev, &out); + const double s_true = entry_oracle_root(&ev, &wt, 150.0, 200.0); + CHECK(fabs(s - s_true) <= 1e-8 && fabs(s - 180.0) <= 1e-8, + "moving-sphere entry at s=180"); + CHECK(evaluations <= 260, "moving-sphere evaluation budget"); +} + +static void test_radius_rate_localize(void) { + /* radius(t) = radius0 + radius_rate * t with radius_rate = -1 and t = -s, so + * R grows as 10 + s; the entry is at s = 45. */ + EntryWorldtube wt = {.radius0 = 10.0, + .radius_rate = -1.0, + .valid_t_min = -1.0e300}; + SpacetimeSource source = entry_source(&wt); + EntryEvaluator ev = {.camera_t = 0.0, + .camera_x = {100.0, 0.0, 0.0}, + .w = {-1.0, 0.0, 0.0}, + .L0 = 0.6, + .L0camera = 0.6}; + AsymptoticRoute out; + unsigned int evaluations = 0; + RayReason reason = RAY_REASON_COUNT; + + CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 30.0, + 60.0, &out, &evaluations, + &reason) == ASYMPTOTIC_OK, + "linear radius-rate localize"); + const double s = path_parameter(&ev, &out); + const double s_true = entry_oracle_root(&ev, &wt, 30.0, 60.0); + CHECK(fabs(s - s_true) <= 1e-8 && fabs(s - 45.0) <= 1e-8, + "linear radius-rate entry at s=45"); + CHECK(evaluations <= 260, "radius-rate evaluation budget"); +} + +static void test_rotated_frame_localize(void) { + /* Camera on a rotated axis: (40,30,0), past direction toward the origin. */ + EntryWorldtube wt = {.radius0 = 10.0, .valid_t_min = -1.0e300}; + SpacetimeSource source = entry_source(&wt); + EntryEvaluator ev = {.camera_t = 0.0, + .camera_x = {40.0, 30.0, 0.0}, + .w = {-0.8, -0.6, 0.0}, + .L0 = 0.2, + .L0camera = 0.1}; + AsymptoticRoute out; + unsigned int evaluations = 0; + RayReason reason = RAY_REASON_COUNT; + + CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 30.0, + 45.0, &out, &evaluations, + &reason) == ASYMPTOTIC_OK, + "rotated flat frame localize"); + const double s = path_parameter(&ev, &out); + const double s_true = entry_oracle_root(&ev, &wt, 30.0, 45.0); + CHECK(fabs(s - s_true) <= 1e-9 && fabs(s - 40.0) <= 1e-9, + "rotated-frame entry at s=40"); + CHECK(fabs(out.x[1] - 6.0) <= 1e-6, "rotated entry position on sphere"); +} + +static void test_grazing_first_entry(void) { + /* Grazing pass: the camera is offset by 9.9 from the sphere axis. The first + * entry at s ~ 48.589 is inside the bracket; the exit at s ~ 51.410 is not. + * Bisection must return the first entry, not the later exit. */ + EntryWorldtube wt = {.radius0 = 10.0, .valid_t_min = -1.0e300}; + SpacetimeSource source = entry_source(&wt); + EntryEvaluator ev = {.camera_t = 0.0, + .camera_x = {50.0, 9.9, 0.0}, + .w = {-1.0, 0.0, 0.0}, + .L0 = 0.0, + .L0camera = 0.0}; + AsymptoticRoute out; + unsigned int evaluations = 0; + RayReason reason = RAY_REASON_COUNT; + + CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 48.0, + 50.0, &out, &evaluations, + &reason) == ASYMPTOTIC_OK, + "grazing first entry localizes"); + const double s = path_parameter(&ev, &out); + const double first = 50.0 - sqrt(100.0 - 9.9 * 9.9); + CHECK(fabs(s - first) <= 1e-8, "grazing entry is the first crossing"); + CHECK(s < 51.4, "not the later exit crossing"); + CHECK(fabs(out.x[1] - 9.9) <= 1e-9, "grazing impact parameter preserved"); + CHECK(evaluations <= 260, "grazing evaluation budget"); +} + +static void test_curved_arc_localize(void) { + /* Circular analytic arc of radius 30 and worldtube centered at (25,0,0) + * radius 8; entry at arc length ~ 87.40. */ + EntryWorldtube wt = {.radius0 = 8.0, + .center0 = {25.0, 0.0, 0.0}, + .valid_t_min = -1.0e300}; + SpacetimeSource source = entry_source(&wt); + EntryEvaluator ev = {.camera_t = 0.0, + .arc_mode = 1, + .arc_center = {0.0, 0.0, 0.0}, + .arc_radius = 30.0, + .arc_theta0 = TEST_PI, + .L0 = 0.9, + .L0camera = 0.8}; + AsymptoticRoute out; + unsigned int evaluations = 0; + RayReason reason = RAY_REASON_COUNT; + + const double lo = 30.0 * (5.9 - TEST_PI); + const double hi = 30.0 * (6.2 - TEST_PI); + const double s_true = entry_oracle_root(&ev, &wt, lo, hi); + CHECK(isfinite(s_true), "curved oracle bracket"); + CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, lo, hi, + &out, &evaluations, + &reason) == ASYMPTOTIC_OK, + "curved arc localize"); + const double s = path_parameter(&ev, &out); + CHECK(s >= s_true - 1e-9 && s - s_true <= 1e-8, + "curved arc entry matches the oracle"); + /* d^2(theta) = 1525 - 1500 cos(theta) = 8^2 on the arc. */ + const double expected = + 30.0 * (2.0 * TEST_PI - acos((1525.0 - 64.0) / 1500.0) - TEST_PI); + CHECK(fabs(s - expected) <= 1e-8, "curved arc entry matches analytic root"); + CHECK(evaluations <= 260, "curved arc evaluation budget"); +} + +static void test_boundary_entry_exact(void) { + EntryWorldtube wt = {.radius0 = 10.0, .valid_t_min = -1.0e300}; + SpacetimeSource source = entry_source(&wt); + EntryEvaluator ev = {.camera_t = 0.0, + .camera_x = {50.0, 0.0, 0.0}, + .w = {-1.0, 0.0, 0.0}, + .L0 = 0.4, + .L0camera = 0.4}; + AsymptoticRoute out; + unsigned int evaluations = 0; + RayReason reason = RAY_REASON_COUNT; + + /* F == 0 exactly at the outside endpoint and strictly inside at 45. */ + CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 40.0, + 45.0, &out, &evaluations, + &reason) == ASYMPTOTIC_OK, + "exact boundary entry localize"); + CHECK(out.kind == ASYMPTOTIC_ROUTE_ENTRY && out.activate_t == -40.0 && + out.x[0] == 10.0, + "boundary endpoint returned directly"); + CHECK(evaluations == 2, "boundary path needs no bisection"); +} + +static void test_unconfirmed_bracket(void) { + EntryWorldtube wt = {.radius0 = 10.0, .valid_t_min = -1.0e300}; + SpacetimeSource source = entry_source(&wt); + EntryEvaluator ev = {.camera_t = 0.0, + .camera_x = {50.0, 0.0, 0.0}, + .w = {-1.0, 0.0, 0.0}}; + AsymptoticRoute out; + unsigned int evaluations = 0; + RayReason reason = RAY_REASON_COUNT; + + /* Both endpoints outside: no strict-inside bracket. */ + CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 10.0, + 20.0, &out, &evaluations, + &reason) == ASYMPTOTIC_INVALID && + reason == RAY_REASON_ENTRY_UNCONFIRMED, + "false candidate outside bracket is unconfirmed, not escaped"); + + /* Both endpoints strictly inside: also no entry bracket. */ + CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 45.0, + 50.0, &out, &evaluations, + &reason) == ASYMPTOTIC_INVALID && + reason == RAY_REASON_ENTRY_UNCONFIRMED, + "both-inside bracket is unconfirmed"); + + /* Reversed bracket ordering. */ + CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 45.0, + 30.0, &out, &evaluations, + &reason) == ASYMPTOTIC_INVALID && + reason == RAY_REASON_INVALID_ARGUMENT, + "reversed bracket rejected"); +} + +static void test_callback_failure_propagation(void) { + EntryEvaluator ev = {.camera_t = 0.0, + .camera_x = {50.0, 0.0, 0.0}, + .w = {-1.0, 0.0, 0.0}}; + AsymptoticRoute out; + unsigned int evaluations = 0; + RayReason reason = RAY_REASON_COUNT; + + EntryWorldtube fail = {.radius0 = 10.0, + .valid_t_min = -1.0e300, + .callback_fails = 1}; + SpacetimeSource source = entry_source(&fail); + CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 30.0, + 45.0, &out, &evaluations, + &reason) == ASYMPTOTIC_INVALID && + reason == RAY_REASON_WORLDTUBE_SAMPLE_FAILED, + "endpoint callback failure propagates"); + + fail = (EntryWorldtube){.radius0 = 10.0, + .valid_t_min = -1.0e300, + .fail_at_call = 2}; + source = entry_source(&fail); + CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 30.0, + 45.0, &out, &evaluations, + &reason) == ASYMPTOTIC_INVALID && + reason == RAY_REASON_WORLDTUBE_SAMPLE_FAILED, + "inside endpoint callback failure propagates"); + + fail = (EntryWorldtube){.radius0 = 10.0, + .valid_t_min = -1.0e300, + .fail_at_call = 3}; + source = entry_source(&fail); + CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 30.0, + 45.0, &out, &evaluations, + &reason) == ASYMPTOTIC_INVALID && + reason == RAY_REASON_WORLDTUBE_SAMPLE_FAILED, + "midpoint callback failure propagates"); +} + +static void test_history_hole_propagation(void) { + /* The inside endpoint falls past the valid history: the driver must report + * TIME_RANGE_EXHAUSTED, never a miss or a fabricated entry. */ + EntryWorldtube wt = {.radius0 = 10.0, .valid_t_min = -40.0}; + SpacetimeSource source = entry_source(&wt); + EntryEvaluator ev = {.camera_t = 0.0, + .camera_x = {50.0, 0.0, 0.0}, + .w = {-1.0, 0.0, 0.0}}; + AsymptoticRoute out; + unsigned int evaluations = 0; + RayReason reason = RAY_REASON_COUNT; + + CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 30.0, + 45.0, &out, &evaluations, + &reason) == ASYMPTOTIC_TIME_RANGE_EXHAUSTED && + reason == RAY_REASON_TIME_RANGE_EXHAUSTED, + "endpoint history hole propagates"); + + /* A midpoint-only history hole: both endpoints are valid, but the first + * bisection midpoint (t = -37.5) falls in an isolated invalid window. */ + wt = (EntryWorldtube){.radius0 = 10.0, + .valid_t_min = -1.0e300, + .hole_center = -37.5, + .hole_halfwidth = 0.5}; + source = entry_source(&wt); + CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 30.0, + 45.0, &out, &evaluations, + &reason) == ASYMPTOTIC_TIME_RANGE_EXHAUSTED && + reason == RAY_REASON_TIME_RANGE_EXHAUSTED, + "midpoint history hole propagates, never a miss"); +} + +static void test_evaluator_failure_propagation(void) { + EntryWorldtube wt = {.radius0 = 10.0, .valid_t_min = -1.0e300}; + SpacetimeSource source = entry_source(&wt); + AsymptoticRoute out; + unsigned int evaluations = 0; + RayReason reason = RAY_REASON_COUNT; + + EntryEvaluator ev = {.camera_t = 0.0, + .camera_x = {50.0, 0.0, 0.0}, + .w = {-1.0, 0.0, 0.0}, + .evaluator_fails_at = 1, + .fail_status = ASYMPTOTIC_TIME_RANGE_EXHAUSTED}; + CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 30.0, + 45.0, &out, &evaluations, + &reason) == ASYMPTOTIC_TIME_RANGE_EXHAUSTED && + reason == RAY_REASON_TIME_RANGE_EXHAUSTED, + "evaluator history failure propagates"); + + ev = (EntryEvaluator){.camera_t = 0.0, + .camera_x = {50.0, 0.0, 0.0}, + .w = {-1.0, 0.0, 0.0}, + .evaluator_fails_at = 1, + .fail_status = ASYMPTOTIC_INVALID}; + CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 30.0, + 45.0, &out, &evaluations, + &reason) == ASYMPTOTIC_INVALID && + reason == RAY_REASON_ENTRY_UNCONFIRMED, + "evaluator invalid failure maps to unconfirmed"); + + ev = (EntryEvaluator){.camera_t = 0.0, + .camera_x = {50.0, 0.0, 0.0}, + .w = {-1.0, 0.0, 0.0}, + .evaluator_fails_at = 2, + .fail_status = ASYMPTOTIC_INVALID}; + CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 30.0, + 45.0, &out, &evaluations, + &reason) == ASYMPTOTIC_INVALID, + "inside endpoint evaluator failure propagates"); +} + +static AsymptoticStatus wide_parameter_path(void *context, double parameter, + AsymptoticRoute *state) { + (void)context; + *state = (AsymptoticRoute){.kind = ASYMPTOTIC_ROUTE_ENTRY, + .end_id = 0, + .activate_t = -parameter, + .x = {2.0 - parameter * 1e-308, 0.0, 0.0}, + .Pi = {1.0, 0.0, 0.0}}; + return ASYMPTOTIC_OK; +} + +static void test_representability_and_state_checks(void) { + EntryWorldtube wt = {.radius0 = 1.0, .valid_t_min = -DBL_MAX}; + SpacetimeSource source = entry_source(&wt); + AsymptoticRoute out; + unsigned int evaluations = 0; + RayReason reason = RAY_REASON_COUNT; + /* Both endpoints are finite, but subtracting them overflows. This must not + * be mistaken for an adjacent bracket and return the far-inside endpoint. */ + CHECK(asymptotic_entry_localize(&source, 0, wide_parameter_path, NULL, + -1.6e308, 1.6e308, &out, &evaluations, + &reason) == ASYMPTOTIC_OK, + "overflow-safe parameter midpoint"); + CHECK(fabs(out.x[0] - 1.0) < 1e-14 && evaluations > 2, + "wide bracket contracts to entry, not initial inside endpoint"); + wt.radius0 = 10.0; + EntryEvaluator ev = {.camera_x = {50.0, 0.0, 0.0}, + .w = {-1.0, 0.0, 0.0}, .nonfinite_at = 1}; + CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 40.0, + 45.0, &out, &evaluations, &reason) == + ASYMPTOTIC_INVALID && reason == RAY_REASON_ENTRY_UNCONFIRMED, + "boundary shortcut rejects nonfinite camera energy reference"); + ev.evaluator_call_count = 0; + ev.nonfinite_at = 3; + CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 35.0, + 45.0, &out, &evaluations, &reason) == + ASYMPTOTIC_INVALID && reason == RAY_REASON_ENTRY_UNCONFIRMED, + "midpoint rejects nonfinite camera energy reference"); + ev.evaluator_call_count = 0; + ev.nonfinite_at = 0; + ev.reversed_time_at = 3; + CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 35.0, + 45.0, &out, &evaluations, &reason) == + ASYMPTOTIC_INVALID && reason == RAY_REASON_ENTRY_UNCONFIRMED, + "midpoint cannot reverse coordinate time"); +} + +int main(void) { + test_geometry_contract(); + test_validate_contract(); + test_fixed_sphere_localize(); + test_too_early_hint(); + test_moving_sphere_localize(); + test_radius_rate_localize(); + test_rotated_frame_localize(); + test_grazing_first_entry(); + test_curved_arc_localize(); + test_boundary_entry_exact(); + test_unconfirmed_bracket(); + test_callback_failure_propagation(); + test_history_hole_propagation(); + test_evaluator_failure_propagation(); + test_representability_and_state_checks(); + if (failures == 0) + puts("asymptotic entry regression passed"); + else + fprintf(stderr, "%d asymptotic entry regression failures\n", failures); + return failures == 0 ? 0 : 1; +} diff --git a/tests/test_asymptotic_quadratic.c b/tests/test_asymptotic_quadratic.c new file mode 100644 index 0000000..2e5a481 --- /dev/null +++ b/tests/test_asymptotic_quadratic.c @@ -0,0 +1,65 @@ +/* Exercise the private numerical kernel directly, including coefficient + * ranges that cannot be represented by a public double worldtube fixture. + * The build rule omits the separately compiled asymptotic.c. */ +#include "../src/asymptotic.c" + +#include + +static int failures; + +#define CHECK(condition, message) do { \ + if (!(condition)) { \ + fprintf(stderr, "FAIL %s:%d: %s\n", __FILE__, __LINE__, message); \ + ++failures; \ + } \ +} while (0) + +static void check_scaled(int exponent) { + const long double scale = scalbnl(1.0L, exponent); + double root = -1.0; + EntryQuadratic k = {scale, -3.0L * scale, 2.0L * scale}; + CHECK(entry_solve(&k, &root) == ENTRY_SOLVE_ENTRY && root == 1.0, + "common scale preserves smallest inward root"); + k = (EntryQuadratic){scale, -scale, scale}; + CHECK(entry_solve(&k, &root) == ENTRY_SOLVE_MISS, + "common scale preserves a clear miss"); + k = (EntryQuadratic){0.0L, -scale, 2.0L * scale}; + CHECK(entry_solve(&k, &root) == ENTRY_SOLVE_ENTRY && root == 2.0, + "common scale preserves linear entry"); + k = (EntryQuadratic){scale, -scale, 0.0L}; + CHECK(entry_solve(&k, &root) == ENTRY_SOLVE_ENTRY && root == 0.0, + "common scale preserves boundary entry"); + k = (EntryQuadratic){-scale, scale, 2.0L * scale}; + CHECK(entry_solve(&k, &root) == ENTRY_SOLVE_ENTRY && root == 2.0, + "common scale preserves concave entry"); +} + +static void test_product_cancellation(void) { + const long double u = scalbnl(1.0L, 1 - LDBL_MANT_DIG); + const EntryQuadratic k = {1.0L + u, -2.0L, 1.0L - 0.5L * u}; + long double scale; + (void)entry_discriminant(&k, &scale); + /* Exact dyadic oracle: 4 - 4(1+u)(1-u/2) = -2u + 2u^2. + * A separately rounded 4*a*c is 4 and loses this nonzero discriminant. */ + double root; + CHECK(entry_solve(&k, &root) == ENTRY_SOLVE_UNCERTAIN, + "product cancellation must remain uncertain, not a proven miss"); +} + +int main(void) { + check_scaled(0); + check_scaled(LDBL_MAX_EXP - 4); + check_scaled(LDBL_MIN_EXP + 4); + check_scaled(LDBL_MIN_EXP - LDBL_MANT_DIG + 2); + test_product_cancellation(); + double root; + EntryQuadratic k = {LDBL_MIN, LDBL_MAX / 8.0L, 1.0L}; + CHECK(entry_solve(&k, &root) == ENTRY_SOLVE_UNCERTAIN, + "scaling cannot silently erase a nonzero coefficient"); + k = (EntryQuadratic){1.0L, 2.0L, -INFINITY}; + CHECK(entry_solve(&k, &root) == ENTRY_SOLVE_UNCERTAIN, + "nonfinite coefficient is not a normal entry"); + if (!failures) + puts("asymptotic quadratic regression passed"); + return failures ? 1 : 0; +} diff --git a/tests/test_asymptotic_schwarzschild.c b/tests/test_asymptotic_schwarzschild.c index 8754aa1..8757f88 100644 --- a/tests/test_asymptotic_schwarzschild.c +++ b/tests/test_asymptotic_schwarzschild.c @@ -6,6 +6,7 @@ #include #include +#include static int failures = 0; #define CHECK(condition, message) \ @@ -189,6 +190,8 @@ static void test_preroute_entry(void) { ASYMPTOTIC_OK && route.kind == ASYMPTOTIC_ROUTE_ENTRY, "outside camera enters"); + CHECK(route.entry_fallback_evaluations == 0, + "analytic entry stays on the fast path"); double value; CHECK(asymptotic_worldtube_value(&source, route.end_id, route.activate_t, route.x, &value) == 0 && @@ -521,6 +524,122 @@ static void test_preroute_branches(void) { spacetime_destroy(&source); } +/* Translated-origin wrapper around the analytic Schwarzschild KS source: the + * inner metric is evaluated at x - origin and the worldtube/end are shifted by + * the same origin. This models a black hole at a large coordinate offset; the + * double reconstruction of the boundary state loses the sub-ULP offset and + * trips the common entry fallback, while the exact inward orbit transfer stays + * valid. It exercises the production fallback path, not a synthetic + * nonlinearity. */ +typedef struct { + SpacetimeSource inner; + double origin[3]; +} ShiftedOriginContext; + +static SpacetimePointStatus shifted_origin_eval(const SpacetimeSource *source, + double t, const double x[3], + MetricData *metric) { + const ShiftedOriginContext *ctx = source->context; + const double local[3] = {x[0] - ctx->origin[0], x[1] - ctx->origin[1], + x[2] - ctx->origin[2]}; + return spacetime_eval(&ctx->inner, t, local, metric); +} + +static SpacetimeRayStatus shifted_origin_classify(const SpacetimeSource *source, + double t, + const double x[3]) { + const ShiftedOriginContext *ctx = source->context; + const double local[3] = {x[0] - ctx->origin[0], x[1] - ctx->origin[1], + x[2] - ctx->origin[2]}; + return spacetime_classify(&ctx->inner, t, local); +} + +static size_t shifted_origin_end_count(const SpacetimeSource *source) { + const ShiftedOriginContext *ctx = source->context; + return spacetime_asymptotic_end_count(&ctx->inner); +} + +static int shifted_origin_end(const SpacetimeSource *source, size_t index, + SpacetimeAsymptoticEnd *out) { + const ShiftedOriginContext *ctx = source->context; + if (spacetime_asymptotic_end(&ctx->inner, index, out)) + return -1; + for (int i = 0; i < 3; ++i) + out->frame_origin[i] = ctx->origin[i]; + return 0; +} + +static int shifted_origin_worldtube(const SpacetimeSource *source, + SpacetimeEndId end_id, double t, + SpacetimeEscapeWorldtubeSample *out) { + const ShiftedOriginContext *ctx = source->context; + if (spacetime_escape_worldtube_sample(&ctx->inner, end_id, t, out)) + return -1; + for (int i = 0; i < 3; ++i) + out->center[i] += ctx->origin[i]; + return 0; +} + +static void shifted_origin_destroy(SpacetimeSource *source) { + ShiftedOriginContext *ctx = source->context; + if (ctx != NULL) { + spacetime_destroy(&ctx->inner); + free(ctx); + } + source->context = NULL; + source->ops = NULL; +} + +static const SpacetimeOps shifted_origin_ops = { + .eval = shifted_origin_eval, + .classify = shifted_origin_classify, + .asymptotic_end_count = shifted_origin_end_count, + .asymptotic_end = shifted_origin_end, + .escape_worldtube_sample = shifted_origin_worldtube, + .destroy = shifted_origin_destroy}; + +static void test_translated_origin_fallback(void) { + ShiftedOriginContext *ctx = malloc(sizeof *ctx); + CHECK(ctx != NULL, "shifted-origin context"); + if (ctx == NULL) + return; + ctx->origin[0] = 1.0e6; + ctx->origin[1] = 2.0e6; + ctx->origin[2] = -3.0e6; + CHECK(spacetime_create_schwarzschild_ks(&ctx->inner, 1.0, 256.0) == 0, + "shifted-origin inner source"); + SpacetimeSource source = {.ops = &shifted_origin_ops, .context = ctx}; + + ObserverCamera cam = {.look_ra_deg = 180.0, .look_dec_deg = 0.0}; + for (int i = 0; i < 3; ++i) + cam.position[i] = ctx->origin[i]; + cam.position[0] += 500.0; + MetricData metric; + CHECK(spacetime_eval(&source, 0.0, cam.position, &metric) == 0, + "shifted-origin camera metric"); + ObserverState observer; + CHECK(observer_from_coordinate_camera(&metric, &cam, &observer, NULL) == + OBSERVER_BUILD_OK, + "shifted-origin camera observer"); + const double direction[3] = {cos(0.3), sin(0.3), 0.0}; + AsymptoticRoute route; + CHECK(asymptotic_route_camera(&source, &observer, direction, &route) == + ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ENTRY, + "shifted-origin entry found"); + CHECK(route.entry_fallback_evaluations > 0, + "shifted-origin entry used the common fallback"); + CHECK(route.failure_reason == RAY_REASON_NONE, + "shifted-origin fallback has no failure reason"); + double value; + CHECK(asymptotic_worldtube_value(&source, route.end_id, route.activate_t, + route.x, &value) == 0 && + value <= 0.0, + "shifted-origin fallback state is inside the worldtube"); + CHECK(route.activate_t < 0.0, "shifted-origin entry is in the past"); + spacetime_destroy(&source); +} + int main(void) { test_round_trip(); test_finish_matches_integration(); @@ -532,6 +651,7 @@ int main(void) { test_time_reference(); test_grazing_reference(); test_preroute_branches(); + test_translated_origin_fallback(); if (failures == 0) puts("asymptotic schwarzschild regression passed"); else diff --git a/tests/test_geodesic.c b/tests/test_geodesic.c index 8089e11..91d6702 100644 --- a/tests/test_geodesic.c +++ b/tests/test_geodesic.c @@ -53,6 +53,8 @@ static int check_reason_names(void) { RAY_REASON_INTEGRATION_ERROR || ray_reason_category(RAY_REASON_INVALID_ESCAPE_DIRECTION) != RAY_REASON_INTEGRATION_ERROR || + ray_reason_category(RAY_REASON_ENTRY_UNCONFIRMED) != + RAY_REASON_INTEGRATION_ERROR || ray_reason_category(RAY_REASON_SLAB_LOAD_FAILED) != RAY_REASON_IO_ERROR) { fputs("detail reasons map to the wrong coarse category\n", stderr); failed = 1;