bound2.py

import argparse, sys, random, math

LOGLEVEL = 0

def line_eq(p1, p2):
    xd = p2[0] - p1[0]
    yd = p2[1] - p1[1]
    return (-yd, xd, yd * p1[0] - xd * p1[1])

def cross(o, a, b):
    """2D cross product of the OA and OB vectors, i.e. the z-component of
    their 3D cross product.

    Returns a positive value if OAB makes a counter-clockwise turn,
    negative for clockwise turn, and zero if the points are collinear.
    """
    return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])

def gmin(a, b):
    if a is None:
        return b
    return min(a, b)

def gmax(a, b):
    if a is None:
        return b
    return max(a, b)

def convex_hull(points):
    """Computes the convex hull of a set of 2D points.

    Input: an iterable sequence of (x, y) pairs representing the points.
    Output: a list of vertices of the convex hull in counter-clockwise order,
      starting from the vertex with the lexicographically smallest coordinates.
    Implements Andrew's monotone chain algorithm. O(n log n) complexity.
    """
    # Boring case: no points.
    if len(points) <= 1:
        return points
    # Sort the points lexicographically (tuples are compared lexicographically).
    points = sorted(points)
    # Build lower hull.
    lower = [points[0]]
    for p in points[1:]:
        while len(lower) >= 2 and cross(lower[-2], lower[-1], p) <= 0:
            lower.pop()
        if lower[-1] != p:
            lower.append(p)
    if (len(lower) == 1):
        return lower
    # Build upper hull.
    upper = [points[-1]]
    for p in reversed(points[:-1]):
        while len(upper) >= 2 and cross(upper[-2], upper[-1], p) <= 0:
            upper.pop()
        if upper[-1] != p:
            upper.append(p)
    # Concatenation of the lower and upper hulls gives the convex hull.
    # The Last point of each list is omitted because it is repeated at the beginning
    # of the other list. 
    return lower[:-1] + upper[:-1]

# DeltaMap: (negate: bool, inc: int)
#     delta -> (1 - 2*negate)*delta + irc
#
# EvenMapConfig:
# - "deltamap": DeltaMap
# - "neg": bool
#     (delta, f, g) -> (deltamap(delta), f, (1 - 2*neg) * (g/2))
#
# OddMapConfig:
# - "deltamap": DeltaMap
# - "swap": bool
# - "neg_f": bool
# - "neg_g": bool
#    (delta, f, g) -> (deltamap(delta), swap ? g : f, ((1 - 2*neg_f)*f + (1 - 2*neg_g)*g)/2)
#
# Config:
# - "even": EvenMapConfig
# - "odd0": OddMapConfig
# - "odd1": OddMapConfig

def apply_deltamap(dm, delta):
    if dm[0]:
        return -delta + dm[1]
    else:
        return delta + dm[1]

def apply_neg(cfg, val):
    if cfg:
        return -val
    else:
        return val

def divstep(state, it, maxmod, cfg):
    ns = {}
    even = cfg["even"]
    even_deltamap = even["deltamap"]
    even_neg = even["neg"]
    odd0 = cfg["odd0"]
    odd0_deltamap = odd0["deltamap"]
    odd0_swap = odd0["swap"]
    odd0_neg_f = odd0["neg_f"]
    odd0_neg_g = odd0["neg_g"]
    odd1 = cfg["odd1"]
    odd1_deltamap = odd1["deltamap"]
    odd1_swap = odd1["swap"]
    odd1_neg_f = odd1["neg_f"]
    odd1_neg_g = odd1["neg_g"]

    for delta, hull in state.items():
        if hull is None:
            continue

        nhull = ns.setdefault(apply_deltamap(even_deltamap, delta), [])
        for f, g in hull:
            nhull.append((f<<1, apply_neg(even_neg, g)))
        if delta >= 0:
            nhull = ns.setdefault(apply_deltamap(odd0_deltamap, delta), [])
            for f, g in hull:
                nhull.append(((g if odd0_swap else f) << 1, apply_neg(odd0_neg_f, f) + apply_neg(odd0_neg_g, g)))
        else:
            nhull = ns.setdefault(apply_deltamap(odd1_deltamap, delta), [])
            for f, g in hull:
                nhull.append(((g if odd1_swap else f) << 1, apply_neg(odd1_neg_f, f) + apply_neg(odd1_neg_g, g)))

    for delta in ns:
        min_g = min(g for _, g in ns[delta])
        max_g = max(g for _, g in ns[delta])
        if min_g*maxmod <= -(1 << it) or max_g*maxmod >= (1 << it):
            ns[delta] = convex_hull(ns[delta])
        else:
            ns[delta] = None

    return ns

def getbound(maxmod, maxit, cfgs):
    STATE = {0: [(0, 0), (1, 0), (1, 1)]}
    it = 0
    while len(STATE) and it < maxit:
        it += 1
        STATE = divstep(STATE, it, maxmod, cfgs[it & 1])
    return it

def print_deltamap(dm):
    return ("-" if dm[0] else "") + "delta" + ("+%i" % dm[1] if dm[1] > 0 else "%i" % dm[1] if dm[1] < 0 else "")

def print_even(cfg):
    return "(%s,f,%sg/2)" % (print_deltamap(cfg["deltamap"]), "-" if cfg["neg"] else "")

def print_odd(cfg):
    return "(%s,%s,(%sf%sg)/2)" % (print_deltamap(cfg["deltamap"]), "g" if cfg["swap"] else "f", "-" if cfg["neg_f"] else "", "-" if cfg["neg_g"] else "+")

def print_config(cfg):
    o0 = print_odd(cfg["odd0"])
    o1 = print_odd(cfg["odd1"])
    e = print_even(cfg["even"])
    if o0 == o1:
        return "%s if g & 1 else %s" % (o0, e)
    else:
        return "(%s if delta >= 0 else %s) if g & 1 else %s" % (o0, o1, e)

DELTAMAPS = [(m, i) for m in [False, True] for i in [-3,-2,-1,0,1,2,3]]
EVENMAPCFGS = [{"deltamap": d, "neg": n} for d in DELTAMAPS for n in [False, True]]
ODDMAPCFGS = [{"deltamap": d, "swap": s, "neg_f": nf, "neg_g": ng} for d in DELTAMAPS for s in [False, True] for nf in [False, True] for ng in [False, True]]

BESTCFG = {"even": {"deltamap": (False, 1), "neg": False}, "odd0": {"deltamap": (True, 0), "neg_f": True, "neg_g": False, "swap": True}, "odd1": {"deltamap": (False, 1), "neg_f": False, "neg_g": False, "swap": False}}
SAFEGCDCFG = {"even": {"deltamap": (False, 1), "neg": False}, "odd0": {"deltamap": (True, -1), "neg_f": True, "neg_g": False, "swap": True}, "odd1": {"deltamap": (False, 1), "neg_f": False, "neg_g": False, "swap": False}}

while True:
    ae = random.choice(EVENMAPCFGS)
    ao0 = random.choice(ODDMAPCFGS)
    ao1 = random.choice(ODDMAPCFGS)
    be = random.choice(EVENMAPCFGS)
    bo0 = random.choice(ODDMAPCFGS)
    bo1 = random.choice(ODDMAPCFGS)

    cfg = [{"even": be, "odd0": bo0, "odd1": bo1}, {"even": ae, "odd0": ao0, "odd1": ao1}]
    bnd100 = getbound(100, 23, cfg)
    if bnd100 < 23:
        bnd1M = getbound(1000000, 60, cfg)
        if bnd1M < 60:
            bnd2p256 = getbound(2**256, 800, cfg)
            if bnd2p256 < 800:
                print("100=%i,1M=%i,2^256=%i: %s; %s" % (bnd100, bnd1M, bnd2p256, print_config(cfg[0]), print_config(cfg[1])))
            else:
                print("100=%i,1M=%i: %s; %s" % (bnd100, bnd1M, print_config(cfg[0]), print_config(cfg[1])))

fgcd.cpp

#include <stdint.h>
#include <assert.h>
#include <stdio.h>
#include <algorithm>
#include <bitset>
#include <set>
#include <vector>
#include <thread>
#include <mutex>
#include <utility>
#include <queue>

typedef signed __int128 int128_t;
typedef unsigned __int128 uint128_t;

#define LOOKAHEAD 21
#define COUNT 1000
#define THREADS 30

#define KEEPBEST (((size_t)1) << (LOOKAHEAD - 5))

namespace {

inline uint64_t rdrand()
{
    uint8_t ok;
    uint64_t r;
    while (1) {
        __asm__ volatile (".byte 0x48, 0x0f, 0xc7, 0xf0; setc %1" : "=a"(r), "=q"(ok) :: "cc"); // rdrand %rax
        if (ok) break;
        __asm__ volatile ("pause" :);
    }
    return r;
}

class RNG {
    uint64_t rng_state[2];
    uint64_t rng_cache = 0;
    int rng_bits = 0;
public:
    RNG() {
        rng_state[0] = rdrand();
        rng_state[1] = rdrand();
    }
    using result_type = uint64_t;
    constexpr uint64_t min() const { return 0; }
    constexpr uint64_t max() const { return 0xFFFFFFFFFFFFFFFF; }
    uint64_t operator()() {
        rng_state[0] += 0x60bee2bee120fc15;
        rng_state[1] += 0xa88615625dc1ee97;
        uint128_t tmp1 = (uint128_t)rng_state[0] * 0xa3b195354a39b70d;
        uint128_t tmp2 = (uint128_t)rng_state[1] * 0xa3b195354a39b70d;
        uint64_t m1 = (tmp1 >> 64) ^ tmp1;
        uint64_t m2 = (tmp2 >> 64) ^ tmp2;
        tmp1 = (uint128_t)m1 * 0x1b03738712fad5c9;
        tmp2 = (uint128_t)m2 * 0x1b03738712fad5c9;
        return tmp1 ^ tmp2 ^ (tmp1 >> 64) ^ (tmp2 >> 64);
    }

    bool bit() {
        if (rng_bits == 0) {
            rng_cache = operator()();
            rng_bits = 64;
        }
        bool ret = rng_cache & 1;
        rng_cache >>= 1;
        rng_bits -= 1;
        return ret;
    }
};

typedef struct {
    int64_t v[5];
} secp256k1_modinv64_signed62;

static int secp256k1_modinv64_divsteps_62_var(uint64_t* peta, uint64_t f0, uint64_t g0, int64_t *t) {
    uint64_t u = 1, v = 0, q = 0, r = 1;
    uint64_t f = f0, g = g0;
    int i;
    int iter = 0;
    uint64_t eta = *peta;
    for (i = 0; i < 62; ++i) {
        if (g & 1) {
            iter = i + 1;
            if ((int64_t)eta < 0) {
                g -= f;
                q -= u;
                r -= v;
                eta = ~eta;
                f += g;
                u += q;
                v += r;
            } else {
                g += f;
                q += u;
                r += v;
                --eta;
            }
        } else {
            --eta;
        }
        g >>= 1;
        u <<= 1;
        v <<= 1;
    }
    t[0] = (int64_t)u;
    t[1] = (int64_t)v;
    t[2] = (int64_t)q;
    t[3] = (int64_t)r;
    *peta = eta;
    return iter;
}

static void secp256k1_modinv64_update_fg_62_var(int len, secp256k1_modinv64_signed62 *f, secp256k1_modinv64_signed62 *g, const int64_t *t) {
    const int64_t M62 = (int64_t)(UINT64_MAX >> 2);
    const int64_t u = t[0], v = t[1], q = t[2], r = t[3];
    int64_t fi, gi;
    int128_t cf, cg;
    int i;

    fi = f->v[0];
    gi = g->v[0];

    cf = (int128_t)u * fi + (int128_t)v * gi;
    cg = (int128_t)q * fi + (int128_t)r * gi;

    cf >>= 62;
    cg >>= 62;

    for (i = 1; i < len; ++i) {
        fi = f->v[i];
        gi = g->v[i];

        cf += (int128_t)u * fi + (int128_t)v * gi;
        cg += (int128_t)q * fi + (int128_t)r * gi;

        f->v[i - 1] = (int64_t)cf & M62; cf >>= 62;
        g->v[i - 1] = (int64_t)cg & M62; cg >>= 62;
    }

    f->v[len - 1] = (int64_t)cf;
    g->v[len - 1] = (int64_t)cg;
}

static int secp256k1_modinv64_iters_var(const secp256k1_modinv64_signed62 *x, const secp256k1_modinv64_signed62* modulus) {
    /* Modular inversion based on the paper "Fast constant-time gcd computation and
     * modular inversion" by Daniel J. Bernstein and Bo-Yin Yang. */

    int64_t t[4];
    secp256k1_modinv64_signed62 f = *modulus;
    secp256k1_modinv64_signed62 g = *x;
    int i, j, len = 5;
    uint64_t eta;
    int64_t cond, fn, gn;
    int iter = 0;

    /* The paper uses 'delta'; eta == -delta (a performance tweak).
     *
     * If g has leading zeros (w.r.t 256 bits), then eta can be set initially to
     * -(1 + clz(g)), and the worst-case divstep count would be only (741 - clz(g)). */
    eta = -(uint64_t)1;

    for (i = 0; i < 12; ++i) {
        iter = secp256k1_modinv64_divsteps_62_var(&eta, f.v[0], g.v[0], t);
        secp256k1_modinv64_update_fg_62_var(len, &f, &g, t);

        if (g.v[0] == 0) {
            cond = 0;
            for (j = 1; j < len; ++j) {
                cond |= g.v[j];
            }
            if (cond == 0) {
                iter += 62 * i;
                break;
            }
        }

        fn = f.v[len - 1];
        gn = g.v[len - 1];

        cond = ((int64_t)len - 2) >> 63;
        cond |= fn ^ (fn >> 63);
        cond |= gn ^ (gn >> 63);

        if (cond == 0) {
            f.v[len - 2] |= fn << 62;
            g.v[len - 2] |= gn << 62;
            --len;
        }
    }

    return iter;
}


static const secp256k1_modinv64_signed62 MODULUS = {{-0x1000003D1LL, 0, 0, 0, 256}};

using num256 = std::bitset<256>;

struct compare_num256{
    bool operator()(const num256& a, const num256& b) const
    {
        for (int i = 0; i < 256; ++i) {
            if (a[i] < b[i]) return true;
            if (a[i] > b[i]) return false;
        }
        return false;
    }
};

using entry = std::pair<double, num256>;

using bestentry = std::pair<int, num256>;

struct compare_entry {
    bool operator()(const entry& a, const entry& b) const
    {
        if (a.first > b.first) return true;
        if (a.first < b.first) return false;
        return false;
    }
};

struct compare_bestentry_full {
    compare_num256 nc;
    bool operator()(const bestentry& a, const bestentry& b) const
    {
        if (a.first > b.first) return true;
        if (a.first < b.first) return false;
        return nc(a.second, b.second);
    }
};

struct compare_bestentry_iter {
    bool operator()(const bestentry& a, const bestentry& b) const
    {
        if (a.first > b.first) return true;
        if (a.first < b.first) return false;
        return false;
    }
};

static void loads62(secp256k1_modinv64_signed62* r, const num256& num) {
    for (int n = 0; n < 5; ++n) {
        r->v[n] = 0;
    }
    for (int p = 0; p < 256; ++p) {
        r->v[p / 62] |= ((uint64_t)num[p]) << (p % 62);
    }
}

num256 masklow(const num256& num, int bits) {
    num256 ret;
    for (int i = 0; i < bits; ++i) ret[i] = num[i];
    return ret;
}

std::string num256_to_string(const num256& num) {
    std::string ret;
    static constexpr const char* HEX = "0123456789abcdef";
    for (int c = 0; c < 64; ++c) {
        uint64_t n = 0;
        for (int bit = 0; bit < 4; ++bit) {
            n = (n << 1) | num[255 - 4*c - bit];
        }
        ret += HEX[n];
    }
    return ret;
}

void analyze(uint64_t count_per_thread, int suffix_len, std::vector<entry>& suffices, std::vector<bestentry>& best) {
    std::vector<std::thread> threads;
    std::mutex mutex;
    for (uint64_t t = 0; t < THREADS; ++t) {
        threads.emplace_back([&](){
            RNG rng;
            std::vector<uint64_t> lsums;
            std::vector<std::pair<int, num256>> lbest;
            lsums.resize(suffices.size());
            lbest.reserve(KEEPBEST * 2);
            int cutoff = 0;
            num256 num;
            for (uint64_t c = 0; c < count_per_thread; ++c) {
                for (int b = suffix_len; b < 256; ++b) {
                    num[b] = rng.bit();
                }
                for (size_t l = 0; l < lsums.size(); ++l) {
                    for (int b = 0; b < suffix_len; ++b) {
                        num[b] = suffices[l].second[b];
                    }
                    secp256k1_modinv64_signed62 s62;
                    loads62(&s62, num);
                    int iter = secp256k1_modinv64_iters_var(&s62, &MODULUS);
//                    printf("  - try %s: %i\n", num256_to_string(num).c_str(), iter);
                    lsums[l] += iter;
                    if (iter >= cutoff) {
                        lbest.emplace_back(iter, num);
                        if (lbest.size() >= 2 * KEEPBEST) {
                            std::partial_sort(lbest.begin(), lbest.begin() + KEEPBEST, lbest.end(), compare_bestentry_iter{});
                            lbest.resize(KEEPBEST);
                            cutoff = lbest.back().first;
                        }
                    }
                }
            }
            std::unique_lock<std::mutex> lock(mutex);
            for (const auto& entry : lbest) {
                 best.push_back(entry);
            }
            for (size_t l = 0; l < lsums.size(); ++l) {
                suffices[l].first += lsums[l];
            }
        });
    }
    for (auto& thread : threads) thread.join();
    if (best.size() > KEEPBEST) {
        std::sort(best.begin(), best.end(), compare_bestentry_full{});
        best.erase(std::unique(best.begin(), best.end()), best.end());
    }
    if (best.size() > KEEPBEST) {
        std::shuffle(best.begin(), best.end(), RNG());
        std::partial_sort(best.begin(), best.begin() + KEEPBEST, best.end(), compare_bestentry_iter{});
        best.resize(KEEPBEST);
    }
}

} // namespace


int main(void) {
    setbuf(stdout, NULL);

    std::vector<bestentry> best;

    std::vector<entry> suffices;
    for (uint64_t l = 0; l >> LOOKAHEAD == 0; ++l) {
        num256 num;
        for (int b = 0; b < LOOKAHEAD; ++b) {
            num[b] = (l >> b) & 1;
        }
        suffices.emplace_back(0, num);
    }

    double weight = 0;
    for (int bits = LOOKAHEAD; bits < 256; ++bits) {
        uint64_t count = COUNT;
        if (bits == LOOKAHEAD) count *= 2;
        if (256 - bits < 28) {
            uint64_t max = (uint64_t)5 << (2 * (256 - bits));
            if (count > max) count = max;
        }
        uint64_t count_per_thread = (count + THREADS - 1) / THREADS;
        analyze(count_per_thread, bits, suffices, best);
        weight += count_per_thread * THREADS;

        std::set<num256, compare_num256> subset;
        for (const auto& entry : best) subset.insert(masklow(entry.second, bits));
        std::shuffle(suffices.begin(), suffices.end(), RNG());
        std::sort(suffices.begin(), suffices.end(), compare_entry{});
        std::vector<entry> nsuffices;
        size_t retained = 0;
        for (size_t pos = 0; pos < suffices.size(); ++pos) {
            const auto& entry = suffices[pos];
            if (subset.count(entry.second)) {
                if ((pos * 2) >> LOOKAHEAD) ++retained;
                num256 num = entry.second;
                nsuffices.emplace_back(entry.first * 0.5, num);
                num[bits] = 1;
                nsuffices.emplace_back(entry.first * 0.5, num);
            }
        }
        for (const auto& entry : suffices) {
            if (!subset.count(entry.second)) {
                num256 num = entry.second;
                nsuffices.emplace_back(entry.first * 0.5, num);
                num[bits] = 1;
                nsuffices.emplace_back(entry.first * 0.5, num);
                if (nsuffices.size() >> LOOKAHEAD) break;
            }
        }
        weight *= 0.5;
        suffices = std::move(nsuffices);
        printf("bit=%i, pat=[%f %s], best=[%i %s], retained=%lu\n", bits, suffices[0].first / weight, num256_to_string(suffices[0].second).c_str(), best[0].first, num256_to_string(best[0].second).c_str(), (unsigned long)retained);
    }

    for (size_t pos = 0; pos < best.size(); ++pos) {
        printf("top#%lu=[%i %s]\n", (unsigned long)pos, best[pos].first, num256_to_string(best[pos].second).c_str());
    }

    return 0;
}

fgcd.py

#/bin/env python3

import sympy
import math

def iter_count(bits):
    return (49*bits+80)//17 if bits < 46 else (49*bits+57)//17

INV256 = [0xFF, 0x55, 0x33, 0x49, 0xC7, 0x5D, 0x3B, 0x11, 0x0F, 0xE5, 0xC3, 0x59, 0xD7, 0xED, 0xCB, 0x21, 0x1F, 0x75, 0x53, 0x69, 0xE7, 0x7D, 0x5B, 0x31, 0x2F, 0x05, 0xE3, 0x79, 0xF7, 0x0D, 0xEB, 0x41, 0x3F, 0x95, 0x73, 0x89, 0x07, 0x9D, 0x7B, 0x51, 0x4F, 0x25, 0x03, 0x99, 0x17, 0x2D, 0x0B, 0x61, 0x5F, 0xB5, 0x93, 0xA9, 0x27, 0xBD, 0x9B, 0x71, 0x6F, 0x45, 0x23, 0xB9, 0x37, 0x4D, 0x2B, 0x81, 0x7F, 0xD5, 0xB3, 0xC9, 0x47, 0xDD, 0xBB, 0x91, 0x8F, 0x65, 0x43, 0xD9, 0x57, 0x6D, 0x4B, 0xA1, 0x9F, 0xF5, 0xD3, 0xE9, 0x67, 0xFD, 0xDB, 0xB1, 0xAF, 0x85, 0x63, 0xF9, 0x77, 0x8D, 0x6B, 0xC1, 0xBF, 0x15, 0xF3, 0x09, 0x87, 0x1D, 0xFB, 0xD1, 0xCF, 0xA5, 0x83, 0x19, 0x97, 0xAD, 0x8B, 0xE1, 0xDF, 0x35, 0x13, 0x29, 0xA7, 0x3D, 0x1B, 0xF1, 0xEF, 0xC5, 0xA3, 0x39, 0xB7, 0xCD, 0xAB, 0x01]

def matl2norm(t):
    u, v, q, r = t
    # 1/2*q^2 + 1/2*r^2 + 1/2*u^2 + 1/2*v^2 + 1/2*sqrt(q^4 + 2*q^2*r^2 + r^4 + u^4 + 8*q*r*u*v + v^4 + 2*(q^2 - r^2)*u^2 - 2*(q^2 - r^2 - u^2)*v^2)
    u2, v2, q2, r2 = u*u, v*v, q*q, r*r
    d = q2*q2 + 2*q2*r2 + r2*r2 + u2*u2 + 8*q*r*u*v + v2*v2 + 2*(q2-r2)*u2 - 2*(q2 - r2 - u2)*v2
    l = q2+r2+u2+v2
    return math.sqrt((l + math.sqrt(d)) * 0.5)

def wyhash(n):
    a = 0xa0761d6478bd642f
    b = n ^ 0xe7037ed1a0b428db
    r = a * b
    a ^= r & 0xffffffffffffffff
    b ^= r >> 64
    a ^= 0xa0761d6478bd642f
    b ^= 0xe7037ed1a0b428db
    r = a * b
    a ^= r & 0xffffffffffffffff
    b ^= r >> 64
    return a ^ b

def naive_modinv(a, n):
    t1, t2 = 0, 1
    r1, r2 = n, a
    while r2 != 0:
        q = r1 // r2
        t1, t2 = t2, t1 - q * t2
        r1, r2 = r2, r1 - q * r2
    if r1 > 1:
        return None
    if t1 < 0:
        t1 += n
    return t1

def divsteps_group(eta, f0, g0, groupsize):
    assert (f0 & 1) == 1
    assert f0 >= 0
    assert g0 >= 0
    f = f0
    g = g0
    u = 1
    v = 0
    q = 0
    r = 1
    i = groupsize
    mask = (1 << (groupsize + 2)) - 1
    while True:
        x = g | ((-1) << i)
        zeros = (x & -x).bit_length() - 1
        g = (g >> zeros) & mask
        u <<= zeros
        v <<= zeros
        eta -= zeros
        i -= zeros
        if i <= 0:
            break
        assert (f & 1) == 1
        assert (g & 1) == 1
#        assert (u * f0 + v * g0) & mask == (f << (groupsize - i)) & mask
#        assert (q * f0 + r * g0) & mask == (g << (groupsize - i)) & mask
        if eta < 0:
            eta, f, g, u, q, v, r = -eta, g, (-f) & mask, q, -u, r, -v
        limit = min(min(eta + 1, i), 8)
        m = (1 << limit) - 1
        w = ((g & m) * INV256[(f >> 1) & 127]) & m
        g = (g + f * w) & mask
        q += u * w
        r += v * w
        assert (g & m) == 0
    assert u >= -(1 << (groupsize + 1)) and u < (1 << (groupsize + 1))
    assert v >= -(1 << (groupsize + 1)) and v < (1 << (groupsize + 1))
    assert q >= -(1 << (groupsize + 1)) and q < (1 << (groupsize + 1))
    assert r >= -(1 << (groupsize + 1)) and r < (1 << (groupsize + 1))
    return eta, (u, v, q, r)

def calc_mont(m, groupsize):
    mm = (1 << groupsize)
    return naive_modinv(mm - (m & (mm -1)), mm)

def update_de(d, e, t, m, mont, groupsize, mask):
    u, v, q, r = t
    cd = u * d + v * e
    ce = q * d + r * e
    md = (mont * (cd & mask)) & mask
    if md >> (groupsize - 1):
        md -= (1 << groupsize)
    me = (mont * (ce & mask)) & mask
    if me >> (groupsize - 1):
        me -= (1 << groupsize)
    cd += m * md
    ce += m * me
    assert (cd & mask) == 0
    assert (ce & mask) == 0
    return cd >> groupsize, ce >> groupsize

def update_fg(f, g, t, groupsize, mask):
    u, v, q, r = t
    cf = u * f + v * g
    assert (cf & mask) == 0
    cg = q * f + r * g
    assert (cg & mask) == 0
    return cf >> groupsize, cg >> groupsize

MAX_GROUPSIZE=16

LOWEST_RATIO=[0 for _ in range(MAX_GROUPSIZE+1)]
HIGHEST_RATIO=[0 for _ in range(MAX_GROUPSIZE+1)]
LOWEST_END_RATIO=[0 for _ in range(MAX_GROUPSIZE+1)]
HIGHEST_END_RATIO=[0 for _ in range(MAX_GROUPSIZE+1)]
LARGEST_T2N=[0 for _ in range(MAX_GROUPSIZE+1)]

def modinv(x, m, mont, groupsize):
    global LOWEST_RATIO, HIGHEST_RATIO, LOWEST_END_RATIO, HIGHEST_END_RATIO
    d = 0
    e = 1
    f = m
    g = x
    eta = -1
    assert (m & 1) == 1
    bits = max(x.bit_length(), m.bit_length())
    max_iter = iter_count(bits)
    max_groups = (max_iter + groupsize - 1) // groupsize
    mask = (1 << groupsize) - 1

    gg = 0
    for group in range(max_groups + 3):
        eta, t = divsteps_group(eta, f & mask, g & mask, groupsize)
        d, e = update_de(d, e, t, m, mont, groupsize, mask)
        f, g = update_fg(f, g, t, groupsize, mask)
        gg += 1

        assert f >= -m and f <= m
        assert g >= -m and f <= m

        if d < 0 or e < 0:
            ratio = min(d / m, e / m)
            if ratio < LOWEST_RATIO[groupsize]:
                LOWEST_RATIO[groupsize] = ratio

        if d > 0 or e > 0:
            ratio = max(d / m, e / m)
            if ratio > HIGHEST_RATIO[groupsize]:
                HIGHEST_RATIO[groupsize] = ratio

        LARGEST_T2N[groupsize] = max(LARGEST_T2N[groupsize], matl2norm(t))

        if g == 0:
            break

    assert (g == 0)
    assert (x == 0) or (f == 1) or (f == -1)

    if d < 0:
        ratio = d / m
        if ratio < LOWEST_END_RATIO[groupsize]:
            LOWEST_END_RATIO[groupsize] = ratio

    if d > 0:
        ratio = d / m
        if ratio > HIGHEST_END_RATIO[groupsize]:
            HIGHEST_END_RATIO[groupsize] = ratio

    if f == -1:
        return (m - (d % m), gg)
    else:
        return (d % m, gg)

import sys

CPU = int(sys.argv[1]) if len(sys.argv) > 1 else 0
CPUS = int(sys.argv[2]) if len(sys.argv) > 1 else 1

LASTLOG = 0
IT = 0
M = 1
while True:
    M += 2
    if not sympy.isprime(M):
        continue
    if CPUS > 1 and wyhash(M) % CPUS != CPU:
        continue
    MONTS = [0 for _ in range(17)]
    for GROUPSIZE in range(1, 17):
        MONTS[GROUPSIZE] = calc_mont(M, GROUPSIZE)
    print("M=%i" % M)
    for X in range(M):
        ITERS = 0
        for GROUPSIZE in range(1, 17):
            I, GG = modinv(X, M, MONTS[GROUPSIZE], GROUPSIZE)
            IT += 1
            if GROUPSIZE == 1:
                ITERS = GG
            else:
                assert GG == (ITERS + GROUPSIZE - 1) // GROUPSIZE
    if IT - LASTLOG > 10000000:
        LASTLOG = IT
        for groupsize in range(1, 17):
                print("groupsize=%i: any=%f..%f end=%f..%f t2n=%f" % (groupsize, LOWEST_RATIO[(groupsize)], HIGHEST_RATIO[(groupsize)], LOWEST_END_RATIO[(groupsize)], HIGHEST_END_RATIO[(groupsize)], LARGEST_T2N[groupsize] / 2**groupsize))

fgcd_recurse.cpp

#include <assert.h>
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <gmp.h>
#include <string>
#include <algorithm>
#include <queue>
#include <thread>
#include <chrono>
#include <atomic>
#include <mutex>
#include <condition_variable>

namespace {


// Data structure size (modulus and range need to fit in this many bits)
static constexpr uint32_t B = 256;

// Modulus to use, as a hex string.
static constexpr const char* MODULUS_STR = "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364141";
//static constexpr const char* MODULUS_STR = "1fffffffe7";

// Number of threads
static constexpr uint32_t THREADS = 32;

// Permitted range on the input x to gcd(MODULUS,x) is [0..2**RANGE_BITS-1]
static constexpr uint32_t RANGE_BITS = 256;

// Prune all subtrees that have a bound up to this value.
static constexpr uint32_t PRUNE_BOUND = 736;

constexpr size_t my_strlen(const char* c) {
    const char* e = c;
    while (*e) ++e;
    return e - c;
}

template<uint32_t LIMBS, bool SIGNED>
class Num {
    static_assert(LIMBS >= 1);
    mp_limb_t data[LIMBS];
    static constexpr uint32_t NEG_FLAG = 0x80000000;
    static constexpr uint32_t ABS_MASK = 0x7FFFFFFF;
    uint32_t size; // (size & ABS_MASK) = length; (size & NEG_FLAG) = negative number

    template<uint32_t LIMBS1, bool SIGNED1> friend class Num;

    uint32_t abs_size() const {
        if constexpr (LIMBS == 1) return 1;
        if constexpr (SIGNED) {
            return size & ABS_MASK;
        } else {
            return size;
        }
    }

    uint32_t neg_flag() const {
        if constexpr (SIGNED) {
            return size & NEG_FLAG;
        } else {
            return 0;
        }
    }

    void normalize() {
        if constexpr (SIGNED) {
            uint32_t a = abs_size();
            while (a && data[a - 1] == 0) --a;
            if (a == 0) {
                size = 1;
            } else {
                size = (size & NEG_FLAG) | a;
            }
        } else {
            while (size > 1 && data[size - 1] == 0) --size;
        }
    }

    int32_t signed_size() const {
        if constexpr (SIGNED) {
            return size & NEG_FLAG ? -(size ^ NEG_FLAG) : size;
        } else {
            return size;
        }
    }

    template<bool SIGNED1, bool SIGNED2>
    void internal_set_sum(const mp_limb_t* u, uint32_t usize, uint32_t neg_u, const mp_limb_t* v, uint32_t vsize, uint32_t neg_v) {
//        printf("internal_set_sum: S1=%i S2=%i usize=%u negu=%x vsize=%u negv=%x\n", (int)SIGNED1, (int)SIGNED2, (unsigned)usize, (unsigned)neg_u, (unsigned)vsize, (unsigned)neg_v);
        assert((usize & NEG_FLAG) == 0);
        assert((vsize & NEG_FLAG) == 0);
        assert(SIGNED1 || (neg_u == 0));
        assert(SIGNED2 || (neg_v == 0));
        if constexpr (SIGNED1 || SIGNED2) {
            if (neg_u != neg_v) {
                if (usize > vsize || ((usize == vsize) && mpn_cmp(u, v, usize) >= 0)) {
                    mp_limb_t borrow = mpn_sub(data, u, usize, v, vsize);
                    assert(borrow == 0);
                    assert(SIGNED || !neg_u);
                    size = SIGNED ? (usize | neg_u) : usize;
                } else {
                    mp_limb_t borrow = mpn_sub(data, v, vsize, u, usize);
                    assert(borrow == 0);
                    assert(SIGNED || !neg_v);
                    size = SIGNED ? (vsize | neg_v) : vsize;
                }
                normalize();
                return;
            }
        }
        if (usize >= vsize) {
            mp_limb_t overflow = mpn_add(data, u, usize, v, vsize);
            assert((SIGNED && SIGNED1) || !neg_u);
            size = (SIGNED && SIGNED1) ? (usize | neg_u) : usize;
            if (overflow) {
                assert(usize < LIMBS);
                data[usize] = overflow;
                size = (SIGNED && SIGNED1) ? ((usize + 1) | neg_u) : usize + 1;
            }
        } else {
            mp_limb_t overflow = mpn_add(data, v, vsize, u, usize);
            assert((SIGNED && SIGNED2) || !neg_v);
            size = (SIGNED && SIGNED2) ? (vsize | neg_v) : vsize;
            if (overflow) {
                assert(vsize < LIMBS);
                data[vsize] = overflow;
                size = (SIGNED && SIGNED2) ? ((vsize + 1) | neg_v) : vsize + 1;
            }
        }
    }

