Minimal cost pairing of 1..p-1 mod p = 2000000011: pair costs are (a*b) mod p; among all pairings of minimal total cost, print the product of the pair costs (exact integer). Theory. By Wilson, the product of the pair costs of any pairing is congruent to -1 (mod p). Costs are at least 1, and replacing a composite cost c = u*v by pairs of costs u and v never raises the total since (u-1)+(v-1) <= uv-1. So the total is at least (p-1)/2 + W, where W is the minimal weight sum e_q*(q-1) over prime multisets q^e_q whose product is -1 mod p. Conversely any such prime multiset q_1..q_r is realised by the chain pairing (1, Q_1), (Q_1^-1, Q_2), ..., (Q_{r-1}^-1, Q_r = p-1) with Q_i = q_1*...*q_i mod p, all other elements paired with their inverses at cost 1. Hence optimal cost multisets are exactly the minimal-weight prime factorisations of integers N == -1 (mod p); for p = 2000000011 the minimiser is unique (verified by the exhaustive search below), so the cost product is that N. Search. Every multiset splits as X * 2^a * 3^b with X over primes >= 5. Hash all residues 2^a 3^b (a + 2b <= 239) to their minimal weight, then DFS over X keeping the running inverse residue and weight; each node looks up -X^-1 in the table. Weight cap 239 is an upper bound reached by the search itself, so the minimum is exact. The unique minimiser N fits in i128 and is printed by repeated division by 10.
# Project Euler 789
# Minimal cost pairing of 1..p-1 mod p = 2000000011: pair costs are
# (a*b) mod p; among all pairings of minimal total cost, print the
# product of the pair costs (exact integer).
#
# Theory. By Wilson, the product of the pair costs of any pairing is
# congruent to -1 (mod p). Costs are at least 1, and replacing a
# composite cost c = u*v by pairs of costs u and v never raises the
# total since (u-1)+(v-1) <= uv-1. So the total is at least
# (p-1)/2 + W, where W is the minimal weight sum e_q*(q-1) over prime
# multisets q^e_q whose product is -1 mod p. Conversely any such prime
# multiset q_1..q_r is realised by the chain pairing (1, Q_1),
# (Q_1^-1, Q_2), ..., (Q_{r-1}^-1, Q_r = p-1) with Q_i = q_1*...*q_i
# mod p, all other elements paired with their inverses at cost 1.
# Hence optimal cost multisets are exactly the minimal-weight prime
# factorisations of integers N == -1 (mod p); for p = 2000000011 the
# minimiser is unique (verified by the exhaustive search below), so
# the cost product is that N.
#
# Search. Every multiset splits as X * 2^a * 3^b with X over primes
# >= 5. Hash all residues 2^a 3^b (a + 2b <= 239) to their minimal
# weight, then DFS over X keeping the running inverse residue and
# weight; each node looks up -X^-1 in the table. Weight cap 239 is an
# upper bound reached by the search itself, so the minimum is exact.
# The unique minimiser N fits in i128 and is printed by repeated
# division by 10.
extern {
function calloc(n: i64, size: i64) -> ptr<void>
function free(p: ptr<void>) -> void
}
const P: i64 = 2000000011
const CAP: i64 = 239
const HB: i64 = 65536
const HMASK: i64 = 65535
const QMAX: i64 = 241
function mulmod(a: i64, b: i64) -> i64 {
# a, b < P < 2^31, so a*b < 2^62 fits i64
return (a * b) % P
}
function modpow(b: i64, e: i64) -> i64 {
let mut r: i64 = 1
let mut bb: i64 = b % P
let mut ee: i64 = e
while ee > 0 {
if ee % 2 == 1 {
r = mulmod(r, bb)
}
bb = mulmod(bb, bb)
ee = ee / 2
}
return r
}
function hslot(tk: ptr<i64>, t: i64) -> i64 {
let mut h: i64 = ((t >> 15) ^ t) & HMASK
while tk[h] != 0 && tk[h] != t {
h = (h + 1) & HMASK
}
return h
}
# DFS over multisets of primes >= 5 (Q ascending, weights Q[j]-1).
# ri = inverse of the running product, w = weight so far, d = depth.
# st: 0 best, 1 ties, 2 saved w, 3 saved target residue, 4 saved depth.
function dfs(i: i64, ri: i64, w: i64, d: i64,
Q: ptr<i64>, QI: ptr<i64>, nq: i64,
tk: ptr<i64>, tw: ptr<i64>,
st: ptr<i64>, stk: ptr<i64>, sol: ptr<i64>) -> i64 {
# target: 2^a 3^b == -1 / X (mod P)
let t: i64 = ((P - 1) * ri) % P
let h: i64 = hslot(tk, t)
if tk[h] == t {
let cand: i64 = w + tw[h]
if cand < st[0] {
st[0] = cand
st[1] = 1
st[2] = w
st[3] = t
st[4] = d
let mut k: i64 = 0
while k < d {
sol[k] = stk[k]
k = k + 1
}
} else {
if cand == st[0] {
st[1] = st[1] + 1
}
}
}
let mut j: i64 = i
while j < nq {
let wj: i64 = Q[j] - 1
if w + wj > st[0] {
j = nq
} else {
stk[d] = Q[j]
dfs(j, mulmod(ri, QI[j]), w + wj, d + 1,
Q, QI, nq, tk, tw, st, stk, sol)
j = j + 1
}
}
return 0
}
function main() -> i32 {
# sieve primes up to QMAX
let comp: ptr<i64> = calloc(QMAX + 1, 8)
let mut i: i64 = 2
while i * i <= QMAX {
if comp[i] == 0 {
let mut j: i64 = i * i
while j <= QMAX {
comp[j] = 1
j = j + i
}
}
i = i + 1
}
let Q: ptr<i64> = calloc(64, 8)
let QI: ptr<i64> = calloc(64, 8)
let mut nq: i64 = 0
i = 5
while i <= QMAX {
if comp[i] == 0 && i - 1 <= CAP && i < P {
Q[nq] = i
QI[nq] = modpow(i, P - 2)
nq = nq + 1
}
i = i + 1
}
# hash residues of 2^a 3^b (a + 2b <= CAP) to minimal weight a + 2b
let tk: ptr<i64> = calloc(HB, 8)
let tw: ptr<i64> = calloc(HB, 8)
let mut b: i64 = 0
let mut r3: i64 = 1
while 2 * b <= CAP {
let mut cur: i64 = r3
let mut a: i64 = 0
while a + 2 * b <= CAP {
let h: i64 = hslot(tk, cur)
if tk[h] == 0 {
tk[h] = cur
tw[h] = a + 2 * b
} else {
if a + 2 * b < tw[h] {
tw[h] = a + 2 * b
}
}
cur = mulmod(cur, 2)
a = a + 1
}
r3 = mulmod(r3, 3)
b = b + 1
}
# exhaustive search up to weight CAP
let st: ptr<i64> = calloc(8, 8)
let stk: ptr<i64> = calloc(64, 8)
let sol: ptr<i64> = calloc(64, 8)
st[0] = CAP + 1
dfs(0, 1, 0, 0, Q, QI, nq, tk, tw, st, stk, sol)
# reconstruct N = (solution primes >= 5) * 2^a 3^b
let mut N: i128 = 1 as i128
let mut k: i64 = 0
while k < st[4] {
N = N * (sol[k] as i128)
k = k + 1
}
let rem: i64 = st[0] - st[2]
let tt: i64 = st[3]
let mut ba: i64 = -1
let mut bb2: i64 = -1
b = 0
r3 = 1
while 2 * b <= rem {
let mut cur: i64 = r3
let mut a: i64 = 0
while a + 2 * b <= rem {
if cur == tt && a + 2 * b == rem {
ba = a
bb2 = b
}
cur = mulmod(cur, 2)
a = a + 1
}
r3 = mulmod(r3, 3)
b = b + 1
}
k = 0
while k < ba {
N = N * (2 as i128)
k = k + 1
}
k = 0
while k < bb2 {
N = N * (3 as i128)
k = k + 1
}
# print N in decimal (exceeds signed 64-bit)
let digs: ptr<i64> = calloc(64, 8)
let mut nd: i64 = 0
let mut x: i128 = N
while x > (0 as i128) {
digs[nd] = (x % (10 as i128)) as i64
x = x / (10 as i128)
nd = nd + 1
}
let mut j2: i64 = nd - 1
while j2 >= 0 {
printf("%c", (48 + digs[j2]) as i32)
j2 = j2 - 1
}
printf("\n")
free(comp)
free(Q)
free(QI)
free(tk)
free(tw)
free(st)
free(stk)
free(sol)
free(digs)
return 0
}
Generated C
#include <stdint.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
/* Flow runtime helpers */
typedef struct flow_temp_node { struct flow_temp_node* next; } flow_temp_node;
static flow_temp_node* flow_temp_head = NULL;
static int flow_temp_atexit_set = 0;
__attribute__((unused)) static void flow_temp_free_all(void) {
while (flow_temp_head) {
flow_temp_node* n = flow_temp_head;
flow_temp_head = n->next;
free(n);
}
}
__attribute__((unused)) static void* flow_temp_alloc(size_t nbytes) {
flow_temp_node* node = (flow_temp_node*)malloc(sizeof(flow_temp_node) + nbytes);
if (!node) return NULL;
node->next = flow_temp_head;
flow_temp_head = node;
if (!flow_temp_atexit_set) {
flow_temp_atexit_set = 1;
atexit(flow_temp_free_all);
}
return (void*)(node + 1);
}
#ifndef FLOW_DIAG
#define FLOW_DIAG(msg) fprintf(stderr, "%s", (msg))
#endif
#ifndef FLOW_LOG
#define FLOW_LOG(fmt, ...) printf(fmt, __VA_ARGS__)
#endif
#ifndef FLOW_LOG_EMPTY
#define FLOW_LOG_EMPTY(fmt) printf(fmt)
#endif
static char* flow_strcat(const char* a, const char* b) {
size_t la = strlen(a ? a : ""), lb = strlen(b ? b : "");
char* r = (char*)flow_temp_alloc(la + lb + 1);
if (!r) return NULL;
if (la) memcpy(r, a, la);
if (lb) memcpy(r + la, b, lb);
r[la + lb] = '\0';
return r;
}
#define __flow_in_arr(arr, val) __extension__ ({ \
int _found = 0; \
size_t _n = sizeof(arr)/sizeof((arr)[0]); \
for (size_t _i = 0; _i < _n; _i++) { \
if ((arr)[_i] == (val)) { _found = 1; break; } \
} _found; })
/* Unified fault handler (MISRA #279) — override with -DFLOW_FAULT_HANDLER=fn */
#ifndef FLOW_FAULT_HANDLER
__attribute__((unused)) static inline void flow_fault_handler(const char* msg) {
fprintf(stderr, "flow: %s\n", msg ? msg : "fault");
abort();
#if defined(__GNUC__) || defined(__clang__)
__builtin_unreachable();
#endif
}
#else
#define flow_fault_handler FLOW_FAULT_HANDLER
#endif
#define flow_div_by_zero_handler() flow_fault_handler("division by zero")
#define flow_shift_ub_handler() flow_fault_handler("invalid shift (amount out of range or left-shift of negative)")
#ifndef FLOW_CHECKED_DIV
#define FLOW_CHECKED_DIV(L, R) (((R) != 0) ? ((L) / (R)) : (flow_div_by_zero_handler(), (L) * 0))
#endif
#ifndef FLOW_CHECKED_MOD
#define FLOW_CHECKED_MOD(L, R) (((R) != 0) ? ((L) % (R)) : (flow_div_by_zero_handler(), (L) * 0))
#endif
#ifndef FLOW_CHECKED_SHL
#define FLOW_CHECKED_SHL(L, R) ((((R) >= 0) && ((unsigned long long)(R) < (sizeof(L) * 8ull)) && ((L) >= 0)) ? ((L) << (R)) : (flow_shift_ub_handler(), (L) * 0))
#endif
#ifndef FLOW_CHECKED_SHR
#define FLOW_CHECKED_SHR(L, R) ((((R) >= 0) && ((unsigned long long)(R) < (sizeof(L) * 8ull))) ? ((L) >> (R)) : (flow_shift_ub_handler(), (L) * 0))
#endif
#include <math.h>
void* _ui_state = NULL;
static inline float i32_to_f32(int32_t v) { return (float)v; }
/* Host stub for @gpu kernels (device codegen replaces this). */
static inline int32_t gpu_thread_id(void) { return 0; }
int64_t mulmod_i64_i64(int64_t a, int64_t b);
int64_t modpow_i64_i64(int64_t b, int64_t e);
int64_t hslot_ptr_i64_i64(int64_t* tk, int64_t t);
int64_t dfs_i64_i64_i64_i64_ptr_i64_ptr_i64_i64_ptr_i64_ptr_i64_ptr_i64_ptr_i64_ptr_i64(int64_t i, int64_t ri, int64_t w, int64_t d, int64_t* Q, int64_t* QI, int64_t nq, int64_t* tk, int64_t* tw, int64_t* st, int64_t* stk, int64_t* sol);
int32_t main(void);
static const int64_t P = 2000000011;
static const int64_t CAP = 239;
static const int64_t HB = 65536;
static const int64_t HMASK = 65535;
static const int64_t QMAX = 241;
int64_t mulmod_i64_i64(int64_t a, int64_t b) {
return FLOW_CHECKED_MOD(((a * b)), (P));
}
int64_t modpow_i64_i64(int64_t b, int64_t e) {
int64_t r = 1;
int64_t bb = FLOW_CHECKED_MOD((b), (P));
int64_t ee = e;
while (ee > 0) {
if (FLOW_CHECKED_MOD((ee), (2)) == 1) {
r = mulmod_i64_i64(r, bb);
}
bb = mulmod_i64_i64(bb, bb);
ee = FLOW_CHECKED_DIV((ee), (2));
}
return r;
}
int64_t hslot_ptr_i64_i64(int64_t* tk, int64_t t) {
int64_t h = ((FLOW_CHECKED_SHR((t), (15)) ^ t) & HMASK);
while ((tk[h] != 0 && tk[h] != t)) {
h = ((h + 1) & HMASK);
}
return h;
}
int64_t dfs_i64_i64_i64_i64_ptr_i64_ptr_i64_i64_ptr_i64_ptr_i64_ptr_i64_ptr_i64_ptr_i64(int64_t i, int64_t ri, int64_t w, int64_t d, int64_t* Q, int64_t* QI, int64_t nq, int64_t* tk, int64_t* tw, int64_t* st, int64_t* stk, int64_t* sol) {
int64_t t = FLOW_CHECKED_MOD((((P - 1) * ri)), (P));
int64_t h = hslot_ptr_i64_i64(tk, t);
if (tk[h] == t) {
int64_t cand = (w + tw[h]);
if (cand < st[0]) {
st[0] = cand;
st[1] = 1;
st[2] = w;
st[3] = t;
st[4] = d;
int64_t k = 0;
while (k < d) {
sol[k] = stk[k];
k = (k + 1);
}
} else {
if (cand == st[0]) {
st[1] = (st[1] + 1);
}
}
}
int64_t j = i;
while (j < nq) {
int64_t wj = (Q[j] - 1);
if ((w + wj) > st[0]) {
j = nq;
} else {
stk[d] = Q[j];
dfs_i64_i64_i64_i64_ptr_i64_ptr_i64_i64_ptr_i64_ptr_i64_ptr_i64_ptr_i64_ptr_i64(j, mulmod_i64_i64(ri, QI[j]), (w + wj), (d + 1), Q, QI, nq, tk, tw, st, stk, sol);
j = (j + 1);
}
}
return 0;
}
int32_t main(void) {
int64_t* comp = (int64_t*)(calloc((QMAX + 1), 8));
int64_t i = 2;
while ((i * i) <= QMAX) {
if (comp[i] == 0) {
int64_t j = (i * i);
while (j <= QMAX) {
comp[j] = 1;
j = (j + i);
}
}
i = (i + 1);
}
int64_t* Q = (int64_t*)(calloc(64, 8));
int64_t* QI = (int64_t*)(calloc(64, 8));
int64_t nq = 0;
i = 5;
while (i <= QMAX) {
if (((comp[i] == 0 && (i - 1) <= CAP) && i < P)) {
Q[nq] = i;
QI[nq] = modpow_i64_i64(i, (P - 2));
nq = (nq + 1);
}
i = (i + 1);
}
int64_t* tk = (int64_t*)(calloc(HB, 8));
int64_t* tw = (int64_t*)(calloc(HB, 8));
int64_t b = 0;
int64_t r3 = 1;
while ((2 * b) <= CAP) {
int64_t cur = r3;
int64_t a = 0;
while ((a + (2 * b)) <= CAP) {
int64_t h = hslot_ptr_i64_i64(tk, cur);
if (tk[h] == 0) {
tk[h] = cur;
tw[h] = (a + (2 * b));
} else {
if ((a + (2 * b)) < tw[h]) {
tw[h] = (a + (2 * b));
}
}
cur = mulmod_i64_i64(cur, 2);
a = (a + 1);
}
r3 = mulmod_i64_i64(r3, 3);
b = (b + 1);
}
int64_t* st = (int64_t*)(calloc(8, 8));
int64_t* stk = (int64_t*)(calloc(64, 8));
int64_t* sol = (int64_t*)(calloc(64, 8));
st[0] = (CAP + 1);
dfs_i64_i64_i64_i64_ptr_i64_ptr_i64_i64_ptr_i64_ptr_i64_ptr_i64_ptr_i64_ptr_i64(0, 1, 0, 0, Q, QI, nq, tk, tw, st, stk, sol);
__int128 N = ((__int128)(1));
int64_t k = 0;
while (k < st[4]) {
N = (N * ((__int128)(sol[k])));
k = (k + 1);
}
int64_t rem = (st[0] - st[2]);
int64_t tt = st[3];
int64_t ba = (-1);
int64_t bb2 = (-1);
b = 0;
r3 = 1;
while ((2 * b) <= rem) {
int64_t cur = r3;
int64_t a = 0;
while ((a + (2 * b)) <= rem) {
if ((cur == tt && (a + (2 * b)) == rem)) {
ba = a;
bb2 = b;
}
cur = mulmod_i64_i64(cur, 2);
a = (a + 1);
}
r3 = mulmod_i64_i64(r3, 3);
b = (b + 1);
}
k = 0;
while (k < ba) {
N = (N * ((__int128)(2)));
k = (k + 1);
}
k = 0;
while (k < bb2) {
N = (N * ((__int128)(3)));
k = (k + 1);
}
int64_t* digs = (int64_t*)(calloc(64, 8));
int64_t nd = 0;
__int128 x = N;
while (x > ((__int128)(0))) {
digs[nd] = ((int64_t)(FLOW_CHECKED_MOD((x), (((__int128)(10))))));
x = FLOW_CHECKED_DIV((x), (((__int128)(10))));
nd = (nd + 1);
}
int64_t j2 = (nd - 1);
while (j2 >= 0) {
printf("%c", ((int32_t)((48 + digs[j2]))));
j2 = (j2 - 1);
}
printf("\n");
free(comp);
free(Q);
free(QI);
free(tk);
free(tw);
free(st);
free(stk);
free(sol);
free(digs);
return 0;
}