Нахождение целых чисел с определенным свойством — Задача Эйлера проекта 221
-
23-08-2019 - |
Вопрос
В последнее время я очень увлекся проектом Эйлера и пытаюсь сделать этот один следующий!Я начал некоторый анализ и уже существенно уменьшил проблему.Вот моя работа:
A = pqr и
1/a = 1/p + 1/q + 1/r, так что pqr/a = pq + pr + qr
И из-за первого уравнения:
пк+пр+qr = 1
Поскольку ровно два из P, Q и R должны быть отрицательными, мы можем упростить уравнение до поиска:
abc, для которого ab = ac+bc+1
Решая для c, получаем:
ab-1 = (a+b)c
с = (ab-1)/(a+b)
Это означает, что нам нужно найти A и B, для которых:
аб = 1 (мод а+б)
И тогда наша ценность с этими a и b есть:
А = 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 решения.Этот код сложен и хотя и быстрее, но недостаточно быстр...
Так что у меня совершенно нет идей!У кого-нибудь есть идеи?
Решение
Как и во многих задачах проекта Эйлера, задача состоит в том, чтобы найти метод, который сводит решение методом грубой силы к чему-то более простому:
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.
Добавляя каждое решение в набор, мы можем остановиться, когда найдем искомое 150-тысячное александрийское целое число.
Другие советы
Эта статья о китайском остатке и быстрой реализации может помочь вам: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 MathStudio.sourceForge.net Интерактивный редактор уравнений и пошаговый решатель.
Библиотека:
Библиотека 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
libapmathhttp://freshmeat.net/projects/libapmathДобро пожаловать на домашнюю страницу проекта APMath.Целью данного проекта является реализация C++-библиотеки произвольной точности, наиболее удобной в использовании, то есть все операции реализованы как перегрузки операторов, именование в основном такое же, как и в .
библиотекаhttp://freshmeat.net/projects/libmatMAT — это библиотека классов математических шаблонов C++.Используйте эту библиотеку для различных матричных операций, поиска корней многочленов, решения уравнений и т. д.Библиотека содержит только заголовочные файлы C++, поэтому компиляция не требуется.
аниматhttp://www.yonsen.bz/animath/animath.htmlAnimath — это библиотека методов конечных элементов, полностью реализованная на C++.Он подходит для моделирования взаимодействия жидкости и структуры и математически основан на тетраэдрических элементах более высокого порядка.
ХОРОШО.Вот еще несколько экспериментов с моим решением китайской теоремы об остатках.Оказывается, 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, похоже, никогда не будет решением).
Я также начал рассматривать основные замены, такие как:
а = и+1 и б = v+1
аб = 1 (мод а+б)
uv+u+v = 0 (по модулю u+v+2)
Но особого улучшения я от этого не вижу...