    template<bool SELF, uint32_t LIMBS1, bool SIGNED1>
    void internal_set_rshift(const Num<LIMBS1, SIGNED1>& in1, unsigned shift) {
        uint32_t a = in1.abs_size();
        if (shift >= GMP_NUMB_BITS * a) {
            if (!SELF) set_zero();
            return;
        }
        unsigned n = shift / GMP_NUMB_BITS;
        a -= n;
        assert(a <= LIMBS);
        if (shift % GMP_NUMB_BITS) {
            mpn_rshift(data, in1.data + n, a, shift % GMP_NUMB_BITS);
            size = in1.neg_flag() | a;
            normalize();
        } else {
            if constexpr (SELF) {
                memmove(data, in1.data + n, sizeof(mp_limb_t) * a);
            } else {
                memcpy(data, in1.data + n, sizeof(mp_limb_t) * a);
            }
            size = in1.neg_flag() | a;
        }
    }

    template<bool SELF, uint32_t LIMBS1, bool SIGNED1>
    void internal_set_lshift(const Num<LIMBS1, SIGNED1>& in1, unsigned shift) {
        if (in1.is_zero()) {
            if (!SELF) set_zero();
            return;
        }
        uint32_t a = in1.abs_size();
        unsigned n = shift / GMP_NUMB_BITS;
        assert(a + n <= LIMBS);
        if (shift % GMP_NUMB_BITS) {
            mp_limb_t overflow = mpn_lshift(data + n, in1.data, a, shift % GMP_NUMB_BITS);
            a += n;
            size = in1.neg_flag() | a;
            if (overflow) {
                assert(a < LIMBS);
                data[a] = overflow;
                size = in1.neg_flag() | (a + 1);
            }
        } else {
            if constexpr (SELF) {
                memmove(data + n, in1.data, sizeof(mp_limb_t) * a);
            } else {
                memcpy(data + n, in1.data, sizeof(mp_limb_t) * a);
            }
            a += n;
            size = in1.neg_flag() | a;
        }
        for (uint32_t i = 0; i < n; ++i) data[i] = 0;
    }


public:
    constexpr Num() : data{0}, size(1) {}

    constexpr uint32_t limbs() const { return LIMBS; }

