Fastest Approximate Common Divisor

Question

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.

---

Your Program

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!

Answer 1

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

Answer 2

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

Answer 3

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("answer = ", res)
        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.

This code is single-threaded.

Answer 4

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