Failed attempt to thoroughly analyze the Stalk Market with Z3

stalk.py

from z3 import *
from ctypes import c_float
import random
import math
import time

def roll(a, b, c, d):
    n = a ^ ((a << 11) & 0xffffffff)
    return n ^ (n >> 8) ^ d ^ (d >> 19)

def rollmany(a, b, c, d, steps):
    for i in range(steps):
        a, b, c, d = b, c, d, roll(a, b, c, d)
    return a, b, c, d

def rollvar(a, b, c, d):
    n = a ^ (a << 11)
    return n ^ LShR(n, 8) ^ d ^ LShR(d, 19)

PAT_RANDOM = 0
PAT_BIGSPIKE = 1
PAT_DECREASING = 2
PAT_SMALLSPIKE = 3

# buyprice: int
# sellprices: a list of int (12x)
# returns a list of steps: (lo, hi) or None
def mksteps(buyprice, sellprices, prevpat, pat):
    def randbool(actual):
        if actual: return 0x80000000, 0x100000000
        else: return 0, 0x80000000

    def randint(min, max, actuallo, actualhi):
        # ((next() * (max - min + 1)) >> 32) + min = [actuallo, actualhi)
        # (next() * (max - min + 1)) >> 32 = [actuallo, actualhi) - min
        # next() * (max - min + 1) = [(actuallo - min) << 32, (actualhi - min) << 32)
        # next() * s = [lo, hi)
        # next() = [ceil(lo / s), ceil(hi / s))
        s = max - min + 1
        lo = (actuallo - min) << 32
        hi = (actualhi - min) << 32
        return (lo + s - 1) // s, (hi + s - 1) // s

    def randfloat(a, b, actuallo, actualhi):
        # assuming a < b: (which is not always the case)
        # a + ((next() >> 9) * 2^-23) * (b - a) = [actuallo, actualhi)
        # (next() >> 9) * 2^-23 = [roundup(actuallo / (b - a)) - a, roundup(actualhi / (b - a)) - a)]
        # XXX not accurate
        actuallo = max(actuallo, min(a, a + (b - a)))
        actualhi = min(actualhi, max(a, a + (b - a)))
        lo = actuallo / (b - a) - a
        hi = actualhi / (b - a) - a
        if lo > hi: lo, hi = hi, lo
        return max(int(lo * 0x800000) << 9, 0), min((int(hi * 0x800000) + 1) << 9, 0x100000000)

    def intceilrate(sellprice, buyprice):
        # int(rate * buyprice + 0.99999f) = sellprice
        # rate * buyprice = [sellprice - 0.99999f, sellprice + (1 - 0.99999f))
        # rate = [roundup((sellprice - 0.99999f) / buyprice),
        #         roundup((sellprice + (1 - 0.99999f)) / buyprice))
        # XXX not accurate, should account for binary32 rounding
        return (sellprice - 0.99999) / buyprice, (sellprice + 0.00001) / buyprice

    steps = []
    steps.append(randint(90, 110, buyprice, buyprice + 1))
    table = {
        PAT_RANDOM: (0, 20, 50, 65, 100),
        PAT_BIGSPIKE: (0, 50, 55, 75, 100),
        PAT_DECREASING: (0, 25, 70, 75, 100),
        PAT_SMALLSPIKE: (0, 45, 70, 85, 100),
    }[prevpat]
    steps.append(randint(0, 99, table[pat], table[pat + 1]))

    if pat == PAT_DECREASING:
        ratelo, ratehi = intceilrate(sellprices[0], buyprice)
        steps.append(randfloat(0, 0.05, 0.9 - ratehi, 0.9 - ratelo))
        for sellprice in sellprices[1:]:
            oldratelo = ratelo
            oldratehi = ratehi
            ratelo, ratehi = intceilrate(sellprice, buyprice)
            steps.append(randfloat(0, 0.02, oldratelo - ratehi - 0.03, oldratehi - ratelo - 0.03))
    else:
        assert False, 'not yet implemented'

    return steps

MAX_SOLUTIONS = 4

