Skip to content

cyclegolf

We need to implement a sorting algorithm using a restricted Python subset. It should sort an array as large as possible, within 1e6 cycles.

Different versions I have tried (the capitalized constants need to be replaced with large numbers, see the Python script below):

Start with basic bubble sort:

Version 1, score 268:

def f(m):
    # simple bubble sort
    m[I] = 1
    while m[I] > 0:
        m[I] = 0
        m[J] = 0
        while m[J] < SCORE - 1:
            if m[m[J]] > m[m[J] + 1]:
                m[TEMP] = m[m[J]]
                m[m[J]] = m[m[J] + 1]
                m[m[J] + 1] = m[TEMP]
                m[I] = 1
            m[J] = m[J] + 1

Version 2, score 323:

def f(m):
    # simple bubble sort
    m[I] = 0

    while m[I] < SCORE:
        m[J] = 0
        while m[J] < SCORE - 1 - m[I]:
            if m[m[J]] > m[m[J] + 1]:
                m[TEMP] = m[m[J]]
                m[m[J]] = m[m[J] + 1]
                m[m[J] + 1] = m[TEMP]
            m[J] = m[J] + 1
        m[I] = m[I] + 1

Version 3, score 332:

def f(m):
    # simple bubble sort
    m[I] = SCORE - 1

    while m[I] > 0:
        m[J] = 0
        while m[J] < m[I]:
            if m[m[J]] > m[m[J] + 1]:
                m[TEMP] = m[m[J]]
                m[m[J]] = m[m[J] + 1]
                m[m[J] + 1] = m[TEMP]
            m[J] = m[J] + 1
        m[I] = m[I] - 1

Attempt some faster sorts:

Version 4, score 571:

def f(m):
    # inplace merge sort
    m[SIZE] = 1
    while m[SIZE] < SCORE:
        m[LEFT] = 0
        while m[LEFT] < SCORE:
            m[MID] = m[LEFT] + m[SIZE] - 1
            if m[MID] >= SCORE - 1:
                break

            m[RIGHT] = m[LEFT] + m[SIZE] * 2 - 1
            if m[RIGHT] >= SCORE:
                m[RIGHT] = SCORE - 1

            m[I] = m[LEFT]
            m[J] = m[MID] + 1

            while ((m[I] <= m[MID]) * (m[J] <= m[RIGHT])) > 0:
                if m[m[I]] <= m[m[J]]:
                    m[I] = m[I] + 1
                else:
                    m[TEMP] = m[m[J]]
                    m[K] = m[J]
                    while m[K] > m[I]:
                        m[m[K]] = m[m[K] - 1]
                        m[K] = m[K] - 1
                    m[m[I]] = m[TEMP]
                    m[I] = m[I] + 1
                    m[MID] = m[MID] + 1
                    m[J] = m[J] + 1

            m[LEFT] = m[LEFT] + m[SIZE] * 2
        m[SIZE] = m[SIZE] * 2

Version 5, score 2559:

def f(m):
    # shell sort
    m[GAP] = SCORE // 2
    while m[GAP] > 0:
        m[I] = m[GAP]
        while m[I] < SCORE:
            m[TEMP] = m[m[I]]
            m[J] = m[I]

            m[K] = m[J] - m[GAP]
            while (m[K] >= 0) & (m[m[K]] > m[TEMP]):
                m[m[J]] = m[m[K]]
                m[J] = m[K]
                m[K] = m[J] - m[GAP]

            m[m[J]] = m[TEMP]
            m[I] = m[I] + 1

        m[GAP] = m[GAP] // 2

Version 6, score 2540:

def f(m):
    # heap sort
    m[I] = 0
    while m[I] < SCORE:
        m[J] = m[I]
        while m[J] > 0:
            m[K] = (m[J] - 1) // 2
            if m[m[J]] > m[m[K]]:
                m[TEMP] = m[m[J]]
                m[m[J]] = m[m[K]]
                m[m[K]] = m[TEMP]
                m[J] = m[K]
            else:
                break
        m[I] = m[I] + 1

    m[END] = SCORE - 1
    while m[END] > 0:
        m[TEMP] = m[0]
        m[0] = m[m[END]]
        m[m[END]] = m[TEMP]

        m[J] = 0
        while 1:
            m[LEFT] = 2 * m[J] + 1
            if m[LEFT] >= m[END]:
                break

            # Find larger child
            # m[RIGHT] = m[LEFT] + 1
            m[MID] = m[LEFT] + (m[LEFT] + 1 < m[END]) * (m[m[LEFT] + 1] > m[m[LEFT]])

            if m[m[J]] < m[m[MID]]:
                m[TEMP] = m[m[J]]
                m[m[J]] = m[m[MID]]
                m[m[MID]] = m[TEMP]
                m[J] = m[MID]
            else:
                break

        m[END] = m[END] - 1

