Schnellster ungefährer gemeinsamer Divisor

13

Überblick

In dieser Herausforderung erhalten Sie zwei Zahlen, die jeweils einen kleinen Versatz größer als ein Vielfaches einer mittelgroßen Zahl sind. Sie müssen eine mittelgroße Zahl ausgeben, die bis auf einen kleinen Versatz fast ein Teiler beider Zahlen ist.

Die Größe der beteiligten Zahlen wird durch einen Schwierigkeitsparameter parametrisiert l. Ihr Ziel ist es, das Problem lin weniger als 1 Minute so schnell wie möglich zu lösen .


Installieren

In einem gegebenen Problem wird es eine Geheimzahl geben p, die eine zufällige l^2( l*l) Bitnummer ist. Es wird zwei Multiplikatoren geben, q1, q2bei denen es sich um zufällige l^3Bitnummern handelt, und es werden zwei Offsets geben, r1, r2bei denen es sich um zufällige lBitnummern handelt.

Die Eingabe für Ihr Programm x1, x2lautet wie folgt :

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

Hier ist ein Programm zum Generieren von Testfällen 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)

Die erste Ausgabezeile ist eine mögliche Lösung, während die zweite und dritte Zeile die Eingabe sind, die Ihrem Programm gegeben wird.


Ihr Programm

Gegeben x1, x2und l, müssen Sie eine l^2Bitnummer finden, p'so dass x1 % p'und x2 % p'beide lBitnummern sind. pwird immer funktionieren, obwohl es andere Möglichkeiten geben kann. Hier ist eine Funktion zum Überprüfen einer Lösung:

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)

Beispiel

Angenommen, es list 3. Das Generatorprogramm wählt eine 9-Bit-Zahl aus p, in diesem Fall ist dies 442. Der Generator wählt zwei 3Bitnummern aus, für r1, r2die 4, 7. Der Generator wählt zwei 27Bitnummern aus, für q1, q2die 117964803, 101808039. Wegen dieser Entscheidungen x1, x2sind 52140442930, 44999153245.

Ihr Programm würde 52140442930, 44999153245als Eingabe angegeben und muss eine 9-Bit-Zahl (im Bereich [256, 511]) ausgeben, sodass 52140442930und 44999153245modulo diese Zahl 3-Bit-Zahlen (im Bereich [4, 7]) ergibt . 442ist der einzige solche Wert in diesem Fall, so müsste Ihr Programm ausgeben 442.


Mehr Beispiele

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

Wertung

Wie oben erwähnt, ist die Punktzahl Ihres Programms die höchste l, die das Programm in weniger als 1 Minute abgeschlossen hat. Genauer gesagt, Ihr Programm wird mit 5 zufälligen Instanzen ausgeführt l, und es muss eine korrekte Antwort auf alle 5 mit einer durchschnittlichen Zeit unter 1 Minute ausgeben. Die Punktzahl eines Programms ist die höchste l, bei der es erfolgreich ist. Tiebreaker wird in dieser Zeit durchschnittlich sein l.

Um Ihnen eine Vorstellung davon zu geben, welche Ergebnisse angestrebt werden sollen, habe ich einen sehr einfachen Brute-Force-Löser geschrieben. Es wurde mit 5 bewertet. Ich habe einen viel schickeren Löser geschrieben. Es hat eine Punktzahl von 12 oder 13, je nach Glück.


Einzelheiten

Um eine perfekte Vergleichbarkeit der Antworten zu gewährleisten, werde ich die Einsendungen auf meinem Laptop zeitlich festlegen, um kanonische Ergebnisse zu erhalten. Ich werde bei allen Einsendungen die gleichen zufällig ausgewählten Instanzen ausführen, um das Glück etwas zu lindern. Mein Laptop hat 4 CPUs, i5-4300U CPU bei 1,9 GHz, 7,5 G RAM.

Es steht Ihnen frei, eine vorläufige Punktzahl zu veröffentlichen, die auf Ihrem eigenen Timing basiert. Machen Sie einfach klar, ob es sich um eine vorläufige oder eine kanonische Punktzahl handelt.


Möge das schnellste Programm gewinnen!