class Sequence:
    def __init__(self):
        self.seq = BitVecs('a b c d', 32)
        self.solver = SolverFor('QF_BV')
        self.nsteps = 0
        self.ncumulbits = 0
        self.lastmodels = None # (initial state, last state)
        self.starttime = time.time()

    # hi is exclusive
    def step(self, lo, hi):
        next = rollvar(*self.seq[-4:])
        nextvar = BitVec(f's{len(self.seq)}', 32)
        self.seq.append(nextvar)
        self.solver.add(next == nextvar)
        self.solver.add(ULE(lo, nextvar), ULE(nextvar, hi - 1)) # hi can overflow!
        self.nsteps += 1
        self.ncumulbits += 32 - math.log(hi - lo, 2)

        # keep models from the last step if they satisfy new constraint
        models = []
        for initial, (a, b, c, d) in self.lastmodels or []:
            a, b, c, d = b, c, d, roll(a, b, c, d)
            if lo <= d < hi:
                models.append((initial, (a, b, c, d)))
        nretained = len(models)

        if nretained < MAX_SOLUTIONS and (not self.lastmodels or len(self.lastmodels) >= MAX_SOLUTIONS):
            # try to find more models unless we know there won't be no more
            self.solver.push()
            vars = self.seq[:4]
            for initial, _ in models:
                self.solver.add(Or(*[i != val for i, val in zip(vars, initial)]))
            while len(models) < MAX_SOLUTIONS:
                ret = self.solver.check()
                if ret == unknown:
                    print(f'\n{self.header}unknown error, most possibly aborted.')
                    raise SystemExit(1)
                if ret != sat: break
                model = self.solver.model()
                initial = [model[i].as_long() if model[i] is not None else 0 for i in vars]
                self.solver.add(Or(*[i != val for i, val in zip(vars, initial)]))
                models.append((initial, rollmany(*initial, self.nsteps)))
                a, b, c, d = initial
                print(
                    f'{self.header}at least {len(models)} solution(s) found ' +
                    f'({nretained} retained), e.g. {a:08X} {b:08X} {c:08X} {d:08X}', end='\r')
            self.solver.pop()

        print(f'\x1b[2K\r{self.header}', end='')
        if len(models) == 0:
            print('no solution possible.')
        else:
            (a, b, c, d), _ = models[0]
            if len(models) >= MAX_SOLUTIONS:
                print('at least ', end='')
            print(
                f'{len(models)} solution(s) found ({nretained} retained), ' +
                f'e.g. {a:08X} {b:08X} {c:08X} {d:08X}')

        self.lastmodels = models
        return len(models)

    @property
    def header(self):
        elapsed = int(time.time() - self.starttime)
        hm, s = divmod(elapsed, 60)
        h, m = divmod(hm, 60)
        return f'[{self.nsteps}/{self.ncumulbits:.2f}/{h}:{m:02d}:{s:02d}] '

def ground_truth():
    a, b, c, d = [random.randint(0, 0xffffffff) for i in range(4)]
    #a, b, c, d = [int(i, 16) for i in '2B71BDF9 73681330 E2A8E331 5E1AF65F'.split()]
    print(f'ground truth: {a:08X} {b:08X} {c:08X} {d:08X}')

    while True:
        a, b, c, d = b, c, d, roll(a, b, c, d)
        yield d

if __name__ == '__main__':
    if False:
        seq = Sequence()
        for value in ground_truth():
            mult = random.randint(3, 20)
            quantized = (value * mult) >> 32
            lo = (quantized << 32) // mult
            hi = ((quantized + 1) << 32) // mult
            nsolutions = seq.step(lo, hi)
            if nsolutions <= 1: break
    else:
        # https://gist.github.com/Treeki/85be14d297c80c8b3c0a76375743325b
        # ./TurnipPrices 3 12
        steps = mksteps(91, [79, 75, 71, 67, 63, 59, 56, 53, 50, 46, 43, 39], PAT_SMALLSPIKE, PAT_DECREASING)
        seq = Sequence()
        for lo, hi in steps:
            nsolutions = seq.step(lo, hi)
            if nsolutions <= 1: break