# Fastest Approximate Common Divisor

Overview

In this challenge, you will be given two numbers which are both a small offset larger than a multiple of a medium-size number. You must output a medium-sized number that is almost a divisor of both of the numbers, except for a small offset.

The size of the numbers involved will be parameterized by a difficulty parameter, l. Your objective is to solve the problem for the largest possible l in under 1 minute.

Setup

In a given problem, there will be a secret number, p, which will be a random l^2 (l*l) bit number. There will be two multipliers, q1, q2, which will be random l^3 bit numbers, and there will be two offsets, r1, r2, which will be random l bit numbers.

The input to your program will be x1, x2, defined as:

x1 = p * q1 + r1
x2 = p * q2 + r2


Here's a program to generate test cases, in Python:

from random import randrange
from sys import argv

l = int(argv[1])

def randbits(bits):
return randrange(2 ** (bits - 1), 2 ** bits)

p = randbits(l ** 2)
print(p)
for i in range(2):
q_i = randbits(l ** 3)
r_i = randbits(l)
print(q_i * p + r_i)


The first line of output is a possible solution, while the second and third lines are the input that your program will be given.

Given x1, x2 and l, you must find a l^2 bit number p' such that x1 % p' and x2 % p' are both l bit numbers. p will always work, though there may be other possibilities. Here's a function to verify a solution:

def is_correct(x1, x2, l, p_prime):
p_prime_is_good = p_prime >> (l**2 - 1) and not p_prime >> l ** 2
x1_is_good = (x1 % p_prime) >> (l-1) and not (x1 % p_prime) >> l
x2_is_good = (x2 % p_prime) >> (l-1) and not (x2 % p_prime) >> l
return bool(p_prime_is_good and x1_is_good and x2_is_good)


Example

Suppose l is 3. The generator program picks a 9-bit number for p, which in this case is 442. The generator picks two 3 bit numbers for r1, r2, which are 4, 7. The generator picks two 27 bit numbers for q1, q2, which are 117964803, 101808039. Because of these choices, x1, x2 are 52140442930, 44999153245.

Your program would be given 52140442930, 44999153245 as input, and must output a 9-bit number (in the range [256, 511]) such that 52140442930 and 44999153245 modulo that number give 3 bit numbers (in the range [4, 7]). 442 is the only such value in this case, so your program would have to output 442.

More examples

l = 2
x1 = 1894
x2 = 2060
p = 11
No other p'.

l = 3
x1 = 56007668599
x2 = 30611458895
p = 424
No other p'.

l = 6
x1 = 4365435975875889219149338064474396898067189178953471159903352227492495111071
x2 = 6466809655659049447127736275529851894657569985804963410176865782113074947167
p = 68101195620
I don't know whether there are other p'.

l = 12
x1 = 132503538560485423319724633262218262792296147003813662398252348727558616998821387759658729802732555377599590456096450977511271450086857949046098328487779612488702544062780731169071526325427862701033062986918854245283037892816922645703778218888876645148150396130125974518827547039720412359298502758101864465267219269598121846675000819173555118275197412936184329860639224312426860362491131729109976241526141192634523046343361089218776687819810873911761177080056675776644326080790638190845283447304699879671516831798277084926941086929776037986892223389603958335825223
x2 = 131643270083452525545713630444392174853686642378302602432151533578354175874660202842105881983788182087244225335788180044756143002547651778418104898394856368040582966040636443591550863800820890232349510212502022967044635049530630094703200089437589000344385691841539471759564428710508659169951391360884974854486267690231936418935298696990496810984630182864946252125857984234200409883080311780173125332191068011865349489020080749633049912518609380810021976861585063983190710264511339441915235691015858985314705640801109163008926275586193293353829677264797719957439635
p = 12920503469397123671484716106535636962543473
I don't know whether there are other p'.

l = 12
x1 = 202682323504122627687421150801262260096036559509855209647629958481910539332845439801686105377638207777951377858833355315514789392768449139095245989465034831121409966815913228535487871119596033570221780568122582453813989896850354963963579404589216380209702064994881800638095974725735826187029705991851861437712496046570494304535548139347915753682466465910703584162857986211423274841044480134909827293577782500978784365107166584993093904666548341384683749686200216537120741867400554787359905811760833689989323176213658734291045194879271258061845641982134589988950037
x2 = 181061672413088057213056735163589264228345385049856782741314216892873615377401934633944987733964053303318802550909800629914413353049208324641813340834741135897326747139541660984388998099026320957569795775586586220775707569049815466134899066365036389427046307790466751981020951925232623622327618223732816807936229082125018442471614910956092251885124883253591153056364654734271407552319665257904066307163047533658914884519547950787163679609742158608089946055315496165960274610016198230291033540306847172592039765417365770579502834927831791804602945514484791644440788
p = 21705376375228755718179424140760701489963164


