Как перечислить x ^ 2 + y ^ 2 = z ^ 2 - 1 (с дополнительными ограничениями)

Ответ 1

Здесь метод, который перечисляет тройки, а не исчерпывающе проверяет их, используя теорию чисел, как описано здесь: https://mathoverflow.net/questions/29644/enumerating-ways-to-decompose-an-integer-into-the- сумма-из-двух квадратов

Поскольку математика заняла у меня некоторое время, чтобы осмыслить, и некоторое время, чтобы реализовать (собрав некоторый код, который заимствован над ней), и так как я не чувствую большой авторитет в этом вопросе, я предоставлю читателю возможность исследовать. Это основано на выражении чисел как гауссовых целочисленных сопряженных. (a + bi)*(a - bi) = a^2 + b^2. Сначала мы делим число z^2 - 1 на простые числа, разлагаем простые числа на гауссовы сопряженные и находим различные выражения, которые мы расширяем и упрощаем, чтобы получить a + bi, которое затем можно повысить, a^2 + b^2.

Приятное чтение о функции суммы квадратов обнаруживает, что мы можем исключить любого кандидата z^2 - 1 который содержит простое число формы 4k + 3 с нечетной степенью. Используя только эту проверку, я смог уменьшить цикл Чернослива на 10 ^ 5 с 214 секунд до 19 секунд (на repl.it), используя приведенный ниже код простого фактора Розетты.

Реализация здесь просто демонстрация. У него нет обработки или оптимизации для ограничения x и y. Скорее, он просто перечисляет по ходу дела. Играть с этим здесь.

Код Python:

# https://math.stackexchange.com/questions/5877/efficiently-finding-two-squares-which-sum-to-a-prime
def mods(a, n):
    if n <= 0:
        return "negative modulus"
    a = a % n
    if (2 * a > n):
        a -= n
    return a

def powmods(a, r, n):
    out = 1
    while r > 0:
        if (r % 2) == 1:
            r -= 1
            out = mods(out * a, n)
        r /= 2
        a = mods(a * a, n)
    return out

def quos(a, n):
    if n <= 0:
        return "negative modulus"
    return (a - mods(a, n))/n

def grem(w, z):
    # remainder in Gaussian integers when dividing w by z
    (w0, w1) = w
    (z0, z1) = z
    n = z0 * z0 + z1 * z1
    if n == 0:
        return "division by zero"
    u0 = quos(w0 * z0 + w1 * z1, n)
    u1 = quos(w1 * z0 - w0 * z1, n)
    return(w0 - z0 * u0 + z1 * u1,
           w1 - z0 * u1 - z1 * u0)

def ggcd(w, z):
    while z != (0,0):
        w, z = z, grem(w, z)
    return w

def root4(p):
    # 4th root of 1 modulo p
    if p <= 1:
        return "too small"
    if (p % 4) != 1:
        return "not congruent to 1"
    k = p/4
    j = 2
    while True:
        a = powmods(j, k, p)
        b = mods(a * a, p)
        if b == -1:
            return a
        if b != 1:
            return "not prime"
        j += 1

def sq2(p):
    if p % 4 != 1:
      return "not congruent to 1 modulo 4"
    a = root4(p)
    return ggcd((p,0),(a,1))

# https://rosettacode.org/wiki/Prime_decomposition#Python:_Using_floating_point
from math import floor, sqrt

