Поиск целых чисел с определенным свойством - задача Эйлера 221 - PullRequest
9 голосов
/ 22 декабря 2008

Недавно я сильно увлекся Project Euler и пытаюсь сделать этот один следующий! Я начал некоторый анализ этого и уже существенно уменьшил проблему. Вот моя работа:

A = pqr и

1 / A = 1 / p + 1 / q + 1 / r, поэтому pqr / A = pq + pr + qr

А из-за первого уравнения:

pq + pr + qr = 1

Поскольку ровно два из p, q и r имеют чтобы быть отрицательным, мы можем упростить уравнение до нахождения:

abc, для которого ab = ac + bc + 1

Решая для c, получаем:

ab-1 = (a + b) c

c = (ab-1) / (a ​​+ b)


Это означает, что нам нужно найти a и b для в том числе:

ab = 1 (mod a + b)

А потом наше значение с этими а и б:

A = abc = ab (ab-1) / (a ​​+ b)

Извините, если много математики! Но теперь все, с чем мы имеем дело, это одно условие и два уравнения. Теперь, поскольку мне нужно найти 150000-е наименьшее целое число, записанное как ab (ab-1) / (a ​​+ b) с ab = 1 (mod a + b), в идеале я хочу найти (a, b), для которого A как можно меньше.

Для простоты я предположил, что Моя первая реализация проста и даже находит достаточно быстро 150 000 решений. Однако поиск 150 000 самых маленьких решений занимает слишком много времени. Вот код в любом случае:

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

Моей следующей мыслью было использовать деревья Штерна-Броко, но это слишком медленно, чтобы найти решение. Мой последний алгоритм состоял в том, чтобы использовать китайскую теорему об остатках, чтобы проверить, дают ли различные значения a + b решения. Этот код сложен и хотя и быстрее, но он недостаточно быстр ...

Так что я абсолютно вне идей! У кого-нибудь есть идеи?

Ответы [ 3 ]

4 голосов
/ 22 декабря 2008

Эта статья о китайском языке, быстрая реализация, может помочь вам: www.codeproject.com/KB/recipes/CRP.aspx

Это дополнительные ссылки для инструментов и библиотек:

Инструменты:

Maxima http://maxima.sourceforge.net/ Maxima - это система для манипулирования символическими и числовыми выражениями, включая дифференцирование, интегрирование, ряды Тейлора, преобразования Лапласа, обыкновенные дифференциальные уравнения, системы линейных уравнений, многочлены и множества, списки, векторы, матрицы и тензоры. Maxima дает числовые результаты с высокой точностью, используя точные дроби, произвольные целые числа точности и числа с плавающей запятой переменной точности. Maxima может отображать функции и данные в двух и трех измерениях.

Mathomatic http://mathomatic.org/math/ Mathomatic - это бесплатное, портативное, универсальное программное обеспечение CAS (Computer Algebra System) и калькулятор, которое может символически решать, упрощать, объединять и сравнивать уравнения, выполнять комплексные числа и полиномиальную арифметику и т. Д. Оно выполняет некоторые вычисления и очень легко использовать.

Scilab www.scilab.org/download/index_download.php Scilab - это система численных расчетов, аналогичная Matlab или Simulink. Scilab включает в себя сотни математических функций, и программы из разных языков (таких как C или Fortran) могут быть добавлены в интерактивном режиме.

mathstudio mathstudio.sourceforge.net Интерактивный редактор формул и пошаговый решатель.

Библиотека:

Armadillo C ++ Library http://arma.sourceforge.net/ Библиотека Armadillo C ++ призвана обеспечить эффективную основу для операций линейной алгебры (матричная и векторная математика), имея простой и удобный интерфейс.

Блиц ++ http://www.oonumerics.org/blitz/ Blitz ++ - библиотека классов C ++ для научных вычислений

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

libapmath http://freshmeat.net/projects/libapmath Добро пожаловать на домашнюю страницу проекта APMath. Целью этого проекта является реализация библиотеки C ++ произвольной точности, которая является наиболее удобной в использовании, это означает, что все операции реализованы как операторные перегрузки, именование в основном такое же, как у.

libmat http://freshmeat.net/projects/libmat MAT - это математическая библиотека классов C ++. Используйте эту библиотеку для различных матричных операций, поиска корней многочленов, решения уравнений и т. Д. Библиотека содержит только заголовочные файлы C ++, поэтому компиляция не требуется.

animath http://www.yonsen.bz/animath/animath.html Animath - это библиотека методов конечных элементов, полностью реализованная на C ++. Он подходит для моделирования взаимодействия жидкости со структурой и математически основан на тетраэдрических элементах высшего порядка.

3 голосов
/ 01 января 2009

Как и во многих задачах Project Euler, уловка заключается в том, чтобы найти метод, который превращает решение грубой силы в нечто более прямолинейное:

A = pqr and 
1/A = 1/p + 1/q + 1/r
1003 * * Итак,
pq + qr + rp = 1  or  -r = (pq - 1)/(p + q)

Без потери общности, 0

Существует k, 1 <= k <= p </p>

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

Но r является целым числом, поэтому k делит p ^ 2 + 1

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

Таким образом, чтобы вычислить A, нам нужно перебрать p, и где k может принимать только те значения, которые являются делителями p в квадрате плюс 1.

Добавляя каждое решение к набору, мы можем остановиться, когда найдем требуемое 150000-е целое число Александрии.

0 голосов
/ 23 декабря 2008

OK. Вот еще игра с моим решением для китайской теоремы об остатках. Оказывается, что a + b не может быть произведением любого простого числа p, если только p = 1 (mod 4) . Это позволяет быстрее вычислять, так как нам нужно только проверить a + b, которые кратны простым, таким как 2, 5, 13, 17, 29, 37 ...

Вот пример возможных значений a + b:

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

А вот полная программа с использованием китайской теоремы об остатках:

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

Это лучше, но я надеюсь улучшить его дальше (например, a + b = 2 ^ n, кажется, никогда не будут решениями).

Я также начал рассматривать основные замены, такие как:

a = u + 1 и b = v + 1

ab = 1 (mod a + b)

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

Однако я не вижу большого улучшения в этом ...

...