Scoring

As mentioned above, your program's score is the highest l that the program completes in under 1 minute. More specifically, your program will be run on 5 random instances with that l, and it must output a correct answer on all 5, with an average time under 1 minute. A program's score will be the highest l that it succeeds on. Tiebreaker will be average time on that l.

To give you an idea of what scores to aim for, I wrote a very simple brute-force solver. It got a score of 5. I wrote a much fancier solver. It got a score of 12 or 13, depending on luck.

Details

For perfect comparability across answers, I will time submissions on my laptop to give canonical scores. I will also run the same randomly chosen instances on all submissions, to alleviate luck somewhat. My laptop has 4 CPUs, i5-4300U CPU @ 1.9 GHz, 7.5G of RAM.

Feel free to post a provisional score based on your own timing, just make it clear whether it's provisional or canonical.

May the fastest program win!

• Does the approximation have to be the closest one? Commented May 30, 2017 at 17:59
• @TuukkaX Any l^2 bit number that's l-bits away from being a factor of both numbers works. There's typically only one, however. Commented May 30, 2017 at 18:08
• Here's a dissertation with some algorithm ideas: tigerprints.clemson.edu/cgi/…. Particularly nice and simple is the one in section 5.1.1 Commented May 31, 2017 at 10:29
• The i5-4300U has 2 cores (4 threads), not 4 cores. Commented May 31, 2017 at 22:11

# Rust, L=13

### src/main.rs

extern crate gmp;
extern crate rayon;

use gmp::mpz::Mpz;
use gmp::rand::RandState;
use rayon::prelude::*;
use std::cmp::max;
use std::env;
use std::mem::swap;

// Solver

fn solve(x1: &Mpz, x2: &Mpz, l: usize) -> Option<Mpz> {
let m = l*l - 1;
let r = -1i64 << l-2 .. 1 << l-2;
let (mut x1, mut x2) = (x1 - (3 << l-2), x2 - (3 << l-2));
let (mut a1, mut a2, mut b1, mut b2) = (Mpz::one(), Mpz::zero(), Mpz::zero(), Mpz::one());
while !x2.is_zero() &&
&(max(a1.abs(), b1.abs()) << l-2) < &x1 &&
&(max(a2.abs(), b2.abs()) << l-2) < &x2
{
let q = &x1/&x2;
swap(&mut x1, &mut x2);
x2 -= &q*&x1;
swap(&mut a1, &mut a2);
a2 -= &q*&a1;
swap(&mut b1, &mut b2);
b2 -= &q*&b1;
}
r.clone().into_par_iter().map(|u| {
let (mut y1, mut y2) = (&x1 - &a1*u + (&b1 << l-2), &x2 - &a2*u + (&b2 << l-2));
for _ in r.clone() {
let d = Mpz::gcd(&y1, &y2);
if d.bit_length() >= m {
let d = d.abs();
let (mut k, k1) = (&d >> l*l, &d >> l*l-1);
while k < k1 {
k += 1;
if (&d%&k).is_zero() {
return Some(&d/&k)
}
}
}
y1 -= &b1;
y2 -= &b2;
}
None
}).find_any(|o| o.is_some()).unwrap_or(None)
}

// Test harness

fn randbits(rnd: &mut RandState, bits: usize) -> Mpz {
rnd.urandom(&(Mpz::one() << bits - 1)) + (Mpz::one() << bits - 1)
}

fn gen(rnd: &mut RandState, l: usize) -> (Mpz, Mpz, Mpz) {
let p = randbits(rnd, l*l);
(randbits(rnd, l*l*l)*&p + randbits(rnd, l),
randbits(rnd, l*l*l)*&p + randbits(rnd, l),
p)
}

fn is_correct(x1: &Mpz, x2: &Mpz, l: usize, p_prime: &Mpz) -> bool {
p_prime.bit_length() == l*l &&
(x1 % p_prime).bit_length() == l &&
(x2 % p_prime).bit_length() == l
}

