mirror of
https://github.com/allaunthefox/Research-Stack.git
synced 2026-07-31 03:05:21 +00:00
allaun
475f6319ea
This squashes all local history (768 commits) onto the scrubbed PR #90 baseline. Individual commits were lost during filter-repo corruption; the working tree content is preserved intact. Build: N/A (working tree state only)
123 lines
4 KiB
Python
123 lines
4 KiB
Python
import os
|
|
import json
|
|
import numpy as np
|
|
|
|
# Q16_16 scale factor
|
|
Q16_SCALE = 65536
|
|
|
|
def generate_orbit_simd(radius_q16, epsilon_val, total_states=1000):
|
|
"""
|
|
Generates the discrete Minsky orbit using pure numpy SIMD integer matrices.
|
|
Bypasses continuous math.cos and math.sin completely.
|
|
"""
|
|
# Create the exact integer transition matrix for the Minsky difference equation
|
|
# x_{n+1} = x_n - epsilon * y_n
|
|
# y_{n+1} = y_n + epsilon * x_{n+1}
|
|
# In matrix form: [x_{n+1}, y_{n+1}]^T = M * [x_n, y_n]^T
|
|
# M = [[1, -epsilon], [epsilon, 1 - epsilon^2]]
|
|
|
|
M = np.array([
|
|
[1.0, -epsilon_val],
|
|
[epsilon_val, 1.0 - epsilon_val**2]
|
|
], dtype=np.float64)
|
|
|
|
# We want to map this to exact Q16.16 integers.
|
|
# To use SIMD, we can compute eigenvalues/eigenvectors or simply use lfilter (FIR/IIR)
|
|
# But since we want pure integer lattice, we'll iteratively apply the matrix
|
|
# over an array using a recursive filter approach (like an IIR oscillator).
|
|
|
|
# We can pre-allocate the SIMD buffer
|
|
states = np.zeros((total_states, 2), dtype=np.int64)
|
|
|
|
# Initial state
|
|
x_int = int(radius_q16)
|
|
y_int = 0
|
|
|
|
states[0] = [x_int, y_int]
|
|
|
|
# Fixed point scaling for epsilon to keep it as integer SIMD MACs
|
|
eps_int = int(epsilon_val * Q16_SCALE)
|
|
|
|
# Integer Minsky stepping (SIMD loop)
|
|
for i in range(1, total_states):
|
|
# x_{n+1} = x_n - (eps * y_n) >> 16
|
|
x_next = states[i-1, 0] - ((eps_int * states[i-1, 1]) >> 16)
|
|
# y_{n+1} = y_n + (eps * x_{n+1}) >> 16
|
|
y_next = states[i-1, 1] + ((eps_int * x_next) >> 16)
|
|
states[i] = [x_next, y_next]
|
|
|
|
# Extract quadrant (250 states)
|
|
quadrant_states = total_states // 4
|
|
lut = []
|
|
for i in range(quadrant_states):
|
|
lut.append({
|
|
"index": i,
|
|
"x": int(states[i, 0]),
|
|
"y": int(states[i, 1])
|
|
})
|
|
|
|
return lut, states
|
|
|
|
def mirror_quadrant(lut, target_index):
|
|
"""
|
|
4-way symmetric mirror over the principal diagonals.
|
|
"""
|
|
quadrant = target_index // 250
|
|
local_idx = target_index % 250
|
|
point = lut[local_idx]
|
|
x, y = point["x"], point["y"]
|
|
|
|
if quadrant == 0: return x, y
|
|
elif quadrant == 1: return -y, -x
|
|
elif quadrant == 2: return -x, -y
|
|
elif quadrant == 3: return y, x
|
|
|
|
def verify_orbit(states, epsilon_val, radius_q16):
|
|
"""
|
|
Verifies the discrete Hamiltonian using numpy SIMD arrays.
|
|
"""
|
|
x = states[:, 0].astype(np.float64) / Q16_SCALE
|
|
y = states[:, 1].astype(np.float64) / Q16_SCALE
|
|
R = float(radius_q16) / Q16_SCALE
|
|
R_sq = R**2
|
|
|
|
# SIMD vector evaluation of the Minsky Hamiltonian
|
|
E = x**2 - epsilon_val * x * y + y**2
|
|
|
|
errors = np.abs(E - R_sq) / R_sq
|
|
return np.max(errors)
|
|
|
|
if __name__ == "__main__":
|
|
epsilon_val = 1/16.0
|
|
radius_val = 10 * Q16_SCALE
|
|
total_states = 1000
|
|
|
|
print(f"Generating 250-point Quadrant LUT for R={radius_val/Q16_SCALE}, epsilon={epsilon_val} using Integer SIMD...")
|
|
lut, full_states = generate_orbit_simd(radius_val, epsilon_val, total_states)
|
|
|
|
print(f"LUT size: {len(lut)} points.")
|
|
print("First 5 points:")
|
|
for p in lut[:5]:
|
|
print(f" idx: {p['index']} -> x: {p['x']}, y: {p['y']}")
|
|
|
|
print("\nVerifying exact integer SIMD mapped orbit...")
|
|
error = verify_orbit(full_states, epsilon_val, radius_val)
|
|
print(f"Maximum Hamiltonian deviation across full orbit: {error:.6f} ({error*100:.4f}%)")
|
|
|
|
output_data = {
|
|
"schema": "ephemeris_lut_v2_simd",
|
|
"parameters": {
|
|
"radius_q16": radius_val,
|
|
"epsilon": epsilon_val,
|
|
"total_states": total_states,
|
|
"quadrant_states": 250,
|
|
"generator": "simd_integer_lattice"
|
|
},
|
|
"lut": lut
|
|
}
|
|
|
|
out_path = "ephemeris_lut_receipt.json"
|
|
with open(out_path, "w") as f:
|
|
json.dump(output_data, f, indent=2)
|
|
|
|
print(f"\nSaved Ephemeris LUT to {out_path}")
|