    constexpr Num(const char* begin, const char* end) : data{0}, size(0) {
        static_assert(GMP_NUMB_BITS % 4 == 0);
        size_t digits = 0;
        uint32_t sign = 0;
        auto it = begin;
        if (it != end && *it == '-') {
            assert(SIGNED);
            sign = NEG_FLAG;
            ++it;
        }
        while (it != end && *it =='0') ++it;
        digits = end - it;
        assert(digits <= LIMBS * (GMP_NUMB_BITS / 4));
        for (size_t i = 0; i < digits; ++i) {
            char c = it[i];
            int v = (c >= '0' && c <= '9') ? (c - '0') : (c >= 'a' && c <= 'f') ? (c - 'a' + 10) : (c >= 'A' && c <= 'F') ? (c - 'A' + 10) : -1;
            assert(v >= 0 && v <= 15);
            size_t bitpos = (digits - 1 - i) * 4;
            data[bitpos / GMP_NUMB_BITS] |= ((mp_limb_t)v) << (bitpos % GMP_NUMB_BITS);
        }
        if (digits == 0) digits = 1;
        size = sign | (((digits * 4) + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS);
    }

    Num(const std::string& str) : Num(str.begin(), str.end()) {}

    constexpr Num(const char* c) : Num(c, c + my_strlen(c)) {}

    std::string to_string() const {
        static_assert(GMP_NUMB_BITS % 4 == 0);
        std::string ret;
        size_t digits = size_t(size & ABS_MASK) * (GMP_NUMB_BITS / 4);
        for (size_t i = 0; i < digits; ++i) {
            size_t bitpos = (digits - 1 - i) * 4;
            unsigned v = (data[bitpos / GMP_NUMB_BITS] >> (bitpos % GMP_NUMB_BITS)) & 0xf;
            if (ret.size() == 0 && v == 0 && i + 1 < digits) continue;
            static constexpr char HEX[16] = {'0', '1', '2', '3', '4', '5', '6', '7', '8', '9', 'a', 'b', 'c', 'd', 'e', 'f'};
            ret += HEX[v];
        }
        if (SIGNED && (size & NEG_FLAG)) return '-' + ret;
        return ret;
    }

    bool is_zero() const { return size == 1 && data[0] == 0; }
    bool is_odd() const { return data[0] & 1; }

    template<uint32_t LIMBS1, bool SIGNED1>
    Num& operator=(const Num<LIMBS1, SIGNED1>& in) {
        assert(in.abs_size() <= LIMBS);
        assert(SIGNED || !in.neg_flag());
        if constexpr (LIMBS1 <= 4 || LIMBS <= 4) {
            memcpy(data, in.data, sizeof(mp_limb_t) * (LIMBS1 < LIMBS ? LIMBS1 : LIMBS));
        } else {
            memcpy(data, in.data, sizeof(mp_limb_t) * in.abs_size());
        }
        size = in.size;
        return *this;
    }

    Num& operator=(const Num& in) {
        assert(in.abs_size() <= LIMBS);
        assert(SIGNED || !in.neg_flag());
        if constexpr (LIMBS <= 4) {
            memcpy(data, in.data, sizeof(mp_limb_t) * LIMBS);
        } else {
            memcpy(data, in.data, sizeof(mp_limb_t) * in.abs_size());
        }
        size = in.size;
        return *this;
    }

    Num(const Num& in) {
        operator=(in);
    }

    template<uint32_t LIMBS1, bool SIGNED1>
    Num(const Num<LIMBS1, SIGNED1>& in) {
        operator=(in);
    }

    template<uint32_t LIMBS1, bool SIGNED1, uint32_t LIMBS2, bool SIGNED2>
    void set_sum(const Num<LIMBS1, SIGNED1>& in1, const Num<LIMBS2, SIGNED2>& in2) {
        internal_set_sum<SIGNED1, SIGNED2>(in1.data, in1.abs_size(), in1.neg_flag(), in2.data, in2.abs_size(), in2.neg_flag());
    }

    template<uint32_t LIMBS1, bool SIGNED1, uint32_t LIMBS2, bool SIGNED2>
    void set_difference(const Num<LIMBS1, SIGNED1>& in1, const Num<LIMBS2, SIGNED2>& in2) {
        internal_set_sum<SIGNED1, true>(in1.data, in1.abs_size(), in1.neg_flag(), in2.data, in2.abs_size(), in2.neg_flag() ^ NEG_FLAG);
    }

    template<uint32_t LIMBS1, bool SIGNED1>
    void set_square(const Num<LIMBS1, SIGNED1>& in1) {
        if constexpr (LIMBS1 == LIMBS && SIGNED1 == SIGNED) assert(&in1 != this);
        mpn_sqr(data, in1.data, in1.abs_size());
        size = in1.abs_size() * 2;
        normalize();
    }

    template<uint32_t LIMBS1, bool SIGNED1, uint32_t LIMBS2, bool SIGNED2>
    void set_product(const Num<LIMBS1, SIGNED1>& in1, const Num<LIMBS2, SIGNED2>& in2) {
        if constexpr (LIMBS1 == LIMBS && SIGNED1 == SIGNED) assert(&in1 != this);
        if constexpr (LIMBS2 == LIMBS && SIGNED2 == SIGNED) assert(&in2 != this);
        mpn_mul(data, in1.data, in1.abs_size(), in2.data, in2.abs_size());
        size = in1.abs_size() + in2.abs_size();
        if constexpr (SIGNED && (SIGNED1 || SIGNED2)) {
            size |= in1.is_neg() ^ in2.is_neg();
        } else {
            assert(in1.is_neg() == in2.is_neg());
        }
        normalize();
    }

    template<uint32_t LIMBS1, bool SIGNED1>
    void set_lshift(const Num<LIMBS1, SIGNED1>& in1, unsigned shift) {
        if constexpr (LIMBS == LIMBS1 && SIGNED == SIGNED1) assert(&in1 != this);
        internal_set_lshift<false>(in1, shift);
    }

    template<uint32_t LIMBS1, bool SIGNED1>
    void set_rshift(const Num<LIMBS1, SIGNED1>& in1, unsigned shift) {
        if constexpr (LIMBS == LIMBS1 && SIGNED == SIGNED1) assert(&in1 != this);
        internal_set_rshift<false>(in1, shift);
    }

    void set_zero() {
        size = 1;
        data[0] = 0;
    }

    template<uint32_t LIMBS1, bool SIGNED1, uint32_t LIMBS2, bool SIGNED2>
    void set_max_abs(const Num<LIMBS1, SIGNED1>& in1, const Num<LIMBS2, SIGNED2>& in2) {
        uint32_t a1 = in1.abs_size();
        uint32_t a2 = in2.abs_size();
        if (a1 > a2 || ((a1 == a2) && mpn_cmp(in1.data, in2.data, a1) >= 0)) {
            assert(a1 <= LIMBS);
            size = a1;
            memcpy(data, in1.data, sizeof(mp_limb_t) * a1);
        } else {
            assert(a2 <= LIMBS);
            size = a2;
            memcpy(data, in2.data, sizeof(mp_limb_t) * a2);
        }
    }

    template<uint32_t LIMBS1, bool SIGNED1>
    Num& operator+=(const Num<LIMBS1, SIGNED1>& in1) {
        set_sum(*this, in1);
        return *this;
    }

    template<uint32_t LIMBS1, bool SIGNED1>
    Num& operator-=(const Num<LIMBS1, SIGNED1>& in1) {
        set_difference(*this, in1);
        return *this;
    }

    Num& operator>>=(unsigned shift) {
        internal_set_rshift<true>(*this, shift);
        return *this;
    }

    Num& operator<<=(unsigned shift) {
        internal_set_lshift<true>(*this, shift);
        return *this;
    }

    template<uint32_t LIMBS1, bool SIGNED1>
    int compare(const Num<LIMBS1, SIGNED1>& in) const {
        int32_t ssize = signed_size();
        int32_t sinsize = in.signed_size();
        if (ssize < sinsize) return -1;
        if (ssize > sinsize) return 1;
        int ret = mpn_cmp(data, in.data, abs_size());
        return (SIGNED && SIGNED1) ? (ssize < 0 ? -ret : ret) : ret;
    }

    void set_bit(uint32_t pos, bool v) {
        assert(pos / GMP_LIMB_BITS < LIMBS);
        if (v) {
            while (abs_size() <= pos / GMP_LIMB_BITS) {
                data[abs_size()] = 0;
                size = (abs_size() + 1) | neg_flag();
            }
            data[pos / GMP_LIMB_BITS] |= ((mp_limb_t)1) << (pos % GMP_LIMB_BITS);
        } else {
            data[pos / GMP_LIMB_BITS] &= ~(((mp_limb_t)1) << (pos % GMP_LIMB_BITS));
            normalize();
        }
    }

    template<uint32_t LIMBS1, bool SIGNED1> bool operator<(const Num<LIMBS1, SIGNED1>& in) const { return compare(in) < 0; }
    template<uint32_t LIMBS1, bool SIGNED1> bool operator>(const Num<LIMBS1, SIGNED1>& in) const { return compare(in) > 0; }
    template<uint32_t LIMBS1, bool SIGNED1> bool operator<=(const Num<LIMBS1, SIGNED1>& in) const { return compare(in) <= 0; }
    template<uint32_t LIMBS1, bool SIGNED1> bool operator>=(const Num<LIMBS1, SIGNED1>& in) const { return compare(in) >= 0; }
    template<uint32_t LIMBS1, bool SIGNED1> bool operator==(const Num<LIMBS1, SIGNED1>& in) const { return compare(in) == 0; }
    template<uint32_t LIMBS1, bool SIGNED1> bool operator!=(const Num<LIMBS1, SIGNED1>& in) const { return compare(in) != 0; }
};

template<uint32_t Bits>
using Int = Num<(Bits + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS, true>;

template<uint32_t Bits>
using UInt = Num<(Bits + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS, false>;


template<uint32_t Bits>
struct Range {
    Int<Bits> cte;
    Int<Bits> var;

    bool known_parity(int bits) const { return bits == 0 || !var.is_odd(); }
    bool is_odd(int bits) const { assert(known_parity(bits)); return cte.is_odd(); }
    bool is_zero(int bits) const { return cte.is_zero() && (bits == 0 || var.is_zero()); }

    template<uint32_t BITS1>
    Range& operator+=(const Range<BITS1>& in) { cte += in.cte; var += in.var; return *this; }

    template<uint32_t BITS1>
    Range& operator-=(const Range<BITS1>& in) { cte -= in.cte; var -= in.var; return *this; }

    Range& operator>>=(int shift) { cte >>= shift; var >>= shift; return *this; }
    Range& operator<<=(int shift) { cte >>= shift; var >>= shift; return *this; }

    template<uint32_t BITS1, uint32_t BITS2>
    void set_sum(const Range<BITS1>& in1, const Range<BITS2>& in2) { cte.set_sum(in1.cte, in2.cte); var.set_sum(in1.var, in2.var); }

    template<uint32_t BITS1, uint32_t BITS2>
    void set_difference(const Range<BITS1>& in1, const Range<BITS2>& in2) { cte.set_difference(in1.cte, in2.cte); var.set_difference(in1.var, in2.var); }

    template<uint32_t BITS1>
    void set_lshift(const Range<BITS1>& in1, unsigned shift) { cte.set_lshift(in1.cte, shift); var.set_lshift(in1.var, shift); }

    template<uint32_t BITS1>
    void set_rshift(const Range<BITS1>& in1, unsigned shift) { cte.set_rshift(in1.cte, shift); var.set_rshift(in1.var, shift); }

    Range() {}

    template<uint32_t BITS1>
    explicit Range(const Int<BITS1>& in) : cte(in), var() {}

    template<uint32_t BITS1>
    explicit Range(const Int<BITS1>& c, const Int<BITS1>& v) : cte(c), var(v) {}

    UInt<Bits> maxabs(int bits) const {
        Int<Bits + 1> s; s.set_lshift(var, bits); 
//        printf("maxabs: var=%s cte=%s s=(var<<%i)=%s", var.to_string().c_str(), cte.to_string().c_str(), bits, s.to_string().c_str());
        s -= var;
//        printf("maxabs: s-=var: %s\n", s.to_string().c_str());
        s += cte;
//        printf("maxabs: s+=cte: %s\n", s.to_string().c_str());
        UInt<Bits> ret; ret.set_max_abs(cte, s);
//        printf("maxabs: ret=%s\n", ret.to_string().c_str());
        return ret;
    }

    std::string to_string() const {
        return "range(" + cte.to_string() + "+" + var.to_string() + "*x)";
    }
};

constexpr UInt<515> BOUNDS[] = {{"1"}, {"5"}, {"5"}, {"5"}, {"11"}, {"11"}, {"41"}, {"41"}, {"9d"}, {"b5"}, {"1a5"}, {"21d"}, {"48d"}, {"5ed"}, {"ae5"}, {"e33"}, {"2065"}, {"2065"}, {"5eb5"}, {"6ae1"}, {"dd45"}, {"135f5"}, {"2077d"}, {"3403d"}, {"65519"}, {"95c1d"}, {"148625"}, {"17af6d"}, {"2bbd5d"}, {"4629fd"}, {"77a055"}, {"bfbc45"}, {"1555ef5"}, {"2225ced"}, {"43b9611"}, {"46bccc5"}, {"c3fd39d"}, {"d8734f5"}, {"1d0e4755"}, {"1fd863b5"}, {"5bfa47cd"}, {"61b867cd"}, {"cf0cee6d"}, {"10c4473e5"}, {"23613be2d"}, {"2f4fd6e55"}, {"6c313dc05"}, {"8804597d5"}, {"131a849a15"}, {"151a815f45"}, {"368949c8f5"}, {"3bde3097fd"}, {"6cb22ba28d"}, {"910bcb385d"}, {"142c5107995"}, {"1b4ecfdf7ed"}, {"3d3ca7347dd"}, {"40000000000"}, {"40000000001"}, {"40000000001"}, {"40000000001"}, {"40000000001"}, {"43b9cb09df6"}, {"6d8e0def9e7"}, {"b1380287676"}, {"11eac71e4038"}, {"1cfbb031a742"}, {"2ee246b1837e"}, {"4bd72fcac731"}, {"7aae7fc05826"}, {"c673f34c89d5"}, {"14105c2a2b4e1"}, {"2074b4d6c38bf"}, {"34805fdd738a7"}, {"54ed7db62f8e9"}, {"8961968b03300"}, {"de3b3800e9c83"}, {"1677cae2203d27"}, {"24583ee1daf718"}, {"3acac8b53d3c04"}, {"5f1a87ca78b399"}, {"99d790dee333b4"}, {"f8dbdb38b94a54"}, {"1928f723d37f12e"}, {"28b311b5e537fde"}, {"41d626a1342e96d"}, {"6a7fb3d6d9bfb41"}, {"ac46751f849deef"}, {"116ad41575c1745e"}, {"1c2cb6d19044fe11"}, {"2d9378346725fb78"}, {"49b998543748f42a"}, {"774268181530fcf4"}, {"c0eac2d472235b2c"}, {"1381147622f88642a"}, {"1f8ceed1c8e1e24f8"}, {"330973a5fe7757644"}, {"528f01ad04e35daee"}, {"858c85f7609374400"}, {"d8083937b1f84551e"}, {"15d7584313b3bc8799"}, {"2354b38822bcfac8b0"}, {"3926f041c6fa159b24"}, {"5c7360af13f0ef2cb8"}, {"958cf37e5a3642466a"}, {"f1eab58eb5c71fc2cb"}, {"18754b01bf20a659f02"}, {"27906d329171087b3f9"}, {"4000000000000000000"}, {"67872ca9d39b5518413"}, {"a77835432cbbf6ceed2"}, {"10ee72c277d72773a90a"}, {"1b63837be79b0397e86f"}, {"2c4e00a1d9d4f617b1cd"}, {"47ab1c7900dfdc526105"}, {"73eec0c69d06927bd769"}, {"bb891ac60df7e9f62c4a"}, {"12f5cbf2b1cc27dbbfdb4"}, {"1eab9ff0160962c74b1b0"}, {"319cfcd322753fa030124"}, {"504170a8ad3822a34ebb9"}, {"81d2d35b0e2fb7462ab81"}, {"d2017f75ce2995d31005a"}, {"153b5f6d8be3a00bb5c81f"}, {"225865a2c0cbfdcd817d95"}, {"378ece003a720be897f1d1"}, {"59df2b8880f49b6691774e"}, {"9160fb876bdc5eaf7bfc6e"}, {"eb2b22d4f4f00e6a02c2ad"}, {"17c6a1f29e2ce6049b50fd6"}, {"2675e437b8ccece6133ca3c"}, {"3e36f6ce2e5294c66ef6cd7"}, {"64a3dc8f4dfc4b725551e72"}, {"a2cc46d794dff774ed9acda"}, {"1f8ceed1c8e1e24f7ee57880"}, {"330973a5fe775764387d6c03"}, {"528f01ad04e35daedf9ae04c"}, {"858c85f76093743ff89a49fe"}, {"d8083937b1f84551d474db7b"}, {"15d7584313b3bc8798a791dac"}, {"2354b38822bcfac8af8eeab52"}, {"3926f041c6fa159b23b619c41"}, {"5c7360af13f0ef2cb726bbbf0"}, {"958cf37e5a364246694f06928"}, {"f1eab58eb5c71fc2cacf0545a"}, {"18754b01bf20a659f0159880c9"}, {"27906d329171087b3f8b2d3aa6"}, {"40000000000000000000000000"}, {"67872ca9d39b5518412bed538a"}, {"a77835432cbbf6ceed18c9033e"}, {"10ee72c277d72773a909c429a76"}, {"1b63837be79b0397e86ea8b52b4"}, {"2c4e00a1d9d4f617b1cc0529fa4"}, {"47ab1c7900dfdc52610401ea587"}, {"73eec0c69d06927bd768087c2bd"}, {"bb891ac60df7e9f62c49d8336d5"}, {"12f5cbf2b1cc27dbbfdb38786e4a"}, {"1eab9ff0160962c74b1af90f95ce"}, {"319cfcd322753fa0301232cb31c5"}, {"504170a8ad3822a34ebb8d12ddc6"}, {"81d2d35b0e2fb7462ab80f78aea8"}, {"d2017f75ce2995d310059e4ebeba"}, {"153b5f6d8be3a00bb5c81e4d70485"}, {"225865a2c0cbfdcd817d9422168af"}, {"378ece003a720be897f1d0d8c9b87"}, {"59df2b8880f49b6691774d35f1b3e"}, {"9160fb876bdc5eaf7bfc6d632f374"}, {"eb2b22d4f4f00e6a02c2acc3baebd"}, {"17c6a1f29e2ce6049b50fd50207136"}, {"2675e437b8ccece6133ca3b21ecaa7"}, {"3e36f6ce2e5294c66ef6cd66dfbe38"}, {"64a3dc8f4dfc4b725551e7131154a2"}, {"a2cc46d794dff774ed9acd979904df"}, {"107589a84d0ba5b0c93fb333d1811de"}, {"1a9fecf5b66fed0012bce7ad42ce0ad"}, {"2b119d47e1277bb933282a1dea337c4"}, {"45ab5055d705d176560bf182e607ddf"}, {"70b2db464113f843922512fbcfcea5b"}, {"b64de0d19c97eddceee7f413c66236a"}, {"126e66150dd23d0a754984d284196848"}, {"1dd09a06054c3f3ce03fb7f4caf6ca8b"}, {"303ab0b51c88d6caff52b2f80f627ffa"}, {"4e0451d88be2190a7ea22d28a815c76a"}, {"7e33bb472387893dfb95e1ff50390fc4"}, {"cc25ce97f9dd5d90e1f5b00b6ed6a593"}, {"14a3c06b4138d76bb7e6b812eeb8d1bbc"}, {"2163217dd824dd0ffe26927f560b1aae4"}, {"36020e4dec7e1154751d36deb205f7ccf"}, {"575d610c4ecef21e629e476ae02b8e418"}, {"8d52ce208af3eb22be3baad45f77b0a35"}, {"e49bc1071be8566c8ed8671034b8bc5b3"}, {"171cd82bc4fc3bcb2dc9aeefbc9b092b53"}, {"25633cdf968d90919a53c1a49c18ddb539"}, {"3c7aad63ad71c7f0b2b3c15164c2826c88"}, {"61d52c06fc829967c056620323d9a96b84"}, {"9e41b4ca45c421ecfe2cb4ea948a319f8c"}, {"10000000000000000000000000000000000"}, {"19e1cb2a74e6d546104afb54e25783c9ff6"}, {"29de0d50cb2efdb3bb463240cf4f12e669a"}, {"43b9cb09df5c9dcea42710a69d5231ac888"}, {"6d8e0def9e6c0e5fa1baa2d4ace3e073358"}, {"b13802876753d85ec73014a7e8edb6e3a52"}, {"11eac71e4037f7149841007a9618f73e45c1"}, {"1cfbb031a741a49ef5da021f0af0874ce2e5"}, {"2ee246b1837dfa7d8b12760cdb50ec7ad2fb"}, {"4bd72fcac7309f6eff6ce1e1b92660ab37b7"}, {"7aae7fc058258b1d2c6be43e5734f4a6b6ee"}, {"c673f34c89d4fe80c048cb2cc71216eeae1f"}, {"14105c2a2b4e08a8d3aee344b77143154b2be"}, {"2074b4d6c38bedd18aae03de2ba9de1315178"}, {"34805fdd738a6574c4016793afae55c587b1b"}, {"54ed7db62f8e802ed7207935c1212ece8b43d"}, {"8961968b032ff73605f650885a2bb9ecf110a"}, {"de3b3800e9c82fa25fc7436326e19d3dbfc45"}, {"1677cae2203d26d9a45dd34d7c6cf6ba5dcd80"}, {"24583ee1daf717abdeff1b58cbcdcd9cde24d3"}, {"3acac8b53d3c039a80b0ab30eebaf2dfeb23c7"}, {"5f1a87ca78b398126d43f54081c4d6de71018d"}, {"99d790dee333b3984cf28ec87b2a9a9415ae48"}, {"f8dbdb38b94a5319bbdb359b7ef8dceaaf2ba5"}, {"1928f723d37f12dc955479c4c4552865611be72"}, {"28b311b5e537fddd3b66b365e64137a7ba68ef9"}, {"41d626a1342e96c324fecccf46047747f0f0b3d"}, {"6a7fb3d6d9bfb4004af39eb50b382b07d5e18b6"}, {"ac46751f849deee4cca0a877a8cdf0c697323f1"}, {"116ad41575c1745d9582fc60b981f77b8cfce563"}, {"1c2cb6d19044fe10e48944bef3f3a96beaebeb33"}, {"2d9378346725fb773bb9fd04f1988da67336914e"}, {"49b998543748f429d526134a1065a11c4092c541"}, {"774268181530fcf380fedfd32bdb2a29d62cf726"}, {"c0eac2d472235b2bfd4acbe03d89ffe496526e1d"}, {"1381147622f886429fa88b4a2a0571da5c8a8577c"}, {"1f8ceed1c8e1e24f7ee5787fd40e43f0c806109bd"}, {"330973a5fe775764387d6c02dbb5a9649339ac2e0"}, {"528f01ad04e35daedf9ae04bbae346eec1a36350a"}, {"858c85f76093743ff89a49fd582c6ab8d9322c8d3"}, {"d8083937b1f84551d474db7ac817df339bbdb4734"}, {"15d7584313b3bc8798a791dab80ae3905ea9be3b23"}, {"2354b38822bcfac8af8eeab517ddec28d39bd74b75"}, {"3926f041c6fa159b23b619c40d2e2f16c87a3622b7"}, {"5c7360af13f0ef2cb726bbbef26c24ad4b96e3080e"}, {"958cf37e5a364246694f0692706376d4e0b7960448"}, {"f1eab58eb5c71fc2cacf0545930a09b21e9cd1b7de"}, {"18754b01bf20a659f0159880c8f66a5ae0f39a3bfdc"}, {"27906d329171087b3f8b2d3aa5228c67e2e0d78cdef"}, {"4000000000000000000000000000000000000000000"}, {"67872ca9d39b5518412bed53895e0f27fd7b80f72db"}, {"a77835432cbbf6ceed18c9033d3c4b99a64cbf2a329"}, {"10ee72c277d72773a909c429a7548c6b221efcd71406"}, {"1b63837be79b0397e86ea8b52b38f81ccd5d830d4626"}, {"2c4e00a1d9d4f617b1cc0529fa3b6db8e9463183e9bc"}, {"47ab1c7900dfdc52610401ea5863dcf91703149049b7"}, {"73eec0c69d06927bd768087c2bc21d338b907bbd9c01"}, {"bb891ac60df7e9f62c49d8336d43b1eb4be8bb6fb1d2"}, {"12f5cbf2b1cc27dbbfdb38786e49982acded90d9f977b"}, {"1eab9ff0160962c74b1af90f95cd3d29adbb640e6554b"}, {"319cfcd322753fa0301232cb31c485bbab878d920fde1"}, {"504170a8ad3822a34ebb8d12ddc50c552caf76434c718"}, {"81d2d35b0e2fb7462ab80f78aea7784c545dd909d7d05"}, {"d2017f75ce2995d310059e4ebeb957161ec68e523fb91"}, {"153b5f6d8be3a00bb5c81e4d70484bb3a2d0f2fb3c325a"}, {"225865a2c0cbfdcd817d9422168aee7b3c442517e82767"}, {"378ece003a720be897f1d0d8c9b8674f6ff110bccb60a7"}, {"59df2b8880f49b6691774d35f1b3dae97735fd0ce0595d"}, {"9160fb876bdc5eaf7bfc6d632f37367378934bd557d921"}, {"eb2b22d4f4f00e6a02c2acc3baebcb7fac8f18eea508b4"}, {"17c6a1f29e2ce6049b50fd50207135b79c40630700711f6"}, {"2675e437b8ccece6133ca3b21ecaa6a5056b91c9f0992b4"}, {"3e36f6ce2e5294c66ef6cd66dfbe373aabcae921d011184"}, {"64a3dc8f4dfc4b725551e7131154a195846f9c59d822eee"}, {"a2cc46d794dff774ed9acd979904de9ee9a3be0197f68a9"}, {"107589a84d0ba5b0c93fb333d1811dd1fc3c2cf37f2c8433"}, {"1a9fecf5b66fed0012bce7ad42ce0ac1f57862d6972ac20f"}, {"2b119d47e1277bb933282a1dea337c31a5cc8fc38f0d978a"}, {"45ab5055d705d176560bf182e607ddee33f3958b76e7a053"}, {"70b2db464113f843922512fbcfcea5afabafacc9e1ffe604"}, {"b64de0d19c97eddceee7f413c6623699ccda453563cf9b59"}, {"126e66150dd23d0a754984d2841968471024b1501c42aad2a"}, {"1dd09a06054c3f3ce03fb7f4caf6ca8a758b3dc974be18adf"}, {"303ab0b51c88d6caff52b2f80f627ff925949b87037959a74"}, {"4e0451d88be2190a7ea22d28a815c769722a15deedd69a320"}, {"7e33bb472387893dfb95e1ff50390fc32018426f0ad2155a1"}, {"cc25ce97f9dd5d90e1f5b00b6ed6a5924ce6b0b7c41bd347f"}, {"14a3c06b4138d76bb7e6b812eeb8d1bbb068d8d427ec5af11f"}, {"2163217dd824dd0ffe26927f560b1aae364c8b234a409461e4"}, {"36020e4dec7e1154751d36deb205f7cce6ef6d1ccd33dfde81"}, {"575d610c4ecef21e629e476ae02b8e417aa6f8ec8a0c608750"}, {"8d52ce208af3eb22be3baad45f77b0a34e6f5d2dd12b12fdca"}, {"e49bc1071be8566c8ed8671034b8bc5b21e8d88ad911737e35"}, {"171cd82bc4fc3bcb2dc9aeefbc9b092b52e5b8c2035e2244d95"}, {"25633cdf968d90919a53c1a49c18ddb5382de58111f5f3b5384"}, {"3c7aad63ad71c7f0b2b3c15164c2826c87a7346df77f8468ed6"}, {"61d52c06fc829967c056620323d9a96b83ce68eff6f23545107"}, {"9e41b4ca45c421ecfe2cb4ea948a319f8b835e337bae642e03c"}, {"1000000000000000000000000000000000000000000000000000"}, {"19e1cb2a74e6d546104afb54e25783c9ff5ee03dcb6911af10bc"}, {"29de0d50cb2efdb3bb463240cf4f12e669932fca8ca13317e9e7"}, {"43b9cb09df5c9dcea42710a69d5231ac887bf35c501736a7b7a3"}, {"6d8e0def9e6c0e5fa1baa2d4ace3e07335760c351897aa5a97e8"}, {"b13802876753d85ec73014a7e8edb6e3a518c60fa6ecd887697f"}, {"11eac71e4037f7149841007a9618f73e45c0c524126d869b879fa"}, {"1cfbb031a741a49ef5da021f0af0874ce2e41eef67000ab96880d"}, {"2ee246b1837dfa7d8b12760cdb50ec7ad2fa2edbec7457ee295c4"}, {"4bd72fcac7309f6eff6ce1e1b92660ab37b64367e5dea6f14a2ea"}, {"7aae7fc058258b1d2c6be43e5734f4a6b6ed90399552b7612de4d"}, {"c673f34c89d4fe80c048cb2cc71216eeae1e36483f7820fddff11"}, {"14105c2a2b4e08a8d3aee344b77143154b2bdd90d31c5ebdecad44"}, {"2074b4d6c38bedd18aae03de2ba9de131517764275f41363de8e0e"}, {"34805fdd738a6574c4016793afae55c587b1a3948fee4127c701f9"}, {"54ed7db62f8e802ed7207935c1212ece8b43cbecf0c964e54a82e3"}, {"8961968b032ff73605f650885a2bb9ecf110945fa09d9bd7acd97d"}, {"de3b3800e9c82fa25fc7436326e19d3dbfc442f32d829a0e90be1a"}, {"1677cae2203d26d9a45dd34d7c6cf6ba5dcd7f433816571c80e0cef"}, {"24583ee1daf717abdeff1b58cbcdcd9cde24d2f555f6482b0faeb69"}, {"3acac8b53d3c039a80b0ab30eebaf2dfeb23c63ba9422ce32b003fb"}, {"5f1a87ca78b398126d43f54081c4d6de71018c1c01c47d7c3fa455b"}, {"99d790dee333b3984cf28ec87b2a9a9415ae4727c264acce0acdc82"}, {"f8dbdb38b94a5319bbdb359b7ef8dceaaf2ba487404460c7c1588fe"}, {"1928f723d37f12dc955479c4c4552865611be7167608bbb602b4e75d"}, {"28b311b5e537fddd3b66b365e64137a7ba68ef8065fda2a09d408812"}, {"41d626a1342e96c324fecccf46047747f0f0b3cdfcb210cb89037e1d"}, {"6a7fb3d6d9bfb4004af39eb50b382b07d5e18b5a5cab08389755207e"}, {"ac46751f849deee4cca0a877a8cdf0c697323f0e3c365e24700dd021"}, {"116ad41575c1745d9582fc60b981f77b8cfce562ddb9e8149026f8fd2"}, {"1c2cb6d19044fe10e48944bef3f3a96beaebeb32787ff980f122e7d51"}, {"2d9378346725fb773bb9fd04f1988da67336914d58f3e6d6285059329"}, {"49b998543748f429d526134a1065a11c4092c540710aab4a4f9fe47e8"}, {"774268181530fcf380fedfd32bdb2a29d62cf725d2f862b79da2be8b3"}, {"c0eac2d472235b2bfd4acbe03d89ffe496526e1c0de5669cedf92852a"}, {"1381147622f886429fa88b4a2a0571da5c8a8577bb75a68c7e63d419c7"}, {"1f8ceed1c8e1e24f7ee5787fd40e43f0c806109bc2b4855681e15ce242"}, {"330973a5fe775764387d6c02dbb5a9649339ac2df106f4d1f9f09adbb0"}, {"528f01ad04e35daedf9ae04bbae346eec1a363509fb16bc47b7c7943f5"}, {"858c85f76093743ff89a49fd582c6ab8d9322c8d290251878f5f099fa0"}, {"d8083937b1f84551d474db7ac817df339bbdb47334cf7f7a03c9f7157a"}, {"15d7584313b3bc8798a791dab80ae3905ea9be3b22831821d3f162057e8"}, {"2354b38822bcfac8af8eeab517ddec28d39bd74b744ac4bf727740296f1"}, {"3926f041c6fa159b23b619c40d2e2f16c87a3622b6445cdf8d1fe12a7d4"}, {"5c7360af13f0ef2cb726bbbef26c24ad4b96e3080d7889136537f6c0289"}, {"958cf37e5a364246694f0692706376d4e0b7960447d7ced4e0cdda8dae0"}, {"f1eab58eb5c71fc2cacf0545930a09b21e9cd1b7ddfe11a3b56a6515ecc"}, {"18754b01bf20a659f0159880c8f66a5ae0f39a3bfdbc8d514418325c72d9"}, {"27906d329171087b3f8b2d3aa5228c67e2e0d78cdeeb990b80efc8a74977"}, {"400000000000000000000000000000000000000000000000000000000000"}, {"67872ca9d39b5518412bed53895e0f27fd7b80f72da446bc42ee7a5da6e5"}, {"a77835432cbbf6ceed18c9033d3c4b99a64cbf2a3284cc5fa79a2b18eaf2"}, {"10ee72c277d72773a909c429a7548c6b221efcd71405cda9ede8ab2e401f8"}, {"1b63837be79b0397e86ea8b52b38f81ccd5d830d4625ea96a5f9cbb3e6555"}, {"2c4e00a1d9d4f617b1cc0529fa3b6db8e9463183e9bb3621da5fa8edcefb0"}, {"47ab1c7900dfdc52610401ea5863dcf91703149049b61a6e1e7e5bda85761"}, {"73eec0c69d06927bd768087c2bc21d338b907bbd9c002ae5a2032b9e2c8b6"}, {"bb891ac60df7e9f62c49d8336d43b1eb4be8bb6fb1d15fb8a570e44d51b86"}, {"12f5cbf2b1cc27dbbfdb38786e49982acded90d9f977a9bc528ba75eabf0c0"}, {"1eab9ff0160962c74b1af90f95cd3d29adbb640e6554add84b7931dfa04e29"}, {"319cfcd322753fa0301232cb31c485bbab878d920fde083f77fc437e8f3c5c"}, {"504170a8ad3822a34ebb8d12ddc50c552caf76434c717af7b2b50c0a833ee4"}, {"81d2d35b0e2fb7462ab80f78aea7784c545dd909d7d04d8f7a38366f5b550d"}, {"d2017f75ce2995d310059e4ebeb957161ec68e523fb9049f1c07e3e4c4a3a2"}, {"153b5f6d8be3a00bb5c81e4d70484bb3a2d0f2fb3c32593952a0b8ac4a34218"}, {"225865a2c0cbfdcd817d9422168aee7b3c442517e82766f5eb365f06f00d073"}, {"378ece003a720be897f1d0d8c9b8674f6ff110bccb60a683a42f8654ee6e8f1"}, {"59df2b8880f49b6691774d35f1b3dae97735fd0ce0595c7203833b9bad3ffe8"}, {"9160fb876bdc5eaf7bfc6d632f37367378934bd557d920ac3ebada0ca32e591"}, {"eb2b22d4f4f00e6a02c2acc3baebcb7fac8f18eea508b38cac00fe833c3ee4d"}, {"17c6a1f29e2ce6049b50fd50207135b79c40630700711f5f0fe9156a80459780"}, {"2675e437b8ccece6133ca3b21ecaa6a5056b91c9f0992b3382b37204dbc9ef6e"}, {"3e36f6ce2e5294c66ef6cd66dfbe373aabcae921d0111831f05623f6d8b1c366"}, {"64a3dc8f4dfc4b725551e7131154a195846f9c59d822eed80ad39d707526889d"}, {"a2cc46d794dff774ed9acd979904de9ee9a3be0197f68a827502204564830609"}, {"107589a84d0ba5b0c93fb333d1811dd1fc3c2cf37f2c8432e240df87159e39267"}, {"1a9fecf5b66fed0012bce7ad42ce0ac1f57862d6972ac20e25d5481f48dfa0645"}, {"2b119d47e1277bb933282a1dea337c31a5cc8fc38f0d97891c03740801533a88c"}, {"45ab5055d705d176560bf182e607ddee33f3958b76e7a052409be3f456f89fa49"}, {"70b2db464113f843922512fbcfcea5afabafacc9e1ffe603c48b9f543e6fb1068"}, {"b64de0d19c97eddceee7f413c6623699ccda453563cf9b58a14164ca16de52908"}, {"126e66150dd23d0a754984d2841968471024b1501c42aad293e7f91f9dd71b59b1"}, {"1dd09a06054c3f3ce03fb7f4caf6ca8a758b3dc974be18ade768afa2c9050ff87d"}, {"303ab0b51c88d6caff52b2f80f627ff925949b87037959a73b7e4a14a4282b1829"}, {"4e0451d88be2190a7ea22d28a815c769722a15deedd69a31f98f50671b46e57520"}, {"7e33bb472387893dfb95e1ff50390fc32018426f0ad2155a0785738907eb0b95b1"}, {"cc25ce97f9dd5d90e1f5b00b6ed6a5924ce6b0b7c41bd347e7c26b6ebd1d02b68a"}, {"14a3c06b4138d76bb7e6b812eeb8d1bbb068d8d427ec5af11edf1e50fd2e63ad3ba"}, {"2163217dd824dd0ffe26927f560b1aae364c8b234a409461e3d7c267e7d3b4178e3"}, {"36020e4dec7e1154751d36deb205f7cce6ef6d1ccd33dfde80f27dc55e47ff1e6cf"}, {"575d610c4ecef21e629e476ae02b8e417aa6f8ec8a0c60874fc58815f9ff3a1b5f8"}, {"8d52ce208af3eb22be3baad45f77b0a34e6f5d2dd12b12fdc9dd00a5bc3ae727ff3"}, {"e49bc1071be8566c8ed8671034b8bc5b21e8d88ad911737e347f84a9f4f4194f507"}, {"171cd82bc4fc3bcb2dc9aeefbc9b092b52e5b8c2035e2244d94dfdb00a23d5b8b95a"}, {"25633cdf968d90919a53c1a49c18ddb5382de58111f5f3b5383376a36b7d7e45e7ea"}, {"3c7aad63ad71c7f0b2b3c15164c2826c87a7346df77f8468ed5a99457b2cb4387f45"}, {"61d52c06fc829967c056620323d9a96b83ce68eff6f235451060c971cb60f0de8528"}, {"9e41b4ca45c421ecfe2cb4ea948a319f8b835e337bae642e03bf229d25db63926a80"}, {"100000000000000000000000000000000000000000000000000000000000000000000"}, {"19e1cb2a74e6d546104afb54e25783c9ff5ee03dcb6911af10bb9e9769b9325cd1d9b"}, {"29de0d50cb2efdb3bb463240cf4f12e669932fca8ca13317e9e68ac63abc6b3dc7aef"}, {"43b9cb09df5c9dcea42710a69d5231ac887bf35c501736a7b7a2acb9007dea8c6aae4"}, {"6d8e0def9e6c0e5fa1baa2d4ace3e07335760c351897aa5a97e72ecf9955110735efa"}, {"b13802876753d85ec73014a7e8edb6e3a518c60fa6ecd887697ea3b73bebfbbf3a3d5"}, {"11eac71e4037f7149841007a9618f73e45c0c524126d869b879f96f6a15d80ce52497d"}, {"1cfbb031a741a49ef5da021f0af0874ce2e41eef67000ab96880cae78b22d549d876dd"}, {"2ee246b1837dfa7d8b12760cdb50ec7ad2fa2edbec7457ee295c3913546e15a694edf6"}, {"4bd72fcac7309f6eff6ce1e1b92660ab37b64367e5dea6f14a2e9d7aafc2fe89b6632f"}, {"7aae7fc058258b1d2c6be43e5734f4a6b6ed90399552b7612de4c77e8138a2dd4683ca"}, {"c673f34c89d4fe80c048cb2cc71216eeae1e36483f7820fddff10dfa3cf16e28365980"}, {"14105c2a2b4e08a8d3aee344b77143154b2bdd90d31c5ebdecad4302a0cfb8d16c3cef0"}, {"2074b4d6c38bedd18aae03de2ba9de131517764275f41363de8e0d9bd6d5432db17d15f"}, {"34805fdd738a6574c4016793afae55c587b1a3948fee4127c701f8f93128e8598e47108"}, {"54ed7db62f8e802ed7207935c1212ece8b43cbecf0c964e54a82e2b128d085dfb43842e"}, {"8961968b032ff73605f650885a2bb9ecf110945fa09d9bd7acd97c1bc0341ca19a44736"}, {"de3b3800e9c82fa25fc7436326e19d3dbfc442f32d829a0e90be1953b9ba3c232f1c957"}, {"1677cae2203d26d9a45dd34d7c6cf6ba5dcd7f433816571c80e0cee6eb4fff9e01c46bef"}, {"24583ee1daf717abdeff1b58cbcdcd9cde24d2f555f6482b0faeb68328cb9642fced4776"}, {"3acac8b53d3c039a80b0ab30eebaf2dfeb23c63ba9422ce32b003fa0cf0fb9327ec5723e"}, {"5f1a87ca78b398126d43f54081c4d6de71018c1c01c47d7c3fa455aa01165dfdbb62db72"}, {"99d790dee333b3984cf28ec87b2a9a9415ae4727c264acce0acdc8136f27bdb408408aa3"}, {"f8dbdb38b94a5319bbdb359b7ef8dceaaf2ba487404460c7c1588fdb62c70d95db851afa"}, {"1928f723d37f12dc955479c4c4552865611be7167608bbb602b4e75c1d49a2271a46c5bfd"}, {"28b311b5e537fddd3b66b365e64137a7ba68ef8065fda2a09d4088115920c1823677985f3"}, {"41d626a1342e96c324fecccf46047747f0f0b3cdfcb210cb89037e1c5678e4999a2b5f401"}, {"6a7fb3d6d9bfb4004af39eb50b382b07d5e18b5a5cab08389755207d237e81913fe12f68c"}, {"ac46751f849deee4cca0a877a8cdf0c697323f0e3c365e24700dd020054cea22f46d98038"}, {"116ad41575c1745d9582fc60b981f77b8cfce562ddb9e8149026f8fd15be27e922b4ced069"}, {"1c2cb6d19044fe10e48944bef3f3a96beaebeb32787ff980f122e7d50f9bec419e80ad482b"}, {"2d9378346725fb773bb9fd04f1988da67336914d58f3e6d62850593285b794a41ce108ca95"}, {"49b998543748f429d526134a1065a11c4092c540710aab4a4f9fe47e775c6d66c2d82f6c0b"}, {"774268181530fcf380fedfd32bdb2a29d62cf725d2f862b79da2be8b24143fe1f2c00adf97"}, {"c0eac2d472235b2bfd4acbe03d89ffe496526e1c0de5669cedf9285290a0ac60a106fa5758"}, {"1381147622f886429fa88b4a2a0571da5c8a8577bb75a68c7e63d419c6d1b95d47f78eb682a"}, {"1f8ceed1c8e1e24f7ee5787fd40e43f0c806109bc2b4855681e15ce241fac2e56c0a274425c"}, {"330973a5fe775764387d6c02dbb5a9649339ac2df106f4d1f9f09adbaf4740ada257b0d0bca"}, {"528f01ad04e35daedf9ae04bbae346eec1a363509fb16bc47b7c7943f4b98eb4ee72f13fffb"}, {"858c85f76093743ff89a49fd582c6ab8d9322c8d290251878f5f099f9f4ed05e38b2f5adb46"}, {"d8083937b1f84551d474db7ac817df339bbdb47334cf7f7a03c9f715791ffc79b38b9ea348a"}, {"15d7584313b3bc8798a791dab80ae3905ea9be3b22831821d3f162057e7fce86d7dee6b0d7b4"}, {"2354b38822bcfac8af8eeab517ddec28d39bd74b744ac4bf727740296f0eb9c9ffc88712f3cb"}, {"3926f041c6fa159b23b619c40d2e2f16c87a3622b6445cdf8d1fe12a7d3d0653d418b1852457"}, {"5c7360af13f0ef2cb726bbbef26c24ad4b96e3080d7889136537f6c0288f56e2e567fb774268"}, {"958cf37e5a364246694f0692706376d4e0b7960447d7ced4e0cdda8dadf5f9179fa4e02746d2"}, {"f1eab58eb5c71fc2cacf0545930a09b21e9cd1b7ddfe11a3b56a6515ecb2d0e1fd12d43e64ad"}, {"18754b01bf20a659f0159880c8f66a5ae0f39a3bfdbc8d514418325c72d83c37a149cbf89bc81"}, {"27906d329171087b3f8b2d3aa5228c67e2e0d78cdeeb990b80efc8a74976d8e49a9fe3c0d8961"}, {"40000000000000000000000000000000000000000000000000000000000000000000000000000"}, {"67872ca9d39b5518412bed53895e0f27fd7b80f72da446bc42ee7a5da6e4c9734766ace7d8375"}, {"a77835432cbbf6ceed18c9033d3c4b99a64cbf2a3284cc5fa79a2b18eaf1acf71ebb8dc69f58e"}, {"10ee72c277d72773a909c429a7548c6b221efcd71405cda9ede8ab2e401f7aa31aab8ffb46ffa5"}, {"1b63837be79b0397e86ea8b52b38f81ccd5d830d4625ea96a5f9cbb3e6554441cd7be7f2a85baa"}, {"2c4e00a1d9d4f617b1cc0529fa3b6db8e9463183e9bb3621da5fa8edcefafeefce8f509fa16fb7"}, {"47ab1c7900dfdc52610401ea5863dcf91703149049b61a6e1e7e5bda857603394925f1781e6e8c"}, {"73eec0c69d06927bd768087c2bc21d338b907bbd9c002ae5a2032b9e2c8b552761db707d783df6"}, {"bb891ac60df7e9f62c49d8336d43b1eb4be8bb6fb1d15fb8a570e44d51b8569a53b7d6dca47cc0"}, {"12f5cbf2b1cc27dbbfdb38786e49982acded90d9f977a9bc528ba75eabf0bfa26d98cbbcc3ad43b"}, {"1eab9ff0160962c74b1af90f95cd3d29adbb640e6554add84b7931dfa04e28b751a0f266a1b33e6"}, {"319cfcd322753fa0301232cb31c485bbab878d920fde083f77fc437e8f3c5b8a0d965fc02ff3b8c"}, {"504170a8ad3822a34ebb8d12ddc50c552caf76434c717af7b2b50c0a833ee345b0f3bbf70b08a05"}, {"81d2d35b0e2fb7462ab80f78aea7784c545dd909d7d04d8f7a38366f5b550cb6c5f457a1cba80ae"}, {"d2017f75ce2995d310059e4ebeb957161ec68e523fb9049f1c07e3e4c4a3a166391c41e1a93deb8"}, {"153b5f6d8be3a00bb5c81e4d70484bb3a2d0f2fb3c32593952a0b8ac4a342177ed0e10b74f1f8ccd"}, {"225865a2c0cbfdcd817d9422168aee7b3c442517e82766f5eb365f06f00d072866911cd5b353680c"}, {"378ece003a720be897f1d0d8c9b8674f6ff110bccb60a683a42f8654ee6e8f08cbc7255bd4789f5e"}, {"59df2b8880f49b6691774d35f1b3dae97735fd0ce0595c7203833b9bad3ffe780711afbb9f1ced05"}, {"9160fb876bdc5eaf7bfc6d632f37367378934bd557d920ac3ebada0ca32e590bf3b51dd4093e7613"}, {"eb2b22d4f4f00e6a02c2acc3baebcb7fac8f18eea508b38cac00fe833c3ee4c9fb15c8f7d9b06c15"}, {"17c6a1f29e2ce6049b50fd50207135b79c40630700711f5f0fe9156a8045977f6ed8b6dc67eb46506"}, {"2675e437b8ccece6133ca3b21ecaa6a5056b91c9f0992b3382b37204dbc9ef6d021022a89b4f4062e"}, {"3e36f6ce2e5294c66ef6cd66dfbe373aabcae921d0111831f05623f6d8b1c36576e146be61bfd1c5b"}, {"64a3dc8f4dfc4b725551e7131154a195846f9c59d822eed80ad39d707526889c691b16ff01d86e40d"}, {"a2cc46d794dff774ed9acd979904de9ee9a3be0197f68a827502204564830608d9de617c8d976e6ae"}, {"107589a84d0ba5b0c93fb333d1811dd1fc3c2cf37f2c8432e240df87159e3926668ad7d001fcc42fbd"}, {"1a9fecf5b66fed0012bce7ad42ce0ac1f57862d6972ac20e25d5481f48dfa0644ff84bda2cbc861d2c"}, {"2b119d47e1277bb933282a1dea337c31a5cc8fc38f0d97891c03740801533a88bd1b6600dc146c0cda"}, {"45ab5055d705d176560bf182e607ddee33f3958b76e7a052409be3f456f89fa48ad33b41a0d0ad73f5"}, {"70b2db464113f843922512fbcfcea5afabafacc9e1ffe603c48b9f543e6fb1067a02b520ab9f6bb1e1"}, {"b64de0d19c97eddceee7f413c6623699ccda453563cf9b58a14164ca16de52907384232a52f3dd3c2e"}, {"126e66150dd23d0a754984d2841968471024b1501c42aad293e7f91f9dd71b59b0b60bdb02b78104288"}, {"1dd09a06054c3f3ce03fb7f4caf6ca8a758b3dc974be18ade768afa2c9050ff87cb002b7e5959b5687c"}, {"303ab0b51c88d6caff52b2f80f627ff925949b87037959a73b7e4a14a4282b182841be95d5c5553cb7a"}, {"4e0451d88be2190a7ea22d28a815c769722a15deedd69a31f98f50671b46e5751fde3ada0a6721c3557"}, {"7e33bb472387893dfb95e1ff50390fc32018426f0ad2155a0785738907eb0b95b0289d1096f48198809"}, {"cc25ce97f9dd5d90e1f5b00b6ed6a5924ce6b0b7c41bd347e7c26b6ebd1d02b6895ec342f25eec1fc58"}, {"14a3c06b4138d76bb7e6b812eeb8d1bbb068d8d427ec5af11edf1e50fd2e63ad3b9cbc4fffea6b1186e1"}, {"2163217dd824dd0ffe26927f560b1aae364c8b234a409461e3d7c267e7d3b4178e2cbd6b6d15b99cc10b"}, {"36020e4dec7e1154751d36deb205f7cce6ef6d1ccd33dfde80f27dc55e47ff1e6ce2e7a8d22687b984d8"}, {"575d610c4ecef21e629e476ae02b8e417aa6f8ec8a0c60874fc58815f9ff3a1b5f7b9ac35ecfedefbb99"}, {"8d52ce208af3eb22be3baad45f77b0a34e6f5d2dd12b12fdc9dd00a5bc3ae727ff221c4bcf2ad73e71b8"}, {"e49bc1071be8566c8ed8671034b8bc5b21e8d88ad911737e347f84a9f4f4194f5062c614915953630e31"}, {"171cd82bc4fc3bcb2dc9aeefbc9b092b52e5b8c2035e2244d94dfdb00a23d5b8b959feddd099c85ff3d9b"}, {"25633cdf968d90919a53c1a49c18ddb5382de58111f5f3b5383376a36b7d7e45e7e93809d1b477d403f04"}, {"3c7aad63ad71c7f0b2b3c15164c2826c87a7346df77f8468ed5a99457b2cb4387f44b50f992b3cbc5e338"}, {"61d52c06fc829967c056620323d9a96b83ce68eff6f235451060c971cb60f0de85272fe26f20058ad9b29"}, {"9e41b4ca45c421ecfe2cb4ea948a319f8b835e337bae642e03bf229d25db63926a7f8f0362583e0fdaf23"}, {"10000000000000000000000000000000000000000000000000000000000000000000000000000000000000"}, {"19e1cb2a74e6d546104afb54e25783c9ff5ee03dcb6911af10bb9e9769b9325cd1d9ab39f60dd00fb0a267"}, {"29de0d50cb2efdb3bb463240cf4f12e669932fca8ca13317e9e68ac63abc6b3dc7aee371a7d634f7ced317"}, {"43b9cb09df5c9dcea42710a69d5231ac887bf35c501736a7b7a2acb9007dea8c6aae3fed1bfe919f98be9f"}, {"6d8e0def9e6c0e5fa1baa2d4ace3e07335760c351897aa5a97e72ecf9955110735ef9fcaa16ea7d54e7db9"}, {"b13802876753d85ec73014a7e8edb6e3a518c60fa6ecd887697ea3b73bebfbbf3a3d427e85beda7d360721"}, {"11eac71e4037f7149841007a9618f73e45c0c524126d869b879f96f6a15d80ce52497c5e079ba2fcf733e75"}, {"1cfbb031a741a49ef5da021f0af0874ce2e41eef67000ab96880cae78b22d549d876dc1f5e0f7d581dee036"}, {"2ee246b1837dfa7d8b12760cdb50ec7ad2fa2edbec7457ee295c3913546e15a694edf5b7291f2fd01ad501a"}, {"4bd72fcac7309f6eff6ce1e1b92660ab37b64367e5dea6f14a2e9d7aafc2fe89b6632ef30eb50eadf34c6c8"}, {"7aae7fc058258b1d2c6be43e5734f4a6b6ed90399552b7612de4c77e8138a2dd4683c99a86ccf97782b35dc"}, {"c673f34c89d4fe80c048cb2cc71216eeae1e36483f7820fddff10dfa3cf16e2836597f00bfcee2cacf5b61b"}, {"14105c2a2b4e08a8d3aee344b77143154b2bdd90d31c5ebdecad4302a0cfb8d16c3ceefdc2c2281201b3af55"}, {"2074b4d6c38bedd18aae03de2ba9de131517764275f41363de8e0d9bd6d5432db17d15e872ea02b4901069d2"}, {"34805fdd738a6574c4016793afae55c587b1a3948fee4127c701f8f93128e8598e4710786a4f7adc92b192e9"}, {"54ed7db62f8e802ed7207935c1212ece8b43cbecf0c964e54a82e2b128d085dfb43842dd3c7e33310ba1ca01"}, {"8961968b032ff73605f650885a2bb9ecf110945fa09d9bd7acd97c1bc0341ca19a447356cd4da02fac68def4"}, {"de3b3800e9c82fa25fc7436326e19d3dbfc442f32d829a0e90be1953b9ba3c232f1c956f51e27d75106c1555"}, {"1677cae2203d26d9a45dd34d7c6cf6ba5dcd7f433816571c80e0cee6eb4fff9e01c46beee7c73b410beb624f6"}, {"24583ee1daf717abdeff1b58cbcdcd9cde24d2f555f6482b0faeb68328cb9642fced4775024f9d84a6faae2cf"}, {"3acac8b53d3c039a80b0ab30eebaf2dfeb23c63ba9422ce32b003fa0cf0fb9327ec5723df66c1b05268590edf"}, {"5f1a87ca78b398126d43f54081c4d6de71018c1c01c47d7c3fa455aa01165dfdbb62db719fad19416e1c9c727"}, {"99d790dee333b3984cf28ec87b2a9a9415ae4727c264acce0acdc8136f27bdb408408aa26d3d018b75f27d3cb"}, {"f8dbdb38b94a5319bbdb359b7ef8dceaaf2ba487404460c7c1588fdb62c70d95db851af986ff4716aff2e0573"}, {"1928f723d37f12dc955479c4c4552865611be7167608bbb602b4e75c1d49a2271a46c5bfc0761b903059b4eeb4"}, {"28b311b5e537fddd3b66b365e64137a7ba68ef8065fda2a09d4088115920c1823677985f2365db9ab5f2ee5308"}, {"41d626a1342e96c324fecccf46047747f0f0b3cdfcb210cb89037e1c5678e4999a2b5f4007f310bef0e3c25fdf"}, {"6a7fb3d6d9bfb4004af39eb50b382b07d5e18b5a5cab08389755207d237e81913fe12f68b2f21874ad6d4c90e9"}, {"ac46751f849deee4cca0a877a8cdf0c697323f0e3c365e24700dd020054cea22f46d98037051b03366c2fddaed"}, {"116ad41575c1745d9582fc60b981f77b8cfce562ddb9e8149026f8fd15be27e922b4ced068342b5cfd2d701a49d"}, {"1c2cb6d19044fe10e48944bef3f3a96beaebeb32787ff980f122e7d50f9bec419e80ad482ae7daec78360a03f9d"}, {"2d9378346725fb773bb9fd04f1988da67336914d58f3e6d62850593285b794a41ce108ca94bcf74f0b4b98321cb"}, {"49b998543748f429d526134a1065a11c4092c540710aab4a4f9fe47e775c6d66c2d82f6c0ade0410a1d8525c4bb"}, {"774268181530fcf380fedfd32bdb2a29d62cf725d2f862b79da2be8b24143fe1f2c00adf96566d5a1ecf14b6375"}, {"c0eac2d472235b2bfd4acbe03d89ffe496526e1c0de5669cedf9285290a0ac60a106fa57571554f2de46d751c32"}, {"1381147622f886429fa88b4a2a0571da5c8a8577bb75a68c7e63d419c6d1b95d47f78eb68299c870d55b3c9f83e4"}, {"1f8ceed1c8e1e24f7ee5787fd40e43f0c806109bc2b4855681e15ce241fac2e56c0a274425bd2066202101724689"}, {"330973a5fe775764387d6c02dbb5a9649339ac2df106f4d1f9f09adbaf4740ada257b0d0bc97bb07f15f35a586fc"}, {"528f01ad04e35daedf9ae04bbae346eec1a363509fb16bc47b7c7943f4b98eb4ee72f13fffa9ac461b8384f5673d"}, {"858c85f76093743ff89a49fd582c6ab8d9322c8d290251878f5f099f9f4ed05e38b2f5adb456e6730429f36257a9"}, {"d8083937b1f84551d474db7ac817df339bbdb47334cf7f7a03c9f715791ffc79b38b9ea3489a1ee6135e925967d2"}, {"15d7584313b3bc8798a791dab80ae3905ea9be3b22831821d3f162057e7fce86d7dee6b0d7b3fb7beee60a5623b8a"}, {"2354b38822bcfac8af8eeab517ddec28d39bd74b744ac4bf727740296f0eb9c9ffc88712f3cab5cf9c6dcdf087130"}, {"3926f041c6fa159b23b619c40d2e2f16c87a3622b6445cdf8d1fe12a7d3d0653d418b185245654d8c38c3443926ea"}, {"5c7360af13f0ef2cb726bbbef26c24ad4b96e3080d7889136537f6c0288f56e2e567fb774267217fcf66976c30e66"}, {"958cf37e5a364246694f0692706376d4e0b7960447d7ced4e0cdda8dadf5f9179fa4e02746d1df500fc0d63441044"}, {"f1eab58eb5c71fc2cacf0545930a09b21e9cd1b7ddfe11a3b56a6515ecb2d0e1fd12d43e64acf2f178cdd4a44891a"}, {"18754b01bf20a659f0159880c8f66a5ae0f39a3bfdbc8d514418325c72d83c37a149cbf89bc80162b66ca191c3282d"}, {"27906d329171087b3f8b2d3aa5228c67e2e0d78cdeeb990b80efc8a74976d8e49a9fe3c0d8960f83f6bc8bd54d0ed6"}, {"4000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"}, {"67872ca9d39b5518412bed53895e0f27fd7b80f72da446bc42ee7a5da6e4c9734766ace7d837403ec2899bf90387a2"}, {"a77835432cbbf6ceed18c9033d3c4b99a64cbf2a3284cc5fa79a2b18eaf1acf71ebb8dc69f58d3df3b4c5af7cff52b"}, {"10ee72c277d72773a909c429a7548c6b221efcd71405cda9ede8ab2e401f7aa31aab8ffb46ffa467e62fa791c0b64b1"}, {"1b63837be79b0397e86ea8b52b38f81ccd5d830d4625ea96a5f9cbb3e6554441cd7be7f2a85ba9f5539f6e2c50e72ef"}, {"2c4e00a1d9d4f617b1cc0529fa3b6db8e9463183e9bb3621da5fa8edcefafeefce8f509fa16fb69f4d81c834cbdeb03"}, {"47ab1c7900dfdc52610401ea5863dcf91703149049b61a6e1e7e5bda857603394925f1781e6e8bf3dccf9d132d04fb0"}, {"73eec0c69d06927bd768087c2bc21d338b907bbd9c002ae5a2032b9e2c8b552761db707d783df56077b80d44217a5a4"}, {"bb891ac60df7e9f62c49d8336d43b1eb4be8bb6fb1d15fb8a570e44d51b8569a53b7d6dca47cbf406b54067d0b85f7e"}, {"12f5cbf2b1cc27dbbfdb38786e49982acded90d9f977a9bc528ba75eabf0bfa26d98cbbcc3ad43ab7cd31b1fe76efc71"}, {"1eab9ff0160962c74b1af90f95cd3d29adbb640e6554add84b7931dfa04e28b751a0f266a1b33e5de0acd76ca8ac4a4a"}, {"319cfcd322753fa0301232cb31c485bbab878d920fde083f77fc437e8f3c5b8a0d965fc02ff3b8b2b3d6d86ab61797ee"}, {"504170a8ad3822a34ebb8d12ddc50c552caf76434c717af7b2b50c0a833ee345b0f3bbf70b08a04806cebd52f6056b3a"}, {"81d2d35b0e2fb7462ab80f78aea7784c545dd909d7d04d8f7a38366f5b550cb6c5f457a1cba80ad24041a7447b133765"}, {"d2017f75ce2995d310059e4ebeb957161ec68e523fb9049f1c07e3e4c4a3a166391c41e1a93deb724ac64ba1d265880f"}, {"153b5f6d8be3a00bb5c81e4d70484bb3a2d0f2fb3c32593952a0b8ac4a342177ed0e10b74f1f8ccc42e8728039c3d7d29"}, {"225865a2c0cbfdcd817d9422168aee7b3c442517e82766f5eb365f06f00d072866911cd5b353680beb1a37bcdc42cf4a4"}, {"378ece003a720be897f1d0d8c9b8674f6ff110bccb60a683a42f8654ee6e8f08cbc7255bd4789f5d441b055529a992d40"}, {"59df2b8880f49b6691774d35f1b3dae97735fd0ce0595c7203833b9bad3ffe780711afbb9f1ced042fad893d56a879d65"}, {"9160fb876bdc5eaf7bfc6d632f37367378934bd557d920ac3ebada0ca32e590bf3b51dd4093e76129beab8b3aaf56965c"}, {"eb2b22d4f4f00e6a02c2acc3baebcb7fac8f18eea508b38cac00fe833c3ee4c9fb15c8f7d9b06c149a1643b78ca4e88ee"}, {"17c6a1f29e2ce6049b50fd50207135b79c40630700711f5f0fe9156a8045977f6ed8b6dc67eb46505b87271c99ab600429"}, {"2675e437b8ccece6133ca3b21ecaa6a5056b91c9f0992b3382b37204dbc9ef6d021022a89b4f4062dd7c9f4f2b19ee72af"}, {"3e36f6ce2e5294c66ef6cd66dfbe373aabcae921d0111831f05623f6d8b1c36576e146be61bfd1c5abfcb815ca50caaddd"}, {"64a3dc8f4dfc4b725551e7131154a195846f9c59d822eed80ad39d707526889c691b16ff01d86e40c166d3baccbbbc3361"}, {"a2cc46d794dff774ed9acd979904de9ee9a3be0197f68a827502204564830608d9de617c8d976e6ad7cbb94c1fe972528b"}, {"107589a84d0ba5b0c93fb333d1811dd1fc3c2cf37f2c8432e240df87159e3926668ad7d001fcc42fbc38f097f78889179b2"}, {"1a9fecf5b66fed0012bce7ad42ce0ac1f57862d6972ac20e25d5481f48dfa0644ff84bda2cbc861d2b5b53243a0f9d9e6de"}, {"2b119d47e1277bb933282a1dea337c31a5cc8fc38f0d97891c03740801533a88bd1b6600dc146c0cd9b0bf76bb17be68b56"}, {"45ab5055d705d176560bf182e607ddee33f3958b76e7a052409be3f456f89fa48ad33b41a0d0ad73f4b5c06927339bce181"}, {"70b2db464113f843922512fbcfcea5afabafacc9e1ffe603c48b9f543e6fb1067a02b520ab9f6bb1e0d8280fe737ad55dc8"}, {"b64de0d19c97eddceee7f413c6623699ccda453563cf9b58a14164ca16de52907384232a52f3dd3c2d2e60c872b1b13d80a"}, {"126e66150dd23d0a754984d2841968471024b1501c42aad293e7f91f9dd71b59b0b60bdb02b781042876149712e98e20c335"}, {"1dd09a06054c3f3ce03fb7f4caf6ca8a758b3dc974be18ade768afa2c9050ff87cb002b7e5959b5687b3c52d8dd270c74f5a"}, {"303ab0b51c88d6caff52b2f80f627ff925949b87037959a73b7e4a14a4282b182841be95d5c5553cb791b5d470c6ed800390"}, {"4e0451d88be2190a7ea22d28a815c769722a15deedd69a31f98f50671b46e5751fde3ada0a6721c3556cf27e0f8cd3e4c250"}, {"7e33bb472387893dfb95e1ff50390fc32018426f0ad2155a0785738907eb0b95b0289d1096f48198808405c91a234619b18b"}, {"cc25ce97f9dd5d90e1f5b00b6ed6a5924ce6b0b7c41bd347e7c26b6ebd1d02b6895ec342f25eec1fc57cd6961bef9f2be545"}, {"14a3c06b4138d76bb7e6b812eeb8d1bbb068d8d427ec5af11edf1e50fd2e63ad3b9cbc4fffea6b1186e0e13d59cf2e1f4ce2f"}, {"2163217dd824dd0ffe26927f560b1aae364c8b234a409461e3d7c267e7d3b4178e2cbd6b6d15b99cc10a7cd895ea11558d5ac"}, {"36020e4dec7e1154751d36deb205f7cce6ef6d1ccd33dfde80f27dc55e47ff1e6ce2e7a8d22687b984d7a49659f46751dfd3e"}, {"575d610c4ecef21e629e476ae02b8e417aa6f8ec8a0c60874fc58815f9ff3a1b5f7b9ac35ecfedefbb9829588ee248f719f51"}, {"8d52ce208af3eb22be3baad45f77b0a34e6f5d2dd12b12fdc9dd00a5bc3ae727ff221c4bcf2ad73e71b737c21c4bd1f7aa802"}, {"e49bc1071be8566c8ed8671034b8bc5b21e8d88ad911737e347f84a9f4f4194f5062c614915953630e30d10e49ba6b695bf7c"}, {"171cd82bc4fc3bcb2dc9aeefbc9b092b52e5b8c2035e2244d94dfdb00a23d5b8b959feddd099c85ff3d9a5db0c3997e427337a"}, {"25633cdf968d90919a53c1a49c18ddb5382de58111f5f3b5383376a36b7d7e45e7e93809d1b477d403f0358d10410dfdafc6de"}, {"3c7aad63ad71c7f0b2b3c15164c2826c87a7346df77f8468ed5a99457b2cb4387f44b50f992b3cbc5e3375291224644d924d15"}, {"61d52c06fc829967c056620323d9a96b83ce68eff6f235451060c971cb60f0de85272fe26f20058ad9b286470ca0b35c4b38c8"}, {"9e41b4ca45c421ecfe2cb4ea948a319f8b835e337bae642e03bf229d25db63926a7f8f0362583e0fdaf22f55343b5635c7edcb"}, {"1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"}, {"19e1cb2a74e6d546104afb54e25783c9ff5ee03dcb6911af10bb9e9769b9325cd1d9ab39f60dd00fb0a266fe40e1e85ca2487c9"}, {"29de0d50cb2efdb3bb463240cf4f12e669932fca8ca13317e9e68ac63abc6b3dc7aee371a7d634f7ced316bdf3fd4aba94866a6"}, {"43b9cb09df5c9dcea42710a69d5231ac887bf35c501736a7b7a2acb9007dea8c6aae3fed1bfe919f98be9e4702d92c0960755c5"}, {"6d8e0def9e6c0e5fa1baa2d4ace3e07335760c351897aa5a97e72ecf9955110735ef9fcaa16ea7d54e7db8b1439cbba73327729"}, {"b13802876753d85ec73014a7e8edb6e3a518c60fa6ecd887697ea3b73bebfbbf3a3d427e85beda7d360720d32f7ac095033fee6"}, {"11eac71e4037f7149841007a9618f73e45c0c524126d869b879f96f6a15d80ce52497c5e079ba2fcf733e744cb413ebf780a6e0e"}, {"1cfbb031a741a49ef5da021f0af0874ce2e41eef67000ab96880cae78b22d549d876dc1f5e0f7d581dee0351085e968f52b322c3"}, {"2ee246b1837dfa7d8b12760cdb50ec7ad2fa2edbec7457ee295c3913546e15a694edf5b7291f2fd01ad5019f42e17df4ef2ca3b9"}, {"4bd72fcac7309f6eff6ce1e1b92660ab37b64367e5dea6f14a2e9d7aafc2fe89b6632ef30eb50eadf34c6c7f9dbbf1c20f7d1a0b"}, {"7aae7fc058258b1d2c6be43e5734f4a6b6ed90399552b7612de4c77e8138a2dd4683c99a86ccf97782b35db2a2b129245cc16295"}, {"c673f34c89d4fe80c048cb2cc71216eeae1e36483f7820fddff10dfa3cf16e2836597f00bfcee2cacf5b61aad85e5fb59fd5b7b3"}, {"14105c2a2b4e08a8d3aee344b77143154b2bdd90d31c5ebdecad4302a0cfb8d16c3ceefdc2c2281201b3af54bd815ace51c2809eb"}, {"2074b4d6c38bedd18aae03de2ba9de131517764275f41363de8e0d9bd6d5432db17d15e872ea02b4901069d11ec4cdd92fc94cf83"}, {"34805fdd738a6574c4016793afae55c587b1a3948fee4127c701f8f93128e8598e4710786a4f7adc92b192e874996203bd48aa7fb"}, {"54ed7db62f8e802ed7207935c1212ece8b43cbecf0c964e54a82e2b128d085dfb43842dd3c7e33310ba1ca00e70f5f4a3d48fe0b5"}, {"8961968b032ff73605f650885a2bb9ecf110945fa09d9bd7acd97c1bc0341ca19a447356cd4da02fac68def3710b3d28e4f2b0332"}, {"de3b3800e9c82fa25fc7436326e19d3dbfc442f32d829a0e90be1953b9ba3c232f1c956f51e27d75106c1554a6a64b4fce81d31d7"}, {"1677cae2203d26d9a45dd34d7c6cf6ba5dcd7f433816571c80e0cee6eb4fff9e01c46beee7c73b410beb624f55aa1e7593a8879645"}, {"24583ee1daf717abdeff1b58cbcdcd9cde24d2f555f6482b0faeb68328cb9642fced4775024f9d84a6faae2ceabd5a596d5d399ea5"}, {"3acac8b53d3c039a80b0ab30eebaf2dfeb23c63ba9422ce32b003fa0cf0fb9327ec5723df66c1b05268590ede3293a23b50f5cb2de"}, {"5f1a87ca78b398126d43f54081c4d6de71018c1c01c47d7c3fa455aa01165dfdbb62db719fad19416e1c9c7266ad8010a1c4d5507f"}, {"99d790dee333b3984cf28ec87b2a9a9415ae4727c264acce0acdc8136f27bdb408408aa26d3d018b75f27d3cac67b9cab85f69d19a"}, {"f8dbdb38b94a5319bbdb359b7ef8dceaaf2ba487404460c7c1588fdb62c70d95db851af986ff4716aff2e05729432ab77093ed6399"}, {"1928f723d37f12dc955479c4c4552865611be7167608bbb602b4e75c1d49a2271a46c5bfc0761b903059b4eeb32eef0cd8049475094"}, {"28b311b5e537fddd3b66b365e64137a7ba68ef8065fda2a09d4088115920c1823677985f2365db9ab5f2ee5307fa5c94a2834d77282"}, {"41d626a1342e96c324fecccf46047747f0f0b3cdfcb210cb89037e1c5678e4999a2b5f4007f310bef0e3c25fde22245e6c596606dda"}, {"6a7fb3d6d9bfb4004af39eb50b382b07d5e18b5a5cab08389755207d237e81913fe12f68b2f21874ad6d4c90e83e7679b76c798effa"}, {"ac46751f849deee4cca0a877a8cdf0c697323f0e3c365e24700dd020054cea22f46d98037051b03366c2fddaec5ef9a2d56fc982cb1"}, {"116ad41575c1745d9582fc60b981f77b8cfce562ddb9e8149026f8fd15be27e922b4ced068342b5cfd2d701a49cce6f386001a221b58"}, {"1c2cb6d19044fe10e48944bef3f3a96beaebeb32787ff980f122e7d50f9bec419e80ad482ae7daec78360a03f9cdeb55771c12defedc"}, {"2d9378346725fb773bb9fd04f1988da67336914d58f3e6d62850593285b794a41ce108ca94bcf74f0b4b98321cac6c4f602701c6d845"}, {"49b998543748f429d526134a1065a11c4092c540710aab4a4f9fe47e775c6d66c2d82f6c0ade0410a1d8525c4ba638830cd27a624f90"}, {"774268181530fcf380fedfd32bdb2a29d62cf725d2f862b79da2be8b24143fe1f2c00adf96566d5a1ecf14b63749c31d3d6612815488"}, {"c0eac2d472235b2bfd4acbe03d89ffe496526e1c0de5669cedf9285290a0ac60a106fa57571554f2de46d751c31bb6000e3e8dda310b"}, {"1381147622f886429fa88b4a2a0571da5c8a8577bb75a68c7e63d419c6d1b95d47f78eb68299c870d55b3c9f83e334f93093ed3c5c3c0"}, {"1f8ceed1c8e1e24f7ee5787fd40e43f0c806109bc2b4855681e15ce241fac2e56c0a274425bd2066202101724688d1866c62ae7fafa28"}, {"330973a5fe775764387d6c02dbb5a9649339ac2df106f4d1f9f09adbaf4740ada257b0d0bc97bb07f15f35a586fbe7caf95130710c993"}, {"528f01ad04e35daedf9ae04bbae346eec1a363509fb16bc47b7c7943f4b98eb4ee72f13fffa9ac461b8384f5673cb87d338b9fd5c854a"}, {"858c85f76093743ff89a49fd582c6ab8d9322c8d290251878f5f099f9f4ed05e38b2f5adb456e6730429f36257a84556356ada3c1e1f7"}, {"d8083937b1f84551d474db7ac817df339bbdb47334cf7f7a03c9f715791ffc79b38b9ea3489a1ee6135e925967d19d477f4f4d4372318"}, {"15d7584313b3bc8798a791dab80ae3905ea9be3b22831821d3f162057e7fce86d7dee6b0d7b3fb7beee60a5623b8923dc67d403999822d"}, {"2354b38822bcfac8af8eeab517ddec28d39bd74b744ac4bf727740296f0eb9c9ffc88712f3cab5cf9c6dcdf08712f47deaa0076ead89be"}, {"3926f041c6fa159b23b619c40d2e2f16c87a3622b6445cdf8d1fe12a7d3d0653d418b185245654d8c38c3443926e9ada56fded1b50d9fc"}, {"5c7360af13f0ef2cb726bbbef26c24ad4b96e3080d7889136537f6c0288f56e2e567fb774267217fcf66976c30e65f909ccde63362d693"}, {"958cf37e5a364246694f0692706376d4e0b7960447d7ced4e0cdda8dadf5f9179fa4e02746d1df500fc0d634410437f6bf1b7442f7acc4"}, {"f1eab58eb5c71fc2cacf0545930a09b21e9cd1b7ddfe11a3b56a6515ecb2d0e1fd12d43e64acf2f178cdd4a448919136493450cc9e3f6d"}, {"18754b01bf20a659f0159880c8f66a5ae0f39a3bfdbc8d514418325c72d83c37a149cbf89bc80162b66ca191c3282cd712ce31f4c2fd0c8"}, {"27906d329171087b3f8b2d3aa5228c67e2e0d78cdeeb990b80efc8a74976d8e49a9fe3c0d8960f83f6bc8bd54d0ed58d71fb72ba71fd939"}, {"400000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"}, {"67872ca9d39b5518412bed53895e0f27fd7b80f72da446bc42ee7a5da6e4c9734766ace7d837403ec2899bf90387a1728921f2197edfa6e"}, {"a77835432cbbf6ceed18c9033d3c4b99a64cbf2a3284cc5fa79a2b18eaf1acf71ebb8dc69f58d3df3b4c5af7cff52aea5219a95795b22a9"}, {"10ee72c277d72773a909c429a7548c6b221efcd71405cda9ede8ab2e401f7aa31aab8ffb46ffa467e62fa791c0b64b02581d571296a74d5d"}, {"1b63837be79b0397e86ea8b52b38f81ccd5d830d4625ea96a5f9cbb3e6554441cd7be7f2a85ba9f5539f6e2c50e72ee9ccc9dca16c5063b9"}, {"2c4e00a1d9d4f617b1cc0529fa3b6db8e9463183e9bb3621da5fa8edcefafeefce8f509fa16fb69f4d81c834cbdeb02540cffb95a41319fa"}, {"47ab1c7900dfdc52610401ea5863dcf91703149049b61a6e1e7e5bda857603394925f1781e6e8bf3dccf9d132d04fafde029b835d8f05154"}, {"73eec0c69d06927bd768087c2bc21d338b907bbd9c002ae5a2032b9e2c8b552761db707d783df56077b80d44217a5a3d4acc8b0be63a4fbb"}, {"bb891ac60df7e9f62c49d8336d43b1eb4be8bb6fb1d15fb8a570e44d51b8569a53b7d6dca47cbf406b54067d0b85f7d3bcb28ee2205d89a0"}, {"12f5cbf2b1cc27dbbfdb38786e49982acded90d9f977a9bc528ba75eabf0bfa26d98cbbcc3ad43ab7cd31b1fe76efc7083df4682b4120cc5a"}, {"1eab9ff0160962c74b1af90f95cd3d29adbb640e6554add84b7931dfa04e28b751a0f266a1b33e5de0acd76ca8ac4a49173058a505ddfc549"}, {"319cfcd322753fa0301232cb31c485bbab878d920fde083f77fc437e8f3c5b8a0d965fc02ff3b8b2b3d6d86ab61797ed67f56deca70ff64ff"}, {"504170a8ad3822a34ebb8d12ddc50c552caf76434c717af7b2b50c0a833ee345b0f3bbf70b08a04806cebd52f6056b39470a027a96ad144e9"}, {"81d2d35b0e2fb7462ab80f78aea7784c545dd909d7d04d8f7a38366f5b550cb6c5f457a1cba80ad24041a7447b133764bf2533e08a8bf27c8"}, {"d2017f75ce2995d310059e4ebeb957161ec68e523fb9049f1c07e3e4c4a3a166391c41e1a93deb724ac64ba1d265880ef522a9fe994bf2506"}, {"153b5f6d8be3a00bb5c81e4d70484bb3a2d0f2fb3c32593952a0b8ac4a342177ed0e10b74f1f8ccc42e8728039c3d7d28f523f82d0cd6aa6c3"}, {"225865a2c0cbfdcd817d9422168aee7b3c442517e82766f5eb365f06f00d072866911cd5b353680beb1a37bcdc42cf4a393cac0cc71b4b40bd"}, {"378ece003a720be897f1d0d8c9b8674f6ff110bccb60a683a42f8654ee6e8f08cbc7255bd4789f5d441b055529a992d3f3a074c75bd77a4bb5"}, {"59df2b8880f49b6691774d35f1b3dae97735fd0ce0595c7203833b9bad3ffe780711afbb9f1ced042fad893d56a879d64ea21e59112263359b"}, {"9160fb876bdc5eaf7bfc6d632f37367378934bd557d920ac3ebada0ca32e590bf3b51dd4093e76129beab8b3aaf56965b574e67a91bf90b73e"}, {"eb2b22d4f4f00e6a02c2acc3baebcb7fac8f18eea508b38cac00fe833c3ee4c9fb15c8f7d9b06c149a1643b78ca4e88ed43d72cb763dd0e2ed"}, {"17c6a1f29e2ce6049b50fd50207135b79c40630700711f5f0fe9156a8045977f6ed8b6dc67eb46505b87271c99ab6004287135541f8314b0fa9"}, {"2675e437b8ccece6133ca3b21ecaa6a5056b91c9f0992b3382b37204dbc9ef6d021022a89b4f4062dd7c9f4f2b19ee72ae17da7466677b148c1"}, {"3e36f6ce2e5294c66ef6cd66dfbe373aabcae921d0111831f05623f6d8b1c36576e146be61bfd1c5abfcb815ca50caaddc24fb58e6314e04c99"}, {"64a3dc8f4dfc4b725551e7131154a195846f9c59d822eed80ad39d707526889c691b16ff01d86e40c166d3baccbbbc33601251d424d4394f63e"}, {"a2cc46d794dff774ed9acd979904de9ee9a3be0197f68a827502204564830608d9de617c8d976e6ad7cbb94c1fe972528a0d35dca07b41e18b0"}, {"107589a84d0ba5b0c93fb333d1811dd1fc3c2cf37f2c8432e240df87159e3926668ad7d001fcc42fbc38f097f78889179b165981b766f1799f5b"}, {"1a9fecf5b66fed0012bce7ad42ce0ac1f57862d6972ac20e25d5481f48dfa0644ff84bda2cbc861d2b5b53243a0f9d9e6ddb1e63bfe7cb94d48c"}, {"2b119d47e1277bb933282a1dea337c31a5cc8fc38f0d97891c03740801533a88bd1b6600dc146c0cd9b0bf76bb17be68b55bf260b2c309c0885a"}, {"45ab5055d705d176560bf182e607ddee33f3958b76e7a052409be3f456f89fa48ad33b41a0d0ad73f4b5c06927339bce180068886d5cf21abab8"}, {"70b2db464113f843922512fbcfcea5afabafacc9e1ffe603c48b9f543e6fb1067a02b520ab9f6bb1e0d8280fe737ad55dc704b7bfb6c6579ec3d"}, {"b64de0d19c97eddceee7f413c6623699ccda453563cf9b58a14164ca16de52907384232a52f3dd3c2d2e60c872b1b13d809c071b611230a65d92"}, {"126e66150dd23d0a754984d2841968471024b1501c42aad293e7f91f9dd71b59b0b60bdb02b781042876149712e98e20c3349e9893e3f67a3588f"}, {"1dd09a06054c3f3ce03fb7f4caf6ca8a758b3dc974be18ade768afa2c9050ff87cb002b7e5959b5687b3c52d8dd270c74f5984a05521d59269aaa"}, {"303ab0b51c88d6caff52b2f80f627ff925949b87037959a73b7e4a14a4282b182841be95d5c5553cb791b5d470c6ed80038fa3768c42904e8b592"}, {"4e0451d88be2190a7ea22d28a815c769722a15deedd69a31f98f50671b46e5751fde3ada0a6721c3556cf27e0f8cd3e4c24fb4f170efd696db300"}, {"7e33bb472387893dfb95e1ff50390fc32018426f0ad2155a0785738907eb0b95b0289d1096f48198808405c91a234619b18ab9febe89e3c9c3667"}, {"cc25ce97f9dd5d90e1f5b00b6ed6a5924ce6b0b7c41bd347e7c26b6ebd1d02b6895ec342f25eec1fc57cd6961bef9f2be544c1c432648756c6675"}, {"14a3c06b4138d76bb7e6b812eeb8d1bbb068d8d427ec5af11edf1e50fd2e63ad3b9cbc4fffea6b1186e0e13d59cf2e1f4ce2e7f5721524304ca66a"}, {"2163217dd824dd0ffe26927f560b1aae364c8b234a409461e3d7c267e7d3b4178e2cbd6b6d15b99cc10a7cd895ea11558d5ab68f0787da2d43b4df"}, {"36020e4dec7e1154751d36deb205f7cce6ef6d1ccd33dfde80f27dc55e47ff1e6ce2e7a8d22687b984d7a49659f46751dfd3d350dc8c5dccb0e2d9"}, {"575d610c4ecef21e629e476ae02b8e417aa6f8ec8a0c60874fc58815f9ff3a1b5f7b9ac35ecfedefbb9829588ee248f719f500e66608b3628cc002"}, {"8d52ce208af3eb22be3baad45f77b0a34e6f5d2dd12b12fdc9dd00a5bc3ae727ff221c4bcf2ad73e71b737c21c4bd1f7aa801dbab626f73df0c3de"}, {"e49bc1071be8566c8ed8671034b8bc5b21e8d88ad911737e347f84a9f4f4194f5062c614915953630e30d10e49ba6b695bf7b46d4367ed725f029f"}, {"171cd82bc4fc3bcb2dc9aeefbc9b092b52e5b8c2035e2244d94dfdb00a23d5b8b959feddd099c85ff3d9a5db0c3997e42733798cd8b5a488d931691"}, {"25633cdf968d90919a53c1a49c18ddb5382de58111f5f3b5383376a36b7d7e45e7e93809d1b477d403f0358d10410dfdafc6dd10bdeb30e6feded5c"}, {"3c7aad63ad71c7f0b2b3c15164c2826c87a7346df77f8468ed5a99457b2cb4387f44b50f992b3cbc5e3375291224644d924d1433278fdb154d82851"}, {"61d52c06fc829967c056620323d9a96b83ce68eff6f235451060c971cb60f0de85272fe26f20058ad9b286470ca0b35c4b38c7d30bf431e7841f9b1"}, {"9e41b4ca45c421ecfe2cb4ea948a319f8b835e337bae642e03bf229d25db63926a7f8f0362583e0fdaf22f55343b5635c7edcae9c7f64e064455616"}, {"100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"}, {"19e1cb2a74e6d546104afb54e25783c9ff5ee03dcb6911af10bb9e9769b9325cd1d9ab39f60dd00fb0a266fe40e1e85ca2487c865fb7e9b7facd7250"}, {"29de0d50cb2efdb3bb463240cf4f12e669932fca8ca13317e9e68ac63abc6b3dc7aee371a7d634f7ced316bdf3fd4aba94866a55e56c8aa16566213b"}, {"43b9cb09df5c9dcea42710a69d5231ac887bf35c501736a7b7a2acb9007dea8c6aae3fed1bfe919f98be9e4702d92c0960755c4a5a9d3572d64a1690"}, {"6d8e0def9e6c0e5fa1baa2d4ace3e07335760c351897aa5a97e72ecf9955110735ef9fcaa16ea7d54e7db8b1439cbba733277285b1418ee19bcb3876"}, {"b13802876753d85ec73014a7e8edb6e3a518c60fa6ecd887697ea3b73bebfbbf3a3d427e85beda7d360720d32f7ac095033fee56904c67e5acac25d1"}, {"11eac71e4037f7149841007a9618f73e45c0c524126d869b879f96f6a15d80ce52497c5e079ba2fcf733e744cb413ebf780a6e0d763c1454e78dd3b79"}, {"1cfbb031a741a49ef5da021f0af0874ce2e41eef67000ab96880cae78b22d549d876dc1f5e0f7d581dee0351085e968f52b322c2f98e93eea6c967f09"}, {"2ee246b1837dfa7d8b12760cdb50ec7ad2fa2edbec7457ee295c3913546e15a694edf5b7291f2fd01ad5019f42e17df4ef2ca3b888176267ddf0392fc"}, {"4bd72fcac7309f6eff6ce1e1b92660ab37b64367e5dea6f14a2e9d7aafc2fe89b6632ef30eb50eadf34c6c7f9dbbf1c20f7d1a0ad048331653bac8c42"}, {"7aae7fc058258b1d2c6be43e5734f4a6b6ed90399552b7612de4c77e8138a2dd4683c99a86ccf97782b35db2a2b129245cc162941777f1521770f4d30"}, {"c673f34c89d4fe80c048cb2cc71216eeae1e36483f7820fddff10dfa3cf16e2836597f00bfcee2cacf5b61aad85e5fb59fd5b7b29c3fd93f9fde85891"}, {"14105c2a2b4e08a8d3aee344b77143154b2bdd90d31c5ebdecad4302a0cfb8d16c3ceefdc2c2281201b3af54bd815ace51c2809ea5ab4513a1e924cc0a"}, {"2074b4d6c38bedd18aae03de2ba9de131517764275f41363de8e0d9bd6d5432db17d15e872ea02b4901069d11ec4cdd92fc94cf822a2fc9f1f8e9ee01b"}, {"34805fdd738a6574c4016793afae55c587b1a3948fee4127c701f8f93128e8598e4710786a4f7adc92b192e874996203bd48aa7fa652fc9415af3bbf78"}, {"54ed7db62f8e802ed7207935c1212ece8b43cbecf0c964e54a82e2b128d085dfb43842dd3c7e33310ba1ca00e70f5f4a3d48fe0b4335aa9b0b85c41b39"}, {"8961968b032ff73605f650885a2bb9ecf110945fa09d9bd7acd97c1bc0341ca19a447356cd4da02fac68def3710b3d28e4f2b0331c6d2d02f1c63f886e"}, {"de3b3800e9c82fa25fc7436326e19d3dbfc442f32d829a0e90be1953b9ba3c232f1c956f51e27d75106c1554a6a64b4fce81d31d6f5de92ed0dea6856a"}, {"1677cae2203d26d9a45dd34d7c6cf6ba5dcd7f433816571c80e0cee6eb4fff9e01c46beee7c73b410beb624f55aa1e7593a88796444898cd66abd789f44"}, {"24583ee1daf717abdeff1b58cbcdcd9cde24d2f555f6482b0faeb68328cb9642fced4775024f9d84a6faae2ceabd5a596d5d399ea46fe42dcf7a4badf54"}, {"3acac8b53d3c039a80b0ab30eebaf2dfeb23c63ba9422ce32b003fa0cf0fb9327ec5723df66c1b05268590ede3293a23b50f5cb2dd8f7438bb063fe1b69"}, {"5f1a87ca78b398126d43f54081c4d6de71018c1c01c47d7c3fa455aa01165dfdbb62db719fad19416e1c9c7266ad8010a1c4d5507e0c52c3ea2b2deab6e"}, {"99d790dee333b3984cf28ec87b2a9a9415ae4727c264acce0acdc8136f27bdb408408aa26d3d018b75f27d3cac67b9cab85f69d1999dec523036c7385e8"}, {"f8dbdb38b94a5319bbdb359b7ef8dceaaf2ba487404460c7c1588fdb62c70d95db851af986ff4716aff2e05729432ab77093ed6398c538132618c03b6f0"}, {"1928f723d37f12dc955479c4c4552865611be7167608bbb602b4e75c1d49a2271a46c5bfc0761b903059b4eeb32eef0cd804947509350e53d8f4d08a612f"}, {"28b311b5e537fddd3b66b365e64137a7ba68ef8065fda2a09d4088115920c1823677985f2365db9ab5f2ee5307fa5c94a2834d77281ed07862bd1aadfaca"}, {"41d626a1342e96c324fecccf46047747f0f0b3cdfcb210cb89037e1c5678e4999a2b5f4007f310bef0e3c25fde22245e6c596606dd9bc5e67d68c0fb20a8"}, {"6a7fb3d6d9bfb4004af39eb50b382b07d5e18b5a5cab08389755207d237e81913fe12f68b2f21874ad6d4c90e83e7679b76c798eff9f2e53522eec765bbd"}, {"ac46751f849deee4cca0a877a8cdf0c697323f0e3c365e24700dd020054cea22f46d98037051b03366c2fddaec5ef9a2d56fc982cb0c27022167a2019041"}, {"116ad41575c1745d9582fc60b981f77b8cfce562ddb9e8149026f8fd15be27e922b4ced068342b5cfd2d701a49cce6f386001a221b573c86aeadc65b3b3b4"}, {"1c2cb6d19044fe10e48944bef3f3a96beaebeb32787ff980f122e7d50f9bec419e80ad482ae7daec78360a03f9cdeb55771c12defedb195e7b0f3074b231b"}, {"2d9378346725fb773bb9fd04f1988da67336914d58f3e6d62850593285b794a41ce108ca94bcf74f0b4b98321cac6c4f602701c6d8448c29976470bf4c84c"}, {"49b998543748f429d526134a1065a11c4092c540710aab4a4f9fe47e775c6d66c2d82f6c0ade0410a1d8525c4ba638830cd27a624f8fd9e8d6238ac345792"}, {"774268181530fcf380fedfd32bdb2a29d62cf725d2f862b79da2be8b24143fe1f2c00adf96566d5a1ecf14b63749c31d3d66128154875649a6aa617cca6ce"}, {"c0eac2d472235b2bfd4acbe03d89ffe496526e1c0de5669cedf9285290a0ac60a106fa57571554f2de46d751c31bb6000e3e8dda310a413a2d6449c124d93"}, {"1381147622f886429fa88b4a2a0571da5c8a8577bb75a68c7e63d419c6d1b95d47f78eb68299c870d55b3c9f83e334f93093ed3c5c3bf5a5b6cbfee4572ea2"}, {"1f8ceed1c8e1e24f7ee5787fd40e43f0c806109bc2b4855681e15ce241fac2e56c0a274425bd2066202101724688d1866c62ae7fafa278f270d998f522320b"}, {"330973a5fe775764387d6c02dbb5a9649339ac2df106f4d1f9f09adbaf4740ada257b0d0bc97bb07f15f35a586fbe7caf95130710c9921d5b199d3a2f4e7e8"}, {"528f01ad04e35daedf9ae04bbae346eec1a363509fb16bc47b7c7943f4b98eb4ee72f13fffa9ac461b8384f5673cb87d338b9fd5c85490c13299a7d7fd2eab"}, {"858c85f76093743ff89a49fd582c6ab8d9322c8d290251878f5f099f9f4ed05e38b2f5adb456e6730429f36257a84556356ada3c1e1f68b50ed37863276c1d"}, {"d8083937b1f84551d474db7ac817df339bbdb47334cf7f7a03c9f715791ffc79b38b9ea3489a1ee6135e925967d19d477f4f4d4372317732c38b63d15194c7"}, {"15d7584313b3bc8798a791dab80ae3905ea9be3b22831821d3f162057e7fce86d7dee6b0d7b3fb7beee60a5623b8923dc67d403999822cd8a330004de7b5196"}, {"2354b38822bcfac8af8eeab517ddec28d39bd74b744ac4bf727740296f0eb9c9ffc88712f3cab5cf9c6dcdf08712f47deaa0076ead89bdcf7c30f767a03f868"}, {"3926f041c6fa159b23b619c40d2e2f16c87a3622b6445cdf8d1fe12a7d3d0653d418b185245654d8c38c3443926e9ada56fded1b50d9fb5c97c0a79d933597e"}, {"5c7360af13f0ef2cb726bbbef26c24ad4b96e3080d7889136537f6c0288f56e2e567fb774267217fcf66976c30e65f909ccde63362d6922364c5a41e0bacb76"}, {"958cf37e5a364246694f0692706376d4e0b7960447d7ced4e0cdda8dadf5f9179fa4e02746d1df500fc0d634410437f6bf1b7442f7acc39bfb7b56c045c804a"}, {"f1eab58eb5c71fc2cacf0545930a09b21e9cd1b7ddfe11a3b56a6515ecb2d0e1fd12d43e64acf2f178cdd4a448919136493450cc9e3f6c55360a140bff1ae01"}, {"18754b01bf20a659f0159880c8f66a5ae0f39a3bfdbc8d514418325c72d83c37a149cbf89bc80162b66ca191c3282cd712ce31f4c2fd0c79e107e6c2bf76d435"}, {"27906d329171087b3f8b2d3aa5228c67e2e0d78cdeeb990b80efc8a74976d8e49a9fe3c0d8960f83f6bc8bd54d0ed58d71fb72ba71fd938191155856cc5aab30"}, {"40000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"}, {"67872ca9d39b5518412bed53895e0f27fd7b80f72da446bc42ee7a5da6e4c9734766ace7d837403ec2899bf90387a1728921f2197edfa6dfeb35c93f09fb4806"}, {"a77835432cbbf6ceed18c9033d3c4b99a64cbf2a3284cc5fa79a2b18eaf1acf71ebb8dc69f58d3df3b4c5af7cff52aea5219a95795b22a85959884eaa82c5726"}, {"10ee72c277d72773a909c429a7548c6b221efcd71405cda9ede8ab2e401f7aa31aab8ffb46ffa467e62fa791c0b64b02581d571296a74d5cb59285a3efc19bf03"}, {"1b63837be79b0397e86ea8b52b38f81ccd5d830d4625ea96a5f9cbb3e6554441cd7be7f2a85ba9f5539f6e2c50e72ee9ccc9dca16c5063b866f2ce1d5b4b379d4"}, {"2c4e00a1d9d4f617b1cc0529fa3b6db8e9463183e9bb3621da5fa8edcefafeefce8f509fa16fb69f4d81c834cbdeb02540cffb95a41319f96b2b097411d2962a6"}, {"47ab1c7900dfdc52610401ea5863dcf91703149049b61a6e1e7e5bda857603394925f1781e6e8bf3dccf9d132d04fafde029b835d8f051539e374ede340c2f275"}, {"73eec0c69d06927bd768087c2bc21d338b907bbd9c002ae5a2032b9e2c8b552761db707d783df56077b80d44217a5a3d4acc8b0be63a4fba9b259fc212e178479"}};

template<uint32_t LIMBS1, uint32_t LIMBS2>
unsigned calc_bound(const Num<LIMBS1, false>& af, const Num<LIMBS2, false>& ag) {
    static constexpr uint32_t LIMBS = LIMBS1 > LIMBS2 ? LIMBS1 : LIMBS2;
    Num<LIMBS*2, false> t; t.set_square(af);
    Num<LIMBS*2+1, false> v; v.set_square(ag); v <<= 2; v += t;
//    printf("calc_bound: af=%s ag=%s t=%s v=%s\n", af.to_string().c_str(), ag.to_string().c_str(), t.to_string().c_str(), v.to_string().c_str());
    return std::upper_bound(std::begin(BOUNDS), std::end(BOUNDS), v) - std::begin(BOUNDS) - 1;
}

template<uint32_t BITS>
unsigned calc_bound(const Range<BITS>& f, const Range<BITS>& g, int bits, int delta) {
    auto af = f.maxabs(bits);
    auto ag = g.maxabs(bits);
    if (delta == 1) {
        return calc_bound(af, ag);
    } else if (delta < 1) {
        UInt<BITS + 2> v; v.set_lshift(af, 1 - delta); v -= af; v += ag; v >>= 1 - delta;
        return calc_bound(af, v) + 1 - delta;
    } else {
        UInt<BITS + 2> v; v.set_sum(af, ag); v >>= delta + 1; v += ag;
        return calc_bound(ag, v) + delta + 1;
    }
}

template<uint32_t BITS>
struct DivStepState {
    Range<BITS> f;
    Range<BITS> g;
    int32_t delta;
    uint32_t iterations_done;

    std::string to_string() const {
        return "DivStepState(f=" + f.to_string() + ", g=" + g.to_string() + ", delta=" + std::to_string(delta) + ", iter=" + std::to_string(iterations_done) + ")";
    }

    DivStepState() {}
};

template<uint32_t BITS>
bool divsteps(DivStepState<BITS>& in, uint32_t range_bits) {
    assert(in.f.known_parity(range_bits) && in.f.is_odd(range_bits));
    while (in.g.known_parity(range_bits) ) {
//        printf("divsteps on %s\n", in.to_string().c_str());
        if (in.g.is_zero(range_bits)) return true;
        in.iterations_done += 1;
        if (in.g.is_odd(range_bits)) {
            if (in.delta > 0) {
                in.delta = 1 - in.delta;
                Range<BITS> t; t.set_difference(in.g, in.f);
                in.f = in.g;
                in.g.set_rshift(t, 1);
            } else {
                in.delta += 1;
                Range<BITS> t; t.set_sum(in.g, in.f);
                in.g.set_rshift(t, 1);
            }
        } else {
            in.delta += 1;
            in.g >>= 1;
        }
    }
//    printf("final: %s\n", in.to_string().c_str());
    return false;
}

template<uint32_t BITS>
struct Entry {
    DivStepState<BITS> divstep_state;
    UInt<BITS> suffix;
    uint32_t best_bound;
    uint32_t predicted_bound;
    uint32_t range_bits;
    uint32_t suffix_len;

    std::string to_string() const {
        return "Entry(suffix=" + suffix.to_string() + "/" + std::to_string(suffix_len) + ", bits=" + std::to_string(range_bits) + ", bound=" + std::to_string(best_bound) + ", state=" + divstep_state.to_string() + ")";
    }

    Entry() {}
    Entry(const DivStepState<BITS>& d, const UInt<BITS>& s, uint32_t b, uint32_t r, uint32_t sl) : divstep_state{d}, suffix{s}, best_bound{b}, range_bits{r}, suffix_len{sl} {}
};

std::mutex lock_best;
uint32_t best_iter[B+1] = {0};
UInt<B> best[B+1];
uint32_t best_complete_iter{0};
UInt<B> best_complete;

template<uint32_t BITS>
void work(std::deque<Entry<BITS>>& queue) {
    Entry<BITS> entry = queue.back();
    queue.pop_back();
    uint32_t local_best[B+1] = {0};
    uint32_t local_best_complete{0};
    for (int b = 0; b < 2; ++b) {
        Entry<BITS> nentry;
        DivStepState<BITS>& state = nentry.divstep_state;
        Range<BITS>& f = state.f;
        Range<BITS>& g = state.g;
        if (b) {
            f.cte.set_sum(entry.divstep_state.f.cte, entry.divstep_state.f.var);
            g.cte.set_sum(entry.divstep_state.g.cte, entry.divstep_state.g.var);
        } else {
            f.cte = entry.divstep_state.f.cte;
            g.cte = entry.divstep_state.g.cte;
        }
        f.var.set_lshift(entry.divstep_state.f.var, 1);
        g.var.set_lshift(entry.divstep_state.g.var, 1);
        state.delta = entry.divstep_state.delta;
        state.iterations_done = entry.divstep_state.iterations_done;
        nentry.suffix = entry.suffix;
        nentry.suffix_len = entry.suffix_len + 1;
        nentry.range_bits = entry.range_bits - 1;
        if (b) nentry.suffix.set_bit(entry.suffix_len, true);
//        printf("split %i from %s\n", b, entry.to_string().c_str());
        bool finished = divsteps(state, entry.range_bits - 1);
        uint32_t new_bound = finished ? state.iterations_done : state.iterations_done + calc_bound(state.f, state.g, nentry.range_bits, state.delta);
        if (new_bound > local_best[nentry.suffix_len]) {
            local_best[nentry.suffix_len] = new_bound;
            std::unique_lock<std::mutex> lock(lock_best);
            if (state.iterations_done > best_iter[nentry.suffix_len]) {
                best_iter[nentry.suffix_len] = new_bound;
                best[nentry.suffix_len] = nentry.suffix;
                printf("iter(%s/%i) %s %u\n", nentry.suffix.to_string().c_str(), (int)nentry.suffix_len, finished ? "=" : "<=", (unsigned)new_bound);
            }
        }
        if (!finished) {
            assert(new_bound <= entry.predicted_bound);
            nentry.best_bound = std::min(entry.best_bound, new_bound);
            if (state.delta == 1) {
                nentry.predicted_bound = 0xffffffff;
            } else {
                nentry.predicted_bound = std::min(entry.predicted_bound, new_bound);
            }
            if (nentry.best_bound <= PRUNE_BOUND) continue;
            queue.emplace_back(nentry);
        } else {
            if (state.iterations_done > entry.best_bound) {
                printf("UGH %s/%i: %i iter with bound %u\n", nentry.suffix.to_string().c_str(), (int)nentry.suffix_len, (int)state.iterations_done, (unsigned)entry.best_bound);
                assert(false);
            }
            if (state.iterations_done > local_best_complete) {
                local_best_complete = state.iterations_done;
                std::unique_lock<std::mutex> lock(lock_best);
                if (state.iterations_done > best_complete_iter) {
                    best_complete_iter = state.iterations_done;
                    best_complete = nentry.suffix;
                }
            }
        }
    }
}


template<typename WorkItem>
class WorkQueue {
    mutable std::mutex mutex;
    std::vector<WorkItem> stealable;
    std::condition_variable have_stealable;
    uint32_t running{0};
    uint64_t total_steps{0};
    uint64_t total_thefts{0};
    uint32_t total_stalls{0};

    void update_stats(uint64_t& steps) {
        total_steps += steps;
        steps = 0;
    }

public:
    template<typename Worker>
    void run(Worker&& worker, uint32_t id) {
        std::deque<WorkItem> queue;

        {
            std::unique_lock<std::mutex> lock(mutex);
            ++running;
        }

        uint64_t steps = 0;
        while (true) {
            while (queue.size()) {
                worker(queue);
                ++steps;
                if ((steps & 0xffff) == 0) {
                    std::unique_lock<std::mutex> lock(mutex);
                    update_stats(steps);
                    if (queue.size() > 1 && stealable.size() < running) {
                        stealable.emplace_back(std::move(queue.front()));
                        queue.pop_front();
                        have_stealable.notify_all();
                    }
                }
            }
            std::unique_lock<std::mutex> lock(mutex);
//            printf("thread %u ran out of work after %lu steps\n", (unsigned)id, (unsigned long)steps);
            update_stats(steps);
            if (stealable.empty()) ++total_stalls;
            while (stealable.empty()) {
                if (--running == 0) {
                    have_stealable.notify_all();
                    return;
                }
                have_stealable.wait(lock);
                ++running;
            }
            queue.emplace_back(std::move(stealable.back()));
            stealable.pop_back();
            ++total_thefts;
//            printf("thread %u stole work\n", (unsigned)id);
        }
    }

    void add_work(WorkItem&& item) {
        std::unique_lock<std::mutex> lock(mutex);
        stealable.emplace_back(std::move(item));
    }

    uint64_t get_steps() const {
        std::unique_lock<std::mutex> lock(mutex);
        return total_steps;
    }

    uint64_t thefts() const {
        std::unique_lock<std::mutex> lock(mutex);
        return total_thefts;
    }

    uint32_t get_stalls() const {
        std::unique_lock<std::mutex> lock(mutex);
        return total_stalls;
    }

    uint32_t running_threads() const {
        std::unique_lock<std::mutex> lock(mutex);
        return running;
    }
};

template<uint32_t BITS>
Entry<BITS> init_entry(const UInt<BITS>& modulus, uint32_t bits) {
    Entry<BITS> ret;
    ret.divstep_state.f.cte = modulus;
    ret.divstep_state.g.var.set_bit(0, true);
    ret.divstep_state.delta = 1;
    ret.divstep_state.iterations_done = 0;
    ret.best_bound = calc_bound(ret.divstep_state.f, ret.divstep_state.g, bits, 1);
    ret.predicted_bound = 0xffffffff;
    ret.range_bits = bits;
    ret.suffix_len = 0;
    return ret;
}



}


void stats_thread(const WorkQueue<Entry<B>>& queue, std::atomic_bool& stop) {
    while (!stop) {
        printf("%lu steps, %lu thefts, %lu stalls, %u threads\n", (unsigned long)queue.get_steps(), (unsigned long)queue.thefts(), (unsigned long)queue.get_stalls(), (unsigned)queue.running_threads());
        std::this_thread::sleep_for(std::chrono::milliseconds{1000});
    }
}

int main(void) {
    WorkQueue<Entry<B>> queue;
    static constexpr UInt<B> MODULUS(MODULUS_STR);
    queue.add_work(init_entry<B>(MODULUS, RANGE_BITS));
    std::vector<std::thread> threads;
    for (uint32_t i = 0; i < THREADS; ++i) {
        threads.emplace_back([&](){queue.run([](std::deque<Entry<B>>& q){work(q);}, i);});
    }
    std::atomic_bool stop{false};
    threads.emplace_back([&](){stats_thread(queue, stop);});
    for (uint32_t i = 0; i < THREADS; ++i) {
        threads[i].join();
    }
    stop = true;
    threads[THREADS].join();
    printf("Completed search after %llu steps\n", (unsigned long long)queue.get_steps());
    for (uint32_t sl = 0; sl <= B; ++sl) {
        if (best_iter[sl]) {
            printf("Highest bound for %u-bit suffix: iter(%s/%u) <= %u\n", (unsigned)sl, best[sl].to_string().c_str(), (unsigned)sl, (unsigned)best_iter[sl]);
        }
    }
    if (best_complete_iter) printf("Best found: %s with %u iterations\n", best_complete.to_string().c_str(), (unsigned)best_complete_iter);
    printf("Proved bound: max_iterations %s %u\n", best_complete_iter ? "=" : "<=", std::max(best_complete_iter, PRUNE_BOUND));
    return 0;
}

fgcd_recurse.py

import math
import bisect
import heapq

TABLE_W = [
    1, 5, 5, 5, 17, 17, 65, 65, 157, 181, 421, 541, 1165, 1517, 2789, 3635, 8293, 8293, 24245, 27361, 56645,
    79349, 132989, 213053, 415001, 613405, 1345061, 1552237, 2866525, 4598269, 7839829, 12565573, 22372085,
    35806445, 71013905, 74173637, 205509533, 226964725, 487475029, 534274997, 1543129037, 1639475149,
    3473731181, 4500780005, 9497198125, 12700184149, 29042662405, 36511782869, 82049276437, 90638999365,
    234231548149, 257130797053, 466845672077, 622968125533, 1386285660565, 1876581808109, 4208169535453
]

assert len(TABLE_W) == 57

def calc_bound(w):
    if w <= 2**(21*2):
        return ((w**19).bit_length() - 1) // 14
    elif w <= 2**(46*2):
        return ((w**49).bit_length() + 45) // 34
    else:
        return ((w**49).bit_length() - 1) // 34

for i in range(len(TABLE_W)):
    assert calc_bound(TABLE_W[i]) >= i

def find_bound(b):
    low = 0
    high = 1
    while calc_bound(high) < b:
        high *= 2
    while True:
        if low + 1 == high:
            return high
        mid = (low + high) // 2
        bnd = calc_bound(mid)
        if bnd < b:
            low = mid
        else:
            high = mid

while True:
    n = find_bound(len(TABLE_W))
    if n.bit_length() > 515:
        break
    TABLE_W.append(n)

maxbits=0

def getbound(f, g):
    global maxbits
    maxbits = max([f.bit_length(), g.bit_length(), maxbits])
    w = (f*f) + ((g*g)<<2)
    return bisect.bisect(TABLE_W, w) - 1

# Entries are tuples:
# - (negation of) Best bound
# - (negation of) Iterations done so far
# - suffix
# - suffixlen
# - delta
# - f: (constant, multiplier)
# - g: (constant, multiplier)

def known_parity(v):
    return v[1] & 1 == 0

def is_odd(v):
    return v[0] & 1

def is_zero(v):
    return v[0] == 0 and v[1] == 0

def shift(a):
    assert(a[0] & 1 == 0)
    assert(a[1] & 1 == 0)
    return (a[0] >> 1, a[1] >> 1)

def addshift(a, b):
    return shift((a[0] + b[0], a[1] + b[1]))

def subshift(a, b):
    return shift((a[0] - b[0], a[1] - b[1]))

def divsteps(delta, f, g, it, bits):
    if bits == 0:
        f = (f[0], 0)
        g = (g[0], 0)
    while known_parity(g):
        if is_zero(g):
            return (None, None, None, it)
        it += 1
        if is_odd(g):
            if delta > 0:
                delta, f, g = 1 - delta, g, subshift(g, f)
            else:
                delta, g = 1 + delta, addshift(g, f)
        else:
            delta, g = 1 + delta, shift(g)
    return (delta, f, g, it)

F=2**256-2**32-977
G=0x84888aea60a71431723f636c02333a3217c923bc6c804a53bdaae9bf95413b69

for bits in range(257):
    delta, f, g, it = divsteps(1, (F, 0), (G & ((1 << bits) - 1), 1 << bits), 0, 256 - bits)
    print("bits=%i: it=%i f=%x+%x*x g=%x+%x*x x=0..2**%i-1 delta=%i" % (bits, it, f[0], f[1], g[0], g[1], 256-bits, delta))

exit()

def newbound(f, g, delta, bits):
    global maxbits
    af = max(abs(f[0]), abs(f[0] + (f[1] << bits) - f[1]))
    ag = max(abs(g[0]), abs(g[0] + (g[1] << bits) - g[1]))
    if delta == 1:
        return getbound(af, ag)
    elif delta < 1:
        rf = af
        t = ((af << (1 - delta)) - af + ag)
#        maxbits=max(maxbits, t.bit_length())
        rg = t >> (1 - delta)
        return 1 - delta + getbound(rf, rg)
    else:
        nd = 1 - delta
        nf = ag
        t = ag + af
        maxbits=max(maxbits, t.bit_length())
        ng = t >> 1
        rf = nf
        rg = nf + (ng >> (1 - nd))
        return 2 - nd + getbound(rf, rg)

WORST=0
COUNTS={-741:1}
steps=0
disc=0

#P = 2**256 - 2**32 - 977
P = 0xfffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364141
Q = [(0, -741, 0, 0, 1, (P, 0), (0, 1), None)]

while len(Q):
    it, bound, suffix, suffixlen, delta, f, g, pbound = heapq.heappop(Q)
    COUNTS[bound] -= 1
    steps += 1
    if steps & 0xffff == 0:
        print("proc %i: %x/%i bound=%i max=%i prog=%f (%s)" % (steps, suffix, suffixlen, -bound, WORST, disc/(1<<256), COUNTS))
        print(maxbits)
    for b in (0, 1):
        nsuffix = suffix + (b << suffixlen)
        nsuffixlen = suffixlen + 1
        nbits = 256 - nsuffixlen
        nd = delta
        nf = (f[0] + (f[1] if b else 0), f[1] * 2)
        ng = (g[0] + (g[1] if b else 0), g[1] * 2)
        nd, nf, ng, nit = divsteps(nd, nf, ng, -it, nbits)
        if nd is None:
            disc += (1 << nbits)
            if nit > WORST:
                print("%x/%i has %i iter" % (nsuffix, nsuffixlen, nit))
                WORST=nit
        else:
            nb = nit + newbound(nf, ng, nd, nbits)
            nbound = max(bound, -nb)
            if pbound is not None:
                assert(nb <= pbound)
            if nd == 1:
                npbound = None
            else:
                npbound = nb
            if nbound <= -1:
                if nbound not in COUNTS:
                    COUNTS[nbound] = 1
                else:
                    COUNTS[nbound] += 1
                heapq.heappush(Q, (-nit, nbound, nsuffix, nsuffixlen, nd, nf, ng, npbound))
            else:
                disc += (1 << nbits)
print("%i steps" % steps)

gen_stable_hull.sage

import functools

# How many bits of precision to use in number field approximations.
CONFIG_PRECISION_BITS = 4096

# What range of delta to use for matrix transactions (1/2-CONFIG_DELTA_RANGE...1/2+CONFIG_DELTA_RANGE)
CONFIG_DELTA_RANGE = 2

# How many bits in the evaluation of a new points fewer than the difference are ignored
# (i.e. if a coordinate x1 is replaced with a coordinate x2, and (x1-x2)/x2 is 1/2^117,
# x2 is assumed to have a relative accuracy of at least 1/2^(117-CONFIG_RADIUS_BITS).
CONFIG_RADIUS_BITS = 10

# How many (assumed to be accurate) bits in an evaluation must match a replacement.
CONFIG_CHECKSUM_BITS = 64

# Transitions to apply in every iteration.
# e = even; + = odd, positive delta; - = odd, negative delta
TRANSITIONS = [
    # Identity.
    "",
    # All 1-step transitions
    "+",
    # All 3-step transitions
    "e+e", "e+-",
    # All 5-step transitions
    "ee+ee", "ee+e-", "ee+-e", "ee+--",
    # All 7-step transitions
    "eee+eee", "eee+ee-", "eee+e-e", "eee+e--", "eee+-ee", "eee+-e-", "eee+--e", "eee+---",
    # All 9-step transitions
    "eeee+eeee", "eeee+eee-", "eeee+ee-e", "eeee+ee--", "eeee+e-ee", "eeee+e-e-", "eeee+e--e", "eeee+e---",
    "eeee+-eee", "eeee+-ee-", "eeee+-e-e", "eeee+-e--", "eeee+--ee", "eeee+--e-", "eeee+---e", "eeee+----",
    # Worst known 54-step transition
    "ee+eee+-+e+ee+-ee+eee+-e+-ee+eee+-e+e+e+-ee+eee+-++e+-",
    # The same, repeated twice.
    "ee+eee+-+e+ee+-ee+eee+-e+-ee+eee+-e+e+e+-ee+eee+-++e+-ee+eee+-+e+ee+-ee+eee+-e+-ee+eee+-e+e+e+-ee+eee+-++e+-",
    # The same, repeated 4 times.
    "ee+eee+-+e+ee+-ee+eee+-e+-ee+eee+-e+e+e+-ee+eee+-++e+-ee+eee+-+e+ee+-ee+eee+-e+-ee+eee+-e+e+e+-ee+eee+-++e+-ee+eee+-+e+ee+-ee+eee+-e+-ee+eee+-e+e+e+-ee+eee+-++e+-ee+eee+-+e+ee+-ee+eee+-e+-ee+eee+-e+e+e+-ee+eee+-++e+-"
]

# Define number field with S = sqrt(D) and A = ((1591853137+3*sqrt(D))/2^55)^(1/54).
print("Defining number field...")
D = 273548757304312537
K.<S> = NumberField(x^2-D)
Kz.<z> = K[]
L.<A> = NumberField(z^54-(1591853137+3*S)/2^55)
S_SYMBOLIC = sqrt(D)
A_SYMBOLIC = ((1591853137+3*S_SYMBOLIC)/2^55)**(1/54)

# Define embedding phi of L into reals with requested precision.
print("Finding embedding...")
RR = RealField(CONFIG_PRECISION_BITS)
for phi in L.embeddings(RR):
    if phi(A) > 0 and phi((A^54)*2^55-1591853137) > 0 and phi(S) > 0:
        break

def convex_hull_matrix(ts):
    """Compute vertices of convex hull of matrices over L."""
    mp = {tuple(int(phi(x) * (1 <<CONFIG_PRECISION_BITS) + 1/2) for x in t.dense_coefficient_list()): t for t in ts}
    vertex_keys = Polyhedron(vertices=mp.keys()).Vrepresentation()
    return [mp[tuple(key)] for key in vertex_keys]

print("Precomputing transition matrices from delta=1/2 to delta=1/2...")
MAT_EVEN = Matrix([[1,0],[0,1/2]])/A
MAT_ODD_POS = Matrix([[0,1],[-1/2,1/2]])/A
MAT_ODD_NEG = Matrix([[1,0],[1/2,1/2]])/A
MATS = []
for seq in TRANSITIONS:
    delta = 1/2
    mat=Matrix([[1,0],[0,1]])
    for step in seq:
        if step == '+':
            assert(delta > 0)
            mat = MAT_ODD_POS * mat
            delta = 1 - delta
        elif step == '-':
            assert(delta < 0)
            mat = MAT_ODD_NEG * mat
            delta = 1 + delta
        elif step == 'e':
            mat = MAT_EVEN * mat
            delta = 1 + delta
        else:
            assert False
    assert delta == 1/2
    MATS.append(mat)
MATS = convex_hull_matrix(MATS)
print(len(MATS))

def cmp_point(a, b):
    """Compare two points with coordinates in L."""
    if a[0] != b[0]:
        d = phi(a[0] - b[0])
        if d > 0:
            return 1
        if d < 0:
            return -1
        assert(false)
    if a[1] != b[1]:
        d = phi(a[1] - b[1])
        if d > 0:
            return 1
        if d < 0:
            return -1
        assert(false)
    return 0

# A key function using cmp_point as a backend.
key_point = functools.cmp_to_key(cmp_point)

def convex_hull(points):
    """Computes the convex hull of a set of 2D points.

    Input: a list of (x, y) pairs representing the points.
    Output: a list of vertices of the convex hull in counter-clockwise order.

    This implements Andrew's monotone chain algorithm, based on
    https://en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain
    """
    # Boring case: no points.
    if len(points) <= 1:
        return points
    # Sort the points lexicographically (tuples are compared lexicographically).
    points = sorted(points, key=key_point)
    # Build lower hull.
    lower = [points[0]]
    for point in points[1:]:
        while len(lower) >= 2 and phi((lower[-1][0] - lower[-2][0]) * (point[1] - lower[-2][1]) - (lower[-1][1] - lower[-2][1]) * (point[0] - lower[-2][0])) <= 0:
            lower.pop()
        if lower[-1] != point:
            lower.append(point)
    if len(lower) == 1:
        return lower
    # Build upper hull.
    upper = [points[-1]]
    for point in reversed(points[:-1]):
        while len(upper) >= 2 and phi((upper[-1][0] - upper[-2][0]) * (point[1] - upper[-2][1]) - (upper[-1][1] - upper[-2][1]) * (point[0] - upper[-2][0])) <= 0:
            upper.pop()
        if upper[-1] != point:
            upper.append(point)
    # Concatenation of the lower and upper hulls gives the convex hull.
    # The last point of each list is omitted because it is repeated at the beginning
    # of the other list.
    return lower[:-1] + upper[:-1]

def approx(s, use_bits, height_bound):
    """Approximate an element of L with a nearby one of the form (x+y*S)*A**z."""
    ls = L(s)
    if s == 0:
        return ls
    rets = set()
    for an in range(54):
        lv = ls/A**an
        val = phi(lv)
        poly = algdep(val, 2, use_bits=use_bits, height_bound=height_bound)
        if poly is None:
            continue
        if poly.degree() == 1:
            rets.add(L(-poly[0] * A**an / poly[1]))
            continue
        if poly.degree() == 2:
            a,b,c = poly[2], poly[1], poly[0]
            disc = b**2 - 4*a*c
            sqf = disc.squarefree_part()
            if sqf == D:
                sq = sqrt(disc / sqf)
                rets.add(L(-b + sq*S) * A**an / (2*a))
                rets.add(L(-b - sq*S) * A**an  / (2*a))
    if rets:
        best = None
        bestd = None
        for ret in rets:
            diff = abs(phi(1 - ret/ls))
            if bestd is None or diff < bestd:
                best = ret
                bestd = diff
        return best
    return None

def stable_hull(h):
    print("Computing stable hull for %s..." % h)
    hull=convex_hull(h)
    iter=0
    while True:
        iter += 1
        newhull = convex_hull([tuple(m*vector(p)) for p in hull for m in MATS])
        if newhull == hull:
            print("- iteration %i: finished" % iter)
            return newhull
        oldhull, hull = hull, newhull
        if len(hull) != len(oldhull):
            print("- iteration %i: len %i -> %i" % (iter, len(oldhull), len(hull)))
        else:
            print("- iteration %i: len %i" % (iter, len(hull)))
            replaced=0
            for i in range(len(hull)):
                if hull[i] == oldhull[i]:
                    continue
                rep = [None, None]
                for j in range(2):
                    if hull[i][j] == oldhull[i][j]:
                        rep[j] = hull[i][j]
                        continue
                    if oldhull[i][j] == 0:
                        continue
                    bits = -log(abs(phi(1-hull[i][j]/oldhull[i][j])),2).n(128)
                    print("  - %s%i changed (%.3f bit difference)" % ("Y" if j else "X", i, bits))
                    if bits > CONFIG_RADIUS_BITS + CONFIG_CHECKSUM_BITS:
                        use_bits = int(bits - CONFIG_RADIUS_BITS)
                        bound_bits = (use_bits - CONFIG_CHECKSUM_BITS - 2) // 3
                        height_bound = 2**bound_bits
                        app = approx(hull[i][j], use_bits, height_bound)
                        if app is not None:
                            rep[j] = app
                if rep[0] is not None and rep[1] is not None:
                    if rep == hull[i]:
                        print("  - point %i: unchanged")
                    else:
                        print("  - point %i: adding limit %s" % (i, rep))
                        hull.append(tuple(rep))
                        replaced += 1
            if replaced > 0:
                hull = convex_hull(hull)
                print("  - %i replaced; new size %i" % (replaced, len(hull)))

with open("hull_trih_out.txt", "w") as f:
    f.write("%s\n" % stable_hull([(0,0),(1,0),(1,1/2)]))
with open("hull_tri_out.txt", "w") as f:
    f.write("%s\n" % stable_hull([(0,0),(1,0),(1,1)]))
with open("hull_sq_out.txt", "w") as f:
    f.write("%s\n" % stable_hull([(0,0),(1,0),(1,1),(0,1)]))

hull-bound-int.py

import argparse, sys, random, math, heapq

LOGLEVEL = 0
MAX_BOUNDING_BOX_COST = 16384
MAX_HORIZONTAL_BAND_COST = 128
MAX_DIRECT_COMPUTE_SET_SIZE = 1024

def line_eq(p1, p2):
    xd = p2[0] - p1[0]
    yd = p2[1] - p1[1]
    return (-yd, xd, yd * p1[0] - xd * p1[1])

def cross(o, a, b):
    """2D cross product of the OA and OB vectors, i.e. the z-component of
    their 3D cross product.

    Returns a positive value if OAB makes a counter-clockwise turn,
    negative for clockwise turn, and zero if the points are collinear.
    """
    return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])

def gmin(a, b):
    if a is None:
        return b
    return min(a, b)

def gmax(a, b):
    if a is None:
        return b
    if b is None:
        return a
    return max(a, b)

def convex_hull(points):
    """Computes the convex hull of a set of 2D points.

    Input: an iterable sequence of (x, y) pairs representing the points.
    Output: a list of vertices of the convex hull in counter-clockwise order,
      starting from the vertex with the lexicographically smallest coordinates.
    Implements Andrew's monotone chain algorithm. O(n log n) complexity.
    """
    # Boring case: no points.
    if len(points) <= 1:
        return points
    # Sort the points lexicographically (tuples are compared lexicographically).
    points = sorted(points)
    # Build lower hull.
    lower = [points[0]]
    for p in points[1:]:
        while len(lower) >= 2 and cross(lower[-2], lower[-1], p) <= 0:
            lower.pop()
        if lower[-1] != p:
            lower.append(p)
    if (len(lower) == 1):
        return lower
    # Build upper hull.
    upper = [points[-1]]
    for p in reversed(points[:-1]):
        while len(upper) >= 2 and cross(upper[-2], upper[-1], p) <= 0:
            upper.pop()
        if upper[-1] != p:
            upper.append(p)
    # Concatenation of the lower and upper hulls gives the convex hull.
    # The Last point of each list is omitted because it is repeated at the beginning
    # of the other list. 
    return lower[:-1] + upper[:-1]

def in_convex_hull(hull, point):
    """Determine if a point is in a convex hull (returned by convex_hull)."""
    if len(hull) == 0:
        # Empty set.
        return False
    elif len(hull) == 1:
        # Single point.
        return hull[0] == point
    elif len(hull) == 2:
        # Two points: test for colinearity, plus whether the point is between the vertices.
        return cross(hull[0], hull[1], point) == 0 and hull[0][0] <= point[0] <= hull[1][0] and (hull[0][0] != hull[1][0] or hull[0][1] <= point[1] <= hull[1][1])
    else:
        # Three or more points: point needs to be left of every edge, when traversing in 
        # counter-clockwise order.
        if cross(hull[-1], hull[0], point) < 0:
            return False
        for n in range(1, len(hull)):
            if cross(hull[n-1], hull[n], point) < 0:
                return False
        return True

def count_divsteps(f, g, delta=1):
    """Compute the number of divsteps safegcd needs with a given input."""
    assert f & 1
    it = 0
    while g != 0:
        it += 1
        if g & 1:
            if delta > 0:
                delta, f, g = 1-delta, g, (g-f)>>1
            else:
                delta, g = 1+delta, (g+f)>>1
        else:
            delta, g = 1+delta, g>>1
    return it

def process_divstep(state, step, cost_factor=1, fixed_bit=None, prove_limit=None):
    new_state = dict()
    max_iter = None
    if LOGLEVEL > 0:
        print("- step %i: %i distinct deltas, avg %f points per delta, max %i points per delta" % (step, len(state), sum(len(x) for x in state.values()) / len(state), max(len(x) for x in state.values())), file=sys.stderr)
    for delta, points in state.items():
#      if step % 12 == 0:
        points = convex_hull(points)
        min_g = min(g for _, g in points)
        max_g = max(g for _, g in points)
        min_g_int = -((-min_g) >> step)
        max_g_int = max_g >> step
        if min_g_int == 0 and max_g_int == 0:
            if max_iter is None or step > max_iter:
                max_iter = step
                if prove_limit is not None and max_iter > prove_limit:
                    return ({}, max_iter)
            continue

