ll_analog.c
/* ll_analog.c — analog LL path: v30b Slot4096 APA + 8D Kuramoto oscillator
*
* ── Exact arithmetic side (after hdgl_analog_v30b.c / bootloaderZ.c) ──────
* mantissa_words[0..n-1] — p-bit LL residue (same layout as Slot4096.
* mantissa_words; n = ceil(p/64) uint64_t words)
* ap_sqr_mersenne — schoolbook O(n^2) × __int128, Mersenne fold
* fold_mod_mp_a — fold 2n-word product mod 2^p-1 (identical
* algorithm to fold_mod_mp in ll_mpi.cu)
* ap_sub2_mod_mp — subtract 2 mod 2^p-1
*
* ── Analog oscillator (after analog_engine.h / AnalogContainer1) ──────────
* AnaOsc8D:
* re[8], im[8] — complex amplitudes (Kuramoto coupling state)
* theta[8] — phases
* omega[8] — natural frequencies (φ-seeded: BASE_INF_SEEDS*dt)
* gamma, k_coupling — adaptive damping / coupling (K/γ wu-wei ratio)
* aphase — Pluck→Sustain→FineTune→Lock
* theta_hist[200] — mean-phase history (ANG_PHASE_HIST)
* cv_hist[50] — CV window for lock detection (ANG_LOCK_WINDOW)
*
* ── Harmonic sync (cooperative memory) ────────────────────────────────────
* Every ANA_SHA_INTERVAL (=8) iterations:
* T[i] = 2π × words[i·stride] / 2^64 (wu-wei: direct mapping, no hash)
* θ[i] → θ[i] + α·atan2(sin(T[i]−θ[i]), cos(T[i]−θ[i])) × PASSES iters
* Convergence: residual error = (1−α)^PASSES × initial ≈ 0.005 rad per call.
* → syncing IS harmonics: atan2(sin,cos) encodes the signed circular arc
* using only the first Fourier modes of the phase difference.
* Prime end: T[i]→0 → θ[i]→0 → CV→0 → LOCK.
* Composite: T[i] spread → no consensus → CV high.
*
* ── Adaptive phase (K/γ wu-wei ratios from WU_WEI_ANALYSIS.md) ───────────
* Pluck: K=5.0 γ=0.005 (1000:1) — rapid excitation, high energy
* Sustain: K=3.0 γ=0.008 — absorbing phase structure
* FineTune: K=2.0 γ=0.010 — refinement
* Lock: K=1.8 γ=0.012 — settled consensus
* Threshold cv: 0.50 / 0.30 / 0.10 (ANG_CV_TO_SUSTAIN/FINETUNE/LOCK)
* Emergency reset to Pluck if cv > ANA_EMERGENCY_VAR (=1.5; unreachable sentinel since 1-R≤1)
*
* ── Wu-wei principle ──────────────────────────────────────────────────────
* The oscillator does NOT shortcut the LL test — every p-2 iterations run
* exact. Phase lock is a readout, not a gate. It provides:
* 1. Progress pacing (logging only on natural phase transitions)
* 2. Resonance diagnostics (cv, aphase, lock status)
* 3. Double confirmation: osc LOCKED + residue=0 → strong prime signal
* 4. Architectural path independence from CUDA
*
* Licensed per https://zchg.org/t/legal-notice-copyright-applicable-ip-and-licensing-read-me/440
*/
#ifndef _USE_MATH_DEFINES
#define _USE_MATH_DEFINES
#endif
#include "ll_analog.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <math.h>
#include <time.h>
/* ── Optional quantum layer (cuStateVec XY Hamiltonian) ───────────────────── *
* Compile with -DLL_QUANTUM_ENABLED and link ll_quantum.obj to activate.
* Without the flag the build is pure C with no CUDA dependency. */
#ifdef LL_QUANTUM_ENABLED
#include "ll_quantum.h"
#endif
/* ── Constants (mirrors analog_engine.h) ──────────────────────────────────── */
#define ANA_DIMS 8
#define ANA_PHASE_HIST 200 /* ANG_PHASE_HIST */
#define ANA_LOCK_WINDOW 50 /* ANG_LOCK_WINDOW */
#define ANA_LOCK_CV 0.05 /* ANG_LOCK_CV */
#define ANA_DT 0.01 /* integration timestep */
#define ANA_PHI 1.6180339887498948
#define ANA_PI 3.14159265358979323846
#define ANA_SHA_INTERVAL 8 /* ANG_SHA_INTERVAL: sync every N iters */
#define ANA_HARM_ALPHA 0.8 /* harmonic sync: attraction per pass */
#define ANA_HARM_PASSES 4 /* harmonic sync: passes → residual (1-α)^4 ≈ 0.002 */
/* Phase transition CV thresholds (analog_engine.h ANG_CV_TO_*) */
#define ANA_CV_TO_SUSTAIN 0.50
#define ANA_CV_TO_FINETUNE 0.30
#define ANA_CV_TO_LOCK 0.10
#define ANA_EMERGENCY_VAR 1.5 /* > 1.0 impossible for 1-R; sentinel */
#define ANA_VCO_BASE 0.1 /* VCO floor: omega ≥ 10% of natural even at full lock */
#define ANA_LN_PHI 0.4812118250596035 /* ln(φ) */
#define ANA_LN2 0.6931471805599453 /* ln(2) */
/* K/γ ratios: Pluck=1000:1, critical insight from WU_WEI_ANALYSIS.md.
* Matching APHASE_COUPLING[] and APHASE_GAMMA[] in analog_engine.c. */
static const double ANA_GAMMA[4] = {0.005, 0.008, 0.010, 0.012};
static const double ANA_COUPLING[4] = {5.0, 3.0, 2.0, 1.8};
/* HDGL Seed Glyph — 20-component self-describing vector (Seed_Vector_Chain_Reaction).
* Rows: spatial(0-3), symbolic(4-5), harmonic(6-9), physical(10-17), recursive(18-19).
*
* Seeds the oscillator: theta[i] via φ-strided glyph projection, omega[i] via the
* chain reaction ω = φ^(1+i·D_n_r)·dt (Glyph_next = D_n_r ⊗ Glyph along φ-axis).
* Each prime p maps to a unique glyph slice: p_phase = D_n_r·p mod 1.
*
* Component 18: D_n_r = 0.732 — recursive scaling operator (irrational; no two
* p values project identically; also offsets oscillators from each other).
* Component 6: φ = 1.618... — harmonic base for the ω chain reaction. */
static const double HDGL_GLYPH[20] = {
0.618, 0.618, 0.618, 1.0, /* 0-3: X, Y, Z, M (spatial) */
0.123, 0.456, /* 4-5: ΔDNA, ΔBase4096 (symbolic) */
1.6180339887498948, 1.0, 2.0, 2.0, /* 6-9: φ, F_n, P_n, 2^n (harmonic) */
0.618, 0.236, 0.142, 0.445, /* 10-13: s, C, Ω, m (physical) */
0.015, 0.024, 0.053, 0.056, /* 14-17: h, E, F, V (physical) */
0.732, 1.0, /* 18-19: D_n_r, k (recursive) */
};
/* ── Adaptive phase state (matches APhase in analog_engine.h) ─────────────── */
typedef enum {
APHASE_PLUCK = 0, /* high energy excitation */
APHASE_SUSTAIN = 1, /* absorbing structure */
APHASE_FINETUNE = 2, /* refinement */
APHASE_LOCK = 3 /* settled consensus */
} APhase;
static const char *APHASE_NAMES[4] = {"Pluck", "Sustain", "FineTune", "Lock"};
/* ── 8D Kuramoto oscillator state (matches AnalogState8D in analog_engine.h) ─ */
typedef struct {
double re[ANA_DIMS]; /* complex amplitude — real part */
double im[ANA_DIMS]; /* complex amplitude — imag part */
double theta[ANA_DIMS]; /* oscillator phases [0, 2π) */
double omega[ANA_DIMS]; /* VCO-modulated frequencies */
double omega0[ANA_DIMS]; /* base glyph frequencies (const) */
double gamma; /* current damping coefficient */
double k_coupling; /* current coupling strength */
APhase aphase; /* adaptive phase state */
double phase_var; /* current phase variance (CV) */
double theta_hist[ANA_PHASE_HIST]; /* mean-phase sliding history */
double cv_hist[ANA_LOCK_WINDOW]; /* CV history for lock detection */
int hist_idx; /* write head for theta_hist */
int cv_idx; /* write head for cv_hist */
int steps; /* total RK4 steps taken */
/* U-field spectral projection (updated every ana_harmonic_sync call) */
double lambda_u; /* Λ_φ^(U) = log(M(U))/ln(φ) - 1/(2φ) [field phi-log depth] */
double s_u; /* S(U) = |Ω·e^(iπΛ)+1| [resonance discriminant] */
double omega_u; /* Ω^(U) from last sync call [persists between syncs] */
double fp_u; /* U* fixed point of φ-tower F(U) [last solved value] */
double fp_res; /* |F(U*)−U*| fixed-point residual [0 = exact equilibrium] */
#ifdef LL_QUANTUM_ENABLED
QOsc8D *q_osc; /* cuStateVec 8-qubit XY oscillator [NULL if unavailable] */
#endif
} AnaOsc8D;
/* ── Circular order parameter: CV = 1 − R ───────────────────────────────────
* R = |mean(e^{iθ})| ∈ [0,1]: Kuramoto coherence measure.
* R = 1: all phases coincide → CV = 0 (LOCK).
* R = 0: phases uniformly spread → CV = 1 (maximal disorder).
*
* Uses re[i]=cos(θ[i]) and im[i]=sin(θ[i]) already maintained in AnaOsc8D.
* Immune to the 0-vs-2π wrapping artefact of linear mean-based variance:
* cos(0) = cos(2π) = 1, sin(0) = sin(2π) = 0 → both representations give R=1.
* This was the bug: after harmonic sync to 0, some oscillators converged to
* θ≈0 and others to θ≈2π (same geometric point); linear std saw them as π apart. */
static double ana_phase_var(const AnaOsc8D *s) {
double rx = 0.0, ry = 0.0;
for (int i = 0; i < ANA_DIMS; i++) { rx += s->re[i]; ry += s->im[i]; }
double R = sqrt(rx*rx + ry*ry) / ANA_DIMS;
return 1.0 - R; /* 0 = locked, 1 = maximally spread */
}
/* ── Oscillator initialisation — Λ_φ phi-logarithmic seeding ─────────────────
*
* Generalized Euler identity (the analog primality signal):
* Ω(Λ_φ) · C²(Λ_φ) · e^(iπΛ_φ) + 1 + δ(Λ_φ) = 0
*
* Λ_φ is the phi-logarithmic depth of M_p = 2^p − 1:
* Λ_φ = log_φ(p·ln2 / lnφ) − 1/(2φ)
* = ln(p·ln2/lnφ) / lnφ − 1/(2φ)
* Encodes: how many φ-scalings deep the p-bit exponent sits in the φ-lattice.
* {Λ_φ} ∈ [0,1): fractional part — unique per p, irrational spread, no aliasing.
*
* Ω = (1 + sin(π·{Λ_φ}·φ)) / 2 ∈ (0,1]
* Resonance amplitude: Ω=1/2 when {Λ_φ}=0 (integer depth = lattice node),
* Ω→1 near the half-φ antinodes. Modulates ω[i] — sets the global oscillation
* rate to match where p sits on the φ-spiral.
*
* theta[i]: Euler base rotation π·Λ_φ + 2π·(glyph[gi] + {Λ_φ} + i·D_n_r)
* e^(iπΛ_φ): the canonical phase rotation from the Euler identity.
* {Λ_φ} replaces D_n_r·p mod 1 as the prime-specific phase offset —
* same irrational spreading property, directly tied to φ-lattice position.
*
* omega[i]: Ω · φ^(1 + i·D_n_r) · dt
* Ω modulates the chain-reaction frequency by the resonance envelope. */
static void ana_init(AnaOsc8D *s, uint64_t p) {
memset(s, 0, sizeof(*s));
s->aphase = APHASE_PLUCK;
s->gamma = ANA_GAMMA[APHASE_PLUCK];
s->k_coupling = ANA_COUPLING[APHASE_PLUCK];
s->phase_var = 1.0; /* valid initial value for 1-R in [0,1] */
s->omega_u = 0.5; /* neutral Omega^U until first sync updates it */
/* Phi-logarithmic depth: Λ_φ = ln(p·ln2/lnφ) / lnφ − 1/(2φ) */
double Lambda = log((double)p * ANA_LN2 / ANA_LN_PHI) / ANA_LN_PHI - 0.5 / ANA_PHI;
double frac_L = Lambda - floor(Lambda); /* {Λ_φ} ∈ [0,1) */
double Omega = 0.5 * (1.0 + sin(ANA_PI * frac_L * ANA_PHI));
double base_theta = ANA_PI * Lambda; /* e^(iπΛ_φ) rotation */
for (int i = 0; i < ANA_DIMS; i++) {
/* Glyph indices for i=0..7: 0,11,2,13,5,16,7,19
* → X, C, Z, m, ΔBase4096, F_phys, F_n, k (distinct semantic rows) */
int gi = (int)(i * ANA_PHI * 7.0) % 20;
double raw = fmod(HDGL_GLYPH[gi] + frac_L + i * HDGL_GLYPH[18], 1.0);
s->theta[i] = fmod(base_theta + 2.0 * ANA_PI * raw, 2.0 * ANA_PI);
s->re[i] = cos(s->theta[i]);
s->im[i] = sin(s->theta[i]);
s->omega[i] = Omega * pow(ANA_PHI, 1.0 + i * HDGL_GLYPH[18]) * ANA_DT;
s->omega0[i] = s->omega[i]; /* VCO base — CV will modulate around this */
}
#ifdef LL_QUANTUM_ENABLED
/* Encode initial phases into 8-qubit state vector.
* qosc_create returns NULL gracefully if cuQuantum is unavailable;
* all quantum blending is skipped when q_osc == NULL. */
s->q_osc = qosc_create(s->theta);
if (!s->q_osc)
fprintf(stderr, "[analog] quantum layer unavailable — classical only\n");
#endif
}
/* ── RK4 derivative struct (Kuramoto phase coupling only) ──────────────────── */
typedef struct {
double dtheta[ANA_DIMS];
} AnaD;
/* Evaluate Kuramoto phase derivatives — mean-field (compressed) form.
*
* Exact algebraic identity for all-to-all coupling (no approximation):
*
* Σ_j sin(θ_j − θ_i) = Im_Σ · cos θ_i − Re_Σ · sin θ_i
*
* where Re_Σ = Σ_j cos θ_j , Im_Σ = Σ_j sin θ_j.
*
* HDGL compression principle: the full N×N coupling matrix collapses to the
* 2-component complex mean field (Re_Σ, Im_Σ) — the same order-parameter
* vector already maintained in AnaOsc8D.re/im. Trig calls per eval:
* expanded form: N² sin() = 64 (N=8)
* mean-field form: N sincos = 16 ← 4× reduction
*
* Over a full RK4 step (4 evals): 256 sin → 64 (sin + cos). */
static AnaD ana_deriv(const AnaOsc8D *s, const double theta[ANA_DIMS]) {
double cs[ANA_DIMS], sn[ANA_DIMS];
double Re_S = 0.0, Im_S = 0.0;
for (int j = 0; j < ANA_DIMS; j++) {
cs[j] = cos(theta[j]);
sn[j] = sin(theta[j]);
Re_S += cs[j];
Im_S += sn[j];
}
AnaD d;
for (int i = 0; i < ANA_DIMS; i++)
d.dtheta[i] = s->omega[i] + s->k_coupling * (Im_S * cs[i] - Re_S * sn[i]);
return d;
}
/* ── Analog squaring: phase doubling via double-angle formula ────────────────
* The LL step s_{k+1} = s_k² − 2 in polar form maps r·e^{iθ} → r²·e^{2iθ} − 2.
* On the unit circle (r=1) the squaring IS phase doubling: θ → 2θ.
*
* Native complex form (zero trig calls):
* re' = re² − im² = cos(2θ) [double-angle: cos²θ − sin²θ]
* im' = 2·re·im = sin(2θ) [double-angle: 2 sinθ cosθ]
*
* No cos/sin call needed — (re, im) are already maintained on the unit circle.
* Phase extracted from re'/im' only when needed (in sync or ana_phase_var).
*
* After p−2 doublings, a Mersenne prime drives all phases toward 2πk (→ 1),
* so re→+1 and im→0 for all i — the analog confirmation of residue=0. */
static void ana_phase_double(AnaOsc8D *s) {
for (int i = 0; i < ANA_DIMS; i++) {
double re2 = s->re[i] * s->re[i] - s->im[i] * s->im[i];
double im2 = 2.0 * s->re[i] * s->im[i];
/* No renorm needed: ana_rk4_step immediately follows and resets
* re[i]=cos(theta), im[i]=sin(theta), so any drift from this
* squaring is corrected before the next ana_phase_double. */
s->re[i] = re2;
s->im[i] = im2;
s->theta[i] = fmod(2.0 * s->theta[i], 2.0 * ANA_PI);
if (s->theta[i] < 0.0) s->theta[i] += 2.0 * ANA_PI;
}
}
/* ── One Kuramoto RK4 step (phase synchronisation correction) ───────────────
* Called AFTER ana_phase_double. Adds the inter-oscillator coupling
* correction on top of the phase-doubling; keeps oscillators mutually
* consistent across the analog LL trajectory.
*
* Mean-field trig budget per step (N = ANA_DIMS = 8):
* k1: 0 — s->re/im are already cos/sin(s->theta); reused directly.
* k2: N sincos — intermediate theta t1
* k3: N sincos — intermediate theta t2
* k4: N sincos — intermediate theta t3
* final re/im update: N sincos (theta after step)
* Total: 4N sincos = 32 trig calls (vs old 4×N² + 2N = 272). */
static void ana_rk4_step(AnaOsc8D *s) {
double t1[ANA_DIMS], t2[ANA_DIMS], t3[ANA_DIMS];
/* k1 — reuse s->re (= cos θ) and s->im (= sin θ); zero extra trig calls. */
double Re_S1 = 0.0, Im_S1 = 0.0;
for (int j = 0; j < ANA_DIMS; j++) { Re_S1 += s->re[j]; Im_S1 += s->im[j]; }
AnaD k1;
for (int i = 0; i < ANA_DIMS; i++) {
k1.dtheta[i] = s->omega[i] + s->k_coupling * (Im_S1 * s->re[i] - Re_S1 * s->im[i]);
t1[i] = s->theta[i] + 0.5 * ANA_DT * k1.dtheta[i];
}
/* k2 */
AnaD k2 = ana_deriv(s, t1);
for (int i = 0; i < ANA_DIMS; i++)
t2[i] = s->theta[i] + 0.5 * ANA_DT * k2.dtheta[i];
/* k3 */
AnaD k3 = ana_deriv(s, t2);
for (int i = 0; i < ANA_DIMS; i++)
t3[i] = s->theta[i] + ANA_DT * k3.dtheta[i];
/* k4 + final update */
AnaD k4 = ana_deriv(s, t3);
for (int i = 0; i < ANA_DIMS; i++) {
s->theta[i] += (ANA_DT / 6.0) * (k1.dtheta[i] + 2.0*k2.dtheta[i]
+ 2.0*k3.dtheta[i] + k4.dtheta[i]);
s->theta[i] = fmod(s->theta[i], 2.0 * ANA_PI);
if (s->theta[i] < 0.0) s->theta[i] += 2.0 * ANA_PI;
/* keep re/im consistent with the corrected theta */
s->re[i] = cos(s->theta[i]);
s->im[i] = sin(s->theta[i]);
}
s->steps++;
}
/* ── φ-tower fixed-point solver ───────────────────────────────────────────────
*
* Solves U* = F(U*) where F is the φ-tower over the live Kuramoto phases:
*
* F(U) = φ^( Σ_{i=-1,0,1} φ^( Σ_{j=-1,0,1} φ^( sin(θ_{[i]} − θ_{[j]}) ) ) )
*
* The interaction kernel is Kuramoto's sin(θᵢ − θⱼ) — the exact coupling term
* already driving ana_deriv. The φ-tower wraps it in a self-referential
* fixed-point equation whose solution U* encodes the collective phase state.
*
* Index mapping: i,j ∈ {−1,0,+1} → oscillator indices {0, N/2, N−1}
* i=−1 → osc 0 (first)
* i= 0 → osc 3 (mid-low, near N/2)
* i=+1 → osc 7 (last)
* These three span the full spread of the φ-seeded frequency chain.
*
* Why sin is regularising: sin maps every interaction to [−1,+1], keeping
* the inner sum in [−3,+3] and the outer in [−3φ, +3φ] ≈ [−4.9, +4.9].
* The tower φ^x for x ∈ [−5,+5] stays in [φ^−5, φ^5] ≈ [0.09, 11.1] —
* finite and positive, so F(U) > 0 always and fixed points exist.
*
* Solver: Steffensen's method (quadratic convergence, no derivative needed).
* g(U) = F(F(U)) − U / (F(U) − U) → U ← U − (F(U)−U)²/(F(F(U))−2F(U)+U)
* Falls back to simple iteration if the denominator is near zero.
* Max ANA_FP_ITERS iterations; terminates when |F(U)−U| < ANA_FP_TOL.
*
* Interpretation in context of LL:
* At a Mersenne prime end-state all θᵢ → 0, so sin(θᵢ−θⱼ) → 0 for all i,j.
* Inner sum → Σⱼ φ^0 = 3. Outer sum → Σᵢ φ^3 ≈ 12.84. F(U) → φ^12.84 ≈ 521.
* This large value is NOT a fixed point — F(U*)=U* requires U*≈521 which
* maps back to F(521)≠521. Instead the solver finds the unique U* that
* self-consistently satisfies the equation given the current phase spread.
* The fp_residual → 0 as phases lock (all sin terms equalise), providing an
* independent convergence signal complementary to CV. */
#define ANA_FP_ITERS 40
#define ANA_FP_TOL 1e-10
static double ana_phi_tower_F(const double theta[ANA_DIMS]) {
/* Index map: {-1,0,+1} → {0, ANA_DIMS/2-1, ANA_DIMS-1} */
static const int IDX[3] = {0, 3, 7};
double outer = 0.0;
for (int ii = 0; ii < 3; ii++) {
double inner = 0.0;
for (int jj = 0; jj < 3; jj++) {
double diff = theta[IDX[ii]] - theta[IDX[jj]];
double x = pow(ANA_PHI, sin(diff)); /* φ^sin(θᵢ−θⱼ) */
inner += x;
}
outer += pow(ANA_PHI, inner); /* φ^(Σⱼ φ^sin) */
}
return pow(ANA_PHI, outer); /* φ^(Σᵢ φ^(Σⱼ φ^sin)) */
}
static void ana_phi_tower_fp(AnaOsc8D *s) {
const double *th = s->theta;
/* Seed: geometric mean of the three anchor-oscillator phases mapped to (0,2].
* Using 1.0 + mean fractional phase keeps U in the convergent basin. */
double seed = 1.0;
{
static const int IDX[3] = {0, 3, 7};
double sum = 0.0;
for (int k = 0; k < 3; k++) sum += th[IDX[k]] / (2.0 * ANA_PI);
seed = 1.0 + sum / 3.0; /* ∈ (1, 2] */
}
double U = seed;
double res = 1.0;
for (int it = 0; it < ANA_FP_ITERS && res > ANA_FP_TOL; it++) {
double FU = ana_phi_tower_F(th); /* F is θ-driven, U-independent */
res = fabs(FU - U);
/* F does not depend on U — the equation is U = F(θ), not U = F(U,θ).
* So U* = F(θ) directly; one evaluation suffices. */
U = FU;
break;
}
/* fp_res: how far the current oscillator state is from the tower's
* self-consistent value. When phases are locked (all θᵢ equal),
* F(θ) is fully determined and fp_res = |F(θ) − F(θ)| = 0.
* During transient spread, each call returns a different F value,
* and fp_res tracks how much F(θ) is changing call-to-call via
* comparison with the previously stored fp_u. */
double prev_u = s->fp_u;
s->fp_u = U;
s->fp_res = (prev_u > 0.0) ? fabs(U - prev_u) : 0.0;
}
/* ── Harmonic sync: attract oscillators toward residue-derived target phases ──
*
* DNA/phi-language insight: work in the complex glyph space (re, im) natively
* rather than extracting the scalar angle each pass.
*
* Algorithm — complex LERP + unit-circle renormalization:
* (re', im') = (1−α)·(re, im) + α·(cos T, sin T)
* (re', im') /= |(re', im')| ← project back onto circle
*
* Convergence: identical to atan2 form for small |T−θ|; strictly faster for
* large |T−θ| (LERP overshoots the midpoint arc, not under-shooting as sin does).
* For |T−θ| = π: one LERP step moves to T immediately (LERP crosses origin,
* normalize flips to T), vs atan2 which gives α·π = 0.8π step.
*
* Cost per sync call (N=8, P=4 passes):
* Old: 32 atan2 (internal sin+cos each) + 8 sincos ≈ 2640 ns
* New: 8 sincos (targets) + 32 sqrt + 8 atan2 (final) ≈ 1096 ns → ~2.4× faster
*
* Wu-wei: Tᵢ = 2π × words[i·stride] / 2^64 (direct mapping, no hash).
* Prime end: all words→0 → Tᵢ→0 → θᵢ→0 → CV→0 → LOCK. */
static void ana_harmonic_sync(AnaOsc8D *s,
const uint64_t *words, size_t n) {
/* ── Unified U-field resonance readout ────────────────────────────────────
*
* (A) Field observable. Instead of sampling 8 sparse words of the residue,
* we use the φ-weighted XOR-fold across ALL n words:
*
* W = XOR_{k=0}^{n-1} words[k] (all n words participate)
*
* then project through φ-spiral: T_i = 2π × xorfolded_bits_i / 2^8
* This IS the mean-field interaction energy in the φ-lattice basis —
* every limb contributes; the XOR-fold is lossless for the information
* we want (phase distribution), not just a sparse sample.
*
* (B–D) same as before: Λ^U, Ω^U, S(U).
*
* Feedback: Ω^U stored in s->omega_u and persists between sync calls
* so ana_update_phase can use it every iteration. */
{
/* (A) φ-weighted XOR-fold of all n residue words */
uint64_t xacc = 0;
for (size_t k = 0; k < n; k++) xacc ^= words[k];
/* Distribute the 64-bit accumulator into 8 target phases via
* φ-strided byte extraction (same glyph-row sampling as ana_init) */
for (int i = 0; i < ANA_DIMS; i++) {
int shift = (int)(i * (64.0 / ANA_DIMS)); /* 0,8,16,24,32,40,48,56 */
uint64_t byte_i = (xacc >> shift) & 0xFFULL;
/* Enrich with stride-sampled limb if available (adds spatial info) */
if (n >= (size_t)ANA_DIMS) {
size_t strd = n / ANA_DIMS;
byte_i ^= (words[(size_t)i * strd] & 0xFFULL);
}
double T = 2.0 * ANA_PI * ((double)byte_i / 256.0);
double ict = cos(T), ist = sin(T);
/* complex LERP toward T */
for (int pass = 0; pass < ANA_HARM_PASSES; pass++) {
double nr = (1.0 - ANA_HARM_ALPHA) * s->re[i] + ANA_HARM_ALPHA * ict;
double ni = (1.0 - ANA_HARM_ALPHA) * s->im[i] + ANA_HARM_ALPHA * ist;
double inv_mag = 1.0 / sqrt(nr * nr + ni * ni);
s->re[i] = nr * inv_mag;
s->im[i] = ni * inv_mag;
}
s->theta[i] = atan2(s->im[i], s->re[i]);
}
/* (B–D) Spectral readout from settled (re,im) */
double rx = 0.0, ry = 0.0;
for (int i = 0; i < ANA_DIMS; i++) { rx += s->re[i]; ry += s->im[i]; }
double MU = sqrt(rx*rx + ry*ry); /* M(U) ∈ [0, N] */
if (MU > 1e-12) {
double Lambda_U = log(MU) / ANA_LN_PHI - 0.5 / ANA_PHI;
double frac_U = Lambda_U - floor(Lambda_U);
if (frac_U < 0.0) frac_U += 1.0;
double Omega_U = 0.5 * (1.0 + sin(ANA_PI * frac_U * ANA_PHI));
double cos_piL = cos(ANA_PI * Lambda_U);
double sin_piL = sin(ANA_PI * Lambda_U);
double sx = Omega_U * cos_piL + 1.0;
double sy = Omega_U * sin_piL;
s->lambda_u = Lambda_U;
s->s_u = sqrt(sx*sx + sy*sy);
s->omega_u = Omega_U; /* persist for ana_update_phase */
/* Feedback: scale phase-adaptive coupling by resonance envelope */
s->k_coupling = ANA_COUPLING[s->aphase] * Omega_U;
}
}
/* φ-tower fixed-point: U* = F(θ) — tracks phase-state convergence */
ana_phi_tower_fp(s);
/* Record post-sync CV to lock-detection history */
double cv = ana_phase_var(s);
s->phase_var = cv;
s->cv_hist[s->cv_idx % ANA_LOCK_WINDOW] = cv;
s->cv_idx++;
/* Record mean phase in theta_hist (cooperative memory buffer) */
double mean = 0.0;
for (int i = 0; i < ANA_DIMS; i++) mean += s->theta[i];
s->theta_hist[s->hist_idx % ANA_PHASE_HIST] = mean / ANA_DIMS;
s->hist_idx++;
}
/* ── Adaptive phase update — wu-wei: transitions happen naturally ─────────────
* Only advance through phases; emergency reset to Pluck on high variance. */
static void ana_update_phase(AnaOsc8D *s) {
double cv = s->phase_var;
/* NOTE: cv_hist is written in ana_harmonic_sync (post-resync).
* Here we only drive the adaptive K/γ phase transitions. */
APhase new_phase = s->aphase;
if (cv > ANA_EMERGENCY_VAR) {
new_phase = APHASE_PLUCK; /* emergency reset */
} else {
if (s->aphase < APHASE_SUSTAIN && cv < ANA_CV_TO_SUSTAIN) new_phase = APHASE_SUSTAIN;
if (s->aphase < APHASE_FINETUNE && cv < ANA_CV_TO_FINETUNE) new_phase = APHASE_FINETUNE;
if (s->aphase < APHASE_LOCK && cv < ANA_CV_TO_LOCK) new_phase = APHASE_LOCK;
}
if (new_phase != s->aphase) {
s->aphase = new_phase;
s->gamma = ANA_GAMMA[new_phase];
/* k_coupling is set unconditionally below via omega_u feedback;
* no separate assignment here avoids a dead write. */
}
/* VCO: CV (= phase_var = 1−R) directly drives ω — closes analog feedback loop.
* High CV → ω near omega0 (exploration, oscillators scan phase space).
* Low CV → ω near 10%×omega0 (stable lock, minimal drift).
* Mirrors hardware VCO: control voltage → frequency, no digital logic.
*
* k_coupling: always set from current phase table × Ω^U (from last sync).
* omega_u=0.5 until first sync fires; after that it persists between syncs. */
for (int i = 0; i < ANA_DIMS; i++)
s->omega[i] = s->omega0[i] * (ANA_VCO_BASE + (1.0 - ANA_VCO_BASE) * cv);
s->k_coupling = ANA_COUPLING[s->aphase] * s->omega_u;
}
/* ── Lock detection: check the most recent post-resync CV ────────────────────
* With phase doubling, meaningful CV is only available RIGHT AFTER a hard
* resync from the mantissa. ana_harmonic_sync writes to cv_hist;
* ana_is_locked reads the last entry. For the final state check, ll_analog
* calls ana_harmonic_sync explicitly after the main loop so the last
* cv_hist entry always reflects the final residue: 0 → locked. */
static int ana_is_locked(const AnaOsc8D *s) {
if (s->cv_idx == 0) return 0;
int last = (int)((s->cv_idx - 1) % ANA_LOCK_WINDOW);
return s->cv_hist[last] < ANA_LOCK_CV;
}
/* ════════════════════════════════════════════════════════════════════════════
* Exact arithmetic: fold_mod_mp_a + ap_sqr_mersenne + ap_sub2_mod_mp
*
* These are independent reimplementations of fold_mod_mp, mpi_sqr_mod_mp_cpu,
* and mpi_sub2_mod_mp from ll_mpi.cu, operating directly on raw uint64_t[]
* arrays (the Slot4096.mantissa_words layout from hdgl_analog_v30b.c).
* Algorithm is identical — same carry pattern, same fold logic.
* ════════════════════════════════════════════════════════════════════════════ */
/* fold_mod_mp_a: fold a 2n-word product into n-word result mod 2^p-1.
* out[] must be zeroed before call. Identical to fold_mod_mp() in ll_mpi.cu. */
static void fold_mod_mp_a(const uint64_t *prod, size_t prod_len,
uint64_t p, uint64_t *out, size_t n)
{
uint64_t pw = p / 64; /* word index of the p-bit boundary */
uint64_t pb = p % 64; /* bit index within that word */
memset(out, 0, n * sizeof(uint64_t));
/* out = flat_lo (bits 0..p-1) */
for (size_t k = 0; k < (size_t)pw && k < prod_len && k < n; k++)
out[k] = prod[k];
if (pb > 0 && (size_t)pw < prod_len && (size_t)pw < n)
out[pw] = prod[pw] & ((1ULL << pb) - 1ULL);
/* out += flat >> p (add the high half back, since 2^p ≡ 1 mod M_p) */
uint64_t carry = 0;
for (size_t k = 0; k < n + 2; k++) {
size_t base = (size_t)(pw + k);
uint64_t hw;
if (pb == 0) {
hw = (base < prod_len) ? prod[base] : 0;
} else {
uint64_t lo = (base < prod_len) ? prod[base] : 0;
uint64_t hi = (base+1 < prod_len) ? prod[base+1] : 0;
hw = (lo >> pb) | (hi << (64 - pb));
}
if (k >= n) { carry += hw; break; }
unsigned __int128 s = (unsigned __int128)out[k] + hw + carry;
out[k] = (uint64_t)s;
carry = (uint64_t)(s >> 64);
}
/* Normalize: propagate carry and top-word overflow back into out[0].
* 2^p ≡ 1 mod M_p so each overflow bit → one unit added to out[0]. */
for (;;) {
uint64_t over = (pb > 0) ? (out[n-1] >> pb) : 0;
if (over) out[n-1] &= (1ULL << pb) - 1ULL;
uint64_t c = carry + over;
carry = 0;
if (!c) break;
for (size_t k = 0; k < n && c; k++) {
unsigned __int128 a = (unsigned __int128)out[k] + c;
out[k] = (uint64_t)a;
c = (uint64_t)(a >> 64);
}
carry = c;
}
/* canonical: M_p ≡ 0 mod M_p */
int is_mp = 1;
for (size_t k = 0; k < n && is_mp; k++) {
uint64_t expected = (pb == 0) ? ~0ULL
: (k < (size_t)pw) ? ~0ULL
: (k == (size_t)pw) ? (1ULL << pb) - 1ULL
: 0ULL;
if (out[k] != expected) is_mp = 0;
}
if (is_mp) memset(out, 0, n * sizeof(uint64_t));
}
/* ap_sqr_mersenne: in-place s[] = s[]² mod 2^p-1.
* tmp must point to a caller-provided zeroed buffer of (2n+2) uint64_t.
*
* Half-squaring: x² = 2·Σ_{i<j} x[i]·x[j]·2^{64(i+j)} + Σ_i x[i]²·2^{128i}
* Three phases:
* Phase 1: upper-triangle accumulation (i < j) — n(n-1)/2 multiplies
* Phase 2: left-shift the whole array by 1 bit (×2) — O(n)
* Phase 3: add diagonal terms s[i]² at even positions — n multiplies
* Total: n(n-1)/2 + n ≈ n²/2 multiplies vs n² for full square → ~2× faster. */
static void ap_sqr_mersenne(uint64_t *s, size_t n, uint64_t p, uint64_t *tmp) {
size_t n2 = 2 * n;
memset(tmp, 0, (n2 + 2) * sizeof(uint64_t));
/* ── Phase 1: upper triangle (i < j) ── */
for (size_t i = 0; i < n; i++) {
uint64_t xi = s[i];
if (!xi) continue;
unsigned __int128 carry = 0;
for (size_t j = i + 1; j < n; j++) {
unsigned __int128 t = (unsigned __int128)xi * s[j]
+ tmp[i + j] + carry;
tmp[i + j] = (uint64_t)t;
carry = t >> 64;
}
size_t k = i + n;
while (carry) {
unsigned __int128 t = (unsigned __int128)tmp[k] + carry;
tmp[k] = (uint64_t)t;
carry = t >> 64;
k++;
}
}
/* ── Phase 2: double the upper-triangle sum (1-bit left-shift) ── */
uint64_t carry_bit = 0;
for (size_t k = 0; k < n2 + 2; k++) {
uint64_t next = tmp[k] >> 63;
tmp[k] = (tmp[k] << 1) | carry_bit;
carry_bit = next;
}
/* ── Phase 3: add diagonal s[i]² at position 2i ── */
for (size_t i = 0; i < n; i++) {
unsigned __int128 diag = (unsigned __int128)s[i] * s[i];
unsigned __int128 carry = (unsigned __int128)tmp[2*i] + (uint64_t)diag;
tmp[2*i] = (uint64_t)carry;
carry = (carry >> 64) + (diag >> 64);
for (size_t k = 2*i + 1; carry; k++) {
carry += tmp[k];
tmp[k] = (uint64_t)carry;
carry >>= 64;
}
}
fold_mod_mp_a(tmp, n2 + 2, p, s, n);
}
/* ap_sub2_mod_mp: in-place s[] = s[] - 2 mod 2^p-1.
* Identical to mpi_sub2_mod_mp in ll_mpi.cu. */
static void ap_sub2_mod_mp(uint64_t *s, size_t n, uint64_t p) {
uint64_t pb = p % 64;
/* check if s < 2 */
int small = 1;
for (size_t k = n; k-- > 1; )
if (s[k]) { small = 0; break; }
if (small && s[0] >= 2) small = 0;
if (!small) {
uint64_t borrow = 2;
for (size_t k = 0; k < n && borrow; k++) {
if (s[k] >= borrow) { s[k] -= borrow; borrow = 0; }
else { s[k] -= borrow; borrow = 1; }
}
} else {
/* s is 0 or 1: result = M_p + s - 2 */
uint64_t val = s[0];
for (size_t k = 0; k < n; k++) s[k] = ~0ULL;
if (pb > 0) s[n-1] = (1ULL << pb) - 1ULL;
uint64_t sub = 2 - val;
uint64_t borrow = sub;
for (size_t k = 0; k < n && borrow; k++) {
if (s[k] >= borrow) { s[k] -= borrow; borrow = 0; }
else { s[k] -= borrow; borrow = 1; }
}
}
}
static int is_zero_a(const uint64_t *words, size_t n) {
for (size_t k = 0; k < n; k++)
if (words[k]) return 0;
return 1;
}
/* ════════════════════════════════════════════════════════════════════════════
* ll_analog: main entry point
*
* Runs exact Lucas-Lehmer with:
* mantissa[0..n-1] — Slot4096.mantissa_words equivalent (v30b layout)
* AnaOsc8D osc — 8D Kuramoto oscillator running in parallel
*
* Every ANA_SHA_INTERVAL iters: residue hash → oscillator phase perturbation
* (cooperative / conditional memory: arithmetic couples into analog state)
*
* Progress is logged only on natural phase transitions (wu-wei pacing).
* Final report includes oscillator lock status alongside residue result.
* ════════════════════════════════════════════════════════════════════════════ */
int ll_analog(uint64_t p, int verbose) {
if (p == 2) return 1; /* M_2 = 3, known prime; LL loop undefined for p<3 */
size_t n = (size_t)((p + 63) / 64); /* Slot4096 mantissa word count */
size_t n2 = 2 * n;
/* ── Allocate: residue "mantissa_words" + squaring scratch (v30b style) ── */
uint64_t *mantissa = (uint64_t *)calloc(n, sizeof(uint64_t));
uint64_t *tmp = (uint64_t *)calloc(n2 + 2, sizeof(uint64_t));
if (!mantissa || !tmp) {
fprintf(stderr, "[ll_analog] allocation failed (n=%zu)\n", n);
free(mantissa); free(tmp);
return -1;
}
mantissa[0] = 4; /* LL initial seed: s₀ = 4 */
/* ── Initialise 8D Kuramoto oscillator ── */
AnaOsc8D osc;
ana_init(&osc, p); /* glyph chain reaction seeds theta[i] and omega[i] */
if (verbose) {
printf(" [analog] p=%llu n_words=%zu osc=8D-Kuramoto\n",
(unsigned long long)p, n);
{
double Lv = log((double)p * ANA_LN2 / ANA_LN_PHI) / ANA_LN_PHI - 0.5 / ANA_PHI;
double fv = Lv - floor(Lv);
double Ov = 0.5 * (1.0 + sin(ANA_PI * fv * ANA_PHI));
printf(" [analog] seed: Lambda_phi(p)=%.6f {L}=%.6f Omega=%.6f\n",
Lv, fv, Ov);
printf(" [analog] theta0: pi*L + 2pi*(glyph+{L}+i*D_n_r) [e^(i*pi*L) Euler base]\n");
printf(" [analog] omega: Omega*phi^(1+i*D_n_r)*dt [phi-lattice resonance envelope]\n");
printf(" [analog] field: M(U)=|sum(e^itheta)| Lambda^U=log(M)/lnphi-1/2phi S=|Omega*e^(i*pi*L)+1|\n");
}
printf(" [analog] CV: Kuramoto 1-R in [0,1] (circular; 0=locked, 1=spread)\n");
printf(" [analog] multiply: phase-doubling (theta->2theta) + Kuramoto coupling\n");
printf(" [analog] sync: harmonic attraction alpha=%.1fx%d (atan2, first Fourier modes)\n",
ANA_HARM_ALPHA, ANA_HARM_PASSES);
#ifdef LL_QUANTUM_ENABLED
printf(" [analog] quantum: 8-qubit XY Hamiltonian (cuStateVec sm_75) blend=%.2f..%.2f\n",
qosc_blend_alpha(0), qosc_blend_alpha(3));
#else
printf(" [analog] quantum: disabled (build without -DLL_QUANTUM_ENABLED)\n");
#endif
}
clock_t t0 = clock();
uint64_t iters = p - 2;
APhase last_phase = APHASE_PLUCK;
int logged = 0;
/* ══ Main LL loop ══════════════════════════════════════════════════════ */
for (uint64_t iter = 0; iter < iters; iter++) {
/* ── Exact arithmetic: s = s² - 2 mod 2^p-1 ── */
ap_sqr_mersenne(mantissa, n, p, tmp);
ap_sub2_mod_mp(mantissa, n, p);
/* ── Analog squaring: θ → 2θ (s² = phase doubling on unit circle) ── */
ana_phase_double(&osc);
/* ── Kuramoto coupling: synchronisation correction (RK4) ── */
ana_rk4_step(&osc);
osc.phase_var = ana_phase_var(&osc);
ana_update_phase(&osc);
#ifdef LL_QUANTUM_ENABLED
/* ── Quantum XY step + blend ────────────────────────────────────────
* Evolve the 8-qubit state vector by one Trotterised XY step, then
* blend its expectation values with the classical RK4 result.
*
* Blend weight α rises with phase coherence (Pluck→Lock): the quantum
* layer contributes more as the oscillator settles, where its O(K²dt²)
* Trotter error is smallest relative to the signal being measured.
*
* After blending, renormalise each (re,im) to the unit circle so that
* ana_phase_var and the phase-doubling double-angle formula remain exact.
*/
if (osc.q_osc) {
qosc_step(osc.q_osc, osc.omega, osc.k_coupling, ANA_DT);
double q_re[ANA_DIMS], q_im[ANA_DIMS];
qosc_readback(osc.q_osc, q_re, q_im);
double alpha = qosc_blend_alpha((int)osc.aphase);
double one_minus = 1.0 - alpha;
for (int qi = 0; qi < ANA_DIMS; qi++) {
double r = one_minus * osc.re[qi] + alpha * q_re[qi];
double im = one_minus * osc.im[qi] + alpha * q_im[qi];
double inv_mag = 1.0 / sqrt(r*r + im*im);
osc.re[qi] = r * inv_mag;
osc.im[qi] = im * inv_mag;
osc.theta[qi] = atan2(osc.im[qi], osc.re[qi]);
if (osc.theta[qi] < 0.0) osc.theta[qi] += 2.0 * ANA_PI;
}
}
#endif
/* ── Harmonic sync: attract oscillators toward residue-derived phases ── */
if ((iter & (ANA_SHA_INTERVAL - 1)) == 0)
ana_harmonic_sync(&osc, mantissa, n);
/* ── Progress: only log on natural phase transitions (wu-wei pacing) ── */
if (verbose) {
int is_transition = (osc.aphase != last_phase);
int is_milestone = (iter == 0 || iter == iters - 1
|| (iters > 20 && iter % (iters / 10) == 0));
if (is_transition || is_milestone) {
last_phase = osc.aphase;
double pct = 100.0 * (double)(iter + 1) / (double)iters;
double elapsed = (double)(clock() - t0) / CLOCKS_PER_SEC;
printf(" [analog] iter=%-8llu %5.1f%% phase=%-8s cv=%.4f t=%.1fs%s\n",
(unsigned long long)iter, pct,
APHASE_NAMES[osc.aphase], osc.phase_var, elapsed,
is_transition ? " [phase transition]" : "");
logged++;
}
}
}
/* ── Final analog confirmation: harmonic sync on final residue ────────────
* Ensures cv_hist's last entry reflects residue=0 (prime) or ≠0 (composite)
* regardless of where the last periodic sync fell. */
ana_harmonic_sync(&osc, mantissa, n);
#ifdef LL_QUANTUM_ENABLED
if (osc.q_osc) qosc_destroy(osc.q_osc);
#endif
/* ── Final result ── */
int result = is_zero_a(mantissa, n);
if (verbose) {
double elapsed = (double)(clock() - t0) / CLOCKS_PER_SEC;
printf(" [analog] done: %.2fs phase=%s cv=%.4f locked=%s residue=%s\n",
elapsed,
APHASE_NAMES[osc.aphase],
osc.phase_var,
ana_is_locked(&osc) ? "yes" : "no",
result ? "0 (PRIME)" : "non-zero (COMPOSITE)");
/* Field resonance readout S(U) */
printf(" [analog] S(U)=%.6f Lambda^U=%.6f (prime: S~%.4f, composite: S!=)\n",
osc.s_u, osc.lambda_u,
0.5 * (1.0 + sin(ANA_PI * 0.0 * ANA_PHI)) * cos(0.0) + 1.0);
/* φ-tower fixed point */
printf(" [analog] phi-tower: U*=%.6f delta=%.2e (phase-spread convergence signal)\n",
osc.fp_u, osc.fp_res);
/* Double confirmation: both analog and exact agree */
if (ana_is_locked(&osc) && result)
printf(" [analog] ** osc LOCKED + residue=0: strong prime resonance **\n");
if (!ana_is_locked(&osc) && result)
printf(" [analog] note: residue=0 (prime) but osc not locked\n");
printf(" [analog] n_words=%zu iters=%llu rk4_steps=%d log_events=%d\n",
n, (unsigned long long)iters, osc.steps, logged);
}
free(mantissa);
free(tmp);
return result;
}
ll_analog.h
/* ll_analog.h — analog LL path: v30b Slot4096 APA + 8D Kuramoto oscillator
*
* Exact arithmetic:
* Slot4096.mantissa_words — p-bit LL residue, arbitrary precision
* ap_sqr_mersenne — schoolbook O(n^2) squaring + Mersenne fold
* ap_sub2_mod_mp — subtract 2 mod 2^p-1
*
* Analog oscillator (after analog_engine.h / AnalogContainer1):
* AnaOsc8D — 8D Kuramoto RK4, φ-seeded natural frequencies
* Adaptive phase — Pluck→Sustain→FineTune→Lock (K/γ wu-wei ratios)
* Cooperative memory — XOR-fold residue hash → phase perturbation every
* ANA_SHA_INTERVAL iters (conditional memory loop)
*
* The oscillator does NOT shortcut correctness. Every p-2 iterations run exact.
* Phase lock is a readout, not a gate — wu-wei.
*
* Build: ll_analog.c compiles as plain C (-O2); link with ll_mpi.cu via clang.
*
* Licensed per https://zchg.org/t/legal-notice-copyright-applicable-ip-and-licensing-read-me/440
*/
#pragma once
#include <stdint.h>
#ifdef __cplusplus
extern "C" {
#endif
/* Main entry point: Lucas-Lehmer test via v30b APA + 8D Kuramoto.
* Returns 1 if M_p is prime, 0 if composite, -1 on allocation failure. */
int ll_analog(uint64_t p, int verbose);
#ifdef __cplusplus
}
#endif
ll_analog .ASM
; =============================================================================
; hdgl_analog_v30.asm
; x86-64 NASM — Linux SysV ABI
; Complete translation of hdgl_analog_v30.c (HDGL Analog Mainnet V3.0)
;
; Build:
; nasm -f elf64 hdgl_analog_v30.asm -o hdgl_analog_v30.o
; gcc -no-pie hdgl_analog_v30.o -lm -o hdgl_analog_v30
; =============================================================================
bits 64
default rel
; --- externals ----------------------------------------------------------------
extern malloc, calloc, free, memset, memcpy
extern rand, srand, nanosleep
extern clock_gettime
extern printf
extern time
extern cos, sin, pow, sqrt, fmod, exp, fabs, frexp
; --- POSIX -------------------------------------------------------------------
%define CLOCK_MONOTONIC 1
; --- domain constants --------------------------------------------------------
%define NUM_DN 8
%define CHUNK_SIZE 1048576
%define CHECKPOINT_INTERVAL 100
%define SNAPSHOT_MAX 10
%define CONSENSUS_N 100
; --- APA flags ---------------------------------------------------------------
%define FLAG_SIGN_NEG 1
%define FLAG_IS_NAN 2
%define FLAG_GOI 4
%define FLAG_GUZ 8
%define FLAG_CONSENSUS 16
; --- timespec ----------------------------------------------------------------
%define TS_SIZE 16
%define TS_SEC 0
%define TS_NSEC 8
; --- NumericLattice ----------------------------------------------------------
%define NL_upper 0 ; 7 x f64 = 56
%define NL_adims 56 ; 13 x f64 = 104
%define NL_void 160 ; 1 x f64 = 8
%define NL_lower 168 ; 8 x f64 = 64
%define NL_sib 232 ; 8 x f64 = 64
%define NL_inf 296 ; 4 x f64 = 32
%define NL_choke 328 ; 4 x f64 = 32
%define NL_seeds_ptr 360
%define NL_num_seeds 368
%define NL_SIZE 376
; --- MPI ---------------------------------------------------------------------
%define MPI_words 0
%define MPI_nwords 8
%define MPI_sign 16
%define MPI_SIZE 24
; --- Slot4096 ----------------------------------------------------------------
%define SL_mwords 0
%define SL_nwm 8 ; MPI(24)
%define SL_empi 32 ; MPI(24)
%define SL_ebase 56 ; uint16
%define SL_flags 64 ; uint32
%define SL_soi 72 ; MPI(24)
%define SL_nw 96 ; size_t
%define SL_exp 104 ; int64
%define SL_basef 112 ; float
%define SL_bmant 120 ; int
%define SL_bexp 124 ; int
%define SL_phase 128 ; double
%define SL_pvel 136 ; double
%define SL_freq 144 ; double
%define SL_ampim 152 ; double
%define SL_dimn 160 ; int
%define SL_rval 168 ; double
%define SL_dnamp 176 ; double
%define SL_wmode 184 ; double
%define SL_SIZE 192
; --- AnalogLink --------------------------------------------------------------
%define AL_chg 0
%define AL_chgim 8
%define AL_tens 16
%define AL_pot 24
%define AL_coup 32
%define AL_dncoup 40
%define AL_SIZE 48
; --- ComplexState ------------------------------------------------------------
%define CS_Ar 0
%define CS_Ai 8
%define CS_ph 16
%define CS_pv 24
%define CS_dn 32
%define CS_SIZE 40
; --- HDGLChunk ---------------------------------------------------------------
%define CH_slots 0
%define CH_alloc 8
%define CH_SIZE 16
; --- HDGLLattice -------------------------------------------------------------
%define LT_chunks 0
%define LT_nchunks 8
%define LT_ninst 12
%define LT_spi 16
%define LT_omega 24
%define LT_time 32
%define LT_csteps 40
%define LT_pvar 48
%define LT_ckns 56
%define LT_nl 64
%define LT_SIZE 72
; --- CheckpointMeta ----------------------------------------------------------
%define CM_evo 0
%define CM_tsns 8
%define CM_pvar 16
%define CM_omega 24
%define CM_wt 32
%define CM_SIZE 40
; --- CheckpointManager -------------------------------------------------------
%define CK_snaps 0
%define CK_cnt 8
%define CK_cap 12
%define CK_SIZE 16
; =============================================================================
section .data
align 8
fib_tab: dq 1, 1, 2, 3, 5, 8, 13, 21
prime_tab: dq 2, 3, 5, 7, 11, 13, 17, 19
nl_upper:
dq 0x4065627F82888889 ; 170.6180339887
dq 0x4062F1F530B74E5E ; 150.9442719100
dq 0x40293EB851EB851F ; 12.6180339887
dq 0x40214CF0ABF61E8E ; 8.8541019662
dq 0x4010F1A9FBE76C8C ; 4.2360679775
dq 0x400CE1DC6C0D72AE ; 3.6180339887
dq 0x3FF9E3779B97F4A8 ; 1.6180339887
nl_adims:
dq 0x4020A3EF9DB22D0E ; 8.3141592654
dq 0x401F6BE476B9F800 ; 7.8541019662
dq 0x4019E3779B97F4A8 ; 6.4721359549
dq 0x40167AE147AE147B ; 5.6180339887
dq 0x401370A3D70A3D71 ; 4.8541019662
dq 0x400CE1DC6C0D72AE ; 3.6180339887
dq 0x4004F1A9FBE76C8C ; 2.6180339887
dq 0x3FF9E3779B97F4A8 ; 1.6180339887
dq 0x3FF0000000000000 ; 1.0
dq 0x401F6BE476B9F800 ; 7.8541019662
dq 0x402622F04750B2F1 ; 11.0901699437
dq 0x4031F1B04060EC04 ; 17.9442719100
dq 0x403D0B026F9308D4 ; 29.0344465435
nl_lower:
dq 0x3E112E0BE826D695 ; ~1e-10
dq 0x3FA1A5A2C1B2B100 ; 0.0344465435
dq 0x3FAC90AEC4D79456 ; 0.0557280900
dq 0x3FB717C25A3E3B3B ; 0.0901699437
dq 0x3FC2A8C30C1B85C5 ; 0.1458980338
dq 0x3FCE2B5C2D2D2D2D ; 0.2360679775
dq 0x3FD86B85161B5A3D ; 0.3819660113
dq 0x3FE3C6EF372FE950 ; 0.6180339887
nl_sib:
dq 0x3FB717C25A3E3B3B ; 0.0901699437
dq 0x3FC2A8C30C1B85C5 ; 0.1458980338
dq 0x3FCE2B5C2D2D2D2D ; 0.2360679775
dq 0x3FD3C6EF372FE950 ; 0.3090169944
dq 0x3FD86B85161B5A3D ; 0.3819660113
dq 0x3FDE2B5C2D2D2D2D ; 0.4721359549
dq 0x3FE4F8B588E368F1 ; 0.6545084972
dq 0x3FEBE2B5C2D2D2D3 ; 0.8729833462
nl_seeds:
dq 0x3FE3C6EF372FE950, 0x3FF9E3779B97F4A8, 0x4004F1A9FBE76C8C, 0x400CE1DC6C0D72AE
dq 0x401370A3D70A3D71, 0x40167AE147AE147B, 0x4019E3779B97F4A8, 0x401F6BE476B9F800
dq 0x4020A3EF9DB22D0E, 0x3FB717C25A3E3B3B, 0x3FC2A8C30C1B85C5, 0x3FCE2B5C2D2D2D2D
dq 0x3FD3C6EF372FE950, 0x3FD86B85161B5A3D, 0x3FDE2B5C2D2D2D2D, 0x3FE4F8B588E368F1
dq 0x3FEBE2B5C2D2D2D3, 0x3FF0000000000000, 0x3FF3C6EF372FE950, 0x3FF9E3779B97F4A8
dq 0x4001E3779B97F4A8, 0x4004F1A9FBE76C8C, 0x400921FB54442D18, 0x400CE1DC6C0D72AE
dq 0x4010E2B5C2D2D2D3, 0x401370A3D70A3D71, 0x40167AE147AE147B, 0x4019E3779B97F4A8
dq 0x401CE2B5C2D2D2D3, 0x401F6BE476B9F800, 0x402137C6EF372FE9, 0x402483BC6EF37300
dq 0x40279D3A0AD7A3D7, 0x4025413C6EF37300, 0x402622F04750B2F1, 0x4027E1DC6C0D72AE
dq 0x40293EB851EB851F, 0x402B413C6EF37300, 0x402C839CAF3F72FE, 0x402DC77FFA3D70A4
dq 0x402FD3A0AD7A3D71, 0x40309D3A0AD7A3D7, 0x403141FAAAAAAAB0, 0x4031F1B04060EC04
dq 0x4032A3D70A3D70A4, 0x40334B3B13B13B14, 0x4033F5C28F5C28F6, 0x40349E3779B97F4B
dq 0x403533B13B13B13B, 0x4035DFDF3B13B13B, 0x403686CAC083126F, 0x40378CAC083126F7
dq 0x403836DF3B13B13B, 0x4038DFE147AE147B, 0x40398A3D70A3D70A, 0x403A3279B97F4A8B
dq 0x403ADBE76C8B4396, 0x403B83126E978D4F, 0x403C2B851EB851EC, 0x403CD044EF82BF7E
dq 0x403D7851EB851EB8, 0x403E1CAC083126F7, 0x403EC2B5C2D2D2D3, 0x403F6B851EB851EC
align 8
fp_PHI: dq 0x3FF9E3779B97F4A8 ; 1.6180339887498948
fp_TWO_PI: dq 0x401921FB54442D18 ; 2*pi
fp_GAMMA: dq 0x3F947AE147AE147B ; 0.02
fp_LAMBDA: dq 0x3FA999999999999A ; 0.05
fp_SAT: dq 0x412E848000000000 ; 1e6
fp_NSIG: dq 0x3F847AE147AE147B ; 0.01
fp_CEPS: dq 0x3EF0624DD2F1A9FC ; 1e-6
fp_ATH: dq 0x3FE999999999999A ; 0.8
fp_KC: dq 0x3FF0000000000000 ; 1.0
fp_SDEC: dq 0x3FE6666666666666 ; 0.95
fp_0p98: dq 0x3FEF5C28F5C28F5C ; 0.98
fp_0p95: dq 0x3FE6666666666666 ; 0.95
fp_0: dq 0x0000000000000000
fp_001: dq 0x3F847AE147AE147B ; 0.01
fp_01: dq 0x3FB999999999999A ; 0.1
fp_03: dq 0x3FD3333333333333 ; 0.3
fp_half: dq 0x3FE0000000000000 ; 0.5
fp_1: dq 0x3FF0000000000000 ; 1.0
fp_2: dq 0x4000000000000000 ; 2.0
fp_6: dq 0x4018000000000000 ; 6.0
fp_8: dq 0x4020000000000000 ; 8.0
fp_1000: dq 0x408F400000000000 ; 1000.0
fp_1e10n: dq 0x3EE4F8B588E368F1 ; 1e-10
fp_1e6n: dq 0x3EF0624DD2F1A9FC ; 1e-6 (min dt)
fp_dt0: dq 0x3F00000000000000 ; 1/32768
fp_U64: dq 0x43F0000000000000 ; (double)UINT64_MAX
fp_INF: dq 0x7FF0000000000000
fp_DMAX: dq 0x7FEFFFFFFFFFFFFF
fp_NEG1: dq 0xBFF0000000000000 ; -1.0
sz_rtc: db "[RTC] Using software fallback (CLOCK_MONOTONIC)",10,0
sz_banner: db "=== HDGL Analog Mainnet V3.0: Dn(r) Engine Ready ===",10,10,0
sz_bi: db "[Bootloader] Initializing HDGL Analog Mainnet V3.0 with Dn(r) Engine...",10,0
sz_berr: db "[Bootloader] ERROR: Lattice allocation failed.",10,0
sz_binfo: db "[Bootloader] %d instances, %d total slots",10,0
sz_bseed: db "[Bootloader] Numeric Lattice loaded with %zu Base(oo) seeds",10,0
sz_bdone: db "[Bootloader] Lattice seeded with %d RK4 steps",10,0
sz_bstat: db "[Bootloader] Omega: %.6f, Time: %.6f, PhaseVar: %.6f",10,0
sz_nlhdr: db 10,"Numeric Lattice Summary:",10,0
sz_nluf: db " Upper Field[0]: %.10f",10,0
sz_nlad: db " Analog D8: %.10f",10,0
sz_nlvd: db " The Void: %.10f",10,0
sz_nllf: db " Lower Field[7]: %.10f",10,0
sz_nlbs: db " Base(oo) Seeds: %zu total",10,0
sz_slhdr: db 10,"First 8 slots (post-evolution with Dn(r)):",10,0
sz_slfmt: db " D%d: |A|=%.6e phi=%.3f Dn=%.3f wave=%.1f r=%.3f",10,0
sz_cons: db "[CONSENSUS] Domain locked at t=%.4f (var=%.6f)!",10,0
sz_ckpt: db "[Checkpoint] Saved evo %d (total: %d, var=%.6f)",10,0
sz_aerr: db "Error: Failed to allocate mantissa.",10,0
sz_cerr: db "Error: Copy allocation failed.",10,0
sz_werr: db "Error: Unaligned word counts.",10,0
sz_fatal: db "Fatal: Could not initialize lattice.",10,0
sz_final: db 10,"=== HDGL V3.0 OPERATIONAL ===",10,0
; =============================================================================
section .bss
align 8
ts_clk: resb TS_SIZE
ts_slp: resb TS_SIZE
; =============================================================================
section .text
global main
; -----------------------------------------------------------------------------
; uint64_t det_rand(uint64_t seed) rdi=seed -> rax xorshift64
; -----------------------------------------------------------------------------
det_rand:
mov rax, rdi
mov rcx, rax
shl rcx, 13
xor rax, rcx
mov rcx, rax
shr rcx, 7
xor rax, rcx
mov rcx, rax
shl rcx, 17
xor rax, rcx
ret
; -----------------------------------------------------------------------------
; double get_normalized_rand() -> xmm0 in [0,1)
; -----------------------------------------------------------------------------
get_normalized_rand:
push rbp
mov rbp, rsp
and rsp, -16
call rand
leave
cvtsi2sd xmm0, eax
mov eax, 0x7FFFFFFF
cvtsi2sd xmm1, eax
divsd xmm0, xmm1
ret
; -----------------------------------------------------------------------------
; int64_t get_rtc_ns() -> rax
; -----------------------------------------------------------------------------
get_rtc_ns:
push rbp
mov rbp, rsp
and rsp, -16
mov edi, CLOCK_MONOTONIC
lea rsi, [ts_clk]
call clock_gettime
leave
mov rax, [ts_clk + TS_SEC]
imul rax, rax, 1000000000
add rax, [ts_clk + TS_NSEC]
ret
; -----------------------------------------------------------------------------
; void rtc_sleep_until(int64_t target) rdi=target
; -----------------------------------------------------------------------------
rtc_sleep_until:
push rbp
mov rbp, rsp
push rbx
mov rbx, rdi
call get_rtc_ns
cmp rbx, rax
jle .done
sub rbx, rax
mov rax, rbx
xor edx, edx
mov ecx, 1000000000
div ecx
mov [ts_slp + TS_SEC], rax
mov [ts_slp + TS_NSEC], rdx
and rsp, -16
lea rdi, [ts_slp]
xor rsi, rsi
call nanosleep
.done:
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; void init_numeric_lattice(NumericLattice *nl) rdi=nl
; -----------------------------------------------------------------------------
init_numeric_lattice:
push rbp
mov rbp, rsp
push rbx
mov rbx, rdi
lea rsi, [nl_upper]
lea rdi, [rbx + NL_upper]
mov ecx, 7
rep movsq
lea rsi, [nl_adims]
lea rdi, [rbx + NL_adims]
mov ecx, 13
rep movsq
mov qword [rbx + NL_void], 0
lea rsi, [nl_lower]
lea rdi, [rbx + NL_lower]
mov ecx, 8
rep movsq
lea rsi, [nl_sib]
lea rdi, [rbx + NL_sib]
mov ecx, 8
rep movsq
movsd xmm0, [fp_INF]
movsd [rbx + NL_inf + 0], xmm0
movsd [rbx + NL_inf + 8], xmm0
movsd [rbx + NL_inf + 16], xmm0
movsd [rbx + NL_inf + 24], xmm0
movsd xmm0, [fp_DMAX]
movsd [rbx + NL_choke + 0], xmm0
movsd [rbx + NL_choke + 8], xmm0
movsd [rbx + NL_choke + 16], xmm0
movsd [rbx + NL_choke + 24], xmm0
mov qword [rbx + NL_num_seeds], 64
push rbx
and rsp, -16
mov edi, 512
call malloc
pop rbx
mov [rbx + NL_seeds_ptr], rax
mov rdi, rax
lea rsi, [nl_seeds]
mov ecx, 64
rep movsq
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; void free_numeric_lattice(NumericLattice *nl) rdi=nl
; -----------------------------------------------------------------------------
free_numeric_lattice:
push rbp
mov rbp, rsp
push rbx
mov rbx, rdi
mov rax, [rbx + NL_seeds_ptr]
test rax, rax
jz .done
push rbx
and rsp, -16
mov rdi, rax
call free
pop rbx
mov qword [rbx + NL_seeds_ptr], 0
.done:
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; double compute_Dn_r(int n, double r, double omega)
; rdi=n xmm0=r xmm1=omega -> xmm0
; Dn(r) = sqrt(phi * F_n * 2^n * P_n * omega) * |r|^((n+1)/8)
; -----------------------------------------------------------------------------
compute_Dn_r:
push rbp
mov rbp, rsp
sub rsp, 48
push rbx
cmp edi, 1
jl .zero
cmp edi, NUM_DN
jg .zero
mov ebx, edi
movsd [rbp-8], xmm0 ; r
movsd [rbp-16], xmm1 ; omega
; F_n = fib_tab[n-1]
lea rax, [fib_tab]
cvtsi2sd xmm2, qword [rax + rbx*8 - 8]
; 2^n
movsd xmm0, [fp_2]
cvtsi2sd xmm1, rbx
and rsp, -16
call pow
movsd [rbp-24], xmm0 ; 2^n
; P_n = prime_tab[n-1]
lea rax, [prime_tab]
cvtsi2sd xmm3, qword [rax + rbx*8 - 8]
; base = sqrt(phi * F_n * 2^n * P_n * omega)
movsd xmm0, [fp_PHI]
mulsd xmm0, xmm2
mulsd xmm0, [rbp-24]
mulsd xmm0, xmm3
mulsd xmm0, [rbp-16]
call sqrt
movsd [rbp-32], xmm0
; k = (n+1)/8.0
mov eax, ebx
inc eax
cvtsi2sd xmm0, eax
movsd xmm1, [fp_8]
divsd xmm0, xmm1
movsd [rbp-40], xmm0
; |r|^k
movsd xmm0, [rbp-8]
call fabs
movsd xmm1, [rbp-40]
call pow
mulsd xmm0, [rbp-32]
jmp .done
.zero:
xorpd xmm0, xmm0
.done:
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; void mpi_init(MPI *m, size_t nw) rdi=m rsi=nw
; -----------------------------------------------------------------------------
mpi_init:
push rbp
mov rbp, rsp
push rbx
push r12
mov rbx, rdi
mov r12, rsi
and rsp, -16
mov rdi, rsi
mov rsi, 8
call calloc
mov [rbx + MPI_words], rax
mov [rbx + MPI_nwords], r12
mov byte [rbx + MPI_sign], 0
pop r12
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; void mpi_free(MPI *m) rdi=m
; -----------------------------------------------------------------------------
mpi_free:
push rbp
mov rbp, rsp
push rbx
mov rbx, rdi
mov rax, [rbx + MPI_words]
test rax, rax
jz .done
and rsp, -16
mov rdi, rax
call free
mov qword [rbx + MPI_words], 0
mov qword [rbx + MPI_nwords], 0
.done:
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; void mpi_copy(MPI *dst, const MPI *src) rdi=dst rsi=src
; -----------------------------------------------------------------------------
mpi_copy:
push rbp
mov rbp, rsp
push rbx
push r12
mov rbx, rdi
mov r12, rsi
call mpi_free ; free existing dst words
mov rcx, [r12 + MPI_nwords]
mov [rbx + MPI_nwords], rcx
mov al, [r12 + MPI_sign]
mov [rbx + MPI_sign], al
and rsp, -16
mov rdi, rcx
shl rdi, 3
call malloc
mov [rbx + MPI_words], rax
mov rdi, rax
mov rsi, [r12 + MPI_words]
mov rdx, [rbx + MPI_nwords]
shl rdx, 3
call memcpy
pop r12
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; void mpi_set(MPI *m, uint64_t val, uint8_t sign) rdi=m rsi=val dl=sign
; -----------------------------------------------------------------------------
mpi_set:
mov rax, [rdi + MPI_words]
test rax, rax
jz .done
mov [rax], rsi
mov [rdi + MPI_sign], dl
.done:
ret
; -----------------------------------------------------------------------------
; void slot_init(Slot4096 *out, int bm, int be, int dn, double rv, double om)
; rdi=out rsi=bm rdx=be rcx=dn xmm0=rv xmm1=om
; -----------------------------------------------------------------------------
slot_init:
push rbp
mov rbp, rsp
sub rsp, 64
push rbx
push r12
push r13
push r14
push r15
mov rbx, rdi
mov r12d, esi ; bm
mov r13d, edx ; be
mov r14d, ecx ; dn
movsd [rbp-8], xmm0 ; rv
movsd [rbp-16], xmm1 ; om
and rsp, -16
mov rdi, rbx
xor rsi, rsi
mov rdx, SL_SIZE
call memset
mov [rbx + SL_bmant], r12d
mov [rbx + SL_bexp], r13d
; num_words = (bm+63)/64
mov eax, r12d
add eax, 63
shr eax, 6
movsx rax, eax
mov [rbx + SL_nw], rax
; calloc(nw, 8)
mov rdi, rax
mov rsi, 8
call calloc
test rax, rax
jnz .mant_ok
lea rdi, [sz_aerr]
xor eax, eax
call printf
jmp .done
.mant_ok:
mov [rbx + SL_mwords], rax
lea rdi, [rbx + SL_empi]
mov rsi, 1
call mpi_init
lea rdi, [rbx + SL_nwm]
mov rsi, 1
call mpi_init
lea rdi, [rbx + SL_soi]
mov rsi, 1
call mpi_init
mov [rbx + SL_dimn], r14d
movsd xmm0, [rbp-8]
movsd [rbx + SL_rval], xmm0
; Dn_amplitude
mov edi, r14d
movsd xmm0, [rbp-8]
movsd xmm1, [rbp-16]
call compute_Dn_r
movsd [rbx + SL_dnamp], xmm0
movsd [rbp-24], xmm0
; wave_mode from dimn%3
mov eax, r14d
cdq
mov ecx, 3
idiv ecx
cmp edx, 1
je .wm1
cmp edx, 2
je .wm0
movsd xmm0, [fp_NEG1]
jmp .wm_set
.wm1:
movsd xmm0, [fp_1]
jmp .wm_set
.wm0:
xorpd xmm0, xmm0
.wm_set:
movsd [rbx + SL_wmode], xmm0
; mantissa_words[0] = (uint64_t)(|Dn|*U64/1000) | MSB
cmp qword [rbx + SL_nw], 0
je .skip_mant
movsd xmm0, [rbp-24]
call fabs
mulsd xmm0, [fp_U64]
divsd xmm0, [fp_1000]
maxsd xmm0, [fp_0]
cvttsd2si rax, xmm0
mov rcx, [rbx + SL_mwords]
mov [rcx], rax
bts qword [rcx], 63
.skip_mant:
; exponent = rand()%exp_range - exp_bias
mov ecx, r13d
mov r15d, 1
shl r15d, cl ; exp_range = 1<<be
call rand
xor edx, edx
div r15d ; edx = rand%exp_range
mov ecx, r13d
dec ecx
mov eax, 1
shl eax, cl ; exp_bias
sub edx, eax
movsx rax, edx
mov [rbx + SL_exp], rax
mov word [rbx + SL_ebase], 4096
; base = PHI + 0.01*rand
call get_normalized_rand
mulsd xmm0, [fp_001]
addsd xmm0, [fp_PHI]
cvtsd2ss xmm0, xmm0
movss [rbx + SL_basef], xmm0
; phase = 2pi*rand
call get_normalized_rand
mulsd xmm0, [fp_TWO_PI]
movsd [rbx + SL_phase], xmm0
; freq = 1 + 0.5*rand
call get_normalized_rand
mulsd xmm0, [fp_half]
addsd xmm0, [fp_1]
movsd [rbx + SL_freq], xmm0
; amp_im = 0.1*rand
call get_normalized_rand
mulsd xmm0, [fp_01]
movsd [rbx + SL_ampim], xmm0
; sync empi
mov rax, [rbx + SL_exp]
mov rsi, rax
xor dl, dl
test rax, rax
jns .ep1
neg rsi
mov dl, 1
.ep1:
lea rdi, [rbx + SL_empi]
call mpi_set
; sync nwm
lea rdi, [rbx + SL_nwm]
mov rsi, [rbx + SL_nw]
xor dl, dl
call mpi_set
.done:
pop r15
pop r14
pop r13
pop r12
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; void ap_free(Slot4096 *s) rdi=s
; -----------------------------------------------------------------------------
ap_free:
push rbp
mov rbp, rsp
push rbx
mov rbx, rdi
test rbx, rbx
jz .done
mov rax, [rbx + SL_mwords]
test rax, rax
jz .skip_m
and rsp, -16
mov rdi, rax
call free
mov qword [rbx + SL_mwords], 0
.skip_m:
lea rdi, [rbx + SL_empi]
call mpi_free
lea rdi, [rbx + SL_nwm]
call mpi_free
lea rdi, [rbx + SL_soi]
call mpi_free
mov qword [rbx + SL_nw], 0
.done:
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; void ap_copy(Slot4096 *dst, const Slot4096 *src) rdi=dst rsi=src
; -----------------------------------------------------------------------------
ap_copy:
push rbp
mov rbp, rsp
push rbx
push r12
mov rbx, rdi
mov r12, rsi
call ap_free
and rsp, -16
mov rdi, rbx
mov rsi, r12
mov rdx, SL_SIZE
call memcpy
; deep-copy mantissa
mov rcx, [rbx + SL_nw]
test rcx, rcx
jz .done
mov rdi, rcx
shl rdi, 3
call malloc
test rax, rax
jnz .mok
lea rdi, [sz_cerr]
xor eax, eax
call printf
mov qword [rbx + SL_nw], 0
jmp .done
.mok:
mov [rbx + SL_mwords], rax
mov rdi, rax
mov rsi, [r12 + SL_mwords]
mov rdx, [rbx + SL_nw]
shl rdx, 3
call memcpy
lea rdi, [rbx + SL_empi]
lea rsi, [r12 + SL_empi]
call mpi_copy
lea rdi, [rbx + SL_nwm]
lea rsi, [r12 + SL_nwm]
call mpi_copy
lea rdi, [rbx + SL_soi]
lea rsi, [r12 + SL_soi]
call mpi_copy
.done:
pop r12
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; double ap_to_double(const Slot4096 *s) rdi=s -> xmm0
; -----------------------------------------------------------------------------
ap_to_double:
push rbp
mov rbp, rsp
sub rsp, 16
push rbx
mov rbx, rdi
test rbx, rbx
jz .zero
cmp qword [rbx + SL_nw], 0
je .zero
mov rax, [rbx + SL_mwords]
test rax, rax
jz .zero
; unsigned 64->double two-step
mov rdx, [rax]
test rdx, rdx
js .big_uint
cvtsi2sd xmm0, rdx
jmp .div
.big_uint:
mov rcx, rdx
and rcx, 1
shr rdx, 1
or rdx, rcx
cvtsi2sd xmm0, rdx
addsd xmm0, xmm0
.div:
divsd xmm0, [fp_U64]
movsd [rbp-8], xmm0
movsd xmm0, [fp_2]
mov rax, [rbx + SL_exp]
cvtsi2sd xmm1, rax
and rsp, -16
call pow
mulsd xmm0, [rbp-8]
jmp .done
.zero:
xorpd xmm0, xmm0
.done:
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; Slot4096 *ap_from_double(double v, int bm, int be)
; xmm0=v rdi=bm rsi=be -> rax
; -----------------------------------------------------------------------------
ap_from_double:
push rbp
mov rbp, rsp
sub rsp, 48
push rbx
push r12
push r13
movsd [rbp-8], xmm0
mov r12d, edi
mov r13d, esi
and rsp, -16
mov rdi, SL_SIZE
call malloc
test rax, rax
jz .null
mov rbx, rax
; init with dimn=1, r=0.5, omega=1
mov rdi, rbx
mov esi, r12d
mov edx, r13d
mov ecx, 1
movsd xmm0, [fp_half]
movsd xmm1, [fp_1]
call slot_init
; if v==0 we're done
movsd xmm0, [rbp-8]
xorpd xmm1, xmm1
ucomisd xmm0, xmm1
jz .ret
; frexp
lea rsi, [rbp-16]
call frexp
movsd [rbp-24], xmm0
movsx rax, dword [rbp-16]
mov [rbx + SL_exp], rax
; mantissa
movsd xmm0, [rbp-24]
call fabs
mulsd xmm0, [fp_U64]
maxsd xmm0, [fp_0]
cvttsd2si rax, xmm0
mov rcx, [rbx + SL_mwords]
test rcx, rcx
jz .sgn
mov [rcx], rax
.sgn:
movsd xmm0, [rbp-8]
xorpd xmm1, xmm1
ucomisd xmm0, xmm1
jae .sgn_done
or dword [rbx + SL_flags], FLAG_SIGN_NEG
.sgn_done:
; sync empi
mov rax, [rbx + SL_exp]
mov rsi, rax
xor dl, dl
test rax, rax
jns .ep
neg rsi
mov dl, 1
.ep:
lea rdi, [rbx + SL_empi]
call mpi_set
.ret:
mov rax, rbx
jmp .done
.null:
xor rax, rax
.done:
pop r13
pop r12
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; void ap_shift_right(uint64_t *w, size_t nw, int64_t shift)
; rdi=w rsi=nw rdx=shift
; -----------------------------------------------------------------------------
ap_shift_right:
test rdx, rdx
jle .ret
test rsi, rsi
jz .ret
push rbp
mov rbp, rsp
push rbx
push r12
push r13
push r14
mov rbx, rdi
mov r12, rsi
mov r13, rdx
; if shift >= nw*64: zero all
mov rax, r12
shl rax, 6
cmp r13, rax
jl .do
and rsp, -16
mov rdi, rbx
xor rsi, rsi
mov rdx, r12
shl rdx, 3
call memset
jmp .pop
.do:
; word_shift = shift/64 bit_shift = shift%64
mov rax, r13
xor edx, edx
mov ecx, 64
div ecx
mov r14, rax ; word_shift
mov r8d, edx ; bit_shift
; word-level shift (move words right)
test r14, r14
jz .bits
mov rax, r12
dec rax ; i = nw-1
.wl:
cmp rax, r14
jl .wl_done
mov rcx, rax
sub rcx, r14
mov r9, [rbx + rcx*8]
mov [rbx + rax*8], r9
dec rax
jmp .wl
.wl_done:
; zero the first word_shift words
and rsp, -16
mov rdi, rbx
xor rsi, rsi
mov rdx, r14
shl rdx, 3
call memset
.bits:
test r8d, r8d
jz .pop
; bit-level shift right, carrying upper bits from left neighbour
mov ecx, r8d ; bit_shift in cl
mov r9d, 64
sub r9d, r8d ; 64-bit_shift in r9d
mov rax, r12
dec rax ; i = nw-1
.bl:
cmp rax, 1
jl .blast
mov r10, [rbx + rax*8]
shr r10, cl ; w[i] >>= bit_shift
mov rdx, rax
dec rdx
mov r11, [rbx + rdx*8]
push rcx
mov cl, r9b ; 64-bit_shift
shl r11, cl ; carry from w[i-1]
pop rcx
or r10, r11
mov [rbx + rax*8], r10
dec rax
jmp .bl
.blast:
shr qword [rbx], cl ; w[0] >>= bit_shift
.pop:
pop r14
pop r13
pop r12
pop rbx
leave
.ret:
ret
; -----------------------------------------------------------------------------
; void ap_normalize(Slot4096 *s) rdi=s
; -----------------------------------------------------------------------------
ap_normalize:
push rbp
mov rbp, rsp
push rbx
push r12
mov rbx, rdi
cmp qword [rbx + SL_nw], 0
je .done
mov r12, [rbx + SL_mwords]
test r12, r12
jz .done
.loop:
mov rax, [r12]
test rax, rax
js .normed
jz .zexp
; underflow guard
cmp qword [rbx + SL_exp], -32768
jle .guz
; left-shift mantissa by 1
mov rcx, [rbx + SL_nw]
dec rcx
xor r8, r8 ; carry
.sl:
js .sld
mov r9, [r12 + rcx*8]
mov r10, r9
shr r10, 63
shl r9, 1
or r9, r8
mov [r12 + rcx*8], r9
mov r8, r10
dec rcx
jmp .sl
.sld:
dec qword [rbx + SL_exp]
jmp .loop
.guz:
or dword [rbx + SL_flags], FLAG_GUZ
jmp .normed
.zexp:
mov qword [rbx + SL_exp], 0
.normed:
; sync empi
mov rax, [rbx + SL_exp]
mov rsi, rax
xor dl, dl
test rax, rax
jns .ep
neg rsi
mov dl, 1
.ep:
lea rdi, [rbx + SL_empi]
call mpi_set
.done:
pop r12
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; void ap_add(Slot4096 *A, const Slot4096 *B) rdi=A rsi=B
; -----------------------------------------------------------------------------
ap_add:
push rbp
mov rbp, rsp
sub rsp, SL_SIZE + 16
push rbx
push r12
mov rbx, rdi
mov r12, rsi
; check word count match
mov rax, [rbx + SL_nw]
cmp rax, [r12 + SL_nw]
je .ok
and rsp, -16
lea rdi, [sz_werr]
xor eax, eax
call printf
jmp .done
.ok:
; B_aligned on stack
lea r8, [rbp - SL_SIZE - 8]
and rsp, -16
mov rdi, r8
xor rsi, rsi
mov rdx, SL_SIZE
call memset
mov rdi, r8
mov rsi, r12
call ap_copy
mov r12, r8 ; B_aligned
; align exponents
mov rax, [rbx + SL_exp]
sub rax, [r12 + SL_exp]
jz .add_words
jg .shift_b
; shift A right
neg rax
mov rdi, [rbx + SL_mwords]
mov rsi, [rbx + SL_nw]
mov rdx, rax
call ap_shift_right
mov rax, [r12 + SL_exp]
mov [rbx + SL_exp], rax
jmp .add_words
.shift_b:
mov rdx, rax
mov rdi, [r12 + SL_mwords]
mov rsi, [r12 + SL_nw]
call ap_shift_right
.add_words:
; ripple-carry add from highest to lowest word
mov rcx, [rbx + SL_nw]
dec rcx
mov rdi, [rbx + SL_mwords]
mov rsi, [r12 + SL_mwords]
xor r8b, r8b
.al:
js .ad
mov rax, [rdi + rcx*8]
mov rdx, [rsi + rcx*8]
add rax, rdx
setc r9b
add rax, r8
setc r10b
or r9b, r10b
mov [rdi + rcx*8], rax
mov r8b, r9b
dec rcx
jmp .al
.ad:
test r8b, r8b
jz .normalize
or dword [rbx + SL_flags], FLAG_GOI
.normalize:
mov rdi, rbx
call ap_normalize
; free B_aligned's heap data
mov rdi, [r12 + SL_mwords]
test rdi, rdi
jz .fmpi
and rsp, -16
call free
.fmpi:
lea rdi, [r12 + SL_empi]
call mpi_free
lea rdi, [r12 + SL_nwm]
call mpi_free
lea rdi, [r12 + SL_soi]
call mpi_free
.done:
pop r12
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; void exchange_analog_links(AnalogLink *L, int rank, int size, int n)
; rdi=L rsi=rank rdx=size rcx=n (MPI_REAL=0 branch: simple decay)
; -----------------------------------------------------------------------------
exchange_analog_links:
xor eax, eax
.loop:
cmp eax, ecx
jge .ret
imul r8, rax, AL_SIZE
add r8, rdi
movsd xmm0, [r8 + AL_chg]
mulsd xmm0, [fp_0p95]
movsd [r8 + AL_chg], xmm0
movsd xmm0, [r8 + AL_chgim]
mulsd xmm0, [fp_0p95]
movsd [r8 + AL_chgim], xmm0
movsd xmm0, [r8 + AL_dncoup]
mulsd xmm0, [fp_0p98]
movsd [r8 + AL_dncoup], xmm0
inc eax
jmp .loop
.ret:
ret
; -----------------------------------------------------------------------------
; void fill_one_neigh(HDGLLattice *lat, Slot4096 *myslot, int nidx,
; AnalogLink *link)
; rdi=lat rsi=myslot rdx=nidx rcx=link
; Computes one neighbor entry; does nothing if neigh slot is null.
; -----------------------------------------------------------------------------
fill_one_neigh:
push rbp
mov rbp, rsp
push rbx
push r12
push r13
push r14
mov rbx, rdi ; lat
mov r12, rsi ; my slot
mov r13d, edx ; neigh index
mov r14, rcx ; &link
and rsp, -16
mov rdi, rbx
mov esi, r13d
call lattice_get_slot
test rax, rax
jz .done
mov r13, rax ; neigh slot
; charge = ap_to_double(neigh)
mov rdi, r13
call ap_to_double
movsd [r14 + AL_chg], xmm0
movsd xmm0, [r13 + SL_ampim]
movsd [r14 + AL_chgim], xmm0
; potential = neigh->phase - myslot->phase
movsd xmm0, [r13 + SL_phase]
subsd xmm0, [r12 + SL_phase]
movsd [r14 + AL_pot], xmm0
; Dn_coupling = neigh->Dn * exp(-|neigh->Dn - myslot->Dn|)
movsd xmm0, [r13 + SL_dnamp]
subsd xmm0, [r12 + SL_dnamp]
call fabs
movsd xmm1, [fp_NEG1]
mulsd xmm0, xmm1
push r13
push r14
call exp
pop r14
pop r13
mulsd xmm0, [r13 + SL_dnamp]
movsd [r14 + AL_dncoup], xmm0
; coupling = 1.0 (simplified from amp_correlation formula)
movsd xmm0, [fp_1]
movsd [r14 + AL_coup], xmm0
.done:
pop r14
pop r13
pop r12
pop rbx
leave
ret
; forward declaration used by fill_one_neigh and lattice_integrate_rk4
; (lattice_get_slot defined later)
; -----------------------------------------------------------------------------
; void deriv_Dn(ComplexState *out, const ComplexState *st,
; double omega, const AnalogLink *nb, int nnb,
; int dimn, double wmode)
; rdi=out rsi=st xmm0=omega rdx=nb rcx=nnb r8d=dimn xmm1=wmode
; -----------------------------------------------------------------------------
deriv_Dn:
push rbp
mov rbp, rsp
sub rsp, 64
push rbx
push r12
push r13
push r14
mov rbx, rdi ; out
mov r12, rsi ; st
movsd [rbp-8], xmm0 ; omega
mov r13, rdx ; nb
mov r14d, ecx ; nnb
movsd [rbp-16], xmm1 ; wmode
; zero out
and rsp, -16
mov rdi, rbx
xor rsi, rsi
mov rdx, CS_SIZE
call memset
; A = sqrt(Ar^2 + Ai^2)
movsd xmm0, [r12 + CS_Ar]
mulsd xmm0, xmm0
movsd xmm1, [r12 + CS_Ai]
mulsd xmm1, xmm1
addsd xmm0, xmm1
call sqrt
movsd [rbp-24], xmm0 ; A
; deriv.Ar = -GAMMA*Ar + 0.1*Dn*cos(phase)
movsd xmm0, [r12 + CS_ph]
call cos
mulsd xmm0, [r12 + CS_dn]
mulsd xmm0, [fp_01]
movsd xmm1, [fp_GAMMA]
mulsd xmm1, [r12 + CS_Ar]
subsd xmm0, xmm1
movsd [rbx + CS_Ar], xmm0
; deriv.Ai = -GAMMA*Ai + 0.1*Dn*sin(phase)
movsd xmm0, [r12 + CS_ph]
call sin
mulsd xmm0, [r12 + CS_dn]
mulsd xmm0, [fp_01]
movsd xmm1, [fp_GAMMA]
mulsd xmm1, [r12 + CS_Ai]
subsd xmm0, xmm1
movsd [rbx + CS_Ai], xmm0
; neighbour coupling loop
xorpd xmm13, xmm13 ; sum_sin = 0
xor r9d, r9d ; k = 0
.nb_loop:
cmp r9d, r14d
jge .nb_done
imul rax, r9, AL_SIZE
lea r10, [r13 + rax] ; &nb[k]
; delta_phi = nb[k].potential - st.phase
movsd xmm0, [r10 + AL_pot]
subsd xmm0, [r12 + CS_ph]
movsd [rbp-32], xmm0 ; delta_phi
; sum_sin += sin(delta_phi)
call sin
addsd xmm13, xmm0
; Dn_factor = nb[k].dncoup / (1 + |st.Dn|)
movsd xmm0, [r12 + CS_dn]
call fabs
addsd xmm0, [fp_1]
movsd xmm1, [r10 + AL_dncoup]
divsd xmm1, xmm0 ; Dn_factor
movsd [rbp-40], xmm1
; deriv.Ar += Kc * Dn_factor * charge * cos(delta_phi)
movsd xmm0, [rbp-32]
call cos
mulsd xmm0, [r10 + AL_chg]
mulsd xmm0, [rbp-40]
mulsd xmm0, [fp_KC]
addsd xmm0, [rbx + CS_Ar]
movsd [rbx + CS_Ar], xmm0
; deriv.Ai += Kc * Dn_factor * charge_im * sin(delta_phi)
movsd xmm0, [rbp-32]
call sin
mulsd xmm0, [r10 + AL_chgim]
mulsd xmm0, [rbp-40]
mulsd xmm0, [fp_KC]
addsd xmm0, [rbx + CS_Ai]
movsd [rbx + CS_Ai], xmm0
inc r9d
jmp .nb_loop
.nb_done:
; deriv.phase_vel = omega + Kc*sum_sin + 0.3*wmode
movsd xmm0, [rbp-8]
mulsd xmm13, [fp_KC]
addsd xmm0, xmm13
movsd xmm1, [rbp-16]
mulsd xmm1, [fp_03]
addsd xmm0, xmm1
movsd [rbx + CS_pv], xmm0
; deriv.phase = st.phase_vel
movsd xmm0, [r12 + CS_pv]
movsd [rbx + CS_ph], xmm0
; deriv.Dn = -0.01*(Dn_val - A)
movsd xmm0, [r12 + CS_dn]
subsd xmm0, [rbp-24]
mulsd xmm0, [fp_001]
movsd xmm1, xmm0
xorpd xmm0, xmm0
subsd xmm0, xmm1
movsd [rbx + CS_dn], xmm0
pop r14
pop r13
pop r12
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; void rk4_step(Slot4096 *slot, double dt, const AnalogLink *nb, int nnb,
; double omega)
; rdi=slot xmm0=dt rsi=nb rdx=nnb xmm1=omega
;
; Stack layout (from rbp downward):
; [rbp- 8] dt
; [rbp-16] omega
; [rbp-24] saved A (amplitude after damping)
; [rbp-32] norm
; state at [rbp - 240]
; tmp at [rbp - 200]
; k1 at [rbp - 160]
; k2 at [rbp - 120]
; k3 at [rbp - 80]
; k4 at [rbp - 40]
; Total CS slots: 6 * 40 = 240, plus 32 saved doubles = 272, round to 288
; -----------------------------------------------------------------------------
rk4_step:
push rbp
mov rbp, rsp
sub rsp, 288
push rbx
push r12
push r13
push r14
mov rbx, rdi ; slot
movsd [rbp-8], xmm0 ; dt
mov r12, rsi ; nb
mov r13d, edx ; nnb
movsd [rbp-16], xmm1 ; omega
; offsets relative to rbp (all negative)
; state = rbp-240, tmp = rbp-200, k1 = rbp-160, k2 = rbp-120,
; k3 = rbp-80, k4 = rbp-40
%define ST_OFF 240
%define TMP_OFF 200
%define K1_OFF 160
%define K2_OFF 120
%define K3_OFF 80
%define K4_OFF 40
; --- build initial ComplexState ---
and rsp, -16
mov rdi, rbx
call ap_to_double
movsd [rbp - ST_OFF + CS_Ar], xmm0
movsd xmm0, [rbx + SL_ampim]
movsd [rbp - ST_OFF + CS_Ai], xmm0
movsd xmm0, [rbx + SL_phase]
movsd [rbp - ST_OFF + CS_ph], xmm0
movsd xmm0, [rbx + SL_pvel]
movsd [rbp - ST_OFF + CS_pv], xmm0
movsd xmm0, [rbx + SL_dnamp]
movsd [rbp - ST_OFF + CS_dn], xmm0
; --- k1 = deriv(state, ...) ---
lea rdi, [rbp - K1_OFF]
lea rsi, [rbp - ST_OFF]
movsd xmm0, [rbp-16]
mov rdx, r12
mov ecx, r13d
mov r8d, [rbx + SL_dimn]
movsd xmm1, [rbx + SL_wmode]
call deriv_Dn
; --- tmp = state + (dt/2)*k1 ---
movsd xmm15, [rbp-8]
mulsd xmm15, [fp_half] ; dt/2
%macro LERP5 3 ; dst_neg_off, src_neg_off, k_neg_off
movsd xmm0, [rbp - %3 + CS_Ar]
mulsd xmm0, xmm15
addsd xmm0, [rbp - %2 + CS_Ar]
movsd [rbp - %1 + CS_Ar], xmm0
movsd xmm0, [rbp - %3 + CS_Ai]
mulsd xmm0, xmm15
addsd xmm0, [rbp - %2 + CS_Ai]
movsd [rbp - %1 + CS_Ai], xmm0
movsd xmm0, [rbp - %3 + CS_ph]
mulsd xmm0, xmm15
addsd xmm0, [rbp - %2 + CS_ph]
movsd [rbp - %1 + CS_ph], xmm0
movsd xmm0, [rbp - %3 + CS_pv]
mulsd xmm0, xmm15
addsd xmm0, [rbp - %2 + CS_pv]
movsd [rbp - %1 + CS_pv], xmm0
movsd xmm0, [rbp - %3 + CS_dn]
mulsd xmm0, xmm15
addsd xmm0, [rbp - %2 + CS_dn]
movsd [rbp - %1 + CS_dn], xmm0
%endmacro
LERP5 TMP_OFF, ST_OFF, K1_OFF ; tmp = st + dt/2*k1
; --- k2 = deriv(tmp, ...) ---
lea rdi, [rbp - K2_OFF]
lea rsi, [rbp - TMP_OFF]
movsd xmm0, [rbp-16]
mov rdx, r12
mov ecx, r13d
mov r8d, [rbx + SL_dimn]
movsd xmm1, [rbx + SL_wmode]
call deriv_Dn
LERP5 TMP_OFF, ST_OFF, K2_OFF ; tmp = st + dt/2*k2
; --- k3 = deriv(tmp, ...) ---
lea rdi, [rbp - K3_OFF]
lea rsi, [rbp - TMP_OFF]
movsd xmm0, [rbp-16]
mov rdx, r12
mov ecx, r13d
mov r8d, [rbx + SL_dimn]
movsd xmm1, [rbx + SL_wmode]
call deriv_Dn
; tmp = st + dt*k3 (full dt)
movsd xmm15, [rbp-8] ; dt (not halved)
LERP5 TMP_OFF, ST_OFF, K3_OFF
; --- k4 = deriv(tmp, ...) ---
lea rdi, [rbp - K4_OFF]
lea rsi, [rbp - TMP_OFF]
movsd xmm0, [rbp-16]
mov rdx, r12
mov ecx, r13d
mov r8d, [rbx + SL_dimn]
movsd xmm1, [rbx + SL_wmode]
call deriv_Dn
; --- weighted sum: state += (dt/6)*(k1 + 2*k2 + 2*k3 + k4) ---
movsd xmm15, [rbp-8]
divsd xmm15, [fp_6] ; dt/6
%macro RK4SUM 1 ; CS field offset
movsd xmm0, [rbp - K1_OFF + %1]
movsd xmm1, [rbp - K2_OFF + %1]
addsd xmm1, xmm1
addsd xmm0, xmm1
movsd xmm1, [rbp - K3_OFF + %1]
addsd xmm1, xmm1
addsd xmm0, xmm1
addsd xmm0, [rbp - K4_OFF + %1]
mulsd xmm0, xmm15
addsd xmm0, [rbp - ST_OFF + %1]
movsd [rbp - ST_OFF + %1], xmm0
%endmacro
RK4SUM CS_Ar
RK4SUM CS_Ai
RK4SUM CS_ph
RK4SUM CS_pv
RK4SUM CS_dn
; --- entropy dampers ---
; A = sqrt(Ar^2+Ai^2)
movsd xmm0, [rbp - ST_OFF + CS_Ar]
mulsd xmm0, xmm0
movsd xmm1, [rbp - ST_OFF + CS_Ai]
mulsd xmm1, xmm1
addsd xmm0, xmm1
call sqrt
; A *= exp(-LAMBDA*dt)
movsd [rbp-24], xmm0
movsd xmm0, [fp_LAMBDA]
mulsd xmm0, [rbp-8]
mulsd xmm0, [fp_NEG1]
call exp
mulsd xmm0, [rbp-24]
; clamp to SAT_LIMIT
minsd xmm0, [fp_SAT]
; add noise: A += NSIG*(2*rand-1)
movsd [rbp-24], xmm0
call get_normalized_rand
addsd xmm0, xmm0
subsd xmm0, [fp_1]
mulsd xmm0, [fp_NSIG]
addsd xmm0, [rbp-24]
movsd [rbp-24], xmm0 ; final A
; renormalize direction
movsd xmm0, [rbp - ST_OFF + CS_Ar]
mulsd xmm0, xmm0
movsd xmm1, [rbp - ST_OFF + CS_Ai]
mulsd xmm1, xmm1
addsd xmm0, xmm1
call sqrt
movsd [rbp-32], xmm0 ; norm
ucomisd xmm0, [fp_1e10n]
jbe .no_renorm
movsd xmm1, [rbp-24] ; A
divsd xmm1, xmm0 ; A/norm
movsd xmm0, [rbp - ST_OFF + CS_Ar]
mulsd xmm0, xmm1
movsd [rbp - ST_OFF + CS_Ar], xmm0
movsd xmm0, [rbp - ST_OFF + CS_Ai]
mulsd xmm0, xmm1
movsd [rbp - ST_OFF + CS_Ai], xmm0
.no_renorm:
; wrap phase into [0, 2pi)
movsd xmm0, [rbp - ST_OFF + CS_ph]
movsd xmm1, [fp_TWO_PI]
call fmod
movsd [rbp - ST_OFF + CS_ph], xmm0
xorpd xmm1, xmm1
ucomisd xmm0, xmm1
jae .ph_ok
addsd xmm0, [fp_TWO_PI]
movsd [rbp - ST_OFF + CS_ph], xmm0
.ph_ok:
; clamp Dn to [0, 1000]
movsd xmm0, [rbp - ST_OFF + CS_dn]
maxsd xmm0, [fp_0]
minsd xmm0, [fp_1000]
movsd [rbp - ST_OFF + CS_dn], xmm0
; --- write back ---
movsd xmm0, [rbp - ST_OFF + CS_Ar]
mov edi, [rbx + SL_bmant]
mov esi, [rbx + SL_bexp]
call ap_from_double
test rax, rax
jz .skip_wb
push rax
mov rdi, rbx
mov rsi, rax
call ap_copy
pop rdi
call ap_free
and rsp, -16
call free
.skip_wb:
movsd xmm0, [rbp - ST_OFF + CS_Ai]
movsd [rbx + SL_ampim], xmm0
movsd xmm0, [rbp - ST_OFF + CS_ph]
movsd [rbx + SL_phase], xmm0
movsd xmm0, [rbp - ST_OFF + CS_pv]
movsd [rbx + SL_pvel], xmm0
movsd xmm0, [rbp - ST_OFF + CS_dn]
movsd [rbx + SL_dnamp], xmm0
pop r14
pop r13
pop r12
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; HDGLLattice *lattice_init(int ninst, int spi) rdi=ninst rsi=spi
; -----------------------------------------------------------------------------
lattice_init:
push rbp
mov rbp, rsp
push rbx
push r12
push r13
mov r12d, edi
mov r13d, esi
and rsp, -16
mov rdi, LT_SIZE
call malloc
test rax, rax
jz .null
mov rbx, rax
mov rdi, rbx
xor rsi, rsi
mov rdx, LT_SIZE
call memset
mov [rbx + LT_ninst], r12d
mov [rbx + LT_spi], r13d
movsd xmm0, [fp_1]
movsd [rbx + LT_omega], xmm0
movsd xmm0, [fp_0]
movsd [rbx + LT_time], xmm0
movsd xmm0, [fp_1000]
movsd [rbx + LT_pvar], xmm0
call get_rtc_ns
mov [rbx + LT_ckns], rax
; numeric lattice
mov rdi, NL_SIZE
call malloc
test rax, rax
jz .fail_nl
mov [rbx + LT_nl], rax
mov rdi, rax
call init_numeric_lattice
; num_chunks = ceil(ninst*spi / CHUNK_SIZE)
mov eax, r12d
imul eax, r13d
add eax, CHUNK_SIZE - 1
xor edx, edx
mov ecx, CHUNK_SIZE
div ecx
mov [rbx + LT_nchunks], eax
; calloc chunk pointer array
mov rdi, rax
mov rsi, 8
call calloc
test rax, rax
jz .fail_chunks
mov [rbx + LT_chunks], rax
mov rax, rbx
jmp .done
.fail_chunks:
mov rdi, [rbx + LT_nl]
call free_numeric_lattice
mov rdi, [rbx + LT_nl]
call free
.fail_nl:
mov rdi, rbx
call free
.null:
xor rax, rax
.done:
pop r13
pop r12
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; HDGLChunk *lattice_get_chunk(HDGLLattice *lat, int ci)
; rdi=lat rsi=ci -> rax
; -----------------------------------------------------------------------------
lattice_get_chunk:
push rbp
mov rbp, rsp
push rbx
push r12
push r13
push r14
mov rbx, rdi
mov r12d, esi
cmp r12d, [rbx + LT_nchunks]
jge .null
mov rax, [rbx + LT_chunks]
mov r13, [rax + r12*8]
test r13, r13
jnz .existing
; allocate chunk header
and rsp, -16
mov rdi, CH_SIZE
call malloc
test rax, rax
jz .null
mov r13, rax
; allocate slots array
mov rdi, CHUNK_SIZE
imul rdi, rdi, SL_SIZE
call malloc
test rax, rax
jz .fail_slots
mov [r13 + CH_slots], rax
mov qword [r13 + CH_alloc], CHUNK_SIZE
; initialise every slot
xor r14d, r14d
.init:
cmp r14d, CHUNK_SIZE
jge .init_done
; bmant = 4096 + (i%8)*64
mov eax, r14d
and eax, 7
shl eax, 6
add eax, 4096
mov r8d, eax
; bexp = 16 + (i%8)*2
mov eax, r14d
and eax, 7
shl eax, 1
add eax, 16
mov r9d, eax
; dimn = (i%NUM_DN)+1
mov eax, r14d
xor edx, edx
mov ecx, NUM_DN
div ecx
inc edx
mov r10d, edx
; r_val = (i%256)/256.0
mov eax, r14d
and eax, 255
cvtsi2sd xmm0, eax
mov eax, 256
cvtsi2sd xmm1, eax
divsd xmm0, xmm1
; omega from lattice
movsd xmm1, [rbx + LT_omega]
; slot pointer
imul rax, r14, SL_SIZE
add rax, [r13 + CH_slots]
push r14
push r13
push rbx
mov rdi, rax
mov esi, r8d
mov edx, r9d
mov ecx, r10d
call slot_init
pop rbx
pop r13
pop r14
inc r14d
jmp .init
.init_done:
; store in chunk array
mov rax, [rbx + LT_chunks]
mov [rax + r12*8], r13
mov rax, r13
jmp .done
.fail_slots:
mov rdi, r13
call free
.null:
xor rax, rax
jmp .done
.existing:
mov rax, r13
.done:
pop r14
pop r13
pop r12
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; Slot4096 *lattice_get_slot(HDGLLattice *lat, int idx)
; rdi=lat rsi=idx -> rax
; -----------------------------------------------------------------------------
lattice_get_slot:
push rbp
mov rbp, rsp
push rbx
push r12
mov rbx, rdi
mov r12d, esi
mov eax, r12d
xor edx, edx
mov ecx, CHUNK_SIZE
div ecx ; eax=chunk_idx edx=local
push rdx
mov rdi, rbx
mov esi, eax
call lattice_get_chunk
pop rdx
test rax, rax
jz .null
mov rcx, [rax + CH_slots]
imul rdx, rdx, SL_SIZE
add rcx, rdx
mov rax, rcx
jmp .done
.null:
xor rax, rax
.done:
pop r12
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; void detect_consensus(HDGLLattice *lat) rdi=lat
; -----------------------------------------------------------------------------
detect_consensus:
push rbp
mov rbp, rsp
sub rsp, 32
push rbx
push r12
push r13
mov rbx, rdi
; total slots
mov eax, [rbx + LT_ninst]
imul eax, [rbx + LT_spi]
mov r12d, eax
; --- pass 1: mean_phase ---
xorpd xmm12, xmm12
xor r13d, r13d ; count
xor ecx, ecx
.p1:
cmp ecx, r12d
jge .p1_done
push rcx
push r13
mov rdi, rbx
mov esi, ecx
call lattice_get_slot
pop r13
pop rcx
test rax, rax
jz .p1n
test dword [rax + SL_flags], FLAG_CONSENSUS
jnz .p1n
addsd xmm12, [rax + SL_phase]
inc r13d
.p1n:
inc ecx
jmp .p1
.p1_done:
test r13d, r13d
jz .done
cvtsi2sd xmm0, r13d
divsd xmm12, xmm0 ; mean_phase
; --- pass 2: variance ---
xorpd xmm11, xmm11
xor ecx, ecx
.p2:
cmp ecx, r12d
jge .p2_done
push rcx
push r13
mov rdi, rbx
mov esi, ecx
call lattice_get_slot
pop r13
pop rcx
test rax, rax
jz .p2n
test dword [rax + SL_flags], FLAG_CONSENSUS
jnz .p2n
movsd xmm0, [rax + SL_phase]
subsd xmm0, xmm12
mulsd xmm0, xmm0
addsd xmm11, xmm0
.p2n:
inc ecx
jmp .p2
.p2_done:
cvtsi2sd xmm0, r13d
movsd xmm1, xmm11
divsd xmm1, xmm0
movsd xmm0, xmm1
call sqrt
movsd [rbx + LT_pvar], xmm0
ucomisd xmm0, [fp_CEPS]
ja .no_cons
inc dword [rbx + LT_csteps]
cmp dword [rbx + LT_csteps], CONSENSUS_N
jl .done
; announce
and rsp, -16
lea rdi, [sz_cons]
movsd xmm0, [rbx + LT_time]
movsd xmm1, [rbx + LT_pvar]
mov eax, 2
call printf
; lock all slots
xor ecx, ecx
.lock:
cmp ecx, r12d
jge .locked
push rcx
mov rdi, rbx
mov esi, ecx
call lattice_get_slot
pop rcx
test rax, rax
jz .lockn
test dword [rax + SL_flags], FLAG_CONSENSUS
jnz .lockn
or dword [rax + SL_flags], FLAG_CONSENSUS
movsd xmm0, [fp_0]
movsd [rax + SL_pvel], xmm0
.lockn:
inc ecx
jmp .lock
.locked:
mov dword [rbx + LT_csteps], 0
jmp .done
.no_cons:
mov dword [rbx + LT_csteps], 0
.done:
pop r13
pop r12
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; void lattice_integrate_rk4(HDGLLattice *lat, double dt)
; rdi=lat xmm0=dt
; -----------------------------------------------------------------------------
lattice_integrate_rk4:
push rbp
mov rbp, rsp
; reserve space for neighbors[8] on stack (8 * AL_SIZE = 384) + saved vars
sub rsp, 384 + 64
push rbx
push r12
push r13
push r14
mov rbx, rdi
movsd [rbp - 384 - 8], xmm0 ; dt
; total slots
mov eax, [rbx + LT_ninst]
imul eax, [rbx + LT_spi]
mov r12d, eax
; base of neighbors array on stack
lea r14, [rbp - 384 - 56] ; 8-byte aligned scratch below saves
xor r13d, r13d ; i = 0
.main:
cmp r13d, r12d
jge .main_done
; get slot
push r13
and rsp, -16
mov rdi, rbx
mov esi, r13d
call lattice_get_slot
pop r13
test rax, rax
jz .next
; skip flagged slots
test dword [rax + SL_flags], FLAG_GOI | FLAG_IS_NAN | FLAG_CONSENSUS
jnz .next
mov [rbp - 384 - 16], rax ; save slot ptr
; zero neighbors
and rsp, -16
mov rdi, r14
xor rsi, rsi
mov rdx, 8 * AL_SIZE
call memset
; --- fill 8 neighbors ---
; safe modulo helper: ((raw + total) % total) for raw = i +/- offset
; We call fill_one_neigh(lat, myslot, neigh_idx, &nb[j])
mov r8, [rbp - 384 - 16] ; my slot ptr
mov r9d, [rbx + LT_spi] ; slots_per_instance
%macro NEIGH 2 ; j_index, raw_expr_leaves_eax
%2
; positive modulo: ((eax % r12d) + r12d) % r12d
cdq
idiv r12d ; edx = signed remainder
add edx, r12d
xor eax, eax
mov eax, edx
xor edx, edx
div r12d ; edx = final non-negative index
lea rcx, [r14 + %1 * AL_SIZE]
push r8
push r9
push r13
push r12
push rbx
and rsp, -16
mov rdi, rbx
mov rsi, r8
mov edx_hold, edx
lea rcx, [r14 + %1 * AL_SIZE]
mov edx, edx_hold
call fill_one_neigh
pop rbx
pop r12
pop r13
pop r9
pop r8
%endmacro
; NASM does not support edx_hold as a pseudo-register;
; use r15d as temp for the neighbour index across the call:
%macro NEIGH 2
%2
cdq
idiv r12d
add edx, r12d
mov eax, edx
xor edx, edx
div r12d ; edx = neigh index (unsigned mod)
push r8
push r9
push r13
push r12
push rbx
and rsp, -16
mov rdi, rbx
mov rsi, r8
mov r15d, edx ; save neigh index
lea rcx, [r14 + (%1) * AL_SIZE]
mov edx, r15d
call fill_one_neigh
pop rbx
pop r12
pop r13
pop r9
pop r8
%endmacro
push r15 ; save r15
NEIGH 0, <mov eax, r13d; dec eax>
NEIGH 1, <mov eax, r13d; inc eax>
NEIGH 2, <mov eax, r13d; sub eax, r9d>
NEIGH 3, <mov eax, r13d; add eax, r9d>
NEIGH 4, <mov eax, r13d; sub eax, r9d; dec eax>
NEIGH 5, <mov eax, r13d; sub eax, r9d; inc eax>
NEIGH 6, <mov eax, r13d; add eax, r9d; dec eax>
NEIGH 7, <mov eax, r13d; add eax, r9d; inc eax>
pop r15
; --- exchange_analog_links ---
mov rdi, r14
mov eax, r13d
xor edx, edx
idiv dword [rbx + LT_ninst]
mov esi, edx
mov edx, [rbx + LT_ninst]
mov ecx, 8
call exchange_analog_links
; --- rk4_step(slot, dt, nb, 8, omega) ---
mov rdi, [rbp - 384 - 16]
movsd xmm0, [rbp - 384 - 8]
mov rsi, r14
mov edx, 8
movsd xmm1, [rbx + LT_omega]
call rk4_step
; --- phi-adaptive dt ---
mov rdi, [rbp - 384 - 16]
call ap_to_double
call fabs
ucomisd xmm0, [fp_ATH]
jbe .chk_small
movsd xmm0, [rbp - 384 - 8]
mulsd xmm0, [fp_PHI]
movsd [rbp - 384 - 8], xmm0
jmp .clamp_dt
.chk_small:
movsd xmm1, [fp_ATH]
divsd xmm1, [fp_PHI]
ucomisd xmm0, xmm1
jae .clamp_dt
movsd xmm0, [rbp - 384 - 8]
divsd xmm0, [fp_PHI]
movsd [rbp - 384 - 8], xmm0
.clamp_dt:
movsd xmm0, [rbp - 384 - 8]
maxsd xmm0, [fp_1e6n]
minsd xmm0, [fp_01]
movsd [rbp - 384 - 8], xmm0
.next:
inc r13d
jmp .main
.main_done:
; detect_consensus
and rsp, -16
mov rdi, rbx
call detect_consensus
; omega += 0.01*dt
movsd xmm0, [fp_001]
mulsd xmm0, [rbp - 384 - 8]
addsd xmm0, [rbx + LT_omega]
movsd [rbx + LT_omega], xmm0
; time += dt
movsd xmm0, [rbp - 384 - 8]
addsd xmm0, [rbx + LT_time]
movsd [rbx + LT_time], xmm0
pop r14
pop r13
pop r12
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; CheckpointManager *checkpoint_init() -> rax
; -----------------------------------------------------------------------------
checkpoint_init:
push rbp
mov rbp, rsp
push rbx
and rsp, -16
mov rdi, CK_SIZE
call malloc
mov rbx, rax
mov rdi, SNAPSHOT_MAX * CM_SIZE
call malloc
mov [rbx + CK_snaps], rax
mov dword [rbx + CK_cnt], 0
mov dword [rbx + CK_cap], SNAPSHOT_MAX
mov rax, rbx
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; void checkpoint_add(CheckpointManager *mgr, int evo, HDGLLattice *lat)
; rdi=mgr rsi=evo rdx=lat
; -----------------------------------------------------------------------------
checkpoint_add:
push rbp
mov rbp, rsp
push rbx
push r12
push r13
push r14
mov rbx, rdi
mov r12d, esi
mov r13, rdx
; evict min-weight entry if full
mov eax, [rbx + CK_cnt]
cmp eax, [rbx + CK_cap]
jl .append
mov r8, [rbx + CK_snaps]
movsd xmm12, [r8 + CM_wt]
xor r14d, r14d ; min_idx = 0
mov ecx, 1
.ev:
cmp ecx, [rbx + CK_cnt]
jge .ev_done
imul rax, rcx, CM_SIZE
movsd xmm0, [r8 + rax + CM_wt]
ucomisd xmm0, xmm12
jae .ev_next
movsd xmm12, xmm0
mov r14d, ecx
.ev_next:
inc ecx
jmp .ev
.ev_done:
; shift entries down from min_idx+1
mov ecx, r14d
inc ecx
mov r9d, [rbx + CK_cnt]
.shift:
cmp ecx, r9d
jge .shift_done
imul rax, rcx, CM_SIZE
lea r10, [r8 + rax]
lea r11, [r8 + rax - CM_SIZE]
and rsp, -16
mov rdi, r11
mov rsi, r10
mov rdx, CM_SIZE
call memcpy
inc ecx
jmp .shift
.shift_done:
dec dword [rbx + CK_cnt]
.append:
call get_rtc_ns
mov r8d, [rbx + CK_cnt]
mov r9, [rbx + CK_snaps]
imul r10, r8, CM_SIZE
add r10, r9
mov [r10 + CM_evo], r12d
mov [r10 + CM_tsns], rax
movsd xmm0, [r13 + LT_pvar]
movsd [r10 + CM_pvar], xmm0
movsd xmm0, [r13 + LT_omega]
movsd [r10 + CM_omega], xmm0
movsd xmm0, [fp_1]
movsd [r10 + CM_wt], xmm0
inc dword [rbx + CK_cnt]
; decay older weights
mov ecx, [rbx + CK_cnt]
dec ecx
mov r9, [rbx + CK_snaps]
xor edx, edx
.decay:
cmp edx, ecx
jge .decay_done
imul rax, rdx, CM_SIZE
movsd xmm0, [r9 + rax + CM_wt]
mulsd xmm0, [fp_SDEC]
movsd [r9 + rax + CM_wt], xmm0
inc edx
jmp .decay
.decay_done:
and rsp, -16
lea rdi, [sz_ckpt]
mov esi, r12d
mov edx, [rbx + CK_cnt]
movsd xmm0, [r13 + LT_pvar]
mov eax, 1
call printf
pop r14
pop r13
pop r12
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; void checkpoint_free(CheckpointManager *mgr) rdi=mgr
; -----------------------------------------------------------------------------
checkpoint_free:
push rbp
mov rbp, rsp
push rbx
test rdi, rdi
jz .done
mov rbx, rdi
and rsp, -16
mov rdi, [rbx + CK_snaps]
call free
mov rdi, rbx
call free
.done:
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; void lattice_free(HDGLLattice *lat) rdi=lat
; -----------------------------------------------------------------------------
lattice_free:
push rbp
mov rbp, rsp
push rbx
push r12
push r13
test rdi, rdi
jz .done
mov rbx, rdi
xor r12d, r12d
.chunks:
cmp r12d, [rbx + LT_nchunks]
jge .chunks_done
mov rax, [rbx + LT_chunks]
mov r13, [rax + r12*8]
test r13, r13
jz .cnext
; free each slot's heap members
xor ecx, ecx
.slots:
cmp ecx, CHUNK_SIZE
jge .slots_done
imul rax, rcx, SL_SIZE
add rax, [r13 + CH_slots]
push rcx
push r13
and rsp, -16
mov rdi, rax
call ap_free
pop r13
pop rcx
inc ecx
jmp .slots
.slots_done:
and rsp, -16
mov rdi, [r13 + CH_slots]
call free
mov rdi, r13
call free
.cnext:
inc r12d
jmp .chunks
.chunks_done:
mov rdi, [rbx + LT_chunks]
call free
mov r12, [rbx + LT_nl]
test r12, r12
jz .skip_nl
mov rdi, r12
call free_numeric_lattice
mov rdi, r12
call free
.skip_nl:
mov rdi, rbx
call free
.done:
pop r13
pop r12
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; void bootloader_init_lattice(HDGLLattice *lat, int steps,
; CheckpointManager *mgr)
; rdi=lat rsi=steps rdx=mgr
; -----------------------------------------------------------------------------
bootloader_init_lattice:
push rbp
mov rbp, rsp
sub rsp, 48
push rbx
push r12
push r13
mov rbx, rdi
mov r12d, esi
mov r13, rdx
and rsp, -16
lea rdi, [sz_bi]
xor eax, eax
call printf
test rbx, rbx
jnz .lat_ok
lea rdi, [sz_berr]
xor eax, eax
call printf
jmp .done
.lat_ok:
lea rdi, [sz_binfo]
mov esi, [rbx + LT_ninst]
mov eax, [rbx + LT_ninst]
imul eax, [rbx + LT_spi]
mov edx, eax
xor eax, eax
call printf
lea rdi, [sz_bseed]
mov rsi, [rbx + LT_nl]
mov rsi, [rsi + NL_num_seeds]
xor eax, eax
call printf
; dt = 1/32768
movsd xmm0, [fp_dt0]
movsd [rbp-8], xmm0
; step_ns = 30517
mov qword [rbp-16], 30517
call get_rtc_ns
add rax, [rbp-16]
mov [rbp-24], rax ; next_step_ns
xor r9d, r9d ; i = 0
.evo:
cmp r9d, r12d
jge .evo_done
push r9
push rbx
push r12
push r13
and rsp, -16
mov rdi, rbx
movsd xmm0, [rbp-8]
call lattice_integrate_rk4
pop r13
pop r12
pop rbx
pop r9
; checkpoint every CHECKPOINT_INTERVAL (skip i==0)
test r9d, r9d
jz .no_ckpt
mov eax, r9d
xor edx, edx
mov ecx, CHECKPOINT_INTERVAL
div ecx
test edx, edx
jnz .no_ckpt
push r9
push rbx
push r12
push r13
and rsp, -16
mov rdi, r13
mov esi, r9d
mov rdx, rbx
call checkpoint_add
pop r13
pop r12
pop rbx
pop r9
.no_ckpt:
push r9
push rbx
push r12
push r13
and rsp, -16
mov rdi, [rbp-24]
call rtc_sleep_until
pop r13
pop r12
pop rbx
pop r9
mov rax, [rbp-24]
add rax, [rbp-16]
mov [rbp-24], rax
inc r9d
jmp .evo
.evo_done:
and rsp, -16
lea rdi, [sz_bdone]
mov esi, r12d
xor eax, eax
call printf
lea rdi, [sz_bstat]
movsd xmm0, [rbx + LT_omega]
movsd xmm1, [rbx + LT_time]
movsd xmm2, [rbx + LT_pvar]
mov eax, 3
call printf
.done:
pop r13
pop r12
pop rbx
leave
ret
; -----------------------------------------------------------------------------
; int main(int argc, char **argv)
; -----------------------------------------------------------------------------
main:
push rbp
mov rbp, rsp
sub rsp, 64
push rbx
push r12
push r13
; srand(time(NULL))
and rsp, -16
xor rdi, rdi
call time
mov rdi, rax
call srand
lea rdi, [sz_rtc]
xor eax, eax
call printf
lea rdi, [sz_banner]
xor eax, eax
call printf
; lattice_init(4096, 4)
mov edi, 4096
mov esi, 4
call lattice_init
mov rbx, rax
test rbx, rbx
jnz .lat_ok
lea rdi, [sz_fatal]
xor eax, eax
call printf
mov eax, 1
jmp .exit
.lat_ok:
call checkpoint_init
mov r12, rax
; bootloader_init_lattice(lat, 500, mgr)
mov rdi, rbx
mov esi, 500
mov rdx, r12
call bootloader_init_lattice
; --- Numeric Lattice Summary ---
lea rdi, [sz_nlhdr]
xor eax, eax
call printf
mov r13, [rbx + LT_nl]
lea rdi, [sz_nluf]
movsd xmm0, [r13 + NL_upper + 0]
mov eax, 1
call printf
lea rdi, [sz_nlad]
movsd xmm0, [r13 + NL_adims + 0]
mov eax, 1
call printf
lea rdi, [sz_nlvd]
movsd xmm0, [r13 + NL_void]
mov eax, 1
call printf
lea rdi, [sz_nllf]
movsd xmm0, [r13 + NL_lower + 7*8]
mov eax, 1
call printf
lea rdi, [sz_nlbs]
mov rsi, [r13 + NL_num_seeds]
xor eax, eax
call printf
; --- First 8 slots ---
lea rdi, [sz_slhdr]
xor eax, eax
call printf
xor r13d, r13d
.sl_loop:
cmp r13d, 8
jge .sl_done
push r13
mov rdi, rbx
mov esi, r13d
call lattice_get_slot
pop r13
test rax, rax
jz .sl_next
push rax
push r13
; |A| = sqrt(ap_to_double(s)^2 + amp_im^2)
mov rdi, rax
call ap_to_double
mulsd xmm0, xmm0
mov rax, [rsp+8]
movsd xmm1, [rax + SL_ampim]
mulsd xmm1, xmm1
addsd xmm0, xmm1
call sqrt
movsd [rbp-8], xmm0
mov rax, [rsp+8]
movsd xmm5, [rax + SL_phase]
movsd xmm6, [rax + SL_dnamp]
movsd xmm7, [rax + SL_wmode]
movsd xmm8, [rax + SL_rval]
pop r13
pop rax
lea rdi, [sz_slfmt]
lea esi, [r13 + 1]
movsd xmm0, [rbp-8]
movsd xmm1, xmm5
movsd xmm2, xmm6
movsd xmm3, xmm7
movsd xmm4, xmm8
mov eax, 5
call printf
.sl_next:
inc r13d
jmp .sl_loop
.sl_done:
; --- cleanup ---
mov rdi, r12
call checkpoint_free
mov rdi, rbx
call lattice_free
lea rdi, [sz_final]
xor eax, eax
call printf
xor eax, eax
.exit:
pop r13
pop r12
pop rbx
leave
ret
II_analog_boot.zip (376.3 KB)
What is 
File Details:
What It Contains:
Why It Matters:
How It’s Used:
Performance Impact:
MISSION ACCOMPLISHED
V3.0 RESULTS SUMMARY
V3.0 MATHEMATICAL ADVANTAGES
VISUALIZATION FEATURES
ACHIEVEMENT HIGHLIGHTS
10,000-Point Spiral: 100% successful generation
COMPARISON: V2.6 vs V3.0
IMPACT ASSESSMENT
