В последнее время я очень увлекся Project Euler, и сейчас я пытаюсь сделать это! Я начал его анализ и уже существенно уменьшил проблему. Вот моя работа:
A = pqr and
1/A = 1/p + 1/q + 1/r so pqr/A = pq + pr + qr
And because of the first equation:
pq+pr+qr = 1
Since exactly two of p, q and r have to be negative, we can simplify the equation down to finding:
abc for which ab = ac+bc+1
Solving for c we get:
ab-1 = (a+b)c
c = (ab-1)/(a+b)
This means we need to find a and b for which:
ab = 1 (mod a+b)
And then our A value with those a and b is:
A = abc = ab(ab-1)/(a+b)
Простите, если это много математики! Но теперь все, что нам нужно, - это одно условие и два уравнения. Теперь, поскольку мне нужно найти 150 000-е наименьшее целое число, записанное как ab (ab-1) / (a + b) с ab = 1 (mod a + b), в идеале я хочу найти (a, b), для которого A как можно меньше.
Для простоты я предположил, что a <b, и я также заметил, что gcd (a, b) = 1.
Моя первая реализация проста и даже достаточно быстро находит 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 решения. Этот код сложен и, хотя и быстрее, он недостаточно быстр ...
Так что у меня совсем нет идей! У кого-нибудь есть идеи?
Эта статья о китайском остатке, быстрой реализации, может вам помочь: www.codeproject.com/KB/recipes/CRP.aspx
Это еще ссылки на инструменты и библиотеки:
Инструменты:
Максима http://maxima.sourceforge.net/ Maxima - это система для управления символьными и числовыми выражениями, включая дифференцирование, интегрирование, ряды Тейлора, преобразования Лапласа, обыкновенные дифференциальные уравнения, системы линейных уравнений, полиномы и множества, списки, векторы, матрицы и тензоры. Maxima дает высокоточные числовые результаты за счет использования точных дробей, целых чисел произвольной точности и чисел с плавающей запятой переменной точности. Maxima может отображать функции и данные в двух и трех измерениях.
Математический http://mathomatic.org/math/ Mathomatic - это бесплатное портативное универсальное программное обеспечение CAS (система компьютерной алгебры) и калькулятора, которое может символически решать, упрощать, комбинировать и сравнивать уравнения, выполнять арифметические операции с комплексными числами и полиномами и т. д. Оно выполняет некоторые вычисления и его очень легко использовать. использовать.
Scilab www.scilab.org/download/index_download.php Scilab - это система численных вычислений, похожая на Matlab или Simulink. Scilab включает сотни математических функций, а программы на различных языках (например, C или Fortran) можно добавлять в интерактивном режиме.
математическая студия mathstudio.sourceforge.net Интерактивный редактор уравнений и пошаговое решение.
Библиотека:
Библиотека Armadillo C++ 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++, поэтому компиляция не требуется.
анимат http://www.yonsen.bz/animath/animath.html Animath - это библиотека методов конечных элементов, полностью реализованная на C++. Он подходит для моделирования взаимодействия жидкости и структуры и математически основан на тетраэдрических элементах более высокого порядка.
Спасибо за интересные библиотеки, но мой вопрос действительно касался конкретной проблемы ... Алгоритм и математика, а не инструменты.
Тем не менее, список ссылок, кажется, тщательно составлен и может заинтересовать некоторых. Поэтому я не понимаю, почему это должно быть -1, поэтому я буду голосовать за. В наши дни люди кажутся слишком беззаботными.
ОК. Вот еще немного игры с моим решением китайской теоремы об остатках. Оказывается, что a + b не может быть произведением любого простого числа p, кроме р = 1 (модуль 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 and b = v+1
ab = 1 (mod a+b)
uv+u+v = 0 (mod u+v+2)
Однако я не вижу в этом особых улучшений ...
Просто из любопытства: ваша программа слишком медленная, чтобы вообще найти решение, или вы просто пытаетесь нарушить правило 1 минуты? (Я еще даже не смотрел на это: P)
Я еще не решил это ... Я могу найти достаточно решений, но они не всегда самые маленькие 150 000.
Взгляните на свой предыдущий вопрос. Я опубликовал гораздо более быструю версию table (), она работает примерно в 3000 раз быстрее, чем та, которую вы использовали в приведенном выше коде.
Как и во многих других проблемах Project Euler, уловка состоит в том, чтобы найти метод, который сводит решение грубой силы к чему-то более прямому:
A = pqr and
1/A = 1/p + 1/q + 1/r
Так,
pq + qr + rp = 1 or -r = (pq - 1)/(p + q)
Без ограничения общности, 0 <p <-q <-r
Существует k, 1 <= k <= 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-е александрийское целое число.
К сожалению, pqr не увеличивается монотонно с увеличением p. Это означает, что нам, вероятно, следует сгенерировать гораздо больше, чем 150000 чисел (сколько?), Затем отсортировать их и, наконец, взять 150000-е из отсортированных чисел. Есть ли способ сгенерировать pqr в монотонно возрастающем порядке?
«К сожалению, pqr не увеличивается монотонно с увеличением p» - не могли бы вы объяснить?
Например: 1806: x = 6, y = 7, z = 43. Но 1428: x = 7, y = 12, z = 17. x1 <x2, но x1y1z1> x2y2z2. Итак, если мы сгенерируем первые 150000 Alexandr. чисел, увеличивая x, может существовать больший x, такой, что xyz меньше 150000-го числа, которое мы вычислили, поэтому действительное 150000-е число - это что-то другое.
Дело в том, что я разложил проблему на части, поэтому вам нужно только перебрать p
Дело в том, что никто не знает, как долго перебирать p - очевидно, что простая генерация 150000 чисел не будет содержать 150000-е число на основе pqr (не на основе p). И проблема требует, чтобы точно производил 150000-е число на основе значения pqr.
В N числах при увеличении p это может / не может иметь 150000-е число pqr -s, независимо от того, насколько велико N, и даже 150000-е число pqr находится среди этих N p-чисел, как мы теперь можем определить, какое именно это число? является? Даже если p-числа отсортированы, 150000-е pqr num может не быть одним из N p-чисел.
умное решение
@PythonPower: пожалуйста, посмотрите мой ответ ниже, в котором показано, как параметризовать A через p и делители p ^ 2 + 1