Divisor comum aproximado mais rápido

13

Visão geral

Nesse desafio, você receberá dois números, ambos com um pequeno deslocamento maior que um múltiplo de um número de tamanho médio. Você deve produzir um número de tamanho médio que seja quase um divisor de ambos os números, exceto por um pequeno deslocamento.

O tamanho dos números envolvidos será parametrizado por um parâmetro de dificuldade l,. Seu objetivo é resolver o problema o maior possível lem menos de 1 minuto.


Configuração

Em um determinado problema, haverá um número secreto p, que será um número de bit aleatório l^2( l*l). Haverá dois multiplicadores, q1, q2que serão l^3números de bits aleatórios , e haverá dois deslocamentos r1, r2, que serão lnúmeros de bits aleatórios .

A entrada para o seu programa será x1, x2definida como:

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

Aqui está um programa para gerar casos de teste, em 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)

A primeira linha de saída é uma solução possível, enquanto a segunda e a terceira linhas são a entrada que seu programa receberá.


O seu programa

Dado x1, x2e l, você deve encontrar um l^2número de bits p'tal que x1 % p'e x2 % p'sejam ambos lnúmeros de bits. psempre funcionará, embora possa haver outras possibilidades. Aqui está uma função para verificar uma solução:

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)

Exemplo

Suponha que lseja 3. O programa gerador escolhe um número de 9 bits para p, que neste caso é 442. O gerador escolhe 3números de dois bits para r1, r2, quais são 4, 7. O gerador escolhe 27números de dois bits para q1, q2, quais são 117964803, 101808039. Por causa dessas escolhas, x1, x2são 52140442930, 44999153245.

Seu programa seria dado 52140442930, 44999153245como entrada e deve gerar um número de 9 bits (no intervalo [256, 511]) tal que, 52140442930e 44999153245esse módulo, forneça números de 3 bits (no intervalo [4, 7]). 442é o único valor nesse caso; portanto, seu programa precisaria produzir 442.


Mais exemplos

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

Pontuação

Como mencionado acima, a pontuação do seu programa é a mais alta possível lem menos de 1 minuto. Mais especificamente, seu programa será executado em 5 instâncias aleatórias com isso le deve gerar uma resposta correta em todas as 5, com um tempo médio inferior a 1 minuto. A pontuação de um programa será a mais alta lpossível. O desempate será o tempo médio para isso l.

Para ter uma idéia de quais pontuações apontar, escrevi um solucionador de força bruta muito simples. Obteve uma pontuação de 5. Escrevi um solucionador muito mais sofisticado. Tem uma pontuação de 12 ou 13, dependendo da sorte.


Detalhes

Para uma perfeita comparabilidade entre as respostas, cronometrarei os envios no meu laptop para fornecer pontuações canônicas. Também executarei as mesmas instâncias escolhidas aleatoriamente em todos os envios, para aliviar um pouco a sorte. Meu laptop possui 4 CPUs, i5-4300U CPU a 1,9 GHz, 7,5G de RAM.

Sinta-se à vontade para postar uma pontuação provisória com base no seu próprio tempo, apenas deixe claro se é provisória ou canônica.


Que ganhe o programa mais rápido!

isaacg
fonte
A aproximação precisa ser a mais próxima?
Yytsi
@TuukkaX Qualquer l^2número de bit que esteja a lmenos de um fator de ambos os números funciona. No entanto, normalmente há apenas um.
isaacg
Aqui está uma dissertação com algumas idéias de algoritmos: tigerprints.clemson.edu/cgi/… . Particularmente agradável e simples é aquele no ponto 5.1.1
isaacg
O i5-4300U possui 2 núcleos (4 threads), não 4 núcleos.
Anders Kaseorg

Respostas:

3

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

Para correr

cargo build --release
target/release/agcd 56007668599 30611458895 3
target/release/agcd --random 13 5
Anders Kaseorg
fonte
O resultado oficial é l = 13, com um tempo médio de 41,53s. l = 14 levou um pouco mais de 2m em média.
Isaacg
2

Mathematica, L = 5

aqui está uma solução rápida 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]]
) &

Entrada
[x1, x2, L]

para que alguém possa testar isso, aqui está o gerador de chaves

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

escolha o valor L executado, o código e você obterá três números.
coloque os dois últimos junto com L como entrada e você deverá obter o primeiro

J42161217
fonte
Eu verifiquei esta solução como obtendo uma pontuação de l = 5. Vou cronometrar se o tempo for necessário como desempate.
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[]]]]]) &

entrada [x1, x2, L]

para que alguém possa testar isso, aqui está o gerador de chaves

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

escolha o valor L, execute o código e você obterá três números.
coloque os dois últimos junto com L como entrada e você deverá obter o primeiro

J42161217
fonte
Quando testei isso l = 12, deu um resultado incorreto. Especificamente, o divisor resultante era um número de 146 bits, que é muito grande. Acho que você precisará apenas de uma pequena alteração para corrigir isso.
Isaacg
Eu adicionei o caso com falha como o exemplo final acima. Seu solucionador deu uma resposta igual a 3 * p, que era um pouco grande demais. Olhando para o seu código, parece que você apenas verifica se sua saída possui pelo menos l^2bits, mas também precisa verificar se ela possui no máximo l^2bits.
Isaacg
No mesmo caso de teste que falhou antes, ainda não está funcionando. Não estou super familiarizado com o Mathematica, mas parecia que ele saiu sem saída. Eu acho que o problema é que ele está encontrando um fator muito grande e, em vez de encontrar um divisor do fator no intervalo desejado, está apenas passando por ele.
Isaacg
A pontuação oficial é L = 12, com um tempo médio de 52,7s. Bem feito!
Isaacg
1

Python, L = 15, tempo médio 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))

Esse código é baseado no fato de que o produto de x1 - r para todos os valores possíveis de r deve ser um múltiplo de p e o produto de x2 - r também. Na maioria das vezes, é gasto o MDC dos dois produtos, cada um com cerca de 60.000.000 de bits. O gcd, que possui apenas cerca de 250.000 bits, é comparado a cada um dos valores xr mais uma vez, para extrair candidatos p '. Se eles são um pouco grandes demais, a divisão de teste é usada para reduzi-los ao intervalo apropriado.

Isso se baseia na biblioteca gmpy2 para Python, que é um invólucro fino para a biblioteca GNU MP, que em particular possui uma rotina gcd realmente boa.

Este código é de thread único.

isaacg
fonte