def fac(n):
    step = lambda x: 1 + (x<<2) - ((x>>1)<<1)
    maxq = long(floor(sqrt(n)))
    d = 1
    q = n % 2 == 0 and 2 or 3 
    while q <= maxq and n % q != 0:
        q = step(d)
        d += 1
    return q <= maxq and [q] + fac(n//q) or [n]

# My code...
# An answer for  https://stackoverflow.com/questions/54110614/

from collections import Counter
from itertools import product
from sympy import I, expand, Add

def valid(ps):
  for (p, e) in ps.items():
    if (p % 4 == 3) and (e & 1):
      return False
  return True

def get_sq2(p, e):
  if p == 2:
    if e & 1:
      return [2**(e / 2), 2**(e / 2)]
    else:
      return [2**(e / 2), 0]
  elif p % 4 == 3:
    return [p, 0]
  else:
    a,b = sq2(p)
    return [abs(a), abs(b)]

def get_terms(cs, e):
  if e == 1:
    return [Add(cs[0], cs[1] * I)]
  res = [Add(cs[0], cs[1] * I)**e]
  for t in xrange(1, e / 2 + 1):
    res.append(
      Add(cs[0] + cs[1]*I)**(e-t) * Add(cs[0] - cs[1]*I)**t)
  return res

def get_lists(ps):
  items = ps.items()
  lists = []
  for (p, e) in items:
    if p == 2:
      a,b = get_sq2(2, e)
      lists.append([Add(a, b*I)])
    elif p % 4 == 3:
      a,b = get_sq2(p, e)
      lists.append([Add(a, b*I)**(e / 2)])
    else:
      lists.append(get_terms(get_sq2(p, e), e))
  return lists


def f(n):
  for z in xrange(2, n / 2):
    zz = (z + 1) * (z - 1)
    ps = Counter(fac(zz))
    is_valid = valid(ps)
    if is_valid:
      print "valid (does not contain a prime of form\n4k + 3 with an odd power)"
      print "z: %s, primes: %s" % (z, dict(ps))
      lists = get_lists(ps)
      cartesian = product(*lists)
      for element in cartesian:
        print "prime square decomposition: %s" % list(element)
        p = 1
        for item in element:
          p *= item
        print "complex conjugates: %s" % p
        vals = p.expand(complex=True, evaluate=True).as_coefficients_dict().values()
        x, y = vals[0], vals[1] if len(vals) > 1 else 0
        print "x, y, z: %s, %s, %s" % (x, y, z)
        print "x^2 + y^2, z^2-1: %s, %s" % (x**2 + y**2, z**2 - 1)
      print ''

if __name__ == "__main__":
  print f(100)

Выход:

valid (does not contain a prime of form
4k + 3 with an odd power)
z: 3, primes: {2: 3}
prime square decomposition: [2 + 2*I]
complex conjugates: 2 + 2*I
x, y, z: 2, 2, 3
x^2 + y^2, z^2-1: 8, 8

valid (does not contain a prime of form
4k + 3 with an odd power)
z: 9, primes: {2: 4, 5: 1}
prime square decomposition: [4, 2 + I]
complex conjugates: 8 + 4*I
x, y, z: 8, 4, 9
x^2 + y^2, z^2-1: 80, 80

valid (does not contain a prime of form
4k + 3 with an odd power)
z: 17, primes: {2: 5, 3: 2}
prime square decomposition: [4 + 4*I, 3]
complex conjugates: 12 + 12*I
x, y, z: 12, 12, 17
x^2 + y^2, z^2-1: 288, 288

valid (does not contain a prime of form
4k + 3 with an odd power)
z: 19, primes: {2: 3, 3: 2, 5: 1}
prime square decomposition: [2 + 2*I, 3, 2 + I]
complex conjugates: (2 + I)*(6 + 6*I)
x, y, z: 6, 18, 19
x^2 + y^2, z^2-1: 360, 360

valid (does not contain a prime of form
4k + 3 with an odd power)
z: 33, primes: {17: 1, 2: 6}
prime square decomposition: [4 + I, 8]
complex conjugates: 32 + 8*I
x, y, z: 32, 8, 33
x^2 + y^2, z^2-1: 1088, 1088

valid (does not contain a prime of form
4k + 3 with an odd power)
z: 35, primes: {17: 1, 2: 3, 3: 2}
prime square decomposition: [4 + I, 2 + 2*I, 3]
complex conjugates: 3*(2 + 2*I)*(4 + I)
x, y, z: 18, 30, 35
x^2 + y^2, z^2-1: 1224, 1224

Ответ 2

Вот простое улучшение в Python (преобразование в более быстрый эквивалент в C-коде оставлено в качестве упражнения для читателя). Чтобы получить точные сроки вычислений, я удалил печать самих решений (после проверки их в предыдущем запуске).

• Используйте внешний цикл для одной свободной переменной (я выбрал z), ограниченный только ее отношением к N
• Используйте внутренний цикл (я выбрал y), ограниченный индексом внешнего цикла.
• Третья переменная вычисляется напрямую в соответствии с требованием 2.

Сроки результаты:

-------------------- 10 
 1 solutions found in 2.3365020751953125e-05  sec.
-------------------- 100 
 6 solutions found in 0.00040078163146972656  sec.
-------------------- 1000 
 55 solutions found in 0.030081748962402344  sec.
-------------------- 10000 
 543 solutions found in 2.2078349590301514  sec.
-------------------- 100000 
 5512 solutions found in 214.93411707878113  sec.

Это 3:35 для большого случая, плюс ваше время, чтобы сопоставить и/или распечатать результаты.

Если вам нужен более быстрый код (это все еще довольно грубая сила), изучите диофантовы уравнения и параметризацию для генерации (y, x) пар, учитывая целевое значение z^2 - 1.

import math
import time

def break3(N):
    """
    10 <= N <= 10^5
    return x, y, z triples such that:
        x <= y <= z
        x^2 + y^2 = z^2 - 1        
        x + y + z <= N
    """

    """
    Observations:
    z <= x + y
    z < N/2
    """

    count = 0
    z_limit = N // 2
    for z in range(3, z_limit):

        # Since y >= x, there a lower bound on y
        target = z*z - 1
        ymin = int(math.sqrt(target/2))
        for y in range(ymin, z):
            # Given y and z, compute x.
            # That a solution iff x is integer.
            x_target = target - y*y
            x = int(math.sqrt(x_target))
            if x*x == x_target and x+y+z <= N:
                # print("solution", x, y, z)
                count += 1

    return count


test = [10, 100, 1000, 10**4, 10**5]
border = "-"*20

for case in test: 
    print(border, case)
    start = time.time()
    print(break3(case), "solutions found in", time.time() - start, "sec.")

Ответ 3

Нет времени, чтобы правильно его протестировать, но, похоже, он дал те же результаты, что и ваш код (при 100 → 6 результатов и при 1000 → 55 результатов).

С N=1000 время 2ms против ваших 144ms также без списка

и N=10000 время 28ms

var N = 1000;
var c = 0;

for (int x = 2; x < N; x+=2)
{
    for (int y = x; y < (N - x); y+=2)
    {
        long z2 = x * x + y * y + 1;
        int z = (int) Math.Sqrt(z2);
        if (x + y + z > N)
            break;
        if (z * z == z2)
            c++;
    }
}

Console.WriteLine(c);

Ответ 4

Границы x и y являются важной частью проблемы. Я лично пошел с этим запросом Wolfram Alpha и проверил точные формы переменных.

Благодаря @Bleep-Bloop и комментариям, была найдена очень элегантная оптимизация границ: x < n и x <= y < n - x. Результаты одинаковы, а времена почти идентичны.

Кроме того, поскольку единственными возможными значениями x и y являются положительные четные целые числа, мы можем уменьшить количество итераций цикла вдвое.

Для дальнейшей оптимизации, так как мы вычисляем верхнюю границу x, мы строим список всех возможных значений для x и проводим параллельные вычисления. Это экономит огромное количество времени при более высоких значениях N но немного медленнее при меньших значениях из-за издержек распараллеливания.

Вот окончательный код:

Непараллельная версия, со значениями int :

List<string> res = new List<string>();
int n2 = n * n;

double maxX = 0.5 * (2.0 * n - Math.Sqrt(2) * Math.Sqrt(n2 + 1));

for (int x = 2; x < maxX; x += 2)
{
    int maxY = (int)Math.Floor((n2 - 2.0 * n * x - 1.0) / (2.0 * n - 2.0 * x));

    for (int y = x; y <= maxY; y += 2)
    {
        int z2 = x * x + y * y + 1;
        int z = (int)Math.Sqrt(z2);

        if (z * z == z2 && x + y + z <= n)
            res.Add(x + "," + y + "," + z);
    }
}

Параллельная версия с long значениями:

using System.Linq;

...

// Use ConcurrentBag for thread safety
ConcurrentBag<string> res = new ConcurrentBag<string>();
long n2 = n * n;

double maxX = 0.5 * (2.0 * n - Math.Sqrt(2) * Math.Sqrt(n2 + 1L));

// Build list to parallelize
int nbX = Convert.ToInt32(maxX);
List<int> xList = new List<int>();
for (int x = 2; x < maxX; x += 2)
    xList.Add(x);

Parallel.ForEach(xList, x =>
{
    int maxY = (int)Math.Floor((n2 - 2.0 * n * x - 1.0) / (2.0 * n - 2.0 * x));

    for (long y = x; y <= maxY; y += 2)
    {
        long z2 = x * x + y * y + 1L;
        long z = (long)Math.Sqrt(z2);

        if (z * z == z2 && x + y + z <= n)
            res.Add(x + "," + y + "," + z);
    }
});

Когда запускается индивидуально на процессоре i5-8400, я получаю следующие результаты:

N: 10; Решения: 1; Прошедшее время: 0,03 мс (не параллельно, int)

N: 100; Решения: 6; Прошедшее время: 0,05 мс (не параллельно, int)

N: 1000; Решения: 55; Прошедшее время: 0,3 мс (не параллельно, int)

N: 10000; Решения: 543; Прошедшее время: 13,1 мс (не параллельно, int)

N: 100000; Решения: 5512; Прошедшее время: 849,4 мс (параллельное, long)

Вы должны использовать long когда N больше 36340, потому что когда оно возведено в квадрат, оно переполняет значение int max. Наконец, параллельная версия начинает становиться лучше, чем простая, когда N составляет около 23000, с int s.

Ответ 5

#include<iostream>
#include<math.h>
int main()
{
    int N = 10000;
    int c = 0;
    for (int x = 2; x < N; x+=2)
    {
        for (int y = x; y < (N - x); y+=2)
        {
            auto z = sqrt(x * x + y * y + 1);
            if(x+y+z>N){
                break;
            }
            if (z - (int) z == 0)
            {
                c++;
            }
        }
    }
    std::cout<<c;
}

Это моё решение. При тестировании предыдущих решений этой проблемы я обнаружил, что x, y всегда четные, а z нечетные. Я не знаю математической природы, стоящей за этим, я сейчас пытаюсь понять это.

Ответ 6

Я хочу, чтобы это было сделано в С#, и оно должно охватывать все контрольные примеры, основанные на условии, указанном в вопросе.

Базовый код, преобразованный в long для обработки верхнего предела N <= 100000, с каждой добавленной оптимизацией, которую я мог. Я использовал альтернативные формы из @Mat (+1) Alpha-запроса Wolfram для максимально возможного предварительного вычисления. Я также провел минимальный тест идеального квадрата, чтобы избежать миллионов вызовов sqrt() на верхнем пределе:

public static void Main()
{
    int c = 0;

    long N = long.Parse(Console.ReadLine());
    long N_squared = N * N;

    double half_N_squared = N_squared / 2.0 - 0.5;
    double x_limit = N - Math.Sqrt(2) / 2.0 * Math.Sqrt(N_squared + 1);

    for (long x = 2; x < x_limit; x += 2)
    {
        long x_squared = x * x + 1;

        double y_limit = (half_N_squared - N * x) / (N - x);

        for (long y = x; y < y_limit; y += 2)
        {
            long z_squared = x_squared + y * y;
            int digit = (int) z_squared % 10;

            if (digit == 3 || digit == 7)
            {
                continue;  // minimalist non-perfect square elimination
            }

            long z = (long) Math.Sqrt(z_squared);

            if (z * z == z_squared)
            {
                c++;
            }
        }
    }

    Console.WriteLine(c);
}

Я следовал за трендом и оставил "вырожденное решение", как подразумевается в коде OP, но явно не указано.