        min_f = min(f for f, _ in points)
        max_f = max(f for f, _ in points)
        min_f_int = -((-min_f) >> step)
        max_f_int = max_f >> step
        min_f_odd = min_f_int + 1 - (min_f_int & 1)
        max_f_odd = max_f_int - 1 + (max_f_int & 1)
        if min_g_int > max_g_int:
            continue
        if min_f_odd > max_f_odd:
            continue
        if (((max_f_odd - min_f_odd) >> 1) + 1) * ((max_g_int - min_g_int) + 1) <= 2048:
            for f in range(min_f_odd, max_f_odd + 1, 2):
                for g in range(min_g_int, max_g_int + 1):
                    if in_convex_hull(points, (f << step, g << step)):
                        it = step + count_divsteps(f, g, delta)
                        if max_iter is None or it > max_iter:
                            max_iter = it
                            if prove_limit is not None and max_iter > prove_limit:
                                return ({}, max_iter)
            continue

        if LOGLEVEL > 1:
            print("  - delta %i: %i points" % (delta, len(points)), file=sys.stderr)

        if fixed_bit is None or not fixed_bit:
            # divsteps^n(delta,f,g) = divsteps^{n-1}(1+delta,f,g//2)
            k = 1+delta
            new_set = [(f << 1, g) for f, g in points]
            if k not in new_state:
                new_state[k] = new_set
            else:
                new_state[k].extend(new_set)

