checksumsets.py: add support for generating animations.

Signed-off-by: Daira Hopwood <daira@jacaranda.org>
This commit is contained in:
Daira Hopwood 2020-06-23 16:27:06 +01:00
parent 5db9b7a1bc
commit 289e616084
2 changed files with 239 additions and 40 deletions

16
animation.sh Executable file
View File

@ -0,0 +1,16 @@
#!/bin/sh
# pip install bintrees Pillow
# apt-get install ffmpeg ffcvt
set -e
./checksumsets.py --animate
if [ -f animation-p.gif ]; then
ffcvt -f animation-p.gif
mv -f animation-p_.webm animation-p.webm
rm -f animation-p.gif
fi
if [ -f animation-q.gif ]; then
ffcvt -f animation-q.gif
mv -f animation-q_.webm animation-q.webm
rm -f animation-q.gif
fi

View File

@ -1,7 +1,15 @@
#!/usr/bin/python3
# -*- coding: utf-8 -*-
# Dependency: <https://pypi.org/project/bintrees/> (pip install bintrees)
# Dependencies:
# <https://pypi.org/project/bintrees/> (pip install bintrees)
# <https://pypi.org/project/Pillow/> (pip install Pillow), if --animate is used
import sys
from dataclasses import dataclass
from typing import Optional
from collections import deque
from math import log
# From the Halo paper:
@ -57,13 +65,25 @@
BRUTEFORCE_THRESHOLD = 100000
def D(q, zeta, mm):
DEBUG = False
@dataclass
class State:
u: Optional[int]
m: Optional[int]
n: Optional[int]
d: Optional[int]
def D(q, zeta, mm, animator=None):
if DEBUG: print("(q, zeta, mm) =", (q, zeta, mm))
Dcheck = [] if BRUTEFORCE_THRESHOLD == 0 else bruteforce_D(q, zeta, min(mm, BRUTEFORCE_THRESHOLD))
(u, m) = (0, 1) # (u + am) : a ∈ Nat is the current arithmetic progression
n = q # the previous min-distance
d = zeta # the current min-distance
# (u + am) : a ∈ Nat is the current arithmetic progression
# n is the previous min-distance
# d is the current min-distance
cur = State(u=0, m=1, n=q, d=zeta)
old = None
while True:
# Consider values of x where D_{q,ζ_q}(x) decreases, i.e. where
@ -77,36 +97,56 @@ def D(q, zeta, mm):
#
# If we set s = floor(n/d), then D_{q,ζ_q}(x) can decrease at a = s
# and potentially also at a = s+1.
assert (m*zeta) % q in (d, q-d)
s = n // d
x0 = u + s*m
d0 = n % d
if DEBUG: print("(x0, d0, u, m, n, d, s) =", (x0, d0, u, m, n, d, s))
assert (cur.m*zeta) % q in (cur.d, q - cur.d)
s = cur.n // cur.d
x0 = cur.u + s*cur.m
d0 = cur.n % cur.d
if DEBUG: print("\n(x0, d0, cur, s) =", (x0, d0, cur, s))
assert dist(0, x0*zeta, q) == d0
if x0-1 < len(Dcheck): assert Dcheck[x0-1] == d
if x0 > mm: return d
if x0-1 < len(Dcheck): assert Dcheck[x0-1] == cur.d
if x0 > mm:
if animator is not None:
animator.render(q, zeta, old, cur, None, s)
return cur.d
if x0 < len(Dcheck): assert Dcheck[x0] == d0
x1 = u + (s+1)*m
d1 = (s+1)*d - n
x1 = cur.u + (s+1)*cur.m
d1 = (s+1)*cur.d - cur.n
if d1 < d0:
if DEBUG: print("(x1, d1, u, m, n, d, s+1) =", (x1, d1, u, m, n, d, s+1))
if DEBUG: print("(x1, d1, cur, s+1) =", (x1, d1, cur, s+1))
assert dist(0, x1*zeta, q) == d1
if x1-1 < len(Dcheck): assert Dcheck[x1-1] == d0
if x1 > mm: return d0
if x1 > mm:
if animator is not None:
animator.render(q, zeta, old, cur, None, s+1)
return d0
if x1 < len(Dcheck): assert Dcheck[x1] == d1
# TODO: This is painfully non-obvious! Need to draw some diagrams to explain it.
(u, m, n, d) = (x0, x1, d0, d1)
# This is the case where the smaller new distance is past zero.
# The next iteration should consider the region of size d0 starting at x = x0
# (i.e. just before we went past zero) and increasing by x1, i.e. dividing
# that region by intervals of d1.
new = State(u=x0, m=x1, n=d0, d=d1)
else:
(u, m, n, d) = (x1, x0, d-d0, d0)
# This is the case where the smaller new distance is short of zero.
# The next iteration should check the region of size cur.d - d0 starting at x = x1
# (i.e. the wraparound past zero) and increasing by x0, i.e. dividing that
# region by intervals of d0.
new = State(u=x1, m=x0, n=cur.d - d0, d=d0)
assert dist(0, new.u*zeta, q) in (new.n, q - new.n)
#if dist(0, new.u*zeta, q) != new.n: print("hmm")
if animator is not None:
animator.render(q, zeta, old, cur, new, s)
(old, cur) = (cur, new)
def bruteforce_D(q, zeta, mm):
# Can't use sortedcontainers because its data structures are backed by
# lists-of-lists, not trees. We must have O(log n) insert, prev, and succ.
from bintrees import RBTree as sortedset
from collections import deque
resD = deque([zeta])
lastd = zeta
@ -120,7 +160,7 @@ def bruteforce_D(q, zeta, mm):
vs = S.succ_key(v)
d = min(v-vp, vs-v)
resD.append(d)
if DEBUG and d < lastd: print((x, d))
#if DEBUG and d < lastd: print((x, d))
lastd = d
return list(resD)
@ -129,29 +169,172 @@ def dist(x, y, q):
z = (x-y+q) % q
return min(z, q-z)
def check_sumset(name, q, zeta, limit):
def signed_mod(x, q):
r = x % q
return r if r <= q//2 else r-q
class Animator:
fontfile = '/usr/share/texlive/texmf-dist/fonts/truetype/google/roboto/Roboto-Regular.ttf'
frame_duration = 20 # ms
pause_frames = 35
zoom_frames = 45
width = 800 # pixels
height = 400 # pixels
oversample = 3
line_halfwidth = 1 # subpixels
ground_colour = '#ffffff' # white
scale_colour = '#0000a0' # blue
midline_colour = '#c00000' # red
old_colour = '#a0a0a0' # grey
cur_colour = '#000000' # black
new_colour = '#008000' # green
final_colour = '#c00000' # red
def __init__(self, name):
# We don't want to depend on PIL unless an Animator is instantiated.
from PIL import Image, ImageDraw, ImageColor, ImageFont
self.Image = Image
self.ImageDraw = ImageDraw
self.ImageColor = ImageColor
self.font = ImageFont.truetype(self.fontfile, 20*self.oversample, index=0, encoding='unic')
self.font_super = ImageFont.truetype(self.fontfile, 12*self.oversample, index=0, encoding='unic')
self.images = deque()
self.name = name
def render(self, q, zeta, old, cur, new, s):
n = min(cur.n, q//2)
for aa in range(1, s+1):
self.render_zoomed(q, zeta, old, cur, None, aa, n)
if new is None:
self.render_zoomed(q, zeta, old, cur, new, s+1, n, final=True)
return
self.render_zoomed( q, zeta, old, cur, new, s+1, n, frames=self.pause_frames)
step = (1.0*n/new.n - 1.0)/self.zoom_frames
for zoom in range(1, self.zoom_frames):
n_scale = int(n/(1.0 + zoom*step))
self.render_zoomed(q, zeta, old, cur, new, s+1, n_scale)
self.render_zoomed( q, zeta, old, cur, new, s+1, new.n, frames=self.pause_frames)
def render_zoomed(self, q, zeta, old, cur, new, aa, n_scale, frames=1, final=False):
px = self.oversample
lx = self.line_halfwidth
(w, h) = (self.width * px, self.height * px)
scale = (w/2)/n_scale
xmid = w//2
ymid = (40*px + h)//2
image = self.Image.new('RGB', (w, h), color=self.ground_colour)
image.convert('P')
draw = self.ImageDraw.Draw(image)
bits = int(log(n_scale, 2))
for tick in range(bits-3, bits+1):
xoff = int(scale*(1<<tick))
draw.text((xmid-xoff-21*px, 7*px), '2', self.scale_colour, font=self.font)
draw.text((xmid-xoff, px), str(tick), self.scale_colour, font=self.font_super)
draw.rectangle((xmid-xoff-lx, 30*px, xmid-xoff+lx, 40*px), fill=self.scale_colour)
draw.text((xmid+xoff-21*px, 7*px), '+2', self.scale_colour, font=self.font)
draw.text((xmid+xoff, px), str(tick), self.scale_colour, font=self.font_super)
draw.rectangle((xmid+xoff-lx, 30*px, xmid+xoff+lx, 40*px), fill=self.scale_colour)
draw.rectangle((xmid-lx, 0, xmid+lx, h), fill=self.midline_colour)
draw.rectangle((0, 40*px-lx, w, 40*px+lx), fill=self.scale_colour)
if old is not None:
old_aa = (old.n // old.d)+1
for a in range(old_aa+1):
x = signed_mod(zeta*(old.u + a*old.m), q)
xpos = w//2 + int(scale*x)
draw.rectangle((xpos-lx, 40*px, xpos+lx, h), fill=self.old_colour)
for a in range(aa+1):
x = signed_mod(zeta*(cur.u + a*cur.m), q)
xpos = w//2 + int(scale*x)
draw.rectangle((xpos-lx, 40*px, xpos+lx, h), fill=self.cur_colour)
if new is not None:
x = signed_mod(zeta*new.u, q)
xpos = w//2 + int(scale*x)
draw.rectangle((xpos, ymid-lx, xmid, ymid+lx), fill=self.new_colour)
draw.rectangle((xpos-lx, 40*px, xpos+lx, h), fill=self.new_colour)
if final:
x = signed_mod(zeta*(cur.u + aa*cur.m), q)
xpos = w//2 + int(scale*x)
draw.rectangle((xpos, ymid-lx, xmid, ymid+lx), fill=self.final_colour)
draw.rectangle((xpos-lx, 40*px, xpos+lx, h), fill=self.final_colour)
image = image.resize((self.width, self.height), self.Image.ANTIALIAS)
for f in range(frames):
self.images.append(image)
sys.stdout.write('.')
sys.stdout.flush()
def save(self):
filename = 'animation-%s.gif' % (self.name,)
print("Saving %s..." % (filename,))
first, *rest = list(self.images)
# Save as animated GIF. We can convert to a more space-efficient format separately.
first.save(fp=filename, format='GIF', append_images=rest,
save_all=True, duration=self.frame_duration, loop=1)
del self.images
def check_sumset(name, q, zeta, limit, animator=None):
print("===== %s =====" % (name,))
Dq = D(q, zeta, limit-1)
print("D_%s = %s" % (name, Dq), end=' ')
Dq = D(q, zeta, limit-1, animator)
print("\nD_%s = %s" % (name, Dq))
print(" %s" % ('' if Dq >= limit else '<'), limit)
if animator is not None:
animator.save()
assert Dq >= limit
print(">=", limit)
halflambda = 64
limit = 3<<halflambda
def main():
args = sys.argv[1:]
if "--help" in args:
print("Usage: checksumsets.py [--animate]")
return
# Tweedledum and Tweedledee
p = (1<<254) + 4707489545178046908921067385359695873
q = (1<<254) + 4707489544292117082687961190295928833
zeta_p = 9513155655832138286304767221959569637168364952810827555227185832555034233288
zeta_q = 24775483399512474214391554062650059912556682109176536098332128018848638018813
halflambda = 64
limit = 3<<halflambda
# Tests
DEBUG = False
assert(D(65537, 6123, 10000) == 3)
assert(D(1299721, 538936, 10000) == 41)
assert(D(179424691, 134938504, 100000) == 121)
# Tweedledum and Tweedledee
p = (1<<254) + 4707489545178046908921067385359695873
q = (1<<254) + 4707489544292117082687961190295928833
zeta_p = 9513155655832138286304767221959569637168364952810827555227185832555034233288
zeta_q = 24775483399512474214391554062650059912556682109176536098332128018848638018813
DEBUG = False
check_sumset("p", p, zeta_p, limit)
check_sumset("q", q, zeta_q, limit)
# Tests
global DEBUG
DEBUG = False
assert(D(65537, 6123, 10000, None) == 3)
assert(D(1299721, 538936, 10000, None) == 41)
assert(D(179424691, 134938504, 100000, None) == 121)
p_params = ("p", p, zeta_p, limit)
q_params = ("q", q, zeta_q, limit)
DEBUG = False
for (name, prime, zeta, limit) in (p_params, q_params):
animator = None
if "--animate" in args:
animator = Animator(name)
check_sumset(name, prime, zeta, limit, animator=animator)
main()