KOKO\Mimi
245 lines
10 KiB
Python
245 lines
10 KiB
Python
"""Generate test-only certified Gauss brackets from exact rational arithmetic.
|
|
|
|
Newton/Decimal produces candidates only. Exact Fraction signs prove n distinct
|
|
Legendre roots; rational interval recurrence encloses the associated weights.
|
|
No runtime solver dependency, reference modification, or libm ULP premise.
|
|
"""
|
|
from decimal import Decimal, localcontext
|
|
from fractions import Fraction as F
|
|
from pathlib import Path
|
|
import math
|
|
import re
|
|
import json
|
|
|
|
|
|
def add(a, b):
|
|
return a[0] + b[0], a[1] + b[1]
|
|
|
|
|
|
def neg(a):
|
|
return -a[1], -a[0]
|
|
|
|
|
|
def mul(a, b):
|
|
v = [x * y for x in a for y in b]
|
|
return min(v), max(v)
|
|
|
|
|
|
def div(a, b):
|
|
assert b[0] > 0 or b[1] < 0
|
|
return mul(a, (1 / b[1], 1 / b[0]))
|
|
|
|
|
|
def scale(a, b):
|
|
return mul(a, (F(b), F(b)))
|
|
|
|
|
|
def legendre(n, x):
|
|
p, q = (F(1), F(1)), x
|
|
for k in range(2, n + 1):
|
|
p, q = q, scale(add(scale(mul(x, q), 2 * k - 1), neg(scale(p, k - 1))), F(1, k))
|
|
derivative = div(scale(add(mul(x, q), neg(p)), n), add(mul(x, x), (-F(1), -F(1))))
|
|
return q, derivative
|
|
|
|
|
|
def bracket(value):
|
|
return math.nextafter(float(value[0]), -math.inf), math.nextafter(float(value[1]), math.inf)
|
|
|
|
|
|
def rule(n):
|
|
result = []
|
|
for i in range(n):
|
|
with localcontext() as context:
|
|
context.prec = 90
|
|
z = Decimal(str(math.cos(math.pi * (i + .75) / (n + .5))))
|
|
for _ in range(30):
|
|
p, q = Decimal(1), z
|
|
for k in range(2, n + 1):
|
|
p, q = q, ((2 * k - 1) * z * q - (k - 1) * p) / k
|
|
d = n * (z * q - p) / (z * z - 1)
|
|
z -= q / d
|
|
lo, hi = bracket((F(z), F(z)))
|
|
pl = legendre(n, (F(lo), F(lo)))[0][0]
|
|
ph = legendre(n, (F(hi), F(hi)))[0][0]
|
|
assert pl * ph < 0
|
|
x = F(lo), F(hi)
|
|
derivative = legendre(n, x)[1]
|
|
w = div((F(2), F(2)), mul(add((F(1), F(1)), neg(mul(x, x))), mul(derivative, derivative)))
|
|
result.append((lo, hi, *bracket(w)))
|
|
result.sort()
|
|
assert all(a[1] < b[0] for a, b in zip(result, result[1:]))
|
|
assert sum(F(r[2]) for r in result) < 2 < sum(F(r[3]) for r in result)
|
|
return result
|
|
|
|
|
|
def slab_fixtures():
|
|
# Approved Appendix A supplies the independently implemented Cholesky /
|
|
# Jacobi spectral candidates; Appendix E supplies certified transcendentals.
|
|
document = Path(__file__).resolve().parents[3] / "docs/shell-transient-heat-transfer/numerical-oracle-scratch.md"
|
|
blocks = re.findall(r"```python\n(.*?)```", document.read_text(encoding="utf-8-sig"), re.S)
|
|
spectral = {}
|
|
exec(blocks[0].split('out={"lambda"')[0], spectral)
|
|
interval_scope = {}
|
|
exec(blocks[-1].split("\nfor t in [.001")[0], interval_scope)
|
|
iv, ia, im, isc = (interval_scope[k] for k in ("interval", "add", "mul", "scale"))
|
|
ex = interval_scope["exneg"]
|
|
n = 4
|
|
|
|
def exact_inverse(matrix):
|
|
a = [list(row) + [F(i == j) for j in range(n)] for i, row in enumerate(matrix)]
|
|
for k in range(n):
|
|
pivot = next(i for i in range(k, n) if a[i][k])
|
|
a[k], a[pivot] = a[pivot], a[k]
|
|
divisor = a[k][k]
|
|
a[k] = [x / divisor for x in a[k]]
|
|
for i in range(n):
|
|
if i != k:
|
|
factor = a[i][k]
|
|
a[i] = [x - factor * y for x, y in zip(a[i], a[k])]
|
|
return [row[n:] for row in a]
|
|
|
|
def exact_mm(a, b):
|
|
return [[sum((x*y for x, y in zip(row, col)), F(0)) for col in zip(*b)] for row in a]
|
|
|
|
def matrix_multiply(a, b):
|
|
return [[sum_intervals([im(x, y) for x, y in zip(row, col)]) for col in zip(*b)] for row in a]
|
|
|
|
def sum_intervals(values):
|
|
total = iv(0)
|
|
for v in values:
|
|
total = ia(total, v)
|
|
return total
|
|
|
|
def apply(a, v):
|
|
return [sum_intervals([im(x, y) for x, y in zip(row, v)]) for row in a]
|
|
|
|
def power(a, count):
|
|
result = [[iv(i == j) for j in range(n)] for i in range(n)]
|
|
while count:
|
|
if count & 1:
|
|
result = matrix_multiply(result, a)
|
|
count >>= 1
|
|
if count:
|
|
a = matrix_multiply(a, a)
|
|
return result
|
|
|
|
c = [[F(0) for _ in range(5)] for _ in range(5)]
|
|
k = [[F(0) for _ in range(5)] for _ in range(5)]
|
|
for offset in (0, 2):
|
|
for i in range(3):
|
|
for j in range(3):
|
|
c[offset+i][offset+j] += F([[4, 2, -1], [2, 16, 2], [-1, 2, 4]][i][j], 60)
|
|
k[offset+i][offset+j] += F(2, 3)*[[7, -8, 1], [-8, 16, -8], [1, -8, 7]][i][j]
|
|
cf, kf = [r[1:] for r in c[1:]], [r[1:] for r in k[1:]]
|
|
inverse = exact_inverse(cf)
|
|
d = exact_mm(inverse, kf)
|
|
y = [-1-sum(inverse[i][j]*c[j+1][0] for j in range(n)) for i in range(n)]
|
|
|
|
def matrix_exponential(t):
|
|
# An independent enclosure of exp(-Cff^-1 Kff t) certifies all spectral
|
|
# eigenpair/coefficient errors together, rather than assuming libm ULPs.
|
|
# For ||B||inf<=1/2, the 64-term exponential remainder is bounded by
|
|
# 2*(1/2)^65/65!. Entrywise intervals preserve subsequent squarings.
|
|
magnitude = max(sum(abs(x*t) for x in row) for row in d)
|
|
shift = 0
|
|
while magnitude > F(1, 2):
|
|
magnitude /= 2
|
|
shift += 1
|
|
a = [[iv(-x*t/F(2)**shift) for x in row] for row in d]
|
|
term = [[iv(i == j) for j in range(n)] for i in range(n)]
|
|
result = [row[:] for row in term]
|
|
for order in range(1, 65):
|
|
term = [[isc(x, F(1, order)) for x in row] for row in matrix_multiply(term, a)]
|
|
result = [[ia(x, z) for x, z in zip(row, terms)] for row, terms in zip(result, term)]
|
|
remainder = 2*F(1, 2)**65/math.factorial(65)
|
|
result = [[ia(x, iv(-remainder, remainder)) for x in row] for row in result]
|
|
for _ in range(shift):
|
|
result = matrix_multiply(result, result)
|
|
return result
|
|
|
|
def outputs(error):
|
|
temperature = [iv(1)] + [ia(iv(1), z) for z in error]
|
|
derivative = [iv(0)] + apply([[iv(-x) for x in row] for row in d], error)
|
|
flux = []
|
|
for offset, s in ((0, -1), (0, 0), (0, 1), (2, -1), (2, 0), (2, 1)):
|
|
deriv = [F(4)*F(s-F(1, 2)), F(-8*s), F(4)*F(s+F(1, 2))]
|
|
flux.append(sum_intervals([isc(temperature[offset+i], -deriv[i]) for i in range(3)]))
|
|
energy = sum_intervals([isc(temperature[j], sum(c[i][j] for i in range(5))) for j in range(5)])
|
|
reaction = sum_intervals([isc(temperature[j], k[0][j]) for j in range(5)] +
|
|
[isc(derivative[j], c[0][j]) for j in range(5)])
|
|
return temperature + flux + [energy, reaction]
|
|
|
|
def spectral_candidate(t, count):
|
|
amplitudes = []
|
|
for eigenvalue, coefficient in zip(spectral["lam"], spectral["coef"]):
|
|
decay = ex(iv(F(eigenvalue)*t)) if count is None else power_scalar(
|
|
interval_scope["div"](iv(1), iv(1+F(eigenvalue)*t/count)), count)
|
|
amplitudes.append(isc(decay, F(coefficient)))
|
|
return outputs(apply([[iv(F(x)) for x in row] for row in spectral["W"]], amplitudes))
|
|
|
|
def power_scalar(x, count):
|
|
result = iv(1)
|
|
while count:
|
|
if count & 1:
|
|
result = im(result, x)
|
|
count >>= 1
|
|
if count:
|
|
x = im(x, x)
|
|
return result
|
|
|
|
table = [[.141, 10.7, .0304, 10.9], [.0153, 1.28, .00468, .824],
|
|
[.00598, .253, .00130, .00632], [.00139, .136, .000666, .00953],
|
|
[.000142, .0284, .0000945, .000526], [.0000414, .00816, .0000159, .0000479]]
|
|
fixtures, metrics = [], []
|
|
for time, space in zip((.001, .01, .05, .1, .5, 1.), table):
|
|
t = F(time)
|
|
temperatures, flux, energy, _, _ = interval_scope["series_interval"](time)
|
|
continuum = temperatures + [flux[i] for i in (0, 1, 2, 2, 3, 4)] + [energy, flux[0]]
|
|
solutions = []
|
|
spectral_enclosure = 0.
|
|
for count in (None, 4096, 8192):
|
|
if count is None:
|
|
operator = matrix_exponential(t)
|
|
else:
|
|
a = [[cf[i][j]+t/count*kf[i][j] for j in range(n)] for i in range(n)]
|
|
operator = power([[iv(x) for x in row] for row in exact_mm(exact_inverse(a), cf)], count)
|
|
exact = outputs(apply(operator, [iv(x) for x in y]))
|
|
candidate = spectral_candidate(t, count)
|
|
enclosed = [(min(a[0], b[0]), max(a[1], b[1])) for a, b in zip(exact, candidate)]
|
|
spectral_enclosure = max(spectral_enclosure, max(float(hi-lo) for lo, hi in enclosed))
|
|
solutions.append(enclosed)
|
|
fixtures.append((time, space, continuum, *solutions))
|
|
metrics.append({"time": time, "spectral_and_matrix_function_enclosure": spectral_enclosure})
|
|
print(json.dumps({"eigenvalues": spectral["lam"], "slab_oracle_metrics": metrics}, indent=2))
|
|
return fixtures
|
|
|
|
|
|
if __name__ == "__main__":
|
|
output = ["// Generated by thermal_oracle_generator.py; exact root/weight certificates.",
|
|
"#ifndef FESA_TESTS_THERMAL_GAUSS_ORACLE_H_", "#define FESA_TESTS_THERMAL_GAUSS_ORACLE_H_",
|
|
"#include <array>", "namespace thermal_oracle {",
|
|
"struct GaussPoint { double x_lo,x_hi,w_lo,w_hi; };"]
|
|
for n in (12, 16):
|
|
rows = rule(n)
|
|
output.append(f"inline constexpr std::array<GaussPoint,{n}> kGauss{n}{{{{")
|
|
output.extend(" {" + ",".join(x.hex() for x in row) + "}," for row in rows)
|
|
output.append("}};")
|
|
output.extend(["}", "#endif", ""])
|
|
destination = Path(__file__).with_name("thermal_gauss_oracle.h")
|
|
destination.write_text("\n".join(output), encoding="utf-8")
|
|
print("Certified Legendre roots/weights: 12 + 16 distinct roots;", destination.name)
|
|
fixtures = slab_fixtures()
|
|
output = ["// Generated by thermal_oracle_generator.py. Outward binary64 brackets.",
|
|
"#ifndef FESA_TESTS_THERMAL_SLAB_ORACLE_H_", "#define FESA_TESTS_THERMAL_SLAB_ORACLE_H_",
|
|
"#include <array>", "namespace thermal_oracle {",
|
|
"struct Bracket { double lo,hi; };",
|
|
"struct SlabCase { double time; std::array<double,4> space; std::array<Bracket,13> fourier,continuous,be4096,be8192; };",
|
|
"inline constexpr std::array<SlabCase,6> kSlab{{"]
|
|
for time, space, *values in fixtures:
|
|
output.append("{" + time.hex() + ",{{" + ",".join(x.hex() for x in space) + "}},")
|
|
for rows in values:
|
|
output.append("{{" + ",".join("{" + ",".join(x.hex() for x in bracket(row)) + "}" for row in rows) + "}},")
|
|
output.append("},")
|
|
output.extend(["}};", "}", "#endif", ""])
|
|
Path(__file__).with_name("thermal_slab_oracle.h").write_text("\n".join(output), encoding="utf-8")
|