fn random_test(l: usize, n: usize) {
let mut rnd = RandState::new();  // deterministic seed
for _ in 0..n {
let (x1, x2, p) = gen(&mut rnd, l);
println!("x1 = {}", x1);
println!("x2 = {}", x2);
println!("l = {}", l);
println!("p = {}", p);
println!("x1 % p = {}", &x1 % &p);
println!("x2 % p = {}", &x2 % &p);
assert!(is_correct(&x1, &x2, l, &p));
let p_prime = solve(&x1, &x2, l).expect("no solution found!");
println!("p' = {}", p_prime);
assert!(is_correct(&x1, &x2, l, &p_prime));
println!("correct");
}
}

// Command line interface

fn main() {
let args = &env::args().collect::<Vec<_>>();
if args.len() == 4 && args[1] == "--random" {
if let (Ok(l), Ok(n)) = (args[2].parse(), args[3].parse()) {
return random_test(l, n);
}
}
if args.len() == 4 {
if let (Ok(x1), Ok(x2), Ok(l)) = (args[1].parse(), args[2].parse(), args[3].parse()) {
match solve(&x1, &x2, l) {
None => println!("no solution found"),
Some(p_prime) => println!("{}", p_prime),
}
return;
}
}
println!("Usage:");
println!("    {} --random L N", args[0]);
println!("    {} X1 X2 L", args[0]);
}


### Cargo.toml

[package]
name = "agcd"
version = "0.1.0"
authors = ["Anders Kaseorg <[email protected]>"]

[dependencies]
rayon = "0.7.1"
rust-gmp = "0.5.0"


### To run

cargo build --release
target/release/agcd 56007668599 30611458895 3
target/release/agcd --random 13 5

• Official result is l=13, with an average time of 41.53s. l=14 took a little over 2m on average. Commented Jun 1, 2017 at 7:58

# Mathematica, L=5

here is a quick 5 solution

