Finding inteiros com uma certa propriedade - Projeto Euler problema 221

StackOverflow https://stackoverflow.com/questions/386862

Pergunta

Eu me tornei muito viciado em Projeto Euler recentemente e estou tentando fazer este um ao lado! Eu comecei algumas análises sobre ele e ter reduzido o problema para baixo substancialmente já. Aqui é o meu trabalho:

A = PQR e

1 / A = 1 / p + 1 / Q + 1 / R de modo PQR / A = pq + pr + qr

E por causa da primeira equação:

pq + PR + qr = 1

Desde exatamente dois de p, q e r têm a ser negativo, podemos simplificar o equação para baixo a encontrar:

ABC para que ab = AC + BC + 1

Resolvendo para c obtemos:

ab-1 = (a + b) c

c = (AB-1) / (a ??+ b)

Isto significa que precisamos de encontrar a e b para que:

AB = 1 (mod a + b)

E então o nosso valor A com aqueles um e b é:

A = ABC = ab (ab-1) / (a ??+ b)

Desculpe se isso é um monte de matemática! Mas agora todos nós temos de lidar com uma condição e duas equações. Agora desde que eu preciso para encontrar o 150000 menor inteiro escrito como ab (ab-1) / (a ??+ b) com ab = 1 (mod a + b), idealmente eu quero procurar (a, b) para que A é tão pequeno quanto possível.

Para facilitar eu assumi a

A minha primeira implementação é para a frente e ainda encontra 150.000 soluções rápido o suficiente. No entanto, leva muito tempo para encontrar as 150.000 menor soluções. Aqui está o código de qualquer maneira: 

n = 150000
seen = set()