isaacg
quelle
Muss die Annäherung die nächste sein?
Yytsi
@TuukkaX Jede l^2Bit-Zahl, die lkein Faktor für beide Zahlen ist, funktioniert. Normalerweise gibt es jedoch nur einen.
isaacg
Hier ist eine Dissertation mit einigen Algorithmus-Ideen: tigerprints.clemson.edu/cgi/… . Besonders schön und einfach ist die in Abschnitt 5.1.1
isaacg
Der i5-4300U verfügt über 2 Kerne (4 Threads) und nicht über 4 Kerne.
Anders Kaseorg

Antworten:

3

Rost, 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"

Laufen

cargo build --release
target/release/agcd 56007668599 30611458895 3
target/release/agcd --random 13 5
Anders Kaseorg
quelle
Offizielles Ergebnis ist l = 13, mit einer durchschnittlichen Zeit von 41,53s. l = 14 dauerte im Durchschnitt etwas mehr als 2 m.
Isaacg
2

Mathematica, L = 5

Hier ist eine schnelle Lösung

(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]]
) &

Eingang
[x1, x2, L]

Damit das jeder testen kann, hier der Schlüsselgenerator

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

Wählen Sie den L-Wert, den Code und Sie erhalten drei Zahlen.
platziere die letzten beiden zusammen mit L als Eingabe, und du solltest die erste bekommen

J42161217
quelle
Ich habe diese Lösung mit l = 5 bewertet. Ich werde es zeitlich festlegen, wenn das Timing als Tiebreaker benötigt wird.
isaacg
1

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[]]]]]) &

Eingang [x1, x2, L]

Damit das jeder testen kann, hier der Schlüsselgenerator

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

Wählen Sie den L-Wert, führen Sie den Code aus und Sie erhalten drei Zahlen.
platziere die letzten beiden zusammen mit L als Eingabe, und du solltest die erste bekommen

J42161217
quelle
Als ich das mit getestet habe l = 12, gab es ein falsches Ergebnis. Insbesondere war der resultierende Divisor eine 146-Bit-Zahl, die zu groß ist. Ich denke, Sie brauchen nur eine kleine Änderung, um dies zu beheben.
isaacg
Ich habe den fehlgeschlagenen Fall als letztes Beispiel oben hinzugefügt. Ihr Löser gab eine Antwort in Höhe von 3 * p, die etwas zu groß war. Wenn Sie Ihren Code betrachten, sehen Sie so aus, als ob Sie nur überprüfen, ob Ihre Ausgabe mindestens l^2Bits enthält, aber Sie müssen auch überprüfen, ob sie höchstens l^2Bits enthält.
isaacg
In demselben Testfall, in dem es zuvor fehlgeschlagen ist, funktioniert es immer noch nicht. Ich bin nicht besonders vertraut mit Mathematica, aber es sah so aus, als würde es ohne Ausgabe beendet. Ich denke, das Problem ist, dass ein Faktor gefunden wird, der zu groß ist. Anstatt einen Teiler des Faktors im gewünschten Bereich zu finden, überspringt man ihn einfach.
isaacg
Lassen Sie uns diese Diskussion im Chat fortsetzen .
Isaacg
Die offizielle Punktzahl ist L = 12 mit einer durchschnittlichen Zeit von 52,7 Sekunden. Gut gemacht!
Isaacg
1

Python, L = 15, durchschnittliche Zeit 39,9 s

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))

Dieser Code basiert auf der Tatsache, dass das Produkt von x1 - r für alle möglichen Werte von r ein Vielfaches von p sein muss, und das Produkt von x2 - r muss es auch sein. Die meiste Zeit wird damit verbracht, den gcd der beiden Produkte zu nehmen, von denen jedes ungefähr 60.000.000 Bits hat. Die gcd, die nur ungefähr 250.000 Bits hat, wird dann noch einmal mit jedem der xr-Werte verglichen, um p'-Kandidaten zu extrahieren. Wenn sie etwas zu groß sind, werden sie mithilfe der Testdivision auf den entsprechenden Bereich reduziert.

Dies basiert auf der gmpy2-Bibliothek für Python, die ein dünner Wrapper für die GNU MP-Bibliothek ist, die insbesondere eine wirklich gute gcd-Routine hat.

Dieser Code ist Singlethread.

isaacg
quelle