Najszybszy przybliżony wspólny dzielnik

13

Przegląd

W tym wyzwaniu otrzymasz dwie liczby, które są małym przesunięciem większym niż wielokrotność liczby średniej wielkości. Musisz wypisać średnią liczbę, która jest prawie dzielnikiem obu liczb, z wyjątkiem niewielkiego przesunięcia.

Wielkość zaangażowanych numery będą programowane przez parametr trudności, l. Twoim celem jest rozwiązanie problemu w jak największym stopniu lw mniej niż 1 minutę.


Ustawiać

W danym problemie będzie tajny numer p, który będzie losową l^2( l*l) liczbą bitów. Będą dwa mnożniki, q1, q2które będą losowymi l^3liczbami bitowymi, i będą dwa przesunięcia r1, r2, które będą losowymi lliczbami bitowymi.

Dane wejściowe do Twojego programu będą x1, x2zdefiniowane jako:

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

Oto program do generowania przypadków testowych w Pythonie:

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)

Pierwszy wiersz wyniku jest możliwym rozwiązaniem, natomiast drugi i trzeci wiersz to dane wejściowe, które otrzyma Twój program.


Twój program

Biorąc pod uwagę x1, x2i ltrzeba znaleźć l^2numer bitu p'takie, że x1 % p'i x2 % p'są zarówno lnumery bitowe. pzawsze będzie działać, choć mogą istnieć inne możliwości. Oto funkcja weryfikacji rozwiązania:

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)

Przykład

Załóżmy, że ljest to 3. Generator wybiera 9-bitową liczbę p, która w tym przypadku jest 442. Generator wybiera dwa 3bitowych liczb dla r1, r2, które są 4, 7. Generator wybiera dwa 27bitowych liczb dla q1, q2, które są 117964803, 101808039. Z powodu tych wyborów x1, x252140442930, 44999153245.

Twój program byłby podawany 52140442930, 44999153245jako dane wejściowe i musi wypisać 9-bitową liczbę (w zakresie [256, 511]), tak aby 52140442930i 44999153245modulo ta liczba dawała 3-bitowe liczby (w zakresie [4, 7]). 442jest jedyną taką wartością w tym przypadku, więc twój program musiałby generować dane wyjściowe 442.


Więcej przykładów

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

Punktacja

Jak wspomniano powyżej, wynik twojego programu jest najwyższy, lktóry program ukończy w czasie krótszym niż 1 minuta. Mówiąc dokładniej, twój program będzie działał na 5 losowych instancjach z tym li musi wypisać poprawną odpowiedź na wszystkich 5, przy średnim czasie poniżej 1 minuty. Wynik programu będzie najwyższy l, jeśli odniesie sukces. Tiebreaker będzie na to średni czas l.

Aby dać ci wyobrażenie o tym, do jakich wyników chcesz dążyć, napisałem bardzo prosty solute brute-force. Otrzymał ocenę 5. Napisałem znacznie bardziej wymyślny solver. Otrzymał wynik 12 lub 13, w zależności od szczęścia.


Detale

Aby uzyskać doskonałą porównywalność między odpowiedziami, będę przesyłać czas na laptopie, aby uzyskać wyniki kanoniczne. Będę również uruchamiał te same losowo wybrane instancje dla wszystkich zgłoszeń, aby nieco złagodzić szczęście. Mój laptop ma 4 procesory, procesor i5-4300U przy 1,9 GHz, 7,5G pamięci RAM.

Możesz opublikować tymczasowy wynik na podstawie własnego harmonogramu, po prostu wyjaśnij, czy jest on tymczasowy, czy kanoniczny.


Niech wygra najszybszy program!

isaacg
źródło
Czy przybliżenie musi być najbliższe?
Yytsi
@TuukkaX Działa dowolna l^2liczba lbitów oddalona od bycia czynnikiem obu liczb. Zazwyczaj jednak jest tylko jeden.
isaacg
Oto rozprawa z kilkoma pomysłami na algorytmy: tigerprints.clemson.edu/cgi/… . Szczególnie ładny i prosty jest ten w sekcji 5.1.1
isaacg
I5-4300U ma 2 rdzenie (4 nici), a nie 4 rdzeni.
Anders Kaseorg,

Odpowiedzi:

3

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

Biegać

cargo build --release
target/release/agcd 56007668599 30611458895 3
target/release/agcd --random 13 5
Anders Kaseorg
źródło
Oficjalny wynik to l = 13, przy średnim czasie 41,53 s. l = 14 zabrał średnio nieco ponad 2 metry.
isaacg
2

Mathematica, L = 5

oto szybkie rozwiązanie 5

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

Wejście
[x1, x2, L]

Aby każdy mógł to przetestować, oto generator kluczy

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

wybierz przebieg L, kod, a otrzymasz trzy liczby.
umieść dwa ostatnie wraz z L jako wejście, a powinieneś dostać pierwszy

J42161217
źródło
Zweryfikowałem to rozwiązanie jako otrzymujące wynik l = 5. Poświęcę czas, jeśli czas jest potrzebny jako remis.
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[]]]]]) &

wejście [x1, x2, L]

Aby każdy mógł to przetestować, oto generator kluczy

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

wybierz wartość L, uruchom kod, a otrzymasz trzy liczby.
umieść dwa ostatnie wraz z L jako wejście, a powinieneś dostać pierwszy

J42161217
źródło
Kiedy to przetestowałem l = 12, dało to niepoprawny wynik. W szczególności wynikowy dzielnik był liczbą 146-bitową, która jest zbyt duża. Myślę, że będziesz potrzebować tylko niewielkiej zmiany, aby to naprawić.
isaacg
Dodałem wadliwy przypadek jako ostatni przykład powyżej. Twój solver dał odpowiedź równą 3 * p, która była nieco za duża. Patrząc na twój kod, wygląda na to, że sprawdzasz tylko, czy twoje wyjście ma co najmniej l^2bity, ale musisz także sprawdzić, czy ma co najmniej l^2bity.
isaacg
W tym samym przypadku testowym, w którym wcześniej się nie udawało, nadal nie działa. Nie jestem zbyt obeznany z Mathematicą, ale wyglądało na to, że został zakończony bez wyjścia. Myślę, że problem polega na tym, że znajduje czynnik, który jest zbyt duży, a następnie zamiast znaleźć dzielnik czynnika w pożądanym zakresie, po prostu przeskakuje obok niego.
isaacg
Oficjalny wynik to L = 12, przy średnim czasie 52,7 s. Dobra robota!
isaacg
1

Python, L = 15, średni czas 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))

Ten kod opiera się na fakcie, że iloczyn x1 - r dla wszystkich możliwych wartości r musi być wielokrotnością p, a iloczyn x2 - r musi być również. Większość czasu spędza na gcd obu produktów, z których każdy ma około 60 000 000 bitów. Gcd, który ma tylko około 250 000 bitów, jest następnie porównywany z każdą z wartości xr ponownie, aby wyodrębnić kandydatów p '. Jeśli są nieco za duże, stosuje się podział próbny, aby zredukować je do odpowiedniego zakresu.

Jest to oparte na bibliotece gmpy2 dla Pythona, która jest cienkim opakowaniem dla biblioteki GNU MP, która w szczególności ma naprawdę dobrą procedurę gcd.

Ten kod jest jednowątkowy.

isaacg
źródło