a = 3
while len(seen) < n:
    for b in range(2, a):
        if (a*b)%(a+b) != 1: continue

        seen.add(a*b*(a*b-1)//(a+b))
        print(len(seen), (a, b), a*b*(a*b-1)//(a+b))

    a += 1

Meu próximo pensamento era usar árvores Stern-Brocot mas que é muito lento para encontrar soluções. Meu algoritmo final era para usar o teorema restante chinês para verificar se diferentes valores de a + b soluções de rendimento. Esse código é complicada e, embora mais rápido, não é rápido o suficiente ...

Então, eu estou absolutamente fora de ideias! Alguém tem alguma idéia?

Foi útil?

Solução

Tal como acontece com muitos dos problemas de projeto Euler, o truque é encontrar uma técnica que reduz a solução de força bruta em algo mais para a frente: 

A = pqr and 
1/A = 1/p + 1/q + 1/r

Assim, 

pq + qr + rp = 1  or  -r = (pq - 1)/(p + q)

Sem perda de generalidade, 0

Existe k, 1 <= k <= p 

-q = p + k
-r = (-p(p + k) – 1) / (p + -p – k)  = (p^2 + 1)/k + p

Mas r é um número inteiro, de modo que divide k p ^ 2 + 1 

pqr = p(p + q)((p^2 + 1)/k + p)

Assim, para calcular Um precisamos iterar p, e onde k só pode assumir valores que são divisores de p quadrado mais 1.

A adição de cada solução a um conjunto, podemos parar quando encontramos o número inteiro de Alexandria exigido 150000.

Outras dicas

Este artigo sobre restante chinês, implementação rápida, pode ajudá-lo: www.codeproject.com/KB/recipes/CRP.aspx

Esta é mais links para ferramentas e bibliotecas:

Ferramentas:

Maxima http://maxima.sourceforge.net/ Maxima é um sistema para a manipulação de expressões simbólicas e numéricas, incluindo diferenciação, integração, série de Taylor, transformadas de Laplace, equações diferenciais ordinárias, sistemas de equações lineares, polinômios, e conjuntos, listas, vetores, matrizes e tensores. Máximos proporciona alta precisão resultados numéricos usando fracções exactas, números inteiros de precisão arbitrária, e números de ponto flutuante de precisão variável. Maxima pode desenhar funções e dados em duas ou três dimensões.

mathomatic http://mathomatic.org/math/ Mathomatic é um, portátil, de uso geral livre CAS (Computer Algebra System) e software da calculadora que podem simbolicamente resolver, simplificar, combinar e comparar equações, executar número complexo e aritmética polinomial, etc. Ele faz alguns cálculos e é muito fácil de usar.

Scilab www.scilab.org/download/index_download.php Scilab é um similar sistema de computação numérica para Matlab ou Simulink. Scilab inclui centenas de funções matemáticas e programas de várias linguagens (como C ou Fortran) podem ser adicionadas de forma interactiva.

mathstudio mathstudio.sourceforge.net Um calculador interactivo editor de equações e passo-a-passo.

Biblioteca:

Biblioteca ++ Armadillo C http://arma.sourceforge.net/ Os objectivos da biblioteca Armadillo C ++ para fornecer uma base eficaz para operações de álgebra linear (matriciais e vectoriais matemática), tendo um simples e fácil de interface de uso.

Blitz ++ http://www.oonumerics.org/blitz/ Blitz ++ é uma biblioteca de classes C ++ para computação científica

BigInteger C # http://msdn.microsoft.com/pt-br/magazine/cc163441. aspx

libapmath http://freshmeat.net/projects/libapmath Bem-vindo à página inicial do APMath-projeto. Objectivo deste projecto é a implementação de uma precisão arbitrária C ++ -. Biblioteca, que é o mais conveniente em uso, isso significa que todas as operações são implementadas como operador-overloadings, nomenclatura é basicamente o mesmo que o da

libmat http://freshmeat.net/projects/libmat MAT é uma biblioteca de classes template C ++ matemática. Use esta biblioteca para várias operações de matriz, encontrar raízes de polinômios, resolver equações, etc. A biblioteca contém apenas os arquivos C ++ cabeçalho, de modo nenhum compilação é necessário.

animath http://www.yonsen.bz/animath/animath.html Animath é uma biblioteca Método dos Elementos Finitos inteiramente implementado em C ++. É adequado para a simulação interacção estrutura-fluido, e baseia-se matematicamente sobre higher-order elementos tetraédricos.

OK. Aqui está um pouco mais de jogo com a minha solução Restante Teorema chinês. Acontece que a + b não pode ser o produto de qualquer procedimento de injeção, p, a menos que p = 1 (mod 4) . Isso permite que a computação mais rápida, já que só tem que verificar a + b que são múltiplos de números primos, tais como 2, 5, 13, 17, 29, 37 ...

Então, aqui está um exemplo de possíveis + b valores de: 

[5, 8, 10, 13, 16, 17, 20, 25, 26, 29, 32, 34, 37, 40, 41, 50, 52, 53, 58, 61, 64, 65, 68, 73, 74, 80, 82, 85, 89, 97, 100]

E aqui está o programa completo usando o teorema chinês do resto: 

cachef = {}
def factors(n):
    if n in cachef: cachef[n]

    i = 2
    while i*i <= n:
        if n%i == 0:
            r = set([i])|factors(n//i)
            cachef[n] = r
            return r

        i += 1

    r = set([n])
    cachef[n] = r
    return r

cachet = {}
def table(n):
    if n == 2: return 1
    if n%4 != 1: return

    if n in cachet: return cachet[n]

    a1 = n-1
    for a in range(1, n//2+1):
        if (a*a)%n == a1:
            cachet[n] = a
            return a

cacheg = {}
def extended(a, b):
    if a%b == 0:
        return (0, 1)
    else:
        if (a, b) in cacheg: return cacheg[(a, b)]

        x, y = extended(b, a%b)
        x, y = y, x-y*(a//b)

        cacheg[(a, b)] = (x, y)
        return (x, y)

def go(n):
    f = [a for a in factors(n)]
    m = [table(a) for a in f]

    N = 1
    for a in f: N *= a

    x = 0
    for i in range(len(f)):
        if not m[i]: return 0

        s, t = extended(f[i], N//f[i])
        x += t*m[i]*N//f[i]

    x %= N
    a = x
    while a < n:
        b = n-a
        if (a*b-1)%(a+b) == 0: return a*b*(a*b-1)//(a+b)
        a += N




li = [5, 8, 10, 13, 16, 17, 20, 25, 26, 29, 32, 34, 37, 40, 41, 50, 52, 53, 58, 61, 64, 65, 68, 73, 74, 80, 82, 85, 89, 97, 100]

r = set([6])
find = 6

for a in li:
    g = go(a)
    if g:
        r.add(g)
        #print(g)
    else:
        pass#print(a)


r = list(r)
r.sort()

print(r)
print(len(r), 'with', len(li), 'iterations')

Esta é melhor, mas espero melhorar ainda mais (por exemplo, a + b = 2 ^ n parecem não ser soluções).

Eu também já começaram a considerar substituições básicas, tais como:

a = L + 1 e b = v + 1

AB = 1 (mod a + b)

uv + u + v = 0 (mod u + v + 2)

No entanto, eu não posso ver muita melhoria com isso ...

Licenciado em: CC-BY-SA com atribuição
Não afiliado a StackOverflow
scroll top