The July 2026 Ponder This challenge is stated as follows.
Continuing the theme of last month, we deal with a movie franchise consisting of n superheroes. They are joined by n supervillains. The producers intend to pair the superheroes and the supervillains to form (hero, villain) pairs where each hero has a unique villain serving as their nemesis.
Each pairing is accepted differently by the audiences. After elaborate work, a method of assigning numerical value f(a,b) to each pairing (a,b) to estimate the audiences’ reaction was developed. The producers wish to find the list of pairings that maximizes the value of the pairing with the minimal value in the list. This minimal value is called the hero-villain value.
The way f(a,b) is computed is as follows: Let p be some prime and define a function $T_{a,b}(x) = x^2 +ax+b ,, $. By setting \(x_0 = 0\) and \(x_{n+1} = T(x_n)\) we obtain a sequence $x_0, x_1, x_2, $ which eventually repeats. Let f(a,b) be the number of steps until the first repeat happens. i.e. if \(x_n\) is the first element in the sequence such that there exists m < n for which \(x_n = x_m\), then f(a,b)=n.
For example, for n=5 and p=101, one possible list of pairings is (1,3),(2,1),(3,4),(4,2),(5,5) which yields the values 14,18,19,22,14 for which the minimum is 14. It turns out that every list of pairings gives a value of at most 14, so 14 is hero-villain value for this case.
Your goal Find the hero-villain value for n=611 and p=14411
A bonus “*” will be given for finding the optimal n in the range 1<n<N for N=1000 which gives the maximal hero-villain value for n and p=17377.
Let’s first create an efficient function for f(a,b)
library(primes)
library(igraph)## Warning: package 'igraph' was built under R version 4.6.1##
## Attaching package: 'igraph'## The following objects are masked from 'package:stats':
##
## decompose, spectrum## The following object is masked from 'package:base':
##
## union# Inputs a, b, and the prime p to be used.
f_ab <- function(a, b, p) {
seen <- logical(p)
x <- 0L
seen[x + 1] <- TRUE
count <- 1
repeat {
x <- (x * x + a * x + b) %% p
if (seen[x + 1]) {
return(count)
}
seen[x + 1] <- TRUE
count <- count + 1
}
}
## Check values
f_ab(1, 3, 101)## [1] 14f_ab(2, 1, 101)## [1] 18f_ab(3, 4, 101)## [1] 19f_ab(4, 2, 101)## [1] 22f_ab(5, 5, 101)## [1] 14We’ll build a n x n matrix that will be a handy reference for all the values of f(a,b) for a given prime p.
f_matrix <- function(n, p) {
if (!is_prime(p)) stop("p must be prime")
F <- matrix(0, nrow = n, ncol = n)
for (a in seq_len(n)) {
for (b in seq_len(n)) {
F[a, b] <- f_ab(a, b, p)
}
## A comforting message that will assure us the function is still at work
## and making progress when n is really big.
if (a %% 100 == 0) {
message("finished a = ", a, " / ", n)
}
}
F
}
f_matrix(5, 101)## [,1] [,2] [,3] [,4] [,5]
## [1,] 13 5 14 8 8
## [2,] 18 9 18 14 12
## [3,] 15 12 20 19 5
## [4,] 18 22 17 18 22
## [5,] 8 17 9 14 14Next up, a function that will have inputs n and p, the number of heroes (and super-villains) to pair with each other and the prime p we’re working with, and output the max-min over all such lists above. To create such a function, a helper function that will take in the matrix produced from the function above and a specific threshold to see if there is a pairing that goes with that threshold.
## For this function, we'll assume a, b, and p are in the global environment
has_perfect_matching <- function(F, threshold) {
n <- nrow(F)
edges <- which(F >= threshold, arr.ind = TRUE)
if (nrow(edges) < n) return(FALSE)
g <- graph_from_edgelist(
cbind(
paste0("a", edges[, 1]),
paste0("b", edges[, 2])
),
directed = FALSE
)
V(g)$type <- grepl("^b", V(g)$name)
m <- max_bipartite_match(g)
m$matching_size == n
}
HVvalue <- function(n, p) {
F <- f_matrix(n, p)
vals <- sort(unique(as.vector(F)))
lo <- 1
hi <- length(vals)
best <- vals[1]
while (lo <= hi) {
mid <- floor((lo + hi) / 2)
candidate <- vals[mid]
message("testing threshold = ", candidate)
if (has_perfect_matching(F, candidate)) {
best <- candidate
lo <- mid + 1
} else {
hi <- mid - 1
}
}
best
}
HVvalue(5, 101)## testing threshold = 14## testing threshold = 18## testing threshold = 15## [1] 14HVvalue(611,14411)## finished a = 100 / 611## finished a = 200 / 611## finished a = 300 / 611## finished a = 400 / 611## finished a = 500 / 611## finished a = 600 / 611## testing threshold = 264## testing threshold = 395## testing threshold = 329## testing threshold = 362## testing threshold = 345## testing threshold = 353## testing threshold = 349## testing threshold = 351## testing threshold = 350## [1] 349While this worked fast to obtain the max-min 14 for n=5, it took some time to find the max-min 349 answer for n = 611 and p = 14411. We’ll need some significant efficiencies in order to optimize over all n for the fixed p = 17377.
One place we can do this is to compute an entire table F[a,b,17377] for all a and b between 2 and 999 inclusive. Then we can search over n using top-left submatrices.
Using the Rcpp package in R, we can use a much faster compiled version.
library(Rcpp)
Rcpp::sourceCpp(code = '
#include <Rcpp.h>
#include <vector>
#include <queue>
#include <functional>
#include <algorithm>
#include <climits>
using namespace Rcpp;
// [[Rcpp::export]]
IntegerMatrix f_table_cpp(int N, int p) {
IntegerMatrix F(N, N);
std::vector<int> seen(p, 0);
int stamp = 1;
for (int a = 1; a <= N; ++a) {
for (int b = 1; b <= N; ++b) {
if (stamp == INT_MAX) {
std::fill(seen.begin(), seen.end(), 0);
stamp = 1;
}
int x = 0;
int count = 1;
seen[0] = stamp;
while (true) {
x = (int)((1LL * x * x + 1LL * a * x + b) % p);
if (seen[x] == stamp) break;
seen[x] = stamp;
++count;
}
F(a - 1, b - 1) = count;
++stamp;
}
if (a % 25 == 0) {
Rcout << "finished a = " << a << " / " << N << "\\n";
}
}
return F;
}
bool has_matching_threshold_inner(const IntegerMatrix& F, int n, int threshold) {
std::vector<int> row_degree(n, 0), col_degree(n, 0);
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n; ++j) {
if (F(i, j) >= threshold) {
++row_degree[i];
++col_degree[j];
}
}
}
for (int i = 0; i < n; ++i) {
if (row_degree[i] == 0 || col_degree[i] == 0) return false;
}
std::vector<int> pairU(n, -1), pairV(n, -1), dist(n);
auto bfs = [&]() {
std::queue<int> q;
bool found_free_col = false;
for (int u = 0; u < n; ++u) {
if (pairU[u] == -1) {
dist[u] = 0;
q.push(u);
} else {
dist[u] = -1;
}
}
while (!q.empty()) {
int u = q.front();
q.pop();
for (int v = 0; v < n; ++v) {
if (F(u, v) < threshold) continue;
int next_u = pairV[v];
if (next_u == -1) {
found_free_col = true;
} else if (dist[next_u] == -1) {
dist[next_u] = dist[u] + 1;
q.push(next_u);
}
}
}
return found_free_col;
};
std::function<bool(int)> dfs = [&](int u) {
for (int v = 0; v < n; ++v) {
if (F(u, v) < threshold) continue;
int next_u = pairV[v];
if (next_u == -1 ||
(dist[next_u] == dist[u] + 1 && dfs(next_u))) {
pairU[u] = v;
pairV[v] = u;
return true;
}
}
dist[u] = -1;
return false;
};
int matching = 0;
while (bfs()) {
for (int u = 0; u < n; ++u) {
if (pairU[u] == -1 && dfs(u)) {
++matching;
}
}
}
return matching == n;
}
// [[Rcpp::export]]
bool has_matching_threshold_cpp(IntegerMatrix F, int n, int threshold) {
return has_matching_threshold_inner(F, n, threshold);
}
// [[Rcpp::export]]
List best_n_cpp(IntegerMatrix F) {
int N = F.nrow();
int maxF = 0;
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
if (F(i, j) > maxF) maxF = F(i, j);
}
}
int best_n = 1;
int best_value = 0;
int matching_checks = 0;
for (int n = 1; n <= N; ++n) {
if (n % 25 == 0) {
Rcout << "checking n = " << n
<< ", current best = " << best_value
<< " at n = " << best_n << "\\n";
}
if (!has_matching_threshold_inner(F, n, best_value + 1)) {
++matching_checks;
continue;
}
++matching_checks;
int lo = best_value + 1;
int hi = maxF;
int local_best = best_value;
while (lo <= hi) {
int mid = lo + (hi - lo) / 2;
++matching_checks;
if (has_matching_threshold_inner(F, n, mid)) {
local_best = mid;
lo = mid + 1;
} else {
hi = mid - 1;
}
}
if (local_best > best_value) {
best_value = local_best;
best_n = n;
Rcout << "new best: HVvalue = " << best_value
<< " at n = " << best_n << "\\n";
}
}
return List::create(
Named("best_n") = best_n,
Named("best_value") = best_value,
Named("matching_checks") = matching_checks
);
}
')We’ll run the code above and save it as an rds file.
p <- 17377
N <- 999
F <- f_table_cpp(N, p)## finished a = 25 / 999
## finished a = 50 / 999
## finished a = 75 / 999
## finished a = 100 / 999
## finished a = 125 / 999
## finished a = 150 / 999
## finished a = 175 / 999
## finished a = 200 / 999
## finished a = 225 / 999
## finished a = 250 / 999
## finished a = 275 / 999
## finished a = 300 / 999
## finished a = 325 / 999
## finished a = 350 / 999
## finished a = 375 / 999
## finished a = 400 / 999
## finished a = 425 / 999
## finished a = 450 / 999
## finished a = 475 / 999
## finished a = 500 / 999
## finished a = 525 / 999
## finished a = 550 / 999
## finished a = 575 / 999
## finished a = 600 / 999
## finished a = 625 / 999
## finished a = 650 / 999
## finished a = 675 / 999
## finished a = 700 / 999
## finished a = 725 / 999
## finished a = 750 / 999
## finished a = 775 / 999
## finished a = 800 / 999
## finished a = 825 / 999
## finished a = 850 / 999
## finished a = 875 / 999
## finished a = 900 / 999
## finished a = 925 / 999
## finished a = 950 / 999
## finished a = 975 / 999saveRDS(F, "F_p17377_N999.rds")
ans <- best_n_cpp(F)## new best: HVvalue = 273 at n = 1
## checking n = 25, current best = 273 at n = 1
## new best: HVvalue = 277 at n = 35
## new best: HVvalue = 285 at n = 45
## checking n = 50, current best = 285 at n = 45
## new best: HVvalue = 303 at n = 58
## new best: HVvalue = 308 at n = 66
## new best: HVvalue = 314 at n = 73
## new best: HVvalue = 317 at n = 74
## checking n = 75, current best = 317 at n = 74
## new best: HVvalue = 323 at n = 78
## checking n = 100, current best = 323 at n = 78
## checking n = 125, current best = 323 at n = 78
## checking n = 150, current best = 323 at n = 78
## new best: HVvalue = 335 at n = 156
## checking n = 175, current best = 335 at n = 156
## new best: HVvalue = 350 at n = 193
## checking n = 200, current best = 350 at n = 193
## checking n = 225, current best = 350 at n = 193
## new best: HVvalue = 363 at n = 237
## checking n = 250, current best = 363 at n = 237
## new best: HVvalue = 364 at n = 259
## checking n = 275, current best = 364 at n = 259
## checking n = 300, current best = 364 at n = 259
## checking n = 325, current best = 364 at n = 259
## checking n = 350, current best = 364 at n = 259
## checking n = 375, current best = 364 at n = 259
## checking n = 400, current best = 364 at n = 259
## new best: HVvalue = 381 at n = 419
## checking n = 425, current best = 381 at n = 419
## checking n = 450, current best = 381 at n = 419
## checking n = 475, current best = 381 at n = 419
## new best: HVvalue = 385 at n = 495
## checking n = 500, current best = 385 at n = 495
## checking n = 525, current best = 385 at n = 495
## checking n = 550, current best = 385 at n = 495
## new best: HVvalue = 389 at n = 574
## checking n = 575, current best = 389 at n = 574
## checking n = 600, current best = 389 at n = 574
## checking n = 625, current best = 389 at n = 574
## checking n = 650, current best = 389 at n = 574
## new best: HVvalue = 400 at n = 671
## checking n = 675, current best = 400 at n = 671
## checking n = 700, current best = 400 at n = 671
## checking n = 725, current best = 400 at n = 671
## checking n = 750, current best = 400 at n = 671
## new best: HVvalue = 405 at n = 753
## checking n = 775, current best = 405 at n = 753
## checking n = 800, current best = 405 at n = 753
## checking n = 825, current best = 405 at n = 753
## checking n = 850, current best = 405 at n = 753
## checking n = 875, current best = 405 at n = 753
## checking n = 900, current best = 405 at n = 753
## new best: HVvalue = 408 at n = 924
## checking n = 925, current best = 408 at n = 924
## checking n = 950, current best = 408 at n = 924
## checking n = 975, current best = 408 at n = 924ans## $best_n
## [1] 924
##
## $best_value
## [1] 408
##
## $matching_checks
## [1] 1196After some computation time, we find the max-min of 408 for n = 924 and p = 17377.