Reproducible source and supporting text for the associated calculation.
View complete source excerpt
#!/usr/bin/env python3
"""Exact exterior charge moments and parity controls; no particle simulation."""
from collections import Counter
import argparse
from datetime import datetime, timezone
from fractions import Fraction as F
from hashlib import sha256
import json
from math import comb, factorial
from pathlib import Path
parser=argparse.ArgumentParser(description=__doc__)
parser.add_argument('--output-dir',type=Path,help='A new directory for paired code and results.')
args=parser.parse_args()
script_bytes=Path(__file__).read_bytes()
checks = 0
def require(condition):
global checks
assert condition
checks += 1
def spectrum(m, parity=None):
values = Counter()
for grade in range(m+1):
if parity is None or grade % 2 == parity:
values[F(grade,m)] += comb(m,grade)
values[F(grade,m)-1] += comb(m,grade)
return values
def moment(values, power):
return sum((mult * charge**power for charge,mult in values.items()), F(0))
def exact(x):
return str(x.numerator) + '/' + str(x.denominator)
def pfaffian(matrix):
if not matrix:
return F(1)
total = F(0)
for j in range(1,len(matrix)):
rest = [i for i in range(len(matrix)) if i not in (0,j)]
minor = [[matrix[a][b] for b in rest] for a in rest]
total += (-1)**(j+1)*matrix[0][j]*pfaffian(minor)
return total
def multiply(a,b):
return [[sum((a[i][k]*b[k][j] for k in range(len(b))),F(0))
for j in range(len(b[0]))] for i in range(len(a))]
def transpose(a):
return list(map(list,zip(*a)))
def rank(matrix):
a = [row[:] for row in matrix]
r = 0
for c in range(len(a[0])):
pivot = next((i for i in range(r,len(a)) if a[i][c]),None)
if pivot is None:
continue
a[r],a[pivot]=a[pivot],a[r]
scale=a[r][c]
a[r]=[v/scale for v in a[r]]
for i in range(len(a)):
if i != r:
scale=a[i][c]
a[i]=[v-scale*w for v,w in zip(a[i],a[r])]
r += 1
if r == len(a):
break
return r
def frame_tensor(x,y):
a,aa=x[0],x[1:]
b,bb=y[0],y[1:]
v=[aa[1]*bb[2]-aa[2]*bb[1],aa[2]*bb[0]-aa[0]*bb[2],aa[0]*bb[1]-aa[1]*bb[0]]
xi=[a*z-b*w for w,z in zip(aa,bb)]
return v,xi
def aggregate_rotation(xs,ys):
omega=[[F(0) for _ in range(6)] for _ in range(6)]
q=F(0)
for x in xs:
for y in ys:
v,xi=frame_tensor(x,y)
require(sum(a*b for a,b in zip(v,xi)) == 0)
q += sum(a*a for a in v+xi)
u=[a+b for a,b in zip(v,xi)]
w=[a-b for a,b in zip(v,xi)]
for i in range(3):
for j in range(3):
omega[2*i][2*j+1] -= u[i]*w[j]
omega[2*i+1][2*j] += w[i]*u[j]
n0=sum(a*a for x in xs for a in x)
n1=sum(a*a for y in ys for a in y)
mixed=[[sum(a*b for a,b in zip(x,y)) for y in ys] for x in xs]
require(q == n0*n1-sum(a*a for row in mixed for a in row))
return omega,q,F(n0*n1)
def gamma_action(state,axis):
plane=axis//2
out={}
for mask,coefficient in state.items():
sign=(-1)**((mask & ((1<<plane)-1)).bit_count())
if axis%2:
sign *= 1j*(-1)**((mask>>plane)&1)
target=mask^(1<<plane)
out[target]=out.get(target,0)+sign*coefficient
return out
def spinor_rotation(m,state):
def gaussian_norm(c):
c=complex(c)
require(c.real==int(c.real) and c.imag==int(c.imag))
return int(c.real)**2+int(c.imag)**2
norm=sum(gaussian_norm(c) for c in state.values())
result=[[F(0) for _ in range(2*m)] for _ in range(2*m)]
for i in range(2*m):
for j in range(i+1,2*m):
acted=gamma_action(gamma_action(state,j),i)
scalar=1j*sum(complex(c).conjugate()*acted.get(mask,0) for mask,c in state.items())
require(scalar.imag == 0 and scalar.real == int(scalar.real))
result[i][j]=F(int(scalar.real),norm)
result[j][i]=-result[i][j]
number=sum(F(mask.bit_count()*gaussian_norm(c),norm) for mask,c in state.items())
return result,number
rows = []
for m in range(1,17):
even,odd,full = spectrum(m,0),spectrum(m,1),spectrum(m)
require(sum(even.values()) == 2**m)
require(sum(full.values()) == 2**(m+1))
require(even+odd == full)
require(all(full[q] == full[-q] for q in full))
if m % 2 == 0:
require(all(even[q] == even[-q] for q in even))
for power in range(0,2*m+4):
e,o,f = moment(even,power),moment(odd,power),moment(full,power)
require(e+o == f)
if power % 2:
require(f == 0)
if m%2 == 0 or power < m:
require(e == 0)
if m%2 and power == m:
require(e == -F(factorial(m),m**m))
# Independent finite-character identity in unnormalized integer weights.
for x in [F(1,2),F(2,3),F(1),F(3,2),F(2)]:
ge = sum(mult*x**int(q*m) for q,mult in even.items())
gi = sum(mult*x**int(-q*m) for q,mult in even.items())
closed = (1+x**(-m))*((1+x)**m+(1-x)**m)/2
supertrace = sum((even[q]-odd[q])*x**int(q*m) for q in full)
require(ge == closed)
require(ge-gi == ((1+x**(-m))*(1-x)**m if m%2 else 0))
require(supertrace == (1+x**(-m))*(1-x)**m)
if m >= 3:
require(moment(even,2)/sum(even.values()) == F(m+1,4*m))
rows.append({'complex_normal_rank':m,'real_normal_directions':2*m,
'orthogonal_complex_structure_orbit_real_dimension':m*(m-1),
'ambient_depth_if_base_four':4+2*m,'even_carrier_rank':sum(even.values()),
'first_nonzero_odd_moment':m if m%2 else None,
'first_odd_moment_trace':exact(-F(factorial(m),m**m)) if m%2 else None,
'normalized_cubic_trace':exact(moment(even,3)),
'full_cubic_trace':exact(moment(full,3)),
'spectrum':[{'q':exact(q),'multiplicity':even[q]} for q in sorted(even)]})
require(moment(spectrum(3,0),3) == -F(2,9))
require(moment(spectrum(3,1),3) == F(2,9))
require(moment(spectrum(3),3) == 0)
require([m for m in range(2,17) if moment(spectrum(m,0),3)] == [3])
require([m for m in range(1,17) if m*(m-1)==2*m] == [3])
# Multivariate parity character with independent plane weights x_i.
# Its lowest polynomial degree is m, with coefficient 2*(-1)^m*x_1...x_m.
for m in range(1,9):
plane_weights = [F(i+1,i+2) for i in range(m)]
degree = m
super_moment = F(0)
for mask in range(1<<m):
grade = mask.bit_count()
w = sum((plane_weights[i] for i in range(m) if mask & (1<<i)),F(0))
determinant_weight = sum(plane_weights,F(0))
super_moment += (-1)**grade * (w**degree + (w-determinant_weight)**degree)
product = F(1)
for w in plane_weights:
product *= w
require(super_moment == 2*(-1)**m*factorial(m)*product)
if m <= 4:
size = 2*m
a = [[F(0) for _ in range(size)] for _ in range(size)]
for i,w in enumerate(plane_weights):
a[2*i][2*i+1] = -w
a[2*i+1][2*i] = w
require(super_moment == 2*factorial(m)*pfaffian(a))
rotation = [[F(int(i==j)) for j in range(size)] for i in range(size)]
if m > 1:
for i,j in [(0,2),(1,3)]:
elementary = [[F(int(x==y)) for y in range(size)] for x in range(size)]
elementary[i][i]=elementary[j][j]=F(3,5)
elementary[i][j]=-F(4,5)
elementary[j][i]=F(4,5)
rotation = multiply(elementary,rotation)
require(multiply(rotation,transpose(rotation)) ==
[[F(int(i==j)) for j in range(size)] for i in range(size)])
rotated = multiply(multiply(rotation,a),transpose(rotation))
require(pfaffian(rotated) == pfaffian(a))
reflection = [[F(int(i==j)) for j in range(size)] for i in range(size)]
reflection[0][0]=-1
reflected=multiply(multiply(reflection,a),transpose(reflection))
require(pfaffian(reflected) == -pfaffian(a))
e=[[F(int(i==j)) for i in range(4)] for j in range(4)]
one=aggregate_rotation([e[0]],[e[1]])
three=aggregate_rotation([e[0],e[1]],[e[0],e[2]])
cancel=aggregate_rotation([e[0],e[1]],[e[2],e[3]])
partial=aggregate_rotation([e[0],e[1]],
[[a+b for a,b in zip(e[0],e[2])],[a+b for a,b in zip(e[1],e[3])]])
require(rank(one[0]) == 2 and pfaffian(one[0]) == 0)
require(rank(three[0]) == 6 and pfaffian(three[0]) == -1)
require(multiply(three[0],three[0]) == [[-F(int(i==j)) for j in range(6)] for i in range(6)])
require(rank(cancel[0]) == 0 and cancel[1] == 4)
require(rank(partial[0]) == 2 and pfaffian(partial[0]) == 0)
require(three[1]/three[2] == partial[1]/partial[2] == F(3,4))
rotation_cases=[]
for name,case in [('single-plane',one),('three-plane',three),('circulation-cancels',cancel),('same-normalized-frame-power',partial)]:
omega,q,n=case
rotation_cases.append({'name':name,'rotation_rank':rank(omega),'pfaffian':exact(pfaffian(omega)),
'frame_power':exact(q),'intensity_product':exact(n),
'normalized_frame_power':exact(q/n),
'volume_efficiency':exact(27*abs(pfaffian(omega))/q**3),
'rotation_matrix':[[exact(a) for a in row] for row in omega]})
require(0 <= 27*abs(pfaffian(omega)) <= q**3)
require(27*abs(pfaffian(three[0]))/three[1]**3 == 1)
spinor_cases=[]
for m,state,name in [(2,{0:1},'rank-two-vacuum'),(2,{3:1},'rank-two-top-singlet'),
(3,{0:1},'rank-three-vacuum'),(3,{3:1},'rank-three-grade-two-basis'),
(3,{3:1,5:2,6:3j},'rank-three-arbitrary-grade-two'),
(3,{0:1,3:1},'rank-three-mixed-charge'),
(4,{3:1},'rank-four-pure-grade-two'),
(4,{3:1,12:1},'rank-four-nonpure-same-charge')]:
j,number=spinor_rotation(m,state)
j0,_=spinor_rotation(m,{0:1})
product=multiply(j,j0)
alignment=-sum(product[i][i] for i in range(2*m))/F(2*m)
require(alignment == 1-2*number/m)
require((1-alignment)/2 == number/m)
is_complex=multiply(j,j)==[[-F(int(i==k)) for k in range(2*m)] for i in range(2*m)]
if m<=3 or len(state)==1:
require(is_complex)
else:
require(not is_complex and rank(j)==0)
spinor_cases.append({'name':name,'complex_normal_rank':m,'grade_expectation':exact(number),
'unshifted_charge_expectation':exact(number/m),'rotation_alignment':exact(alignment),
'rotation_is_orthogonal_complex_structure':is_complex,'rotation_rank':rank(j)})
# The same rank-four counterexample embeds at every higher rank tested.
for m in range(5,9):
j,number=spinor_rotation(m,{3:1,12:1})
require(number==2 and rank(j)==2*(m-4))
require(multiply(j,j)!=[[-F(int(i==k)) for k in range(2*m)] for i in range(2*m)])
result = {'status':'passed','checks':checks,'comparison_ranks':[1,16],
'scope':'Exact pointwise characters, not selected particles or physical anomaly coefficients',
'rows':rows,'aggregate_rotation_cases':rotation_cases,
'spinor_rotation_cases':spinor_cases,
'sources':{'experiments/normal-charge-character-checks.py':sha256(script_bytes).hexdigest()}}
output_dir=args.output_dir or Path(__file__).parent/('run-'+datetime.now(timezone.utc).strftime('%Y%m%dT%H%M%S.%fZ'))
output_dir.mkdir(parents=True,exist_ok=False)
(output_dir/'executed-source.py').write_bytes(script_bytes)
output = output_dir/'normal-charge-character-checks.json'
output.write_text(json.dumps(result,indent=2)+'\n')
print(json.dumps({'status':result['status'],'checks':checks,'output':str(output),
'rank_three_cubic_trace':'-2/9','full_completion_cubic_trace':'0'},indent=2))