(l = #3;
a = #1 - Range[2^(l - 1), 2^l];
b = #2 - Range[2^(l - 1), 2^l];
Last@Intersection[
Flatten[Table[Select[Divisors@a[[i]], # < 2^l^2 &], {i, Length@a}],
1],
Flatten[Table[Select[Divisors@b[[i]], # < 2^l^2 &], {i, Length@b}],
1]]
) &


Input
[x1,x2,L]

in order that anybody can test this, here is the key generator

l = 5;
a = RandomInteger[{2^(l^2 - 1), 2^(l^2)}]
b = RandomInteger[{2^(l^3 - 1), 2^(l^3)}];
c = RandomInteger[{2^(l - 1), 2^l - 1}];
f = RandomInteger[{2^(l^3 - 1), 2^(l^3)}];
g = RandomInteger[{2^(l - 1), 2^l - 1}];
x = a*b + c
y = a*f + g


choose the L value run, the code and you will get three numbers.
place the last two along with L as input, and you should get the first one

• I've verified this solution as getting a score of l=5. I'll time it if timing is needed as a tiebreaker. Commented May 31, 2017 at 11:08

# Python, L=15, average time 39.9s

from sys import argv
from random import seed, randrange

from gmpy2 import gcd, mpz
import gmpy2

def mult_buckets(buckets):
if len(buckets) == 1:
return buckets[0]
new_buckets = []
for i in range(0, len(buckets), 2):
if i == len(buckets) - 1:
new_buckets.append(buckets[i])
else:
new_buckets.append(buckets[i] * buckets[i+1])
return mult_buckets(new_buckets)

def get_products(xs, l):
num_buckets = 1000
lower_r = 1 << l - 1
upper_r = 1 << l
products = []
bucket_size = (upper_r - lower_r)//num_buckets + 1
for x in xs:
buckets = []
for bucket_num in range(num_buckets):
product = mpz(1)
for r in range(lower_r + bucket_num * bucket_size,
min(upper_r, lower_r + (bucket_num + 1) * bucket_size)):
product *= x - mpz(r)
buckets.append(product)
products.append(mult_buckets(buckets))
return products

def solve(xs, l):
lower_r = 2**(l - 1)
upper_r = 2**l
lower_p = 2**(l**2 - 1)
upper_p = 2**(l**2)

products = get_products(xs, l)
overlap = gcd(*products)
candidate_lists = []
for x in xs:
candidates = []
candidate_lists.append(candidates)
for r in range(lower_r, upper_r):
common_divisor = gcd(x-r, overlap)
if common_divisor >= lower_p:
candidates.append(common_divisor)
to_reduce = []
for l_candidate in candidate_lists[0]:
for r_candidate in candidate_lists[1]:
best_divisor = gcd(l_candidate, r_candidate)
if lower_p <= best_divisor < upper_p:
return best_divisor
elif best_divisor > upper_p:
to_reduce.append(best_divisor)
for divisor in to_reduce:
cutoff = divisor // lower_p
non_divisors = []
for sub_divisor in range(2, cutoff + 1):
if any(sub_divisor % non_divisor == 0 for non_divisor in non_divisors):
continue
quotient, remainder = gmpy2.c_divmod(divisor, sub_divisor)
if remainder == 0 and lower_p <= quotient < upper_p:
return quotient
if quotient < lower_p:
break
if remainder != 0:
non_divisors.append(sub_divisor)

def is_correct(x1, x2, l, p_prime):
p_prime_is_good = p_prime >> (l**2 - 1) and not p_prime >> l ** 2
x1_is_good = (x1 % p_prime) >> (l-1) and not (x1 % p_prime) >> l
x2_is_good = (x2 % p_prime) >> (l-1) and not (x2 % p_prime) >> l
return bool(p_prime_is_good and x1_is_good and x2_is_good)

if __name__ == '__main__':
seed(0)

l = int(argv[1])
reps = int(argv[2])

def randbits(bits):
return randrange(2 ** (bits - 1), 2 ** bits)

for _ in range(reps):
p = randbits(l ** 2)
print("p = ", p)
xs = []
for i in range(2):
q_i = randbits(l ** 3)
print("q", i, "=", q_i)
r_i = randbits(l)
print("r", i, "=", r_i)
x_i = q_i * p + r_i
print("x", i, "=", x_i)
xs.append(x_i)

res = solve(xs, l)
print("Correct?", is_correct(xs[0], xs[1], l, res))


This code is based on the fact that the product of x1 - r for all possible values of r must be a multiple of p, and the product of x2 - r must be as well. Most of the time is spent taking the gcd of the two products, each of which has about 60,000,000 bits. The gcd, which only has about 250,000 bits, is then compared to each of the x-r values once again, to extract p' candidates. If they're a bit too large, trial division is used to reduce them to the appropriate range.

This is based on the gmpy2 library for Python, which is a thin wrapper for the GNU MP library, which in particular has a really good gcd routine.

# Mathematica, L=12

ClearAll[l, a, b, k, m];
(l = #3;
a = #1 - Range[2^(l - 1), 2^l];
b = #2 - Range[2^(l - 1), 2^l];
k = Length@a;
m = Length@b;
For[i = 1, i <= k, i++,
For[j = 1, j <= m, j++, If[2^(l^2 - 1) < GCD[a[[i]], b[[j]]],
If[GCD[a[[i]], b[[j]]] > 2^l^2,
Print@Max@
Select[Divisors[GCD[a[[i]], b[[j]]]],
2^(l^2 - 1) < # < 2^l^2 &]; Abort[],
Print[GCD[a[[i]], b[[j]]]];
Abort[]]]]]) &


input [x1,x2,L]

in order that anybody can test this, here is the key generator

l = 12;
a = RandomInteger[{2^(l^2 - 1), 2^(l^2)}]
b = RandomInteger[{2^(l^3 - 1), 2^(l^3)}];
c = RandomInteger[{2^(l - 1), 2^l - 1}];
f = RandomInteger[{2^(l^3 - 1), 2^(l^3)}];
g = RandomInteger[{2^(l - 1), 2^l - 1}];
x = a*b + c
y = a*f + g


choose the L value, run the code and you will get three numbers.
place the last two along with L as input, and you should get the first one

• When I tested this with l = 12, it gave an incorrect result. Specifically, the resulting divisor was a 146-bit number, which is too large. I think you'll only need a small change to fix this. Commented May 31, 2017 at 22:13
• I added the failing case as the final example above. Your solver gave an answer equal to 3*p, which was a little bit too big. Looking at your code, it looks like you only check that your output has at least l^2 bits, but you also need to check that it has at most l^2 bits. Commented May 31, 2017 at 22:47
• On the same test case that it failed before, it's still not working. I'm not super familiar with Mathematica, but it looked like it exited with no output. I think the problem is that it's finding a factor that's too large, and then instead of finding a divisor of the factor in the desired range, it's just skipping past it. Commented May 31, 2017 at 23:18
• Commented Jun 1, 2017 at 0:02
• Official score is L=12, with an average time of 52.7s. Well done! Commented Jun 1, 2017 at 14:33