Diviseur commun approximatif le plus rapide

13

Aperçu

Dans ce défi, vous recevrez deux nombres qui sont tous deux un petit décalage plus grand qu'un multiple d'un nombre de taille moyenne. Vous devez sortir un nombre de taille moyenne qui est presque un diviseur des deux nombres, à l'exception d'un petit décalage.

La taille des effectifs concernés sera paramétré par un paramètre de difficulté, l. Votre objectif est de résoudre le problème le plus largement possible len moins d'une minute.


Installer

Dans un problème donné, il y aura un nombre secret p, qui sera un nombre aléatoire l^2( l*l). Il y aura deux multiplicateurs, q1, q2qui seront des l^3nombres de bits aléatoires , et il y aura deux décalages r1, r2, qui seront des lnombres de bits aléatoires .

L'entrée dans votre programme sera x1, x2définie comme suit:

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

Voici un programme pour générer des cas de test, en 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)

La première ligne de sortie est une solution possible, tandis que les deuxième et troisième lignes sont l'entrée que votre programme recevra.


Votre programme

Étant donné x1, x2et l, vous devez trouver un l^2nombre de bits p'tel que x1 % p'et x2 % p'sont tous deux ldes nombres de bits. pfonctionnera toujours, bien qu'il puisse y avoir d'autres possibilités. Voici une fonction pour vérifier une 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)

Exemple

Supposons que l3. Le programme générateur choisit un nombre de 9 bits pour p, qui dans ce cas est 442. Le générateur sélectionne deux 3nombres de bits pour r1, r2, qui sont 4, 7. Le générateur sélectionne deux 27nombres de bits pour q1, q2, qui sont 117964803, 101808039. En raison de ces choix, le x1, x2sont 52140442930, 44999153245.

Votre programme serait donné 52140442930, 44999153245en entrée et doit sortir un nombre de 9 bits (dans la plage [256, 511]) tel que 52140442930et 44999153245modulo ce nombre donne des nombres de 3 bits (dans la plage [4, 7]). 442est la seule valeur de ce type dans ce cas, donc votre programme devrait sortir 442.


Plus d'exemples

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

Notation

Comme mentionné ci-dessus, le score de votre programme est le plus élevé lque le programme termine en moins d'une minute. Plus précisément, votre programme sera exécuté sur 5 instances aléatoires avec cela l, et il doit produire une réponse correcte sur les 5, avec un temps moyen inférieur à 1 minute. Le score d'un programme sera le plus élevé lauquel il réussira. Tiebreaker sera le temps moyen à ce sujet l.

Pour vous donner une idée des scores à viser, j'ai écrit un solveur de force brute très simple. Il a obtenu un score de 5. J'ai écrit un solveur beaucoup plus sophistiqué. Il a obtenu un score de 12 ou 13, selon la chance.


Détails

Pour une parfaite comparabilité entre les réponses, je chronomètre les soumissions sur mon ordinateur portable pour donner des scores canoniques. Je vais également exécuter les mêmes instances choisies au hasard sur toutes les soumissions, pour alléger quelque peu la chance. Mon ordinateur portable possède 4 processeurs, un processeur i5-4300U à 1,9 GHz, 7,5 G de RAM.

N'hésitez pas à publier un score provisoire basé sur votre propre timing, indiquez simplement s'il est provisoire ou canonique.


Que le programme le plus rapide gagne!

isaacg
la source
L'approximation doit-elle être la plus proche?
Yytsi
@TuukkaX Tout l^2nombre de bits qui est à l-bits d'être un facteur des deux nombres fonctionne. Cependant, il n'y en a généralement qu'un.
isaacg
Voici une dissertation avec quelques idées d'algorithmes: tigerprints.clemson.edu/cgi/… . Particulièrement agréable et simple est celui de la section 5.1.1
isaacg
Le i5-4300U a 2 cœurs (4 threads), pas 4 cœurs.
Anders Kaseorg

Réponses:

3

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

Courir

cargo build --release
target/release/agcd 56007668599 30611458895 3
target/release/agcd --random 13 5
Anders Kaseorg
la source
Le résultat officiel est l = 13, avec un temps moyen de 41,53s. l = 14 a pris un peu plus de 2m en moyenne.
isaacg
2

Mathematica, L = 5

voici une solution rapide à 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]]
) &

Entrée
[x1, x2, L]

afin que tout le monde puisse le tester, voici le générateur de clés

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

choisissez la valeur L exécutée, le code et vous obtiendrez trois nombres.
placez les deux derniers avec L comme entrée, et vous devriez obtenir le premier

J42161217
la source
J'ai vérifié que cette solution obtient un score de l = 5. Je vais le chronométrer si le timing est nécessaire comme bris d'égalité.
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[]]]]]) &

contribution [x1, x2, L]

afin que tout le monde puisse le tester, voici le générateur de clés

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

choisissez la valeur L, exécutez le code et vous obtiendrez trois nombres.
placez les deux derniers avec L comme entrée, et vous devriez obtenir le premier

J42161217
la source
Lorsque j'ai testé cela avec l = 12, cela a donné un résultat incorrect. Plus précisément, le diviseur résultant était un nombre de 146 bits, ce qui est trop grand. Je pense que vous n'aurez besoin que d'un petit changement pour résoudre ce problème.
isaacg
J'ai ajouté le cas défaillant comme dernier exemple ci-dessus. Votre solveur a donné une réponse égale à 3 * p, ce qui était un peu trop grand. En regardant votre code, il semble que vous ne vérifiez que votre sortie comporte au moins des l^2bits, mais vous devez également vérifier qu'elle contient au plus de l^2bits.
isaacg
Sur le même cas de test qu'il a échoué auparavant, il ne fonctionne toujours pas. Je ne suis pas très familier avec Mathematica, mais il semble qu'il soit sorti sans sortie. Je pense que le problème est qu'il trouve un facteur trop grand, puis au lieu de trouver un diviseur du facteur dans la plage souhaitée, il ne fait que le dépasser.
isaacg
Continuons cette discussion dans le chat .
isaacg
Le score officiel est L = 12, avec un temps moyen de 52,7 s. Bien joué!
isaacg
1

Python, L = 15, temps moyen 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))

Ce code est basé sur le fait que le produit de x1 - r pour toutes les valeurs possibles de r doit être un multiple de p, et le produit de x2 - r doit également l'être. La plupart du temps est consacré à prendre le gcd des deux produits, chacun ayant environ 60 000 000 bits. Le pgcd, qui n'a que 250 000 bits environ, est ensuite à nouveau comparé à chacune des valeurs xr, pour extraire p 'candidats. S'ils sont un peu trop grands, la division d'essai est utilisée pour les réduire à la plage appropriée.

Ceci est basé sur la bibliothèque gmpy2 pour Python, qui est un wrapper mince pour la bibliothèque GNU MP, qui a en particulier une très bonne routine gcd.

Ce code est monothread.

isaacg
la source