Based on the uniform distribution, use flash sort:

Version 7, score 11878:

def f(m):
    # flash sort

    # count elements in each bucket
    m[J] = 0
    while m[J] < SCORE:
        m[BUFFER + m[m[J]] // DIST] = m[BUFFER + m[m[J]] // DIST] + 1
        m[J] = m[J] + 1

    # compute prefix sum
    m[J] = 1
    while m[J] < BUCKETS:
        m[BUFFER + m[J]] = m[BUFFER + m[J]] + m[BUFFER + m[J] - 1]
        m[J] = m[J] + 1

    # permutation
    m[I] = 0
    m[J] = 0
    m[K] = BUCKETS - 1
    while m[I] < SCORE - 1:
        while m[J] > m[BUFFER + m[K]] - 1:
            m[J] = m[J] + 1
            m[K] = m[m[J]] // DIST

        # found a non-full bucket
        m[TEMP] = m[m[J]]
        while m[J] != m[BUFFER + m[K]]:
            m[K] = m[TEMP] // DIST
            m[BUFFER + m[K]] = m[BUFFER + m[K]] - 1

            m[TEMP2] = m[m[BUFFER + m[K]]]
            m[m[BUFFER + m[K]]] = m[TEMP]
            m[TEMP] = m[TEMP2]
            m[I] = m[I] + 1

    # insertion sort step
    m[J] = 1
    while m[J] < SCORE:
        m[TEMP] = m[m[J]]
        m[I] = m[J] - 1
        while (m[I] >= 0) & (m[m[I]] > m[TEMP]):
            m[m[I] + 1] = m[m[I]]
            m[I] = m[I] - 1
        m[m[I] + 1] = m[TEMP]
        m[J] = m[J] + 1

Version 8, score 13534, with minor opts:

def f(m):
    # flash sort

    # count elements in each bucket
    m[J] = 0
    while m[J] < SCORE:
        m[TEMP] = BUFFER + m[m[J]] // DIST
        m[m[TEMP]] = m[m[TEMP]] + 1
        m[J] = m[J] + 1

    # compute prefix sum
    m[J] = 1
    while m[J] < BUCKETS:
        m[TEMP] = BUFFER + m[J]
        m[m[TEMP]] = m[m[TEMP]] + m[m[TEMP] - 1]
        m[J] = m[J] + 1

    # permutation
    m[I] = 0
    m[J] = 0
    m[K] = BUFFER + BUCKETS - 1
    while m[I] < SCORE - 1:
        while m[J] > m[m[K]] - 1:
            m[J] = m[J] + 1
            m[K] = BUFFER + m[m[J]] // DIST

        # found a non-full bucket
        m[TEMP] = m[m[J]]
        while 1:
            m[K] = BUFFER + m[TEMP] // DIST
            m[TEMP3] = m[m[K]] - 1
            m[m[K]] = m[TEMP3]

            m[TEMP2] = m[m[TEMP3]]
            m[m[TEMP3]] = m[TEMP]
            m[TEMP] = m[TEMP2]
            m[I] = m[I] + 1
            if m[J] == m[TEMP3]:
                break

    # insertion sort step
    m[J] = 1
    while m[J] < SCORE:
        if m[m[J] - 1] <= m[m[J]]:
            # fast path
            m[J] = m[J] + 1
            continue

        m[TEMP] = m[m[J]]
        m[I] = m[J] - 1
        while m[m[I]] > m[TEMP]:
            m[m[I] + 1] = m[m[I]]
            m[I] = m[I] - 1
        m[m[I] + 1] = m[TEMP]
        m[J] = m[J] + 1

Test harness:

from pwn import *
import json

if args.PROG:
    file = args.PROG
else:
    file = "sort.py"

left = 100
right = 15000
buckets = 4096
while left + 1 < right:
    mid = (left + right) // 2

    score = mid
    code = (
        open(file)
        .read()
        .replace("SCORE", str(score))
        .replace("BUCKETS", str(buckets))
        .replace("DIST", str((2**64) // buckets))
        .replace("SIZE", "999995")
        .replace("MID", "999993")
        .replace("RIGHT", "999992")
        .replace("END", "999990")
        .replace("I", "999999")
        .replace("J", "999998")
        .replace("TEMP3", "999988")
        .replace("TEMP2", "999989")
        .replace("TEMP", "999997")
        .replace("K", "999996")
        .replace("LEFT", "999994")
        .replace("GAP", "999991")
        .replace("BUFFER", "900000")
    )
    print(code)

    req = json.dumps(
        {
            "code": code,
            "score": score,
        }
    )

    p = process(["python3", "main.py"])
    p.sendline(req.encode())
    res = json.loads(p.recvall().decode())
    print(res, mid, left, right)
    if res["ok"]:
        left = mid
    else:
        right = mid

Idea from pitust on Discord after competition: since the keys are uniformly distributed 64-bit random numbers, use each value's own top bits as its address in a large sparse table, and on collision do hash-table linear probing with a sorted swap (the smaller value swaps forward), so every probe run stays sorted and the whole table ends up sorted, then a final scan copies out the non-zero cells. This pushes the score up to ~26k.


AI-generated and AI-solved continuation (documented after the competition)

The following sections document the AI-generated and AI-solved attempts that pushed the score from ~13k to 65,536.

Unlike the human attempts above — which ran under the original platform's default jail (JAIL_CPU=100, i.e. only 10% of one core, on a shared CPU) — for the AI part the jail's CPU limit was deliberately changed from the default (JAIL_CPU=100, 10%) to JAIL_CPU=1000 (100% of one core) for local testing on a much stronger CPU. The 20 s wall-time limit still applied, but because the local CPU was far faster and used at full 100% rather than 10%, the wall time per cycle shrank to the point where the binding constraint became the algorithmic one — the 1e6-cycle budget — rather than CPU throttling; this was done deliberately to see the algorithmic limit. The submissions still had to satisfy the platform's constraints (1e6-cycle budget, 4,096-byte code limit, and the 20 s wall-time timeout on the default jail), but during development the strong local CPU meant wall time was never the bottleneck, so the AI could iterate purely on the algorithmic (cycle) limit.

In the competition itself, the human competitor reached ~27k, while the AI reached ~42k; the push to 65,536 documented below was completed afterwards with local testing. The final score of 65,536 is the maximum the platform allows with a 1e6-cycle budget, because the bitonic-SIMD sort requires the array size to be a power of two, and the next power of two (131,072) needs more than 1e6 cycles for packing/unpacking alone.

Phase 1: bucket/counting sorts (13,000 → 20,448)

Starting from the flash sort (version 8, 13,534), the biggest improvement was replacing the in-place flash permutation with a linked-list bucket sort. Each value's top bits select a bucket; every value is stored into a per-bucket singly linked list threaded through the input area. The pointer trick: a bucket's head pointer lives at BUFFER + bucket, and the "next" pointer of element i lives at N + i. This is a pure counting sort with no comparisons at all, so the build is ~23 cycles per element. A final pass walks every bucket in order, copies the values out in list order, and runs an insertion sort on any bucket whose run exceeds one element (buckets are small, so this is nearly free).

# score 13,000 (bucket counting sort with linked lists)
# memory layout:
#   m[0 .. 12999]        input values
#   m[13000 + i]         "next" pointer of element i (index + 1, 0 = end of list)
#   m[26000 .. 38999]    output buffer (sorted values are written here)
#   m[39000 + b]         head of bucket b (index + 1, 0 = empty)
# bucket = top 13 bits of the value (2^13 = 8192 buckets, one per 2^51 range)
def f(m):
    # build: push every value onto the head of its bucket's linked list
    m[-1] = 0
    while m[-1] < 13000:                      # for i in range(N):
        m[-5] = m[m[-1]] // 2251799813685248  #   b = value >> 51
        m[13000 + m[-1]] = m[39000 + m[-5]]   #   next[i] = head[b]
        m[39000 + m[-5]] = m[-1] + 1          #   head[b] = i + 1
        m[-1] = m[-1] + 1

    # walk the buckets in order, copying each run out, and insertion-sort
    # the runs of length >= 2 (buckets hold ~N/8192 = 1.6 elements on average,
    # so the insertion sorts are nearly free)
    m[-9] = 26000                             # w = output write pointer
    m[-1] = 0
    while m[-1] < 8192:                       # for b in range(2**13):
        m[-10] = m[39000 + m[-1]]             #   node = head[b]
        m[-6] = m[-9]                         #   s = w  (start of the run)
        while m[-10] > 0:                     #   while node: walk the list
            m[m[-9]] = m[m[-10] - 1]          #     out[w] = a[node - 1]
            m[-10] = m[13000 + m[-10] - 1]    #     node = next[node - 1]
            m[-9] = m[-9] + 1                 #     w = w + 1
        m[-7] = m[-9]                         #   e = w  (end of the run)
        m[-8] = m[-7] - m[-6]                 #   length = e - s
        if m[-8] >= 2:                        #   if length >= 2: insertion sort
            m[-2] = m[-6] + 1
            while m[-2] < m[-7]:
                m[-4] = m[m[-2]]              #     v = out[j]
                m[-3] = m[-2] - 1
                while m[-3] >= m[-6]:         #     shift larger values right
                    if m[m[-3]] <= m[-4]:
                        break
                    m[m[-3] + 1] = m[m[-3]]
                    m[-3] = m[-3] - 1
                m[m[-3] + 1] = m[-4]          #     place v
                m[-2] = m[-2] + 1
        m[-1] = m[-1] + 1

    # copy the output buffer back over the input area
    m[-1] = 0
    while m[-1] < 13000:
        m[m[-1]] = m[26000 + m[-1]]
        m[-1] = m[-1] + 1

The submissions between 13,000 and 20,448 are refinements of this scheme (the one outlier, 16,384, already belongs to the bigint phase below):

  • 8-cell unrolled builds — process 8 input values per loop iteration (the loop overhead amortizes 8×), which alone took the build cost from ~23 to ~9.5 cycles/element.
  • k = 14 buckets instead of 13 — a finer split shortens the insertion-sort runs inside each bucket.
  • Walk shortcuts — for bucket runs of length 1/2/3/4, emit specialized copy-out code that avoids the generic list-walk overhead.
  • Write-in-sorted-order for runs of length 2 and 3 — a length-2 bucket can be emitted directly in sorted order with one comparison instead of walking a list and re-sorting.

Best of this family: score 20,448 at 999,973 cycles. The ceiling here is set by the per-element cost of the unrolled build plus the bucket walk and the copy-out; with the 1e6 budget this approach tops out around 20k–21k.

Phase 2: sorted-swap linear probing (23,000 → 26,670)

After the competition, pitust on Discord pointed out that because the keys are uniformly distributed 64-bit random numbers, we can use each value's own top bits as its address in a large sparse table, and on collision do hash-table linear probing with a sorted swap: when a new value lands on an occupied slot, if the new value is smaller, swap it with the occupant and continue probing with the (larger) displaced value. Every probe run therefore stays sorted, and when all values are inserted, scanning the table left to right yields the sorted order.

My implementation uses table[i] = value // RADIX + N with the table size and RADIX tuned so the probe index never overflows the table. The submitted source unrolls the insert block 8× per loop iteration (batch indices m[-2]..m[-9]) and the copy-out scan 12×; the single-value form below is the equivalent, readable version:

# score 23,000 (sorted-swap probing)
# sparse table of 65,538 cells at m[23000 .. 88537]:
#   slot = 23000 + value // 2^48   (the top 16 bits of the value pick the slot;
#                                   the table has 2^16 slots plus two slack
#                                   cells, so linear probing never runs off
#                                   the end)
# registers: m[-12] = value being inserted, m[-10] = slot address,
#            m[-11] = current occupant,  m[-14] = copy-out write pointer
def f(m):
    m[-2] = 0
    while m[-2] < 23000:
        # insert a[m[-2]] with sorted-swap linear probing:
        m[-12] = m[m[-2]]                            # v = a[i]
        m[-10] = 23000 + m[-12] // 281474976710656   # slot = N + top16(v)
        m[-11] = m[m[-10]]                           # occ = table[slot]
        while m[-11] != 0:                           # probe while the slot is full
            if m[-12] < m[-11]:                      #   if v < occ:
                m[m[-10]] = m[-12]                   #     table[slot] = v
                m[-12] = m[-11]                      #     v = occ
            m[-10] = m[-10] + 1                      #   slot = slot + 1
            m[-11] = m[m[-10]]                       #   occ = table[slot]
        m[m[-10]] = m[-12]                           # empty slot found: table[slot] = v
        m[-2] = m[-2] + 1

    # every probe run stayed sorted, so a left-to-right scan of the table
    # yields the sorted array: copy the non-zero cells back to the front
    m[-14] = 0
    m[-2] = 23000
    while m[-2] < 88538:
        if m[m[-2]] != 0:
            m[m[-14]] = m[m[-2]]
            m[-14] = m[-14] + 1
        m[-2] = m[-2] + 1

(Reformatted for readability; the 8× unrolling of the insert block and the 12× unrolling of the scan are the only differences from the submitted code.)

The score growth 23,000 → 23,570 → 23,760 → 23,850 → 23,870 → 24,040 → 24,048 → 25,800 → 26,200 → 26,278 → 26,358 → 26,412 → 26,416 → 26,440 → 26,520 → 26,560 → 26,640 → 26,670 came from:

  • shrinking the probe distance by choosing the table size / RADIX split so that the expected probe length is minimal,
  • unrolling the 8-value insertion block,
  • unrolling the final copy-out scan (checking several table cells per iteration),
  • micro-optimizing the swap (only read m[-15] when actually swapping).

The hard limit of this family is the per-value cost — roughly 13 cycles for the address computation and probe loop, plus the copy-out scan and loop overhead — so ~26.7k is the ceiling.

Phase 3: SIMD bigint bitonic sort (16,384 → 65,536)

The decisive step was to stop doing per-element work at all — an idea also suggested on Discord, just like the sorted-swap probing in Phase 2. In the restricted subset we cannot use shifts (<</>>), but * and // by powers of two are allowed and work on arbitrarily large integers. So we can pack the whole array into one giant Python bigint, with each value in its own 65-bit lane (values < 2^64, plus one spare high bit per lane), and run a bitonic sorting network using only bigint-wide SWAR-style operations. Each compare-swap stage is O(1) bigint ops, so the whole sort is O(log² N) operations instead of O(N log N) per-element work — at the cost of doing arithmetic on multi-megabit integers.

The core compare-swap for a stage at distance d (with lane width 65, C = 2^65, and a high-bit mask H with bit 64 set in every lane):

# lane-wise compare-swap of two halves of the packed bigint x, distance d
hi   = (x & NOTMASK) // C ** d            # upper-lane values, shifted down
lo   = x & MASK                           # lower-lane values
ge   = ((lo + H - hi) & H) - ((lo + H - hi) & H) // 2**64   # full 64-bit lane mask where lo >= hi
mn   = (lo & ge) | (hi & NOTGE)           # lane-wise minima
mx   = (lo & NOTGE) | (hi & GE)           # lane-wise maxima
# the stage's direction then picks which blend occupies the lower lanes;
# the other blend is shifted back up by d lanes and OR'd in.

The first working version sorted 16,384 values (16 bigints of 1,024 lanes) in 727,755 cycles — comparable per value to the probing approach, but with much better scaling, since the sort adds only O(log² N) stages: by 32,768 the all-in cost was down to ~23 cycles per value and it kept dropping.

The wall-clock problem. The first attempt to scale up — a single 32,768-lane bigint (2.1M-bit integers) — ran the whole sort in ~745k cycles, but every operation on a 2.1M-bit integer is dozens of times slower in wall time even though it only counts as one cycle, so the run exceeded the platform's 20 s wall-time limit (default jail, JAIL_CPU=100, i.e. 10% of a core) and timed out. The architectural fix was to split the array into 16 smaller bigints of 2,048 lanes each (133k-bit integers), so every bigint op is ~16× cheaper in wall time. That version (score 32,768) fit the platform at 743,095 cycles — with its unpack already being the dominant cost, since it extracted each 64-bit lane with one bigint division per lane (x // 2**65), i.e. 32,768 slow divisions per run.

The second fix was to batch the unpack: mask off the low 16 lanes with one &, shift the bigint down by 16 lanes with one //, then split the small (≤ 16-lane) chunk into individual lanes with cheap small-integer divisions. That single change cut the total wall time from ~11.6 s to ~1.4 s and gave plenty of headroom for the push to 65,536.

Pushing past 32,768. The next power of two is 65,536, and there the cycle budget (1e6) becomes the binding constraint, not wall time. The per-value cost of pack + unpack had to fit in ~14.4 cycles. That forced five more ideas:

  1. Fewer, larger bigints. The sorting network's cycle cost is proportional to (number of bigints) × (number of stages), so going from 16×2,048 to 4×16,384 lanes roughly halves the sort's cycle count (the per-op cost is counted once regardless of operand size).

  2. Expression-chain packing. Pack several values per statement by summing v0 + v1*C + v2*C**2 + ... inside one assignment (only one accumulator read per batch). With 16 values per statement and the powers C**1..C**16 kept in memory cells (built once), packing costs ~6 cycles per value.

  3. Zero-padding instead of full padding. For N < 65,536 the leftover lanes are left as zero, and zeros sort to the front automatically. The pack then only packs the N real values, and the unpack only extracts the real lanes (skipping the leading zero lanes with one bigint division per bigint). This makes the cost proportional to N, not to the padded power of two.

  4. Pattern-doubling mask construction. Building the SWAR masks (the per-lane direction patterns) with the straightforward loop costs O(N) cycles (~230k at 16,384-lane bigints). Instead, each mask is derived from the previous one by mask[k] = mask[k-1] | mask[k-1] * C**(2**(k-1)) — doubling the pattern — so all masks are built in ~150 cycles total.

  5. Batch alignment. The pack/unpack process a fixed number of lanes per iteration, so N must be aligned to the batch size to keep the per-bigint real counts exact (8 lanes in the earlier versions, 16 in the later ones; this is why the final N is 65,536 and not, say, 65,532).

The intermediate milestones in this push were 57,600 (999,343), 57,632 (999,879), 59,000 (999,433) and 61,440 (984,626), each a different N/P split that had to obey the batch alignment (8 or 16 lanes, depending on the version) and stay under 1e6 cycles.

Final submission: score 65,536 at 997,852 cycles, 4,081 bytes. Four bigints of 16,384 lanes each; the constants (masks, powers) live in cells just above the input, and the multiplier/divisor cells use negative indices to keep the code under the 4,096-byte cap. It runs in ~5.7 s locally.

For reference, here is the final code, reformatted from the 4,081-byte submission; comments are restored and the pack statement chain is spread over its 16 lanes:

# score 65,536 (final submission), reformatted for readability:
# four bigints of 16,384 lanes live in m[-30 .. -27]; every lane holds one
# 65-bit value (64 value bits + 1 spare high bit).  C = 2^65; the masks and
# lane powers are built once into m[80002 .. 80031] and the negative-index
# register cells.
def f(m):
    # ---- constants ----
    # M = C^16384 - 1: a mask with all 65 bits set in all 16,384 lanes
    m[-39] = 0x20000000000000000               # C
    m[-11] = 0
    while m[-11] < 14:                         # square 14x: C, C^2, C^4, ..., C^16384
        m[-39] = m[-39] * m[-39]
        m[-11] = m[-11] + 1
    m[-39] = m[-39] - 1                        # M (all-lanes mask)
    m[-40] = (m[-39] // 0x1ffffffffffffffff) * 0x10000000000000000
                                               # H: bit 64 set in every lane (sign bit)

    m[80002] = 0x20000000000000000             # m[80002 + k] = C^(2^k), k = 0..14
    m[-11] = 1
    while m[-11] < 15:
        m[80002 + m[-11]] = m[80002 + m[-11] - 1] * m[80002 + m[-11] - 1]
        m[-11] = m[-11] + 1

    m[-12] = 0x1ffffffffffffffff               # C - 1 (the 65 low bits of lane 0)
    m[80017] = (m[-39] // (m[80002 + 1] - 1)) * 0x1ffffffffffffffff
    m[-11] = 1
    while m[-11] < 14:                         # m[80017 + k] = half-block mask k:
        m[-12] = m[-12] | m[-12] * m[80002 + m[-11] - 1]     # pattern doubling
        m[80017 + m[-11]] = (m[-39] // (m[80002 + m[-11] + 1] - 1)) * m[-12]
        m[-11] = m[-11] + 1
    m[80017 + 14] = m[-39]                     # m[80017 + 14] = M

    m[-76] = 0x20000000000000000               # m[-76 + k] = C^(k+1), k = 0..15
    m[-11] = 1
    while m[-11] < 16:
        m[-76 + m[-11]] = m[-76 + m[-11] - 1] * 0x20000000000000000
        m[-11] = m[-11] + 1
    m[-60] = m[-76 + 15] - 1                   # C^16 - 1: mask of the low 16 lanes
    m[-59] = 0x10000000000000000               # 2^64: normalizes the sign borrow

    # ---- pack: 65,536 values -> 4 bigints x 16,384 lanes, 16 lanes/statement ----
    m[-3] = 0                                  # bigint index
    m[-11] = 0                                 # input index
    while m[-3] < 4:
        m[-6] = 0                              # accumulator (the bigint)
        m[-5] = 0
        while m[-5] < 16384:
            # one statement packs 16 consecutive values into the low lanes;
            # the pre-computed lane powers m[-76 + k] = C^(k+1) do the shifting
            m[-6] = (
                m[-6] * m[-76 + 15]
                + m[m[-11]]
                + m[m[-11] + 1] * m[-76]
                + m[m[-11] + 2] * m[-75]
                + m[m[-11] + 3] * m[-74]
                + m[m[-11] + 4] * m[-73]
                + m[m[-11] + 5] * m[-72]
                + m[m[-11] + 6] * m[-71]
                + m[m[-11] + 7] * m[-70]
                + m[m[-11] + 8] * m[-69]
                + m[m[-11] + 9] * m[-68]
                + m[m[-11] + 10] * m[-67]
                + m[m[-11] + 11] * m[-66]
                + m[m[-11] + 12] * m[-65]
                + m[m[-11] + 13] * m[-64]
                + m[m[-11] + 14] * m[-63]
                + m[m[-11] + 15] * m[-62]
            )
            m[-5] = m[-5] + 16
            m[-11] = m[-11] + 16
        m[-30 + m[-3]] = m[-6]                 # bigint b stored at m[-30 + b]
        m[-3] = m[-3] + 1

    # ---- bitonic sorting network ----
    # m[-1] = subarray size (2, 4, ..., 65536), m[-2] = compare-swap distance,
    # m[-16] = log2(size), m[-14] = index of the mask for this distance
    m[-1] = 2; m[-16] = 1
    while m[-1] <= 65536:
        m[-2] = m[-1] // 2
        m[-14] = m[-16] - 1
        while m[-2] > 0:
            if m[-2] >= 16384:
                # distance spans whole bigints: compare-swap pairs of bigints.
                # m[-5] = partner bigint; direction comes from (m[-4] & m[-1]).
                m[-3] = 0; m[-4] = 0
                while m[-3] < 4:
                    m[-5] = m[-3] ^ (m[-2] // 16384)
                    if m[-5] > m[-3]:
                        m[-6] = m[-30 + m[-3]]
                        m[-7] = m[-30 + m[-5]]
                        # lane-wise sign of (lo - hi), one bit per lane:
                        m[-8] = m[-6] - m[-7]      # borrows from the spare bit
                        m[-8] = m[-8] & m[-40]     # keep only the sign bits
                        m[-8] = m[-8] - (m[-8] // m[-59])  # 0x1..1 where lo < hi
                        m[-13] = m[-39] ^ m[-8]    # NOT of the sign mask
                        m[-9] = (m[-6] & m[-8]) | (m[-7] & m[-13])   # lane-wise min
                        m[-10] = (m[-6] & m[-13]) | (m[-7] & m[-8])  # lane-wise max
                        if (m[-4] & m[-1]) == 0:
                            m[-30 + m[-3]] = m[-9]
                            m[-30 + m[-5]] = m[-10]
                        else:
                            m[-30 + m[-3]] = m[-10]
                            m[-30 + m[-5]] = m[-9]
                    m[-4] = m[-4] + 16384
                    m[-3] = m[-3] + 1
            else:
                # intra-bigint compare-swap at distance m[-2] (< 16384).
                # m[-13] = half-block mask, m[-15] = C^d, m[-17] = direction
                m[-13] = m[80017 + m[-14]]
                m[-15] = m[80002 + m[-14]]
                m[-17] = m[80017 + m[-14]] & m[80017 + m[-16]]
                m[-18] = m[-39] ^ m[-13]
                m[-3] = 0; m[-4] = 0
                while m[-3] < 4:
                    if m[-1] >= 16384:
                        # subarray spans bigints: direction alternates per block
                        m[-17] = m[-13] * (1 - (m[-4] & m[-1]) // m[-1])
                    m[-6] = m[-30 + m[-3]]
                    m[-7] = m[-6] & m[-18]       # hi half (upper d lanes of each block)
                    m[-7] = m[-7] // m[-15]      # shift it down by d lanes
                    m[-6] = m[-6] & m[-13]       # lo half
                    m[-8] = m[-6] - m[-7]        # lane-wise sign of lo - hi
                    m[-8] = m[-8] & m[-40]
                    m[-8] = m[-8] - (m[-8] // m[-59])
                    m[-9] = m[-39] ^ m[-8]
                    m[-10] = (m[-6] & m[-8]) | (m[-7] & m[-9])      # min
                    m[-9] = (m[-6] & m[-9]) | (m[-7] & m[-8])       # max
                    m[-12] = m[-13] ^ m[-17]     # blend mask for this direction
                    m[-8] = (m[-10] & m[-17]) | (m[-9] & m[-12])
                    m[-9] = (m[-9] & m[-17]) | (m[-10] & m[-12])
                    m[-10] = m[-9] * m[-15]      # shift the upper half back up
                    m[-30 + m[-3]] = m[-8] | m[-10]
                    m[-4] = m[-4] + 16384
                    m[-3] = m[-3] + 1
            m[-2] = m[-2] // 2
            m[-14] = m[-14] - 1
        m[-1] = m[-1] * 2
        m[-16] = m[-16] + 1

    # ---- unpack: 16 lanes per iteration ----
    m[-3] = 0; m[-11] = 0
    while m[-3] < 4:
        m[-6] = m[-30 + m[-3]]
        m[-5] = 0
        while m[-5] < 16384:
            m[-7] = m[-6] & m[-60]               # peel off the low 16 lanes
            m[-6] = m[-6] // m[-76 + 15]         # shift the bigint down 16 lanes
            m[m[-11]] = m[-7] & 0xffffffffffffffff
            m[m[-11] + 1] = (m[-7] // m[-76]) & 0xffffffffffffffff
            m[m[-11] + 2] = (m[-7] // m[-75]) & 0xffffffffffffffff
            m[m[-11] + 3] = (m[-7] // m[-74]) & 0xffffffffffffffff
            m[m[-11] + 4] = (m[-7] // m[-73]) & 0xffffffffffffffff
            m[m[-11] + 5] = (m[-7] // m[-72]) & 0xffffffffffffffff
            m[m[-11] + 6] = (m[-7] // m[-71]) & 0xffffffffffffffff
            m[m[-11] + 7] = (m[-7] // m[-70]) & 0xffffffffffffffff
            m[m[-11] + 8] = (m[-7] // m[-69]) & 0xffffffffffffffff
            m[m[-11] + 9] = (m[-7] // m[-68]) & 0xffffffffffffffff
            m[m[-11] + 10] = (m[-7] // m[-67]) & 0xffffffffffffffff
            m[m[-11] + 11] = (m[-7] // m[-66]) & 0xffffffffffffffff
            m[m[-11] + 12] = (m[-7] // m[-65]) & 0xffffffffffffffff
            m[m[-11] + 13] = (m[-7] // m[-64]) & 0xffffffffffffffff
            m[m[-11] + 14] = (m[-7] // m[-63]) & 0xffffffffffffffff
            m[m[-11] + 15] = (m[-7] // m[-62]) & 0xffffffffffffffff
            m[-11] = m[-11] + 16
            m[-5] = m[-5] + 16
        m[-3] = m[-3] + 1

This is the end of the road for this particular platform: the score is capped at 65,536 because the sort needs a power-of-two size and the next power of two would exceed the 1e6-cycle budget. Reaching it required combining the two Discord-suggested ideas (sorted-swap probing, and the SIMD bigint bitonic sort) with batching to fix wall-clock time and the cycle-level micro-optimizations above; all of the implementation and verification was done by the AI agent.

Comments