        if fixed_bit is None or fixed_bit:
            if delta > 0:
                # divsteps^n(delta,f,g) = divsteps^{n-1}(1-delta,g,(g-f)//2)
                k = 1-delta
                new_set = [(g << 1, g-f) for f, g in points]
                if k not in new_state:
                    new_state[k] = new_set
                else:
                    new_state[k].extend(new_set)
            else:
                # divsteps^n(delta,f,g) = divsteps^{n-1}(1+delta,f,(f+g)//2)
                k = 1+delta
                new_set = [(f << 1, f+g) for f, g in points]
                if k not in new_state:
                    new_state[k] = new_set
                else:
                    new_state[k].extend(new_set)


    return (new_state, max_iter)


def countsteps(minmod, maxmod, half_range, cost_factor=1, fixed_bits=None, prove_limit=None):
    assert minmod & 1
    assert maxmod & 1
    if half_range:
        state = {1: [(minmod, 0), (minmod, minmod//2), (maxmod, 0), (maxmod, maxmod//2)]}
    else:
        state = {1: [(minmod, 0), (minmod, minmod-1), (maxmod, 0), (maxmod, maxmod-1)]}
    steps = 0
    max_iter = 0
    while len(state):
        fixed_bit = None
        if fixed_bits is not None and len(fixed_bits) > steps:
            fixed_bit = fixed_bits[steps]
        state, new_max_iter = process_divstep(state=state, step=steps, cost_factor=cost_factor, fixed_bit=fixed_bit, prove_limit=prove_limit)
        max_iter = gmax(max_iter, new_max_iter)
        if prove_limit is not None and max_iter > prove_limit:
            return max_iter
        steps += 1
    return max_iter


parser = argparse.ArgumentParser(description="Analyse the maximum number of divsteps needed in the safegcd algorithm.")
parser.add_argument('--verbose', '-v', action='count', default=0, help="Increase verbosity (can be repeated up to 3 times)")
parser.add_argument('--threshold', '-t', action='store_true', default=False, help="Find the largest modulus that can be processed with a given number of divsteps")
parser.add_argument('--exact-modulus', '-e', action='store_true', default=False, help="Make computations for the exact modulus specified, rather than assuming any modulus up to that number")
parser.add_argument('--half-range', action='store_true', default=False, help="Assume inputs in range 0..MODULUS/2 rather than 0..MODULUS-1")
parser.add_argument('--prove', type=int, nargs=1, default=None, help="Recursive split the range in order to prove the specified bound")
parser.add_argument('--self-test', action='store_true', default=False, help="Run a selftest")
parser.add_argument('--cost-factor', type=int, default=None, action='append', help="Multiply the maximum work per step thresholds by these number. Can be provided multiple times in prove mode.")
parser.add_argument('number', type=int, nargs=1, help="The maximal modulus to consider, or the number of iterations in --threshold mode")
args = parser.parse_args()

cost_factors = sorted(args.cost_factor) if args.cost_factor is not None else [1]
#print("COST FACTORS: %s" % cost_factors)

if args.self_test:
    real_half_cum=0
    real_full_cum=0
    for M in range(3, 30000, 2):
        real_half = max(count_divsteps(M, x) for x in range(0, M//2))
        real_full = max(real_half, max(count_divsteps(M, x) for x in range(M//2, M)))
        real_half_cum = max(real_half_cum, real_half)
        real_full_cum = max(real_full_cum, real_full)
        calc_half = countsteps(minmod=M, maxmod=M, half_range=True)
        calc_full = countsteps(minmod=M, maxmod=M, half_range=False)
        calc_half_cum = countsteps(minmod=1, maxmod=M, half_range=True)
        calc_full_cum = countsteps(minmid=1, maxmod=M, half_range=False)
        print("M=%i half=%i(%i) halfcum=%i(%i) full=%i(%i) fullcum=%i(%i)" % (M, real_half, calc_half, real_half_cum, calc_half_cum, real_full, calc_full, real_full_cum, calc_full_cum))
        assert real_half <= calc_half
        assert real_full <= calc_full
        assert real_half_cum <= calc_half_cum
        assert real_full_cum <= calc_full_cum

LOGLEVEL = args.verbose
comment = "  # "
if args.exact_modulus:
    comment += "exact-modulus "
if args.half_range:
    comment += "half-range "
if comment == "  # ":
    comment = ""
if args.threshold:
    assert args.prove is None
    assert not args.exact_modulus
    assert len(cost_factors) == 1
    goal = args.number[0]
    bits = (goal + 1.3) / 2.828339631
    shift = int(max(0, bits - 64))
    bits -= shift
    proc = (int(2.0**bits) << shift) | 1
    high, highsteps = None, None
    low, lowsteps = None, None
    while True:
        steps = countsteps(minmod=1, maxmod=proc, half_range=args.half_range, cost_factor=cost_factors[0])
        if steps <= goal:
            low, lowsteps = proc, steps
        else:
            high, highsteps = proc, steps
        if low is None:
            proc = min(high-2, ((high*2) // 3) | 1)
        elif high is None:
            proc = max(low+2, ((low*3) // 2) | 1)
        elif low + 2 == high:
            print("%i %i %s" % (low, lowsteps, comment))
#            print("%i %i %s" % (high, highsteps, comment))
            exit()
        else:
            proc = ((low + high) // 2) | 1
elif args.prove is not None:
    assert len(cost_factors) >= 1
    modulus = args.number[0]
    todo = [[]]
    prog = 1 << modulus.bit_length()
    while len(todo):
        fixed_bits = todo.pop(0)
        prune=False
        for cost_factor in cost_factors:
            steps = countsteps(minmod=modulus if args.exact_modulus else 1, maxmod=modulus, half_range=args.half_range, fixed_bits=fixed_bits, cost_factor=cost_factor, prove_limit=args.prove[0])
            print("%i %i %s  # fixed_bits=%s progress=%f cost_factor=%i queue=%i" % (modulus, steps, comment, fixed_bits, modulus.bit_length() - math.log2(prog), cost_factor, len(todo)))
            if steps <= args.prove[0]:
                prune=True
                prog -= (1 << (modulus.bit_length() - len(fixed_bits)))
                break
        if not prune:
            if len(fixed_bits) >= modulus.bit_length():
                print("Disproven")
                exit()
            bit = random.getrandbits(1)
            todo.append(fixed_bits + [bit])
            todo.append(fixed_bits + [bit ^ 1])
    print("Proven")
else:
    assert len(cost_factors) == 1
    modulus = args.number[0]
    steps = countsteps(minmod=modulus if args.exact_modulus else 1, maxmod=modulus, half_range=args.half_range, cost_factor=cost_factors[0])
    print("%i %i %s" % (modulus, steps, comment))

hull-bound.py

import argparse, sys, random, math

LOGLEVEL = 0
MAX_BOUNDING_BOX_COST = 16384
MAX_HORIZONTAL_BAND_COST = 128
MAX_DIRECT_COMPUTE_SET_SIZE = 1048576

def line_eq(p1, p2):
    xd = p2[0] - p1[0]
    yd = p2[1] - p1[1]
    return (-yd, xd, yd * p1[0] - xd * p1[1])

def cross(o, a, b):
    """2D cross product of the OA and OB vectors, i.e. the z-component of
    their 3D cross product.

    Returns a positive value if OAB makes a counter-clockwise turn,
    negative for clockwise turn, and zero if the points are collinear.
    """
    return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])

def gmin(a, b):
    if a is None:
        return b
    return min(a, b)

def gmax(a, b):
    if a is None:
        return b
    return max(a, b)

def convex_hull(points):
    """Computes the convex hull of a set of 2D points.

    Input: an iterable sequence of (x, y) pairs representing the points.
    Output: a list of vertices of the convex hull in counter-clockwise order,
      starting from the vertex with the lexicographically smallest coordinates.
    Implements Andrew's monotone chain algorithm. O(n log n) complexity.
    """
    # Boring case: no points.
    if len(points) <= 1:
        return points
    # Sort the points lexicographically (tuples are compared lexicographically).
    points = sorted(points)
    # Build lower hull.
    lower = [points[0]]
    for p in points[1:]:
        while len(lower) >= 2 and cross(lower[-2], lower[-1], p) <= 0:
            lower.pop()
        if lower[-1] != p:
            lower.append(p)
    if (len(lower) == 1):
        return lower
    # Build upper hull.
    upper = [points[-1]]
    for p in reversed(points[:-1]):
        while len(upper) >= 2 and cross(upper[-2], upper[-1], p) <= 0:
            upper.pop()
        if upper[-1] != p:
            upper.append(p)
    # Concatenation of the lower and upper hulls gives the convex hull.
    # The Last point of each list is omitted because it is repeated at the beginning
    # of the other list. 
    return lower[:-1] + upper[:-1]

def in_convex_hull(hull, point):
    """Determine if a point is in a convex hull (returned by convex_hull)."""
    if len(hull) == 0:
        # Empty set.
        return False
    elif len(hull) == 1:
        # Single point.
        return hull[0] == point
    elif len(hull) == 2:
        # Two points: test for colinearity, plus whether the point is between the vertices.
        return cross(hull[0], hull[1], point) == 0 and hull[0][0] <= point[0] <= hull[1][0] and (hull[0][0] != hull[1][0] or hull[0][1] <= point[1] <= hull[1][1])
    else:
        # Three or more points: point needs to be left of every edge, when traversing in 
        # counter-clockwise order.
        if cross(hull[-1], hull[0], point) < 0:
            return False
        for n in range(1, len(hull)):
            if cross(hull[n-1], hull[n], point) < 0:
                return False
        return True

def convex_hull_intersect_y_line(hull, y):
    if len(hull) == 0:
        return []
    elif len(hull) == 1:
        if hull[0][1] == y:
            return hull
        else:
            return []
    elif len(hull) == 2:
        if hull[0][1] == y and hull[1][1] == y:
            return hull
        elif hull[0][1] == hull[1][1]:
            return []
        else:
            if hull[0][0] == hull[1][0] and (y < hull[0][1] or y > hull[1][1]):
                return []
            a, b, c = line_eq(hull[0], hull[1])
            n = -c - b*y
            r = n // a
            if r & 1 and r * a == n and hull[0][0] <= r <= hull[1][0]:
                return [(r, y)]
            else:
                return []
    else:
        min_x = None
        max_x = None
        for i in range(len(hull)):
            i2 = 0 if i + 1 == len(hull) else i + 1
            p1 = hull[i]
            p2 = hull[i2]
            if p1[1] == p2[1]:
                if p1[1] == y:
                    min_x = gmax(min_x, min(hull[i][0], hull[i2][0]))
                    max_x = gmin(max_x, max(hull[i][0], hull[i2][0]))
                elif (p1[0] < p2[0] and y < p1[1]) or (p1[0] > p2[0] and y > p1[1]):
                    return []
                else:
                    pass
            else:
                a, b, c = line_eq(p1, p2)
                n = -c - b*y
                if a > 0:
                    min_x = gmax(min_x, (n + a - 1) // a)
                else:
                    max_x = gmin(max_x, n // a)
        if min_x & 1 != 1:
            min_x += 1
        if max_x & 1 != 1:
            max_x -= 1
        if min_x > max_x:
            return []
        elif min_x == max_x:
            return [(min_x, y)]
        else:
            return [(min_x, y), (max_x, y)]

#for N in range(10000):
#    lst = []
#    for _ in range(random.randrange(0, 17)):
#        lst.append((2*random.randrange(-4, 4) + 1, random.randrange(-8, 8)))
#    hull = convex_hull(lst)
#    for p in lst:
#        i = in_convex_hull(hull, p)
#        if not i:
#            print(lst, hull, p)
#        assert in_convex_hull(hull, p)
#    for g in range(-9, 9):
#        intersection = convex_hull_intersect_y_line(hull, g)
#        assert all(f & 1 for f, _ in intersection)
#        assert all(gg == g for _, gg in intersection)
#        assert len(intersection) <= 2
#        for f in range(-11, 12, 2):
#             expect = len(intersection) and f >= intersection[0][0] and f <= intersection[-1][0]
#             real = in_convex_hull(hull, (f, g))
#             assert expect == in_convex_hull(hull, (f, g))
#
#h1 = convex_hull([(0,0),(3,0),(0,3),(3,3)])
#assert convex_hull_intersect_y_line(h1, -1) == []
#assert convex_hull_intersect_y_line(h1, 0) == [(0,0), (3,0)]
#assert convex_hull_intersect_y_line(h1, 1) == [(0,1), (3,1)]
#assert convex_hull_intersect_y_line(h1, 3) == [(0,3), (3,3)]
#assert convex_hull_intersect_y_line(h1, 4) == []
#h2 = convex_hull([(1,0), (0,1), (1,2), (2,1)])
#assert convex_hull_intersect_y_line(h2, -1) == []
#assert convex_hull_intersect_y_line(h2, 0) == [(1,0)]
#assert convex_hull_intersect_y_line(h2, 1) == [(0,1), (2,1)]
#assert convex_hull_intersect_y_line(h2, 2) == [(1,2)]
#assert convex_hull_intersect_y_line(h2, 3) == []
#h3 = convex_hull([(0,2), (1,0), (2,1), (3,3)])
#assert convex_hull_intersect_y_line(h3, -1) == []
#assert convex_hull_intersect_y_line(h3, 0) == [(1,0)]
#assert convex_hull_intersect_y_line(h3, 1) == [(1,1), (2,1)]
#assert convex_hull_intersect_y_line(h3, 2) == [(0,2), (2,2)]
#assert convex_hull_intersect_y_line(h3, 3) == [(3,3)]
#assert convex_hull_intersect_y_line(h3, 4) == []
#h4 = convex_hull([(0,2), (1,0), (3,3)])
#assert convex_hull_intersect_y_line(h4, -1) == []
#assert convex_hull_intersect_y_line(h4, 0) == [(1,0)]
#assert convex_hull_intersect_y_line(h4, 1) == [(1,1)]
#assert convex_hull_intersect_y_line(h4, 2) == [(0,2), (2,2)]
#assert convex_hull_intersect_y_line(h4, 3) == [(3,3)]
#assert convex_hull_intersect_y_line(h4, 4) == []

def count_divsteps(f, g, delta=0):
    """Compute the number of divsteps safegcd needs with a given input."""
    assert f & 1
    it = 0
    while g != 0:
        it += 1
        if g & 1:
            if delta >= 2:
                delta, f, g = -delta, g, (g-f)>>1
            else:
                delta, g = 1+delta, (g+f)>>1
        else:
            delta, g = 1+delta, g>>1
    return it

def divstep_parity_subset(points_hull, delta, g_parity, cost_factor=1):
    """Minimize a set of points defined by points_hull to one that includes all its odd/even g points.
    """
    assert(len(points_hull))
    # Compute lower and upper bound of any g in the set.
    min_g = min(g for _, g in points_hull)
    max_g = max(g for _, g in points_hull)
    # Correct the lower and upper bound to be one that matches the expected g parity.
    corrected_min_g = min_g + ((min_g & 1) != g_parity)
    corrected_max_g = max_g - ((max_g & 1) != g_parity)
    # If the set contains no points with the expected g parity, there is nothing to be done.
    if corrected_min_g > corrected_max_g:
        return (None, None)
    # Compute the lower and upper bound of any f in the set.
    min_f = min(f for f, _ in points_hull)
    max_f = max(f for f, _ in points_hull)

    if len(points_hull) == 1:
        # Fast path for a single point: either it's good or not.
        if points_hull[0][1] & 1 == g_parity:
            out = points_hull
        else:
            return (None, None)
    elif len(points_hull) == 2:
        # In case the hull is a line, use exact logic to trim the vertices with incorrect parity.
        p1 = points_hull[0]
        p2 = points_hull[1]
        # Count the number of times the line goes through points with integer coordinates and odd f.
        m = math.gcd((p2[0] - p1[0]) >> 1, p2[1] - p1[1])
        # Measure the length in f and g dimension of each integer-point crossing segment.
        df = (p2[0] - p1[0]) // m
        dg = (p2[1] - p1[1]) // m
        # Chop off the segment ending in p1 if needed.
        if (p1[1] & 1) != g_parity:
            if (dg & 1) == 0:
                # If the integer-crossing segment has even g length there are no points on the line
                # with correct g parity.
                return (None, 0)
            p1 = (p1[0] + df, p1[1] + dg)
            m -= 1
        # Chop off the segment ending in p2 if needed.
        if (p2[1] & 1) != g_parity:
            p2 = (p2[0] - df, p2[1] - dg)
            m -= 1
        # If the resulting number of odd-f correct-g-parity points on the line is sufficiently small, evaluate
        # directly in all of them.
        if dg & 1:
            m >>= 1
            df <<= 1
            dg <<= 1
        if m + 1 <= MAX_DIRECT_COMPUTE_SET_SIZE * cost_factor:
            iter_pruned = count_divsteps(p2[0], p2[1], delta)
            while p1 != p2:
                iter_pruned = max(iter_pruned, count_divsteps(p1[0], p1[1], delta))
                p1 = (p1[0] + df, p1[1] + dg)
            return (None, iter_pruned)
        # Otherwise continue with just the (corrected) vertices.
        out = [p1, p2]
    elif (((corrected_max_g - corrected_min_g) >> 1) + 1) * len(points_hull) <= MAX_HORIZONTAL_BAND_COST * cost_factor:
        intersections = []
        num_points = 0
        for y in range(corrected_min_g, corrected_max_g + 1, 2):
            intersection = convex_hull_intersect_y_line(points_hull, y)
            if len(intersection) == 1:
                intersections.append(intersection)
                num_points += 1
            elif len(intersection) == 2:
                intersections.append(intersection)
                num_points += (intersection[1][0] - intersection[0][0] + 1)
        if num_points <= MAX_DIRECT_COMPUTE_SET_SIZE * cost_factor:
            iter_pruned = 0
            for intersection in intersections:
                if len(intersection) == 1:
                    iter_pruned = max(iter_pruned, count_divsteps(intersection[0][0], intersection[0][1], delta))
                else:
                    for f in range(intersection[0][0], intersection[1][0] + 1, 2):
                        iter_pruned = max(iter_pruned, count_divsteps(f, intersection[0][1], delta))
            return (None, iter_pruned)
        out = sum(intersections, start=[])
    elif (((corrected_max_g - corrected_min_g)>>1) + 1) * (((max_f - min_f)>>1) + 1) * len(points_hull) <= MAX_BOUNDING_BOX_COST * cost_factor:
        # If the rectangular bounding box (min_f..max_f, corrected_min_g..corrected_max_g) is sufficiently
        # small, just iterate over all its points with the expected g parity (and odd f) that are in the hull.
        out = []
        for f in range(min_f, max_f + 1, 2):
            for g in range(corrected_min_g, corrected_max_g + 1, 2):
                if in_convex_hull(points_hull, (f, g)):
                    out.append((f, g))
        if len(out) == 0:
            return (None, None)
        # If the resulting set is small enough, just count the divsteps needed explicitly for all of them.
        # Note that this must happen before turning the set back into a hull, as the worst-case step
        # count may be in the middle of the hull somewhere.
        if len(out) <= MAX_DIRECT_COMPUTE_SET_SIZE * cost_factor:
            iter_pruned = max(count_divsteps(f, g, delta) for f, g in out)
            return (None, iter_pruned)
    else:
        # The bounding rectangle of the hull is too large; just replace all vertices with mismatching
        # g parity with both (f, g-1) and (f, g+1). This is a strict overestimate, but very simple to
        # compute.
        out = []
        for f, g in points_hull:
            if (g & 1) == g_parity:
                out.append((f, g))
            else:
                out.append((f, g-1))
                out.append((f, g+1))

    # Minimize the computed set again.
    out = convex_hull(out)
    if len(out) <= 2 and all(g == 0 for _, g in out):
        # If the entire set has g=0, we're done.
        return (None, 0)
    else:
        return (out, None)

def process_divstep(state, step, cost_factor=1, fixed_bit=None, prove_limit=None):
    new_state = dict()
    max_iter = 0
    if LOGLEVEL > 0:
        print("- step %i: %i distinct deltas, %i points in total, max %i points per delta" % (step, len(state), sum(len(x) for x in state.values()), max(len(x) for x in state.values())), file=sys.stderr)
    for delta, points in state.items():
        point_hull = convex_hull(points)
        if fixed_bit is None or not fixed_bit:
            if LOGLEVEL > 1:
                print("  - delta %i: %i points -> %i point hull" % (delta, len(points), len(point_hull)), file=sys.stderr)
            even_points, iter_pruned = divstep_parity_subset(points_hull=point_hull, delta=delta, g_parity=0, cost_factor=cost_factor)
            if iter_pruned is not None:
                if LOGLEVEL > 2:
                    print("    - even points were pruned (needing %i steps)" % iter_pruned, file=sys.stderr)
                max_iter = max(max_iter, iter_pruned + step)
                if prove_limit is not None and max_iter > prove_limit:
                    return (None, max_iter)
            if even_points is not None:
                assert len(even_points)
                if LOGLEVEL > 2:
                    print("    - processing %i point even-subset hull" % len(even_points), file=sys.stderr)
                # divsteps^n(delta,f,g) = divsteps^{n-1}(1+delta,f,g//2)
                k = 1+delta
                if k not in new_state:
                    new_state[k] = []
                new_set = new_state[k]
                for f, g in even_points:
                    assert (f & 1) == 1
                    assert (g & 1) == 0
                    new_set.append((f, g>>1))

        if fixed_bit is None or fixed_bit:
            odd_points, iter_pruned = divstep_parity_subset(points_hull=point_hull, delta=delta, g_parity=1, cost_factor=cost_factor)
            if iter_pruned is not None:
                if LOGLEVEL > 2:
                    print("    - odd points were pruned (needing %i steps)" % iter_pruned, file=sys.stderr)
                max_iter = max(max_iter, iter_pruned + step)
            if odd_points is not None:
                assert len(odd_points)
                if LOGLEVEL > 2:
                    print("    - processing %i point odd-subset hull" % len(odd_points), file=sys.stderr)
                if delta >= 2:
                    # divsteps^n(delta,f,g) = divsteps^{n-1}(1-delta,g,(g-f)//2)
                    k = -delta
                    if k not in new_state:
                        new_state[k] = []
                    new_set = new_state[k]
                    for f, g in odd_points:
                        assert (f & 1) == 1
                        assert (g & 1) == 1
                        new_set.append((g, (g-f)>>1))
                else:
                    # divsteps^n(delta,f,g) = divsteps^{n-1}(1+delta,f,(f+g)//2)
                    k = 1+delta
                    if k not in new_state:
                        new_state[k] = []
                    new_set = new_state[k]
                    for f, g in odd_points:
                        assert (f & 1) == 1
                        assert (g & 1) == 1
                        new_set.append((f, (f+g)>>1))

        if prove_limit is not None and max_iter > prove_limit:
            return ({}, max_iter)
    
    return (new_state, max_iter)


def countsteps(minmod, maxmod, half_range, cost_factor=1, fixed_bits=None, prove_limit=None):
    assert minmod & 1
    assert maxmod & 1
    if half_range:
        state = {0: [(minmod, 0), (minmod, minmod//2), (maxmod, 0), (maxmod, maxmod//2)]}
    else:
        state = {0: [(minmod, 0), (minmod, minmod-1), (maxmod, 0), (maxmod, maxmod-1)]}
    steps = 0
    max_iter = 0
    while len(state):
        fixed_bit = None
        if fixed_bits is not None and len(fixed_bits) > steps:
            fixed_bit = fixed_bits[steps]
        state, new_max_iter = process_divstep(state=state, step=steps, cost_factor=cost_factor, fixed_bit=fixed_bit, prove_limit=prove_limit)
        max_iter = max(max_iter, new_max_iter)
        if prove_limit is not None and max_iter > prove_limit:
            return max_iter
        steps += 1
    return max_iter


parser = argparse.ArgumentParser(description="Analyse the maximum number of divsteps needed in the safegcd algorithm.")
parser.add_argument('--verbose', '-v', action='count', default=0, help="Increase verbosity (can be repeated up to 3 times)")
parser.add_argument('--threshold', '-t', action='store_true', default=False, help="Find the largest modulus that can be processed with a given number of divsteps")
parser.add_argument('--exact-modulus', '-e', action='store_true', default=False, help="Make computations for the exact modulus specified, rather than assuming any modulus up to that number")
parser.add_argument('--half-range', action='store_true', default=False, help="Assume inputs in range 0..MODULUS/2 rather than 0..MODULUS-1")
parser.add_argument('--prove', type=int, nargs=1, default=None, help="Recursive split the range in order to prove the specified bound")
parser.add_argument('--self-test', action='store_true', default=False, help="Run a selftest")
parser.add_argument('--cost-factor', type=int, default=None, action='append', help="Multiply the maximum work per step thresholds by these number. Can be provided multiple times in prove mode.")
parser.add_argument('number', type=int, nargs=1, help="The maximal modulus to consider, or the number of iterations in --threshold mode")
args = parser.parse_args()

cost_factors = sorted(args.cost_factor) if args.cost_factor is not None else [1]
print("COST FACTORS: %s" % cost_factors)

if args.self_test:
    real_half_cum=0
    real_full_cum=0
    for M in range(3, 30000, 2):
        real_half = max(count_divsteps(M, x) for x in range(0, M//2))
        real_full = max(real_half, max(count_divsteps(M, x) for x in range(M//2, M)))
        real_half_cum = max(real_half_cum, real_half)
        real_full_cum = max(real_full_cum, real_full)
        calc_half = countsteps(minmod=M, maxmod=M, half_range=True)
        calc_full = countsteps(minmod=M, maxmod=M, half_range=False)
        calc_half_cum = countsteps(minmod=1, maxmod=M, half_range=True)
        calc_full_cum = countsteps(minmid=1, maxmod=M, half_range=False)
        print("M=%i half=%i(%i) halfcum=%i(%i) full=%i(%i) fullcum=%i(%i)" % (M, real_half, calc_half, real_half_cum, calc_half_cum, real_full, calc_full, real_full_cum, calc_full_cum))
        assert real_half <= calc_half
        assert real_full <= calc_full
        assert real_half_cum <= calc_half_cum
        assert real_full_cum <= calc_full_cum

LOGLEVEL = args.verbose
comment = "  # "
if args.exact_modulus:
    comment += "exact-modulus "
if args.half_range:
    comment += "half-range "
if comment == "  # ":
    comment = ""
if args.threshold:
    assert args.prove is None
    assert not args.exact_modulus
    assert len(cost_factors) == 1
    goal = args.number[0]
    bits = (goal + 1.3) / 2.828339631
    shift = int(max(0, bits - 64))
    bits -= shift
    proc = (int(2.0**bits) << shift) | 1
    high, highsteps = None, None
    low, lowsteps = None, None
    while True:
        steps = countsteps(minmod=1, maxmod=proc, half_range=args.half_range, cost_factor=cost_factors[0])
        if steps <= goal:
            low, lowsteps = proc, steps
        else:
            high, highsteps = proc, steps
        if low is None:
            proc = min(high-2, ((high*2) // 3) | 1)
        elif high is None:
            proc = max(low+2, ((low*3) // 2) | 1)
        elif low + 2 == high:
            print("%i %i %s" % (low, lowsteps, comment))
            print("%i %i %s" % (high, highsteps, comment))
            exit()
        else:
            proc = ((low + high) // 2) | 1
elif args.prove is not None:
    assert len(cost_factors) >= 1
    modulus = args.number[0]
    todo = [[]]
    prog = 1 << modulus.bit_length()
    while len(todo):
        fixed_bits = todo.pop(0)
        prune=False
        for cost_factor in cost_factors:
            steps = countsteps(minmod=modulus if args.exact_modulus else 1, maxmod=modulus, half_range=args.half_range, fixed_bits=fixed_bits, cost_factor=cost_factor, prove_limit=args.prove[0])
            print("%i %i %s  # fixed_bits=%s progress=%f cost_factor=%i queue=%i" % (modulus, steps, comment, fixed_bits, modulus.bit_length() - math.log2(prog), cost_factor, len(todo)))
            if steps <= args.prove[0]:
                prune=True
                prog -= (1 << (modulus.bit_length() - len(fixed_bits)))
                break
        if not prune:
            if len(fixed_bits) >= modulus.bit_length():
                print("Disproven")
                exit()
            bit = random.getrandbits(1)
            todo.append(fixed_bits + [bit])
            todo.append(fixed_bits + [bit ^ 1])
    print("Proven")
else:
    assert len(cost_factors) == 1
    modulus = args.number[0]
    steps = countsteps(minmod=modulus if args.exact_modulus else 1, maxmod=modulus, half_range=args.half_range, cost_factor=cost_factors[0])
    print("%i %i %s" % (modulus, steps, comment))

hull.py

import functools

def cmp(a, b):
    return (a > b) - (a < b)

def turn(p, q, r):
    return cmp((q[0] - p[0])*(r[1] - p[1]) - (r[0] - p[0])*(q[1] - p[1]), 0)

def _keep_left(hull, r):
    while len(hull) > 1 and turn(hull[-2], hull[-1], r) <= 0:
        hull.pop()
    if not len(hull) or hull[-1] != r:
        hull.append(r)
    return hull


def convex_hull(points):
    '''
    Returns points on convex hull in CCW order according to Graham's scan algorithm. 
    By Tom Switzer <thomas.switzer@gmail.com>.
    '''

    points = sorted(points)
    l = functools.reduce(_keep_left, points, [])
    u = functools.reduce(_keep_left, reversed(points), [])
    return l.extend(u[i] for i in range(1, len(u) - 1)) or l


def divstep(d):
    nd = dict()
    for (modf, modg, delta), de in d.items():
        # g is odd
        if delta > 0:
            # divsteps^n(delta,f,g) = divsteps^{n-1}(1-delta,g,(g-f)//2)
            k = (modg, True, 1-delta)
            if k not in nd:
                nd[k] = []
            nde = nd[k]
            for x, y in de:
                nde.append((y<<1, y-x))
        else:
            # divsteps^n(delta,f,g) = divsteps^{n-1}(1+delta,f,(f+g)//2)
            k = (modf, True, 1+delta)
            if k not in nd:
                nd[k] = []
            nde = nd[k]
            for x, y in de:
                nde.append((x<<1, x+y))
        
        # g is even
        # divsteps^n(delta,f,g) = divsteps^{n-1}(1+delta,f,g//2)
        k = (modf, True, 1+delta)
        if k not in nd:
            nd[k] = []
        nde = nd[k]
        for x, y in de:
            nde.append((x<<1,y))

    return nd

def shiftstep(d, modulus, bits):
    nd = dict()
    for (modf, modg, delta), de in d.items():
        k = (False, False, delta)
        if k not in nd:
            nd[k] = []
        nde = nd[k]
        for x, y in de:
            if modf:
                low_x = (x - ((1 << (bits-1))*modulus)) >> bits
                high_x = (x + (((1 << (bits-1)) - 1)*modulus) + (1 << bits) - 1) >> bits
            else:
                assert x == ((x >> bits) << bits)
                low_x = x >> bits
                high_x = x >> bits
            if modg:
                low_y = (y - ((1 << (bits-1))*modulus)) >> bits
                high_y = (y + (((1 << (bits-1)) - 1)*modulus) + (1 << bits) - 1) >> bits
            else:
                assert y == ((y >> bits) << bits)
                low_y = y >> bits
                high_y = y >> bits
            nde.append((low_x, low_y))
            nde.append((high_x, low_y))
            nde.append((low_x, high_y))
            nde.append((high_x, high_y))

    return nd

def range_x(d):
    low = None
    high = None
    for de in d.values():
        for x, _ in de:
            if low is None:
                low = x
                high = x
            else:
                low = min(low, x)
                high = max(high, x)
    return (low, high)


def compact(d):
    return {k: convex_hull(v) for k, v in d.items()}

x = {(False, False, 1): [(0, 1)]}
for j in range(120):
    for _ in range(4):
        x = compact(divstep(x))
    x = compact(shiftstep(x, 2**256-1, 4))
    rn = range_x(x)
    print("big step %i: 0x%x..0x%x" % (j, rn[0], rn[1]))

hull_bound.cpp

#include <gmpxx.h>
#include <stdio.h>
#include <algorithm>
#include <map>
#include <mutex>
#include <thread>
#include <random>
#include <cmath>

static constexpr int THREADS = 14;

namespace {

struct Point {
    mpz_class x;
    mpz_class y;

    template<typename F, typename G>
    explicit Point(F&& ff, G&& gg) : x(std::forward<F>(ff)), y(std::forward<G>(gg)) {}

    Point() : x(0), y(0) {}

    Point(const Point& a) : x(a.x), y(a.y) {}
    Point(Point&& a) : x(std::move(a.x)), y(std::move(a.y)) {}

    Point& operator=(const Point& a) { x = a.x; y = a.y; return *this; }
    Point& operator=(Point&& a) { x = std::move(a.x); y = std::move(a.y); return *this; }

    friend int cmp(const Point& a, const Point& b) {
        int fc = cmp(a.x, b.x);
        if (fc) return fc;
        return cmp(a.y, b.y);
    }

    friend bool operator<(const Point& a, const Point& b) { return cmp(a, b) < 0; }
    friend bool operator>(const Point& a, const Point& b) { return cmp(a, b) > 0; }
    friend bool operator<=(const Point& a, const Point& b) { return cmp(a, b) <= 0; }
    friend bool operator>=(const Point& a, const Point& b) { return cmp(a, b) >= 0; }
    friend bool operator==(const Point& a, const Point& b) { return a.x == b.x && a.y == b.y; }
    friend bool operator!=(const Point& a, const Point& b) { return a.x != b.x || a.y != b.y; }

    std::string to_string() const { return "(" + x.get_str() + ", " + y.get_str() + ")"; }
};

int Turn(const Point& o, const Point& a, const Point& b, mpz_class& tmp1, mpz_class& tmp2, mpz_class& tmp3) {
    tmp1 = a.x - o.x;
    tmp2 = b.y - o.y;
    tmp1 *= tmp2;
    tmp2 = a.y - o.y;
    tmp3 = b.x - o.x;
    tmp2 *= tmp3;
    return cmp(tmp1, tmp2);
}

void MakeHull(std::vector<Point>& points, mpz_class& tmp1, mpz_class& tmp2, mpz_class& tmp3) {
    if (points.size() <= 1) return;
    std::sort(points.begin(), points.end());
    std::vector<Point> ret;
    ret.reserve(points.size() + 1);
    ret.push_back(points[0]);
    for (auto it = points.begin() + 1; it != points.end(); ++it) {
        while (ret.size() >= 2 && Turn(*(ret.end() - 2), *(ret.end() - 1), *it, tmp1, tmp2, tmp3) <= 0) {
            ret.pop_back();
        }
        if (*it != ret.back()) {
            ret.push_back(*it);
        }
    }
    if (ret.size() == 1) {
        points = std::move(ret);
        return;
    }
    ret.pop_back();
    size_t lower_size = ret.size();
    ret.push_back(*points.rbegin());
    for (auto it = points.rbegin() + 1; it != points.rend(); ++it) {
        while (ret.size() >= 2 + lower_size && Turn(*(ret.end() - 2), *(ret.end() - 1), *it, tmp1, tmp2, tmp3) <= 0) {
            ret.pop_back();
        }
        if (*it != ret.back()) {
            ret.push_back(*it);
        }
    }
    ret.pop_back();
    points = std::move(ret);
}

struct HullData {
    std::vector<Point> points;
    mpz_class max_abs_g;

    std::string to_string() const {
        std::string ret = "[[";
        for (size_t i = 0; i < points.size(); ++i) {
            if (i) ret += ", ";
            ret += points[i].to_string();
        }
        ret += "], " + max_abs_g.get_str() + "]";
        return ret;
    }
};

bool ProcessHull(HullData& out, const std::map<int, HullData>& in, int delta, mpz_class& tmp1, mpz_class& tmp2, mpz_class& tmp3) {
    out.points.clear();
    const mpz_class* max_parent_abs_g_ptr = nullptr;
    std::map<int, HullData>::const_iterator it;
    if (1 - delta >= 0 && (it = in.find(1 - delta), it != in.end())) {
        max_parent_abs_g_ptr = &it->second.max_abs_g;
        for (const Point& point : it->second.points) {
            out.points.emplace_back(point.y << 1, point.y - point.x);
        }
    }
    if ((it = in.find(delta - 1), it != in.end())) {
        if (max_parent_abs_g_ptr == nullptr || it->second.max_abs_g > *max_parent_abs_g_ptr) {
            max_parent_abs_g_ptr = &it->second.max_abs_g;
        }
        for (const Point& point : it->second.points) {
            out.points.emplace_back(point.x << 1, point.y);
        }
        if (delta - 1 < 0) {
            for (const Point& point : it->second.points) {
                out.points.emplace_back(point.x << 1, point.x + point.y);
            }
        }
    }
    if (out.points.size() == 0) return false;
    MakeHull(out.points, tmp1, tmp2, tmp3);
    mpz_class max_parent_abs_g = (*max_parent_abs_g_ptr) << 1;
    const mpz_class* min_g_ptr = nullptr;
    const mpz_class* max_g_ptr = nullptr;
    for (const auto& entry : out.points) {
        if (min_g_ptr == nullptr || entry.y < *min_g_ptr) min_g_ptr = &entry.y;
        if (max_g_ptr == nullptr || entry.y > *max_g_ptr) max_g_ptr = &entry.y;
    }
    mpz_class abs_min_g = abs(*min_g_ptr), abs_max_g = abs(*max_g_ptr);
    if (abs_min_g > abs_max_g) abs_max_g = std::move(abs_min_g);
    if (max_parent_abs_g < abs_max_g) {
        out.max_abs_g = std::move(max_parent_abs_g);
    } else {
        out.max_abs_g = std::move(abs_max_g);
    }
    return true;
}

void ProcessDivStep(std::map<int, HullData>& new_state, const std::map<int, HullData>& old_state, int step) {
    std::vector<int> new_deltas;
    int min_new_delta = std::min({1 + old_state.begin()->first, 1 - old_state.begin()->first, 1 + old_state.rbegin()->first, 1 - old_state.rbegin()->first});
    int max_new_delta = std::max({1 + old_state.begin()->first, 1 - old_state.begin()->first, 1 + old_state.rbegin()->first, 1 - old_state.rbegin()->first});
    for (int new_delta = min_new_delta; new_delta <= max_new_delta; new_delta += 2) {
        new_deltas.push_back(new_delta);
    }
    new_state.clear();
    std::random_device rng;
    std::mutex mutex;
    std::vector<std::thread> threads;
    for (int t = 0; t < THREADS; ++t) {
        threads.emplace_back([&](){
            mpz_class tmp1, tmp2, tmp3;
            while (true) {
                std::vector<int> deltas;
                std::map<int, HullData> new_hulls;
                {
                    std::unique_lock<std::mutex> lock(mutex);
                    if (new_deltas.empty()) return;
                    int count = std::max<int>(new_deltas.size() / (2 * THREADS), 1);
                    while (count--) {
                        std::swap(new_deltas.back(), new_deltas[std::uniform_int_distribution<size_t>(0, new_deltas.size() - 1)(rng)]);
                        deltas.push_back(new_deltas.back());
                        new_deltas.pop_back();
                    }
                }
                for (int delta : deltas) {
                    HullData new_hull;
                    if (ProcessHull(new_hull, old_state, delta, tmp1, tmp2, tmp3)) {
                        new_hulls[delta] = std::move(new_hull);
                    }
                }
                {
                    std::unique_lock<std::mutex> lock(mutex);
                    for (auto& item : new_hulls) {
                        new_state[item.first] = std::move(item.second);
                    }
                }
            }
        });
    }
    for (auto& thread : threads) thread.join();
}

mpz_class MaxAbsG(const std::map<int, HullData>& state) {
    const mpz_class* max_abs_g_ptr = nullptr;
    for (const auto& item : state) {
        if (max_abs_g_ptr == nullptr || (item.second.max_abs_g > *max_abs_g_ptr)) max_abs_g_ptr = &item.second.max_abs_g;
    }
    return *max_abs_g_ptr;
}

}

int main() {
    setbuf(stdout, nullptr);
    std::map<int, HullData> state;
    state[1].points = std::vector<Point>{Point{0,0},Point{1,0},Point{1,1}};
    state[1].max_abs_g = 1;
    int steps = 0;
    while (true) {
        std::map<int, HullData> new_state;
        ProcessDivStep(new_state, state, steps);
        state = std::move(new_state);
        ++steps;
        mpz_class max_abs_g = MaxAbsG(state);
        mpz_class alpha = ((mpz_class{1} << steps) - 1) / max_abs_g;
        int shift = std::max<int>(0, mpz_sizeinbase(alpha.get_mpz_t(), 2) - 128);
        mpz_class shifted_alpha = alpha >> shift;
        printf("%i 0x%s (2**%f)\n", steps, alpha.get_str(16).c_str(), log2(shifted_alpha.get_d()) + shift);
    }
    return 0;
}

hull_bound_minimal.py

#!/bin/env python

#import matplotlib.pyplot as plt

"""Compute ranges of moduli that can be handled with given numbers of divsteps.

Output will include lines like these:

  ...
  721 0x68c99daf6eb26776b8a7a2ace253df30d84c7c4328ac408aa128dc8ca448a777 (2**254.711324)
  722 0x9e4ccd3a0a07119e8a58e4ff434f504ad45434b4af37b913db51296f40fa5a5b (2**255.306518)
  723 0xaa040650e6d9c377b515b9e3d71fbfacc6c7f06f115dd0693166a8a64511ef26 (2**255.409524)
  724 0x1030596cf6d817d1357f908ef70cdb00b38d047fbba852139babb6c8646fb15b2 (2**256.016930)
  ...

indicating that e.g. for all moduli M up to 0x9e4c..5a5b (and 0 <= input <= M), 722
divsteps is sufficient to guarantee g=0.

This code runs significantly faster in Pypy.
"""

import math


def convex_hull(points):
    """Computes the convex hull of a set of 2D points.

    Input: a list of (x, y) pairs representing the points.
    Output: a list of vertices of the convex hull in counter-clockwise order.

    This implements Andrew's monotone chain algorithm, based on
    https://en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain
    """
    # Boring case: no points.
    if len(points) <= 1:
        return points
    # Sort the points lexicographically (tuples are compared lexicographically).
    points = sorted(points)
    # Build lower hull.
    lower = [points[0]]
    for point in points[1:]:
        while len(lower) >= 2 and ((lower[-1][0] - lower[-2][0]) * (point[1] - lower[-2][1]) <=
                                   (lower[-1][1] - lower[-2][1]) * (point[0] - lower[-2][0])):
            lower.pop()
        if lower[-1] != point:
            lower.append(point)
    if len(lower) == 1:
        return lower
    # Build upper hull.
    upper = [points[-1]]
    for point in reversed(points[:-1]):
        while len(upper) >= 2 and ((upper[-1][0] - upper[-2][0]) * (point[1] - upper[-2][1]) <=
                                   (upper[-1][1] - upper[-2][1]) * (point[0] - upper[-2][0])):
            upper.pop()
        if upper[-1] != point:
            upper.append(point)
    # Concatenation of the lower and upper hulls gives the convex hull.
    # The last point of each list is omitted because it is repeated at the beginning
    # of the other list.
    return lower[:-1] + upper[:-1]

SQ=math.sqrt(2.0)
QQ=0.6388567272696165831006765275289

def process_divstep(state, step):
    """Apply one divstep to state.

    state is a dict of the form delta -> [points, abs_g], where points is a set of (f, g)
    coordinates (multiplied by 2**step so that they remain integral) that defines a convex
    hull of possible combinations of f and g. abs_g is the largest abs(g) in points, or
    the largest such value in any of the contributing hulls if that is smaller.

    Every invocation of this function will compute the new state based on the previous state,
    taking into account that not every transaction is possible for every delta, but always
    exploring both the "odd g" and "even g" variants.
    """

    new_state = dict()
    for delta, (points, abs_g) in state.items():
        # Odd g:
        if delta > 0:
            # divsteps^n(delta,f,g) = divsteps^{n-1}(1-delta,g,(g-f)/2)
            new_state.setdefault(1 - delta, [[], 0])
            new_state[1 - delta][0].extend((g << 1, g - f) for f, g in points)
            new_state[1 - delta][1] = max(new_state[1 - delta][1], abs_g << 1)
        else:
            # divsteps^n(delta,f,g) = divsteps^{n-1}(1+delta,f,(f+g)/2)
            new_state.setdefault(1 + delta, [[], 0])
            new_state[1 + delta][0].extend((f << 1, g + f) for f, g in points)
            new_state[1 + delta][1] = max(new_state[1 + delta][1], abs_g << 1)

        # Even g:
        # divsteps^n(delta,f,g) = divsteps^{n-1}(1+delta,f,g/2)
        new_state.setdefault(1 + delta, [[], 0])
        new_state[1 + delta][0].extend((f << 1, g) for f, g in points)
        new_state[1 + delta][1] = max(new_state[1 + delta][1], abs_g << 1)

    for delta in sorted(new_state.keys()):

        # Minimize the set of points by computing a new convex hull over them.
        new_state[delta][0] = convex_hull(new_state[delta][0])
        # Update abs_g values for every hull.
        new_state[delta][1] = min(new_state[delta][1], max(abs(g) for _, g in new_state[delta][0]))

#        print("step %i, delta %i: %s" % (step, delta, new_state[delta]))


    return new_state



# Define the initial state: a single hull for delta=1, defining the triangle f=0..1 g=0..f.
STEP = 0
STATE = {1: [[(0, 0), (1, 0), (1, 1)], 1]}
while True:
#    if STEP % 140 == 0 and 1 in STATE:
#        pp = STATE[1][0] + [STATE[1][0][0]]
#        plt.plot([f * (QQ**STEP) for f, _ in pp], [g * (QQ**STEP) for _, g in pp], '.-', label="delta=1, %i steps" % STEP)
#        print(len(pp)-1)
#        plt.legend(loc='upper left', fontsize='large')
#        plt.grid()
#        plt.show()
    # Compute next state.
    STATE = process_divstep(STATE, STEP)
    STEP += 1
    # Compute the maximum of all abs_g values over all hulls at the current state.
    # For any input f,g, and any valid sequence of divsteps applied to it, abs(g)
    # must at some point have been less than or equal to that value (divided by 2**STEP).
    max_abs_g = max(abs_g for _, abs_g in STATE.values())
    # Thus, if the input triangle would have been alpha=(2**STEP / abs_g)-epsilon times
    # larger, abs(g) would at some point have been < 1. In the actual algorithm, g is an
    # integer, and thus abs(g)<1 implies g=0, which means the computation would have
    # been complete for any input.
    alpha = ((1 << STEP) - 1) // max_abs_g
    # That scaled triangle is f=0..alpha g=0..f, covering all moduli up alpha, and all
    # inputs up to the modulus.
    print("%i 0x%x (2**%f)" % (STEP, alpha, math.log(alpha, 2)))

verify_stable_hull.sage

import functools

# How many bits of precision to use in number field approximations.
CONFIG_PRECISION_BITS = 4096

# What denominator to use in approximate formula
CONFIG_DENOMINATOR = 19929

# Define number field with S = sqrt(D) and A = ((1591853137+3*sqrt(D))/2^55)^(1/54).
print("Defining number field...")
D = 273548757304312537
K.<S> = NumberField(x^2-D)
Kz.<z> = K[]
L.<A> = NumberField(z^54-(1591853137+3*S)/2^55)
S_SYMBOLIC = sqrt(D)
A_SYMBOLIC = ((1591853137+3*S_SYMBOLIC)/2^55)**(1/54)

# Define embedding phi of L into reals with requested precision.
print("Finding embedding...")
RR = RealField(CONFIG_PRECISION_BITS)
for phi in L.embeddings(RR):
    if phi(A) > 0 and phi((A^54)*2^55-1591853137) > 0 and phi(S) > 0:
        break
NUMERATOR = ceil(log(1/2,phi(A)) * CONFIG_DENOMINATOR)

def cmp_point(a, b):
    """Compare two points with coordinates in L."""
    if a[0] != b[0]:
        d = phi(a[0] - b[0])
        if d > 0:
            return 1
        if d < 0:
            return -1
        assert(false)
    if a[1] != b[1]:
        d = phi(a[1] - b[1])
        if d > 0:
            return 1
        if d < 0:
            return -1
        assert(false)
    return 0

# A key function using cmp_point as a backend.
key_point = functools.cmp_to_key(cmp_point)

def convex_hull(points):
    """Computes the convex hull of a set of 2D points.

    Input: a list of (x, y) pairs representing the points.
    Output: a list of vertices of the convex hull in counter-clockwise order.

    This implements Andrew's monotone chain algorithm, based on
    https://en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain
    """
    # Boring case: no points.
    if len(points) <= 1:
        return points
    # Sort the points lexicographically (tuples are compared lexicographically).
    points = sorted(points, key=key_point)
    # Build lower hull.
    lower = [points[0]]
    for point in points[1:]:
        while len(lower) >= 2 and phi((lower[-1][0] - lower[-2][0]) * (point[1] - lower[-2][1]) - (lower[-1][1] - lower[-2][1]) * (point[0] - lower[-2][0])) <= 0:
            lower.pop()
        if lower[-1] != point:
            lower.append(point)
    if len(lower) == 1:
        return lower
    # Build upper hull.
    upper = [points[-1]]
    for point in reversed(points[:-1]):
        while len(upper) >= 2 and phi((upper[-1][0] - upper[-2][0]) * (point[1] - upper[-2][1]) - (upper[-1][1] - upper[-2][1]) * (point[0] - upper[-2][0])) <= 0:
            upper.pop()
        if upper[-1] != point:
            upper.append(point)
    # Concatenation of the lower and upper hulls gives the convex hull.
    # The last point of each list is omitted because it is repeated at the beginning
    # of the other list.
    return lower[:-1] + upper[:-1]

MAT_EVEN = Matrix([[1,0],[0,1/2]]) / A
MAT_ODD_POS = Matrix([[0,1],[-1/2,1/2]]) / A
MAT_ODD_NEG = Matrix([[1,0],[1/2,1/2]]) / A

DELTA_MATS = {}
DELTA_MATS[1/2] = [Matrix([[1,0],[0,1]])]
for delta2 in range(3, 11, 2):
    DELTA_MATS[delta2/2] = [Matrix([[1,0],[0,2**(1/2-delta2/2)]]) * A**(1/2-delta2/2)]
for delta2 in range(-9, -1, 2):
    DELTA_MATS[delta2/2] = [Matrix([[0,2**(delta2/2-1/2)],[-1/2,2**(delta2/2-3/2)]]) * A**(delta2/2-3/2)]
DELTA_MATS[-1/2] = [MAT_EVEN*DELTA_MATS[-3/2][0], MAT_ODD_NEG*DELTA_MATS[-3/2][0], MAT_ODD_POS*DELTA_MATS[3/2][0]]

def verify(points, require, name):
    print("* Analysis for %s:" % name)

    # Construct hulls for delta=-9/2 through delta=9/2
    hulls = {}
    for delta in DELTA_MATS:
        hulls[delta] = convex_hull([tuple(mat*vector(p)) for mat in DELTA_MATS[delta] for p in points])
    print("  * Hull size: %i for delta=-1/2; %i for other deltas" % (len(hulls[-1/2]), len(hulls[1/2])))

    assert convex_hull(hulls[1/2] + require) == hulls[1/2]
    print("  * Verified: contains %s" % require)

    max_abs_g = max([(abs(phi(g)), g*sgn(phi(g))) for _, g in hulls[1/2]])[1]
    print("  * Max |g|: %s (%f)" % (max_abs_g, phi(max_abs_g)))

    for delta in hulls:
        if 1+delta in hulls:
            assert convex_hull([tuple(MAT_EVEN*vector(p)) for p in hulls[delta]] + hulls[delta+1]) == hulls[delta+1]
        if delta > 0 and 1-delta in hulls:
            assert convex_hull([tuple(MAT_ODD_POS*vector(p)) for p in hulls[delta]] + hulls[1-delta]) == hulls[1-delta]
        if delta < 0 and 1+delta in hulls:
            assert convex_hull([tuple(MAT_ODD_NEG*vector(p)) for p in hulls[delta]] + hulls[1+delta]) == hulls[1+delta]
    print("  * Verified stability for delta=-9/2...9/2")

    cterm = ceil(log(phi(1/max_abs_g),phi(A)) * CONFIG_DENOMINATOR)
    print("  * Formula: floor((%i*log2(M) + %i) / %i)" % (NUMERATOR, cterm + CONFIG_DENOMINATOR - 1, CONFIG_DENOMINATOR))

POINTS_TRIH = [((14602479/1099511627776*S - 7748334001375541/1099511627776)*A^15, (5248155/2199023255552*S - 2784764000070745/2199023255552)*A^15), ((-1718383502173803/17927291358695426424832*S - 3328015/65536)*A^15, (-809666612337607/35854582717390852849664*S - 834107/131072)*A^15), ((6590337/274877906944*S - 3496949542445723/274877906944)*A^17, (3294357/549755813888*S - 1748044174949303/549755813888)*A^17), ((-93793201277/547097514608625074*S - 185/2)*A^17, (-14968274485/273548757304312537*S - 17)*A^17), ((-93793201277/1094195029217250148*S - 185/4)*A^15, (-213539397157/4376780116869000592*S - 321/16)*A^15), (2349441345/18014398509481984*S - 1246655191745249755/18014398509481984, 4860195819/72057594037927936*S - 2578905986969811401/72057594037927936), (-513578441/547097514608625074*S - 1/2, -549425803/1094195029217250148*S - 1/4), ((1257108893252316260373/590295810358705651712*S - 657560095225890870616184601143/590295810358705651712)*A^39, (3224422004397901489371/2361183241434822606848*S - 1686609053234009001585839547161/2361183241434822606848)*A^39), ((-38161344812487/547097514608625074*S - 158779/2)*A^39, (-33495837853102/273548757304312537*S - 10262)*A^39), ((441890701100964968991957/19342813113834066795298816*S - 231141226551696595856763237472887/19342813113834066795298816)*A^24, (1194248705831577428763675/77371252455336267181195264*S - 624679609654461465106321162014425/77371252455336267181195264)*A^24), ((-18248880445865/17507120467476002368*S - 44661/64)*A^24, (-64802864674791/70028481869904009472*S - 91099/256)*A^24), ((169815994083840534566115/4835703278458516698824704*S - 88826212144405858523285231252465/4835703278458516698824704)*A^26, (648253709382430664311701/19342813113834066795298816*S - 339084205958673171958256968491191/19342813113834066795298816)*A^26), ((-816134148489/547097514608625074*S - 2261/2)*A^26, (-546824990194/273548757304312537*S - 778)*A^26), ((58010606597823904440661455/158456325028528675187087900672*S - 30343799334589979591382321997704405/158456325028528675187087900672)*A^11, (240457877026750695243912981/633825300114114700748351602688*S - 125777094859665540269489404107599671/633825300114114700748351602688)*A^11), ((-9853115849/547097514608625074*S - 21/2)*A^11, (-5565952754/273548757304312537*S - 10)*A^11), ((21/256*S - 11142971959/256)*A^46, (111/1024*S - 58898566069/1024)*A^46), ((307802625912117/547097514608625074*S - 1817839/2)*A^46, (-306128047637303/273548757304312537*S - 240875)*A^46), ((6741/8388608*S - 3576893998839/8388608)*A^31, (40623/33554432*S - 21555283328117/33554432)*A^31), ((-1338339117871/547097514608625074*S - 9603/2)*A^31, (-5190734070041/547097514608625074*S - 8485/2)*A^31), ((2338581/274877906944*S - 1240892500326199/274877906944)*A^16, (14602479/1099511627776*S - 7748334001375541/1099511627776)*A^16), ((-227179222459049/4481822839673856606208*S - 623477/16384)*A^16, (-1718383502173803/17927291358695426424832*S - 3328015/65536)*A^16), ((823995/68719476736*S - 437226341874105/68719476736)*A^18, (6590337/274877906944*S - 3496949542445723/274877906944)*A^18), ((-33920103337/547097514608625074*S - 117/2)*A^18, (-93793201277/547097514608625074*S - 185/2)*A^18), ((273516027/2251799813685248*S - 145132448533242233/2251799813685248)*A^3, (2349441345/9007199254740992*S - 1246655191745249755/9007199254740992)*A^3), ((-441883717/547097514608625074*S - 1/2)*A^3, (-513578441/273548757304312537*S - 1)*A^3), ((136726168839761822937/147573952589676412928*S - 71517808110915902565678564067/147573952589676412928)*A^42, (1257108893252316260373/295147905179352825856*S - 657560095225890870616184601143/295147905179352825856)*A^42), ((153482668387355/547097514608625074*S - 394241/2)*A^42, (-38161344812487/273548757304312537*S - 158779)*A^42), ((32855849367829369553049/4835703278458516698824704*S - 17186017500157080615992137601059/4835703278458516698824704)*A^27, (441890701100964968991957/9671406556917033397649408*S - 231141226551696595856763237472887/9671406556917033397649408)*A^27), ((2514055834299/4376780116869000592*S - 10721/16)*A^27, (-18248880445865/8753560233738001184*S - 44661/32)*A^27), ((-34701431782727265153339/1208925819614629174706176*S + 18151392381363899097100318683449/1208925819614629174706176)*A^29, (169815994083840534566115/2417851639229258349412352*S - 88826212144405858523285231252465/2417851639229258349412352)*A^29), ((1926197476085/547097514608625074*S - 559/2)*A^29, (-816134148489/273548757304312537*S - 2261)*A^29), ((-8303257154159872740241077/19807040628566084398385987584*S + 4343212106986950186917804764310807/19807040628566084398385987584)*A^14, (58010606597823904440661455/79228162514264337593543950336*S - 30343799334589979591382321997704405/79228162514264337593543950336)*A^14), ((14968274485/547097514608625074*S + 17/2)*A^14, (-9853115849/273548757304312537*S - 21)*A^14), ((-3/16*S + 1591853137/16)*A^49, (21/128*S - 11142971959/128)*A^49), ((3372432258834775/547097514608625074*S - 3526517/2)*A^49, (307802625912117/273548757304312537*S - 1817839)*A^49), ((-1275/524288*S + 676537583225/524288)*A^34, (6741/4194304*S - 3576893998839/4194304)*A^34), ((16747918926551/547097514608625074*S + 5131/2)*A^34, (-1338339117871/273548757304312537*S - 9603)*A^34), ((-474171/17179869184*S + 251603531274809/17179869184)*A^19, (2338581/137438953472*S - 1240892500326199/137438953472)*A^19), ((64802864674791/280113927479616037888*S + 91099/1024)*A^19, (-227179222459049/2240911419836928303104*S - 623477/8192)*A^19), ((-257397/4294967296*S + 136579407301463/4294967296)*A^21, (823995/34359738368*S - 437226341874105/34359738368)*A^21), ((273412495097/547097514608625074*S + 389/2)*A^21, (-33920103337/273548757304312537*S - 117)*A^21), ((-95555829/140737488355328*S + 50703615384095191/140737488355328)*A^6, (273516027/1125899906842624*S - 145132448533242233/1125899906842624)*A^6), ((2782976377/547097514608625074*S + 5/2)*A^6, (-441883717/273548757304312537*S - 1)*A^6), ((-95555829/70368744177664*S + 50703615384095191/70368744177664)*A^8, (88980099/281474976710656*S - 47214416574573521/281474976710656)*A^8), ((2782976377/273548757304312537*S + 5)*A^8, (-4287163095/2188390058434500296*S - 11/8)*A^8), ((-280095681103138609359/36893488147419103232*S + 146510571778743742012626509269/36893488147419103232)*A^43, (136726168839761822937/147573952589676412928*S - 71517808110915902565678564067/147573952589676412928)*A^43), ((306128047637303/547097514608625074*S + 240875/2)*A^43, (153482668387355/547097514608625074*S - 394241/2)*A^43), ((-102258712933283899859727/1208925819614629174706176*S + 53488802262884878810192774967957/1208925819614629174706176)*A^28, (32855849367829369553049/4835703278458516698824704*S - 17186017500157080615992137601059/4835703278458516698824704)*A^28), ((5190734070041/1094195029217250148*S + 8485/4)*A^28, (2514055834299/4376780116869000592*S - 10721/16)*A^28), ((-102258712933283899859727/604462909807314587353088*S + 53488802262884878810192774967957/604462909807314587353088)*A^30, (-34701431782727265153339/1208925819614629174706176*S + 18151392381363899097100318683449/1208925819614629174706176)*A^30), ((5190734070041/547097514608625074*S + 8485/2)*A^30, (1926197476085/547097514608625074*S - 559/2)*A^30), ((-36735986828656825085541711/19807040628566084398385987584*S + 19215613799989435647506810504886101/19807040628566084398385987584)*A^15, (-13231976781521376785111163/39614081257132168796771975168*S + 6921293739298699373099247954679033/39614081257132168796771975168)*A^15), ((1718383502173803/17927291358695426424832*S + 3328015/65536)*A^15, (809666612337607/35854582717390852849664*S + 834107/131072)*A^15), ((-16578465937995944295225633/4951760157141521099596496896*S + 8671752860394232444575031690503803/4951760157141521099596496896)*A^17, (-8303257154159872740241077/9903520314283042199192993792*S + 4343212106986950186917804764310807/9903520314283042199192993792)*A^17), ((93793201277/547097514608625074*S + 185/2)*A^17, (14968274485/273548757304312537*S + 17)*A^17), (1, 1/2), (513578441/547097514608625074*S + 1/2, 549425803/1094195029217250148*S + 1/4), ((-501/32768*S + 265839473879/32768)*A^39, (-1275/131072*S + 676537583225/131072)*A^39), ((38161344812487/547097514608625074*S + 158779/2)*A^39, (33495837853102/273548757304312537*S + 10262)*A^39), ((-175797/1073741824*S + 93281001975063/1073741824)*A^24, (-474171/4294967296*S + 251603531274809/4294967296)*A^24), ((18248880445865/17507120467476002368*S + 44661/64)*A^24, (64802864674791/70028481869904009472*S + 91099/256)*A^24), ((-67587/268435456*S + 35862859323473/268435456)*A^26, (-257397/1073741824*S + 136579407301463/1073741824)*A^26), ((816134148489/547097514608625074*S + 2261/2)*A^26, (546824990194/273548757304312537*S + 778)*A^26), ((-23066991/8796093022208*S + 12239753994833589/8796093022208)*A^11, (-95555829/35184372088832*S + 50703615384095191/35184372088832)*A^11), ((9853115849/547097514608625074*S + 21/2)*A^11, (5565952754/273548757304312537*S + 10)*A^11), ((-52102731242862554037/4611686018427387904*S + 27253547486207455572288134167/4611686018427387904)*A^46, (-280095681103138609359/18446744073709551616*S + 146510571778743742012626509269/18446744073709551616)*A^46), ((-307802625912117/547097514608625074*S + 1817839/2)*A^46, (306128047637303/273548757304312537*S + 240875)*A^46), ((-16889320287639158676597/151115727451828646838272*S + 8834352470380244928273114071127/151115727451828646838272)*A^31, (-102258712933283899859727/604462909807314587353088*S + 53488802262884878810192774967957/604462909807314587353088)*A^31), ((1338339117871/547097514608625074*S + 9603/2)*A^31, (5190734070041/547097514608625074*S + 8485/2)*A^31), ((-5876002511783862075107637/4951760157141521099596496896*S + 3073580015172684068601890637551767/4951760157141521099596496896)*A^16, (-36735986828656825085541711/19807040628566084398385987584*S + 19215613799989435647506810504886101/19807040628566084398385987584)*A^16), ((227179222459049/4481822839673856606208*S + 623477/16384)*A^16, (1718383502173803/17927291358695426424832*S + 3328015/65536)*A^16), ((-2068802195959017888746139/1237940039285380274899124224*S + 1082135188351820564414306731548249/1237940039285380274899124224)*A^18, (-16578465937995944295225633/4951760157141521099596496896*S + 8671752860394232444575031690503803/4951760157141521099596496896)*A^18), ((33920103337/547097514608625074*S + 117/2)*A^18, (93793201277/547097514608625074*S + 185/2)*A^18), ((-687711828538431775806670299/40564819207303340847894502572032*S + 359723694493774133192627747855670809/40564819207303340847894502572032)*A^3, (-5910461518043306451322618593/162259276829213363391578010288128*S + 3091604601235971043889713709288607163/162259276829213363391578010288128)*A^3), ((441883717/547097514608625074*S + 1/2)*A^3, (513578441/273548757304312537*S + 1)*A^3), ((-57/8192*S + 30245209603/8192)*A^42, (-501/16384*S + 265839473879/16384)*A^42), ((-153482668387355/547097514608625074*S + 394241/2)*A^42, (38161344812487/273548757304312537*S + 158779)*A^42), ((-13305/268435456*S + 7059868662595/268435456)*A^27, (-175797/536870912*S + 93281001975063/536870912)*A^27), ((-2514055834299/4376780116869000592*S + 10721/16)*A^27, (18248880445865/8753560233738001184*S + 44661/32)*A^27), ((13659/67108864*S - 7247707332761/67108864)*A^29, (-67587/134217728*S + 35862859323473/134217728)*A^29), ((-1926197476085/547097514608625074*S + 559/2)*A^29, (816134148489/273548757304312537*S + 2261)*A^29), ((3294357/1099511627776*S - 1748044174949303/1099511627776)*A^14, (-23066991/4398046511104*S + 12239753994833589/4398046511104)*A^14), ((-14968274485/547097514608625074*S - 17/2)*A^14, (9853115849/273548757304312537*S + 21)*A^14), ((7736717960909434203/288230376151711744*S - 4046870582507585955985131673/288230376151711744)*A^49, (-52102731242862554037/2305843009213693952*S + 27253547486207455572288134167/2305843009213693952)*A^49), ((-3372432258834775/547097514608625074*S + 3526517/2)*A^49, (-307802625912117/273548757304312537*S + 1817839)*A^49), ((3224422004397901489371/9444732965739290427392*S - 1686609053234009001585839547161/9444732965739290427392)*A^34, (-16889320287639158676597/75557863725914323419136*S + 8834352470380244928273114071127/75557863725914323419136)*A^34), ((-16747918926551/547097514608625074*S - 5131/2)*A^34, (1338339117871/273548757304312537*S + 9603)*A^34), ((1194248705831577428763675/309485009821345068724781056*S - 624679609654461465106321162014425/309485009821345068724781056)*A^19, (-5876002511783862075107637/2475880078570760549798248448*S + 3073580015172684068601890637551767/2475880078570760549798248448)*A^19), ((-64802864674791/280113927479616037888*S - 91099/1024)*A^19, (227179222459049/2240911419836928303104*S + 623477/8192)*A^19), ((648253709382430664311701/77371252455336267181195264*S - 339084205958673171958256968491191/77371252455336267181195264)*A^21, (-2068802195959017888746139/618970019642690137449562112*S + 1082135188351820564414306731548249/618970019642690137449562112)*A^21), ((-273412495097/547097514608625074*S - 389/2)*A^21, (33920103337/273548757304312537*S + 117)*A^21), ((240457877026750695243912981/2535301200456458802993406410752*S - 125777094859665540269489404107599671/2535301200456458802993406410752)*A^6, (-687711828538431775806670299/20282409603651670423947251286016*S + 359723694493774133192627747855670809/20282409603651670423947251286016)*A^6), ((-2782976377/547097514608625074*S - 5/2)*A^6, (441883717/273548757304312537*S + 1)*A^6), ((240457877026750695243912981/1267650600228229401496703205376*S - 125777094859665540269489404107599671/1267650600228229401496703205376)*A^8, (-223626975755840540281378659/5070602400912917605986812821504*S + 116973299817054296461569171874035569/5070602400912917605986812821504)*A^8), ((-2782976377/273548757304312537*S - 5)*A^8, (4287163095/2188390058434500296*S + 11/8)*A^8), ((111/2048*S - 58898566069/2048)*A^43, (-57/8192*S + 30245209603/8192)*A^43), ((-306128047637303/547097514608625074*S - 240875/2)*A^43, (-153482668387355/547097514608625074*S + 394241/2)*A^43), ((40623/67108864*S - 21555283328117/67108864)*A^28, (-13305/268435456*S + 7059868662595/268435456)*A^28), ((-5190734070041/1094195029217250148*S - 8485/4)*A^28, (-2514055834299/4376780116869000592*S + 10721/16)*A^28), ((40623/33554432*S - 21555283328117/33554432)*A^30, (13659/67108864*S - 7247707332761/67108864)*A^30), ((-5190734070041/547097514608625074*S - 8485/2)*A^30, (-1926197476085/547097514608625074*S + 559/2)*A^30)]
POINTS_TRI =  [((1846473450626721320152174965/20282409603651670423947251286016*S - 965840958204476561172733424304297815/20282409603651670423947251286016)*A^4, (665089515621596580624079257/40564819207303340847894502572032*S - 347890566659208443392825736261137187/40564819207303340847894502572032)*A^4), ((-84664131269339/17927291358695426424832*S - 162751/65536)*A^4, (-34422267523319/35854582717390852849664*S - 51115/131072)*A^4), ((833288864836361180826811995/5070602400912917605986812821504*S - 435871155039625342008769097650459545/5070602400912917605986812821504)*A^6, (417352740546254942175642231/10141204801825835211973625643008*S - 218306074588639278069576333232161421/10141204801825835211973625643008)*A^6), ((-4646117901/547097514608625074*S - 9/2)*A^6, (-648136613/273548757304312537*S - 1)*A^6), ((75/512*S - 39796328425/512)*A^45, (39/1024*S - 20694090781/1024)*A^45), ((-4646117901/1094195029217250148*S - 9/4)*A^4, (-9831210805/4376780116869000592*S - 17/16)*A^4), ((75/1024*S - 39796328425/1024)*A^43, (153/4096*S - 81184509987/4096)*A^43), ((-136385605326647/547097514608625074*S - 848107/2)*A^43, (-182411716069134/273548757304312537*S + 61962)*A^43), ((25143/33554432*S - 13341321141197/33554432)*A^28, (64281/134217728*S - 34108637166499/134217728)*A^28), ((-2512664496887/547097514608625074*S - 6571/2)*A^28, (-1243371010574/273548757304312537*S - 1270)*A^28), ((8831607/1099511627776*S - 4686207102567053/1099511627776)*A^13, (23848665/4398046511104*S - 12654524064504035/4398046511104)*A^13), ((-968555423609/17507120467476002368*S - 2053/64)*A^13, (-2936420841239/70028481869904009472*S - 4939/256)*A^13), ((3394545/274877906944*S - 1801205702312555/274877906944)*A^15, (12945591/1099511627776*S - 6869159881222989/1099511627776)*A^15), ((-45280900441/547097514608625074*S - 101/2)*A^15, (-24847514514/273548757304312537*S - 42)*A^15), (1159161045/9007199254740992*S - 615071381923816055/9007199254740992, 4803581751/36028797018963968*S - 2548865559721767629/36028797018963968), (-499747661/547097514608625074*S - 1/2, -536844899/547097514608625074*S - 1/2), ((2618703102965742784887/4722366482869645213696*S - 1369773669566168251844320581517/4722366482869645213696)*A^35, (14078724832947026855925/18889465931478580854784*S - 7364205035461198383123768809175/18889465931478580854784)*A^35), ((-2823734205755/547097514608625074*S - 54943/2)*A^35, (-12459632581479/273548757304312537*S - 16731)*A^35), ((848898991156486958702007/154742504910672534362390528*S - 444036395302140825660159507407437/154742504910672534362390528)*A^20, (5139877741542827385672885/618970019642690137449562112*S - 2688532803577882204434149400096535/618970019642690137449562112)*A^20), ((-118186186975/547097514608625074*S - 371/2)*A^20, (-244061016553/547097514608625074*S - 437/2)*A^20), ((295345983751281184882023927/5070602400912917605986812821504*S - 154487597886317029444976922010790157/5070602400912917605986812821504)*A^5, (1846473450626721320152174965/20282409603651670423947251286016*S - 965840958204476561172733424304297815/20282409603651670423947251286016)*A^5), ((-12560465936505/4481822839673856606208*S - 27909/16384)*A^5, (-84664131269339/17927291358695426424832*S - 162751/65536)*A^5), ((103984031072526559662792441/1267650600228229401496703205376*S - 54391270112746515984798191104574531/1267650600228229401496703205376)*A^7, (833288864836361180826811995/5070602400912917605986812821504*S - 435871155039625342008769097650459545/5070602400912917605986812821504)*A^7), ((-2053571449/547097514608625074*S - 5/2)*A^7, (-4646117901/547097514608625074*S - 9/2)*A^7), ((9/128*S - 4775559411/128)*A^46, (75/512*S - 39796328425/512)*A^46), ((1050136912573131/547097514608625074*S - 3040017/2)*A^46, (-136385605326647/273548757304312537*S - 848107)*A^46), ((2787/8388608*S - 1478831564273/8388608)*A^31, (25143/16777216*S - 13341321141197/16777216)*A^31), ((2408974593931/547097514608625074*S - 9553/2)*A^31, (-2512664496887/273548757304312537*S - 6571)*A^31), ((661539/274877906944*S - 351024310799281/274877906944)*A^16, (8831607/549755813888*S - 4686207102567053/549755813888)*A^16), ((7688642603/4376780116869000592*S - 305/16)*A^16, (-968555423609/8753560233738001184*S - 2053/32)*A^16), ((-690489/68719476736*S + 366385693571331/68719476736)*A^18, (3394545/137438953472*S - 1801205702312555/137438953472)*A^18), ((62937414789/547097514608625074*S + 33/2)*A^18, (-45280900441/273548757304312537*S - 101)*A^18), ((-165762327/1125899906842624*S + 87956426743789933/1125899906842624)*A^3, (1159161045/4503599627370496*S - 615071381923816055/4503599627370496)*A^3), ((648136613/547097514608625074*S + 1/2)*A^3, (-499747661/273548757304312537*S - 1)*A^3), ((-388913470253112406329/295147905179352825856*S + 203430251672668351724425441539/295147905179352825856)*A^38, (2618703102965742784887/2361183241434822606848*S - 1369773669566168251844320581517/2361183241434822606848)*A^38), ((91205858034567/547097514608625074*S - 30981/2)*A^38, (-2823734205755/273548757304312537*S - 54943)*A^38), ((-162073798004585406847929/9671406556917033397649408*S + 84776476104466232965854429867139/9671406556917033397649408)*A^23, (848898991156486958702007/77371252455336267181195264*S - 444036395302140825660159507407437/77371252455336267181195264)*A^23), ((621685505287/547097514608625074*S + 635/2)*A^23, (-118186186975/273548757304312537*S - 371)*A^23), ((-60027218710804860344131449/316912650057057350374175801344*S + 31398635284095342052362666141995459/316912650057057350374175801344)*A^8, (295345983751281184882023927/2535301200456458802993406410752*S - 154487597886317029444976922010790157/2535301200456458802993406410752)*A^8), ((2936420841239/280113927479616037888*S + 4939/1024)*A^8, (-12560465936505/2240911419836928303104*S - 27909/8192)*A^8), ((-32583548226173843864902167/79228162514264337593543950336*S + 17043584043836612128398407771045997/79228162514264337593543950336)*A^10, (103984031072526559662792441/633825300114114700748351602688*S - 54391270112746515984798191104574531/633825300114114700748351602688)*A^10), ((12423757257/547097514608625074*S + 21/2)*A^10, (-2053571449/273548757304312537*S - 5)*A^10), ((-3/8*S + 1591853137/8)*A^49, (9/64*S - 4775559411/64)*A^49), ((4241495580332569/547097514608625074*S - 2335195/2)*A^49, (1050136912573131/273548757304312537*S - 3040017)*A^49), ((-3/4*S + 1591853137/4)*A^51, (3/16*S - 1591853137/16)*A^51), ((4241495580332569/273548757304312537*S - 2335195)*A^51, (21043686181502665/2188390058434500296*S - 50975467/8)*A^51), ((-5589/2097152*S + 2965622394231/2097152)*A^32, (2787/8388608*S - 1478831564273/8388608)*A^32), ((12459632581479/547097514608625074*S + 16731/2)*A^32, (2408974593931/547097514608625074*S - 9553/2)*A^32), ((-2042517/68719476736*S + 1083795697941943/68719476736)*A^17, (661539/274877906944*S - 351024310799281/274877906944)*A^17), ((244061016553/1094195029217250148*S + 437/4)*A^17, (7688642603/4376780116869000592*S - 305/16)*A^17), ((-2042517/34359738368*S + 1083795697941943/34359738368)*A^19, (-690489/68719476736*S + 366385693571331/68719476736)*A^19), ((244061016553/547097514608625074*S + 437/2)*A^19, (62937414789/547097514608625074*S + 33/2)*A^19), ((-733949781/1125899906842624*S + 389446753761770999/1125899906842624)*A^4, (-264121593/2251799813685248*S + 140147595455495747/2251799813685248)*A^4), ((84664131269339/17927291358695426424832*S + 162751/65536)*A^4, (34422267523319/35854582717390852849664*S + 51115/131072)*A^4), ((-331230843/281474976710656*S + 175756952166901497/281474976710656)*A^6, (-165762327/562949953421312*S + 87956426743789933/562949953421312)*A^6), ((4646117901/547097514608625074*S + 9/2)*A^6, (648136613/273548757304312537*S + 1)*A^6), ((4646117901/1094195029217250148*S + 9/4)*A^4, (9831210805/4376780116869000592*S + 17/16)*A^4), ((-187976035826178449451/18446744073709551616*S + 98325245077427287723046626441/18446744073709551616)*A^43, (-388913470253112406329/73786976294838206464*S + 203430251672668351724425441539/73786976294838206464)*A^43), ((136385605326647/547097514608625074*S + 848107/2)*A^43, (182411716069134/273548757304312537*S - 61962)*A^43), ((-63185799322567022846871/604462909807314587353088*S + 33050804462912941164125871079661/604462909807314587353088)*A^28, (-162073798004585406847929/2417851639229258349412352*S + 84776476104466232965854429867139/2417851639229258349412352)*A^28), ((2512664496887/547097514608625074*S + 6571/2)*A^28, (1243371010574/273548757304312537*S + 1270)*A^28), ((-22210825153880377826634711/19807040628566084398385987584*S + 11617889573150773218583724259549101/19807040628566084398385987584)*A^13, (-60027218710804860344131449/79228162514264337593543950336*S + 31398635284095342052362666141995459/79228162514264337593543950336)*A^13), ((968555423609/17507120467476002368*S + 2053/64)*A^13, (2936420841239/70028481869904009472*S + 4939/256)*A^13), ((-8535473706168775220480913/4951760157141521099596496896*S + 4464678384786445507074787429726283/4951760157141521099596496896)*A^15, (-32583548226173843864902167/19807040628566084398385987584*S + 17043584043836612128398407771045997/19807040628566084398385987584)*A^15), ((45280900441/547097514608625074*S + 101/2)*A^15, (24847514514/273548757304312537*S + 42)*A^15), (1, 1), (499747661/547097514608625074*S + 1/2, 536844899/547097514608625074*S + 1/2), ((-1047/262144*S + 555556744813/262144)*A^35, (-5589/1048576*S + 2965622394231/1048576)*A^35), ((2823734205755/547097514608625074*S + 54943/2)*A^35, (12459632581479/273548757304312537*S + 16731)*A^35), ((-338007/8589934592*S + 179352501092653/8589934592)*A^20, (-2042517/34359738368*S + 1083795697941943/34359738368)*A^20), ((118186186975/547097514608625074*S + 371/2)*A^20, (244061016553/547097514608625074*S + 437/2)*A^20), ((-117457047/281474976710656*S + 62324789576568813/281474976710656)*A^5, (-733949781/1125899906842624*S + 389446753761770999/1125899906842624)*A^5), ((12560465936505/4481822839673856606208*S + 27909/16384)*A^5, (84664131269339/17927291358695426424832*S + 162751/65536)*A^5), ((-41367129/70368744177664*S + 21950131355777891/70368744177664)*A^7, (-331230843/281474976710656*S + 175756952166901497/281474976710656)*A^7), ((2053571449/547097514608625074*S + 5/2)*A^7, (4646117901/547097514608625074*S + 9/2)*A^7), ((-21876829653177867753/2305843009213693952*S + 11443185444951688930589304723/2305843009213693952)*A^46, (-187976035826178449451/9223372036854775808*S + 98325245077427287723046626441/9223372036854775808)*A^46), ((-1050136912573131/547097514608625074*S + 3040017/2)*A^46, (136385605326647/273548757304312537*S + 848107)*A^46), ((-6870899990778915423171/151115727451828646838272*S + 3593984321068147631630795842961/151115727451828646838272)*A^31, (-63185799322567022846871/302231454903657293676544*S + 33050804462912941164125871079661/302231454903657293676544)*A^31), ((-2408974593931/547097514608625074*S + 9553/2)*A^31, (2512664496887/273548757304312537*S + 6571)*A^31), ((-1651314187709068283943171/4951760157141521099596496896*S + 863758358839244400847126659162961/4951760157141521099596496896)*A^16, (-22210825153880377826634711/9903520314283042199192993792*S + 11617889573150773218583724259549101/9903520314283042199192993792)*A^16), ((-7688642603/4376780116869000592*S + 305/16)*A^16, (968555423609/8753560233738001184*S + 2053/32)*A^16), ((1744281776916879550864857/1237940039285380274899124224*S - 912387222369318901793511370466787/1237940039285380274899124224)*A^18, (-8535473706168775220480913/2475880078570760549798248448*S + 4464678384786445507074787429726283/2475880078570760549798248448)*A^18), ((-62937414789/547097514608625074*S - 33/2)*A^18, (45280900441/273548757304312537*S + 101)*A^18), ((417352740546254942175642231/20282409603651670423947251286016*S - 218306074588639278069576333232161421/20282409603651670423947251286016)*A^3, (-2915802718799189781131605749/81129638414606681695789005144064*S + 1525178545569862089965500057369676759/81129638414606681695789005144064)*A^3), ((-648136613/547097514608625074*S - 1/2)*A^3, (499747661/273548757304312537*S + 1)*A^3), ((153/16384*S - 81184509987/16384)*A^38, (-1047/131072*S + 555556744813/131072)*A^38), ((-91205858034567/547097514608625074*S + 30981/2)*A^38, (2823734205755/273548757304312537*S + 54943)*A^38), ((64281/536870912*S - 34108637166499/536870912)*A^23, (-338007/4294967296*S + 179352501092653/4294967296)*A^23), ((-621685505287/547097514608625074*S - 635/2)*A^23, (118186186975/273548757304312537*S + 371)*A^23), ((23848665/17592186044416*S - 12654524064504035/17592186044416)*A^8, (-117457047/140737488355328*S + 62324789576568813/140737488355328)*A^8), ((-2936420841239/280113927479616037888*S - 4939/1024)*A^8, (12560465936505/2240911419836928303104*S + 27909/8192)*A^8), ((12945591/4398046511104*S - 6869159881222989/4398046511104)*A^10, (-41367129/35184372088832*S + 21950131355777891/35184372088832)*A^10), ((-12423757257/547097514608625074*S - 21/2)*A^10, (2053571449/273548757304312537*S + 5)*A^10), ((7646596679165302887/144115188075855872*S - 3999730546410763808204919517/144115188075855872)*A^49, (-21876829653177867753/1152921504606846976*S + 11443185444951688930589304723/1152921504606846976)*A^49), ((-4241495580332569/547097514608625074*S + 2335195/2)*A^49, (-1050136912573131/273548757304312537*S + 3040017)*A^49), ((7646596679165302887/72057594037927936*S - 3999730546410763808204919517/72057594037927936)*A^51, (-7115116487006282433/288230376151711744*S + 3721727449270462561192192603/288230376151711744)*A^51), ((-4241495580332569/273548757304312537*S + 2335195)*A^51, (-21043686181502665/2188390058434500296*S + 50975467/8)*A^51), ((14078724832947026855925/37778931862957161709568*S - 7364205035461198383123768809175/37778931862957161709568)*A^32, (-6870899990778915423171/151115727451828646838272*S + 3593984321068147631630795842961/151115727451828646838272)*A^32), ((-12459632581479/547097514608625074*S - 16731/2)*A^32, (-2408974593931/547097514608625074*S + 9553/2)*A^32), ((5139877741542827385672885/1237940039285380274899124224*S - 2688532803577882204434149400096535/1237940039285380274899124224)*A^17, (-1651314187709068283943171/4951760157141521099596496896*S + 863758358839244400847126659162961/4951760157141521099596496896)*A^17), ((-244061016553/1094195029217250148*S - 437/4)*A^17, (-7688642603/4376780116869000592*S + 305/16)*A^17), ((5139877741542827385672885/618970019642690137449562112*S - 2688532803577882204434149400096535/618970019642690137449562112)*A^19, (1744281776916879550864857/1237940039285380274899124224*S - 912387222369318901793511370466787/1237940039285380274899124224)*A^19), ((-244061016553/547097514608625074*S - 437/2)*A^19, (-62937414789/547097514608625074*S - 33/2)*A^19)]
POINTS_SQ =   [((-860416572691442404083/293720741620865866544447488*S - 1476704240119/1073741824)*A^25, (445608937309069749521/587441483241731733088894976*S - 1976456843843/2147483648)*A^25), ((261/4096*S - 138491222919/4096)*A^42, (105/8192*S - 55714859795/8192)*A^42), ((-50847064456796597/8963645679347713212416*S - 74757521/32768)*A^27, (6406140982921043/4481822839673856606208*S - 33782057/16384)*A^27), ((92613/134217728*S - 49142098192327/134217728)*A^27, (47145/268435456*S - 25015972047955/268435456)*A^27), ((92613/268435456*S - 49142098192327/268435456)*A^25, (186903/1073741824*S - 99174042288237/1073741824)*A^25), ((-263353894716287/8963645679347713212416*S - 432019/32768)*A^10, (-25100876787023/2240911419836928303104*S - 72995/8192)*A^10), ((33030789/8796093022208*S - 17526721695745031/8796093022208)*A^10, (68485719/35184372088832*S - 36339735543283501/35184372088832)*A^10), ((-117512179890065308911/8963645679347713212416*S + 115290442685/32768)*A^49, (23272812518801081473/2240911419836928303104*S - 62034341779/8192)*A^49), ((17681400309253310889/288230376151711744*S - 9248668379872891391411340499/288230376151711744)*A^49, (45376161456752145159/1152921504606846976*S - 23735058441352222932079387069/1152921504606846976)*A^49), ((-19941404219504323841/286836661739126822797312*S + 576139539/1048576)*A^34, (7443909988058856945/1147346646956507291189248*S - 115719918307/4194304)*A^34), ((6216001889276217061737/9444732965739290427392*S - 3251424611007345369358345614867/9444732965739290427392)*A^34, (16801555131234771894087/37778931862957161709568*S - 8788444860536297589118076398717/37778931862957161709568)*A^34), ((-1197519426877436385/8963645679347713212416*S + 485789107/32768)*A^36, (26026627338255919/2240911419836928303104*S - 480142205/8192)*A^36), ((2388702287747140125999/2361183241434822606848*S - 1249466384517933842226784619509/2361183241434822606848)*A^36, (9120076222508354351913/9444732965739290427392*S - 4770468351253887637011426387283/9444732965739290427392)*A^36), ((-7156650679964705/8963645679347713212416*S - 5121933/32768)*A^21, (-941488191060817/2240911419836928303104*S - 3071741/8192)*A^21), ((816052743639197548524747/77371252455336267181195264*S - 426855400273616306478571124412777/77371252455336267181195264)*A^21, (3382734704347094901306537/309485009821345068724781056*S - 1769418199373082106264238930611667/309485009821345068724781056)*A^21), ((-14890813620579/8963645679347713212416*S - 22439/32768)*A^2, (-8721664706095/4481822839673856606208*S - 17539/16384)*A^2), ((1158291950616558835600962201/40564819207303340847894502572032*S - 605871591104825681580427953386886691/40564819207303340847894502572032)*A^2, (6227601851890326445007737659/162259276829213363391578010288128*S - 3257492241713080563110505743830304569/162259276829213363391578010288128)*A^2), ((-8799062025554149399/8963645679347713212416*S + 11082799029/32768)*A^41, (-781093389465341681/8963645679347713212416*S - 7196486173/32768)*A^41), ((593481615204152757534461503994537803305/1329227995784915872903807060280344576*S - 310402327363246185733860430307463122530991443539/1329227995784915872903807060280344576)*A^41, (3593446380892697563637967782478275993835/5316911983139663491615228241121378304*S - 1879441740584388850180644842189426980826897938633/5316911983139663491615228241121378304)*A^41), ((-326506377500128038401/73430185405216466636111872*S + 124938150931/268435456)*A^26, (-860416572691442404083/293720741620865866544447488*S - 1476704240119/1073741824)*A^26), ((39/1024*S - 20694090781/1024)*A^43, (261/4096*S - 138491222919/4096)*A^43), ((-76471628388480769/8963645679347713212416*S + 60370707/32768)*A^28, (-50847064456796597/8963645679347713212416*S - 74757521/32768)*A^28), ((11367/33554432*S - 6031531536093/33554432)*A^28, (92613/134217728*S - 49142098192327/134217728)*A^28), ((-388447655556493/8963645679347713212416*S - 128137/32768)*A^13, (-263353894716287/4481822839673856606208*S - 432019/16384)*A^13), ((3825831/1099511627776*S - 2030053692993949/1099511627776)*A^13, (33030789/4398046511104*S - 17526721695745031/4398046511104)*A^13), ((-724901539971013230301/8963645679347713212416*S + 1338420796519/32768)*A^52, (-117512179890065308911/4481822839673856606208*S + 115290442685/16384)*A^52), ((1917009867751946877/72057594037927936*S - 1002736674566612810538658607/72057594037927936)*A^52, (17681400309253310889/144115188075855872*S - 9248668379872891391411340499/144115188075855872)*A^52), ((-16817030661642957117/71709165434781705699328*S + 29362084231/262144)*A^37, (-19941404219504323841/143418330869563411398656*S + 576139539/524288)*A^37), ((461612634148469822781/2361183241434822606848*S - 241457243121434629739240111471/2361183241434822606848)*A^37, (6216001889276217061737/4722366482869645213696*S - 3251424611007345369358345614867/4722366482869645213696)*A^37), ((-4008984318044403859/8963645679347713212416*S + 9139642601/32768)*A^39, (-1197519426877436385/4481822839673856606208*S + 485789107/16384)*A^39), ((-488492339816733493479/590295810358705651712*S + 255517299425021527582768132189/590295810358705651712)*A^39, (2388702287747140125999/1180591620717411303424*S - 1249466384517933842226784619509/1180591620717411303424)*A^39), ((-6406140982921043/8963645679347713212416*S + 33782057/32768)*A^24, (-7156650679964705/4481822839673856606208*S - 5121933/16384)*A^24), ((-116822059178687781966537/9671406556917033397649408*S + 61106499819029148353565694671667/9671406556917033397649408)*A^24, (816052743639197548524747/38685626227668133590597632*S - 426855400273616306478571124412777/38685626227668133590597632)*A^24), ((25100876787023/8963645679347713212416*S + 72995/32768)*A^5, (-14890813620579/4481822839673856606208*S - 22439/16384)*A^5), ((-172045375002540621137803191/2535301200456458802993406410752*S + 89992341774912719898076367729352781/2535301200456458802993406410752)*A^5, (1158291950616558835600962201/20282409603651670423947251286016*S - 605871591104825681580427953386886691/20282409603651670423947251286016)*A^5), ((-23272812518801081473/8963645679347713212416*S + 62034341779/32768)*A^44, (-8799062025554149399/4481822839673856606208*S + 11082799029/16384)*A^44), ((-113312595955014955689661454405916411495/83076749736557242056487941267521536*S + 59264672405915643311191471954189850827120225501/83076749736557242056487941267521536)*A^44, (593481615204152757534461503994537803305/664613997892457936451903530140172288*S - 310402327363246185733860430307463122530991443539/664613997892457936451903530140172288)*A^44), ((-7443909988058856945/4589386587826029164756992*S + 115719918307/16777216)*A^29, (-326506377500128038401/36715092702608233318055936*S + 124938150931/134217728)*A^29), ((-9/64*S + 4775559411/64)*A^46, (39/512*S - 20694090781/512)*A^46), ((-26026627338255919/8963645679347713212416*S + 480142205/32768)*A^31, (-76471628388480769/4481822839673856606208*S + 60370707/16384)*A^31), ((-3657/2097152*S + 1940468974003/2097152)*A^31, (11367/16777216*S - 6031531536093/16777216)*A^31), ((941488191060817/8963645679347713212416*S + 3071741/32768)*A^16, (-388447655556493/4481822839673856606208*S - 128137/16384)*A^16), ((-1347081/68719476736*S + 714785038547699/68719476736)*A^16, (3825831/549755813888*S - 2030053692993949/549755813888)*A^16), ((941488191060817/4481822839673856606208*S + 3071741/16384)*A^18, (-5273674297843071/35854582717390852849664*S + 1021549/131072)*A^18), ((-1347081/34359738368*S + 714785038547699/34359738368)*A^18, (1239375/137438953472*S - 657634327223125/137438953472)*A^18), ((-6875163/549755813888*S + 3648083262978777/549755813888)*A^14, (3484077/2199023255552*S - 1848712967333183/2199023255552)*A^14), ((-254852820410751994657/8963645679347713212416*S + 877259025779/32768)*A^53, (-724901539971013230301/8963645679347713212416*S + 1338420796519/32768)*A^53), ((-3941097610375341003/18014398509481984*S + 2061482926326569645218170473/18014398509481984)*A^53, (1917009867751946877/72057594037927936*S - 1002736674566612810538658607/72057594037927936)*A^53), ((781093389465341681/17927291358695426424832*S + 7196486173/65536)*A^38, (-16817030661642957117/71709165434781705699328*S + 29362084231/262144)*A^38), ((-1438597313781936809739/590295810358705651712*S + 752491841971477684904776375849/590295810358705651712)*A^38, (461612634148469822781/2361183241434822606848*S - 241457243121434629739240111471/2361183241434822606848)*A^38), ((781093389465341681/8963645679347713212416*S + 7196486173/32768)*A^40, (-4008984318044403859/8963645679347713212416*S + 9139642601/32768)*A^40), ((-1438597313781936809739/295147905179352825856*S + 752491841971477684904776375849/295147905179352825856)*A^40, (-488492339816733493479/590295810358705651712*S + 255517299425021527582768132189/590295810358705651712)*A^40), ((860416572691442404083/293720741620865866544447488*S + 1476704240119/1073741824)*A^25, (-445608937309069749521/587441483241731733088894976*S + 1976456843843/2147483648)*A^25), ((-516788307391310453586507/9671406556917033397649408*S + 270318164515324446387735582696937/9671406556917033397649408)*A^25, (-186170407002571649211687/19342813113834066795298816*S + 97380768852999533105273768940317/19342813113834066795298816)*A^25), ((50847064456796597/8963645679347713212416*S + 74757521/32768)*A^27, (-6406140982921043/4481822839673856606208*S + 33782057/16384)*A^27), ((-233218700704471332622821/2417851639229258349412352*S + 121990475023161363708034204771111/2417851639229258349412352)*A^27, (-116822059178687781966537/4835703278458516698824704*S + 61106499819029148353565694671667/4835703278458516698824704)*A^27), ((-233218700704471332622821/4835703278458516698824704*S + 121990475023161363708034204771111/4835703278458516698824704)*A^25, (-466862819061846896555895/19342813113834066795298816*S + 244203474661219660415165594114445/19342813113834066795298816)*A^25), ((263353894716287/8963645679347713212416*S + 432019/32768)*A^10, (25100876787023/2240911419836928303104*S + 72995/8192)*A^10), ((-83146082851193716046172837/158456325028528675187087900672*S + 43491495804983650092406520069764967/158456325028528675187087900672)*A^10, (-172045375002540621137803191/633825300114114700748351602688*S + 89992341774912719898076367729352781/633825300114114700748351602688)*A^10), ((117512179890065308911/8963645679347713212416*S - 115290442685/32768)*A^49, (-23272812518801081473/2240911419836928303104*S + 62034341779/8192)*A^49), ((-44174638197447982076507684900028388425/5192296858534827628530496329220096*S + 23104187485572614315315743891353310834881979315/5192296858534827628530496329220096)*A^49, (-113312595955014955689661454405916411495/20769187434139310514121985316880384*S + 59264672405915643311191471954189850827120225501/20769187434139310514121985316880384)*A^49), ((19941404219504323841/286836661739126822797312*S - 576139539/1048576)*A^34, (-7443909988058856945/1147346646956507291189248*S + 115719918307/4194304)*A^34), ((-3/4*S + 1591853137/4)*A^51, (-9/16*S + 4775559411/16)*A^51), ((1197519426877436385/8963645679347713212416*S - 485789107/32768)*A^36, (-26026627338255919/2240911419836928303104*S + 480142205/8192)*A^36), ((-939/131072*S + 498250031881/131072)*A^36, (-3657/524288*S + 1940468974003/524288)*A^36), ((7156650679964705/8963645679347713212416*S + 5121933/32768)*A^21, (941488191060817/2240911419836928303104*S + 3071741/8192)*A^21), ((-323307/4294967296*S + 171552420721353/4294967296)*A^21, (-1347081/17179869184*S + 714785038547699/17179869184)*A^21), ((-1294905/68719476736*S + 687099528788995/68719476736)*A^17, (-6875163/274877906944*S + 3648083262978777/274877906944)*A^17), ((14890813620579/8963645679347713212416*S + 22439/32768)*A^2, (8721664706095/4481822839673856606208*S + 17539/16384)*A^2), ((-460006905/2251799813685248*S + 244087811588636995/2251799813685248)*A^2, (-2475792219/9007199254740992*S + 1313699203458447001/9007199254740992)*A^2), ((8799062025554149399/8963645679347713212416*S - 11082799029/32768)*A^41, (781093389465341681/8963645679347713212416*S + 7196486173/32768)*A^41), ((-237526243491300829065/73786976294838206464*S + 124243635636614039330502060915/73786976294838206464)*A^41, (-1438597313781936809739/295147905179352825856*S + 752491841971477684904776375849/295147905179352825856)*A^41), ((326506377500128038401/73430185405216466636111872*S - 124938150931/268435456)*A^26, (860416572691442404083/293720741620865866544447488*S + 1476704240119/1073741824)*A^26), ((-82654475097184701093705/2417851639229258349412352*S + 43234348915581228320615453439155/2417851639229258349412352)*A^26, (-516788307391310453586507/9671406556917033397649408*S + 270318164515324446387735582696937/9671406556917033397649408)*A^26), ((76471628388480769/8963645679347713212416*S - 60370707/32768)*A^28, (50847064456796597/8963645679347713212416*S + 74757521/32768)*A^28), ((-29099160381445887664071/604462909807314587353088*S + 15220993801033053838617127524861/604462909807314587353088)*A^28, (-233218700704471332622821/2417851639229258349412352*S + 121990475023161363708034204771111/2417851639229258349412352)*A^28), ((388447655556493/8963645679347713212416*S + 128137/32768)*A^13, (263353894716287/4481822839673856606208*S + 432019/16384)*A^13), ((-9674109193880065875089415/19807040628566084398385987584*S + 5060268205004778797392899059992765/19807040628566084398385987584)*A^13, (-83146082851193716046172837/79228162514264337593543950336*S + 43491495804983650092406520069764967/79228162514264337593543950336)*A^13), ((724901539971013230301/8963645679347713212416*S - 1338420796519/32768)*A^52, (117512179890065308911/4481822839673856606208*S - 115290442685/16384)*A^52), ((-4802829659332247634965400073542188445/1298074214633706907132624082305024*S + 2511972512700549908688939929967520419381428111/1298074214633706907132624082305024)*A^52, (-44174638197447982076507684900028388425/2596148429267413814265248164610048*S + 23104187485572614315315743891353310834881979315/2596148429267413814265248164610048)*A^52), ((16817030661642957117/71709165434781705699328*S - 29362084231/262144)*A^37, (19941404219504323841/143418330869563411398656*S - 576139539/524288)*A^37), (0, 1), ((4008984318044403859/8963645679347713212416*S - 9139642601/32768)*A^39, (1197519426877436385/4481822839673856606208*S - 485789107/16384)*A^39), ((105/16384*S - 55714859795/16384)*A^39, (-939/65536*S + 498250031881/65536)*A^39), ((6406140982921043/8963645679347713212416*S - 33782057/32768)*A^24, (7156650679964705/4481822839673856606208*S + 5121933/16384)*A^24), ((47145/536870912*S - 25015972047955/536870912)*A^24, (-323307/2147483648*S + 171552420721353/2147483648)*A^24), ((186903/4294967296*S - 99174042288237/4294967296)*A^20, (-1294905/34359738368*S + 687099528788995/34359738368)*A^20), ((-25100876787023/8963645679347713212416*S - 72995/32768)*A^5, (14890813620579/4481822839673856606208*S + 22439/16384)*A^5), ((68485719/140737488355328*S - 36339735543283501/140737488355328)*A^5, (-460006905/1125899906842624*S + 244087811588636995/1125899906842624)*A^5), ((23272812518801081473/8963645679347713212416*S - 62034341779/32768)*A^44, (8799062025554149399/4481822839673856606208*S - 11082799029/16384)*A^44), ((45376161456752145159/4611686018427387904*S - 23735058441352222932079387069/4611686018427387904)*A^44, (-237526243491300829065/36893488147419103232*S + 124243635636614039330502060915/36893488147419103232)*A^44), ((7443909988058856945/4589386587826029164756992*S - 115719918307/16777216)*A^29, (326506377500128038401/36715092702608233318055936*S - 124938150931/134217728)*A^29), ((16801555131234771894087/151115727451828646838272*S - 8788444860536297589118076398717/151115727451828646838272)*A^29, (-82654475097184701093705/1208925819614629174706176*S + 43234348915581228320615453439155/1208925819614629174706176)*A^29), ((26026627338255919/8963645679347713212416*S - 480142205/32768)*A^31, (76471628388480769/4481822839673856606208*S - 60370707/16384)*A^31), ((9120076222508354351913/37778931862957161709568*S - 4770468351253887637011426387283/37778931862957161709568)*A^31, (-29099160381445887664071/302231454903657293676544*S + 15220993801033053838617127524861/302231454903657293676544)*A^31), ((-941488191060817/8963645679347713212416*S - 3071741/32768)*A^16, (388447655556493/4481822839673856606208*S + 128137/16384)*A^16), ((3382734704347094901306537/1237940039285380274899124224*S - 1769418199373082106264238930611667/1237940039285380274899124224)*A^16, (-9674109193880065875089415/9903520314283042199192993792*S + 5060268205004778797392899059992765/9903520314283042199192993792)*A^16), ((-941488191060817/4481822839673856606208*S - 3071741/16384)*A^18, (5273674297843071/35854582717390852849664*S - 1021549/131072)*A^18), ((3382734704347094901306537/618970019642690137449562112*S - 1769418199373082106264238930611667/618970019642690137449562112)*A^18, (-3145687244766485486891439/2475880078570760549798248448*S + 1645425002815848345564330064690549/2475880078570760549798248448)*A^18), ((254852820410751994657/8963645679347713212416*S - 877259025779/32768)*A^53, (724901539971013230301/8963645679347713212416*S - 1338420796519/32768)*A^53), ((9842952134528933610385571206621549995/324518553658426726783156020576256*S - 5148053743218016101656700990346447603875137801/324518553658426726783156020576256)*A^53, (-4802829659332247634965400073542188445/1298074214633706907132624082305024*S + 2511972512700549908688939929967520419381428111/1298074214633706907132624082305024)*A^53), ((-781093389465341681/17927291358695426424832*S - 7196486173/65536)*A^38, (16817030661642957117/71709165434781705699328*S - 29362084231/262144)*A^38), ((3593446380892697563637967782478275993835/10633823966279326983230456482242756608*S - 1879441740584388850180644842189426980826897938633/10633823966279326983230456482242756608)*A^38, (-1154406540740524496637724249478026432605/42535295865117307932921825928971026432*S + 603776878321580635690238600270277999421033609679/42535295865117307932921825928971026432)*A^38), ((-781093389465341681/8963645679347713212416*S - 7196486173/32768)*A^40, (4008984318044403859/8963645679347713212416*S - 9139642601/32768)*A^40), ((3593446380892697563637967782478275993835/5316911983139663491615228241121378304*S - 1879441740584388850180644842189426980826897938633/5316911983139663491615228241121378304)*A^40, (1219519920076086533500121766500124780615/10633823966279326983230456482242756608*S - 637832431131404107245203120959574490702932164477/10633823966279326983230456482242756608)*A^40)]
verify(POINTS_TRIH, [(0,0), (1,0), (1,1/2)], "trih")
verify(POINTS_TRI, [(0,0), (1,0), (1,1)], "tri")
verify(POINTS_SQ, [(0,0), (1,0), (1,1), (0,1)], "sq")