Генерация уникальных, упорядоченных пифагорейских триплетов - PullRequest
27 голосов
/ 22 февраля 2009

Это программа, которую я написал для вычисления пифагорейских триплетов. Когда я запускаю программу, она печатает каждый набор триплетов дважды из-за оператора if. Можно ли как-то сказать программе, чтобы новый набор триплетов печатался только один раз? Спасибо.

import math

def main():
    for x in range (1, 1000):
        for y in range (1, 1000):
            for z in range(1, 1000):
                if x*x == y*y + z*z:
                    print y, z, x
                    print '-'*50

if __name__ == '__main__':
    main()  

Ответы [ 16 ]

73 голосов
/ 23 февраля 2009

Пифагорейские тройки являются хорошим примером для утверждения, что «for циклы считаются вредными », потому что for циклы соблазняют нас задуматься о подсчете, часто самой несущественной части задачи.

(я собираюсь придерживаться псевдокода, чтобы избежать смещения языка и чтобы упростить псевдокод, я не буду оптимизировать множественные вычисления, например, x * x и y * y.)

Версия 1 :

for x in 1..N {
    for y in 1..N {
        for z in 1..N {
            if x * x + y * y == z * z then {
                // use x, y, z
            }
        }
    }
}

- худшее решение. Он генерирует дубликаты и пересекает части пространства, которые бесполезны (например, всякий раз, когда z < y). Его сложность по времени является кубической на N.

Версия 2 , первое улучшение, связано с требованием удержания x < y < z, например:

for x in 1..N {
    for y in x+1..N {
        for z in y+1..N {
            if x * x + y * y == z * z then {
                // use x, y, z
            }
        }
    }
}

, что сокращает время выполнения и устраняет дублированные решения. Тем не менее, он все еще кубический на N; улучшение - это всего лишь уменьшение коэффициента N -куб.

Бессмысленно продолжать изучать растущие значения z после того, как z * z < x * x + y * y больше не удерживается. Этот факт мотивирует Версия 3 , первый шаг от грубой итерации за z:

for x in 1..N {
    for y in x+1..N {
        z = y + 1
        while z * z < x * x + y * y {
            z = z + 1
        }
        if z * z == x * x + y * y and z <= N then {
            // use x, y, z
        }
    }
}

Для N из 1000 это примерно в 5 раз быстрее, чем в версии 2, но это все еще куб. При N.

Следующее понимание заключается в том, что x и y являются единственными независимыми переменными; z зависит от их значений, а последнее z значение, рассмотренное для предыдущего значения y, является хорошим начальным поисковым значением для следующего значения y. Это приводит к версии 4 :

for x in 1..N {
    y = x+1
    z = y+1
    while z <= N {
        while z * z < x * x + y * y {
            z = z + 1
        }
        if z * z == x * x + y * y and z <= N then {
            // use x, y, z
        }
        y = y + 1
    }
}

, что позволяет y и z "подметать" значения выше x только один раз. Он не только более чем в 100 раз быстрее для N из 1000, он также квадратичен для N, поэтому ускорение увеличивается с ростом N.

Я встречал такого рода улучшения достаточно часто, чтобы не доверять «счетным циклам» для любых, кроме самых тривиальных целей (например, обход массива).

Обновление: Очевидно, я должен был указать на некоторые вещи в V4, которые легко упустить из виду.

  1. Оба из while циклов управляются значением z (один напрямую, другой косвенно через квадрат z). Внутреннее while фактически ускоряет внешнее while, а не ортогонально ему. Важно посмотреть, что делают циклы, а не просто подсчитать, сколько существует циклов.

  2. Все вычисления в V4 являются строго целочисленной арифметикой. Преобразование в / из плавающей запятой, а также вычисления с плавающей запятой дорогостоящие в сравнении.

  3. V4 работает в постоянной памяти, требуя только трех целочисленных переменных. Нет массивов или хеш-таблиц, которые можно было бы выделить и инициализировать (и, возможно, вызвать ошибку нехватки памяти).

  4. Исходный вопрос позволял всем x, y и x варьироваться в одном и том же диапазоне. V1..V4 следовал этому образцу.

Ниже приведен не очень научный набор моментов времени (использующий Java под Eclipse на моем старом ноутбуке с запущенным другим ...), где «use x, y, z» было реализовано путем создания экземпляра объекта Triple с три значения и положить его в ArrayList. (Для этих прогонов N было установлено на 10000, что в каждом случае составляло 12 471 тройку.)

Version 4:           46 sec.
using square root:  134 sec.
array and map:      400 sec.

Алгоритм "массива и карты" по существу :

squares = array of i*i for i in 1 .. N
roots = map of i*i -> i for i in 1 .. N
for x in 1 .. N
    for y in x+1 .. N
        z = roots[squares[x] + squares[y]]
        if z exists use x, y, z

Алгоритм "использования квадратного корня" , по существу :

for x in 1 .. N
    for y in x+1 .. N
        z = (int) sqrt(x * x + y * y)
        if z * z == x * x + y * y then use x, y, z

Фактический код для V4:

public Collection<Triple> byBetterWhileLoop() {
    Collection<Triple> result = new ArrayList<Triple>(limit);
    for (int x = 1; x < limit; ++x) {
        int xx = x * x;
        int y = x + 1;
        int z = y + 1;
        while (z <= limit) {
            int zz = xx + y * y;
            while (z * z < zz) {++z;}
            if (z * z == zz && z <= limit) {
                result.add(new Triple(x, y, z));
            }
            ++y;
        }
    }
    return result;
}

Обратите внимание, что x * x - это , рассчитанный во внешнем цикле (хотя я не удосужился кешировать z * z); аналогичные оптимизации выполняются в других вариантах.

Я буду рад предоставить исходный код Java по запросу для других вариантов, которые я рассчитал, на случай, если я что-то неправильно реализовал.

40 голосов
/ 25 ноября 2011

Значительно быстрее, чем любое другое решение. Находит тройки через тройное дерево.

Вольфрам говорит:

Холл (1970) и Робертс (1977) доказывают, что это примитивная пифагорейская тройка тогда и только тогда, когда

(a,b,c)=(3,4,5)M

где M - конечное произведение матриц U, A, D.

И там у нас есть формула для генерации каждой примитивной тройки.

В приведенной выше формуле гипотенуза постоянно растет, поэтому довольно легко проверить максимальную длину.

В Python:

import numpy as np

def gen_prim_pyth_trips(limit=None):
    u = np.mat(' 1  2  2; -2 -1 -2; 2 2 3')
    a = np.mat(' 1  2  2;  2  1  2; 2 2 3')
    d = np.mat('-1 -2 -2;  2  1  2; 2 2 3')
    uad = np.array([u, a, d])
    m = np.array([3, 4, 5])
    while m.size:
        m = m.reshape(-1, 3)
        if limit:
            m = m[m[:, 2] <= limit]
        yield from m
        m = np.dot(m, uad)

Если вам нужны все тройки, а не только примитивы:

def gen_all_pyth_trips(limit):
    for prim in gen_prim_pyth_trips(limit):
        i = prim
        for _ in range(limit//prim[2]):
            yield i
            i = i + prim

list(gen_prim_pyth_trips(10**4)) потребовалось 2,81 миллисекунды, чтобы вернуться с 1593 элементами, в то время как list(gen_all_pyth_trips(10**4)) потребовалось 19,8 миллисекунд, чтобы вернуться с 12471 элементами

Для справки: принятый ответ (в python) заняло 38 секунд для 12471 элементов.

Просто для удовольствия, установив верхний предел в один миллион list(gen_all_pyth_trips(10**6)) возвращается за 2,66 секунды с 1980642 элементами (почти 2 миллиона в три раза за 3 секунды). list(gen_all_pyth_trips(10**7)) ставит мой компьютер на колени, так как список становится настолько большим, что он потребляет каждый последний бит оперативной памяти. Выполнение чего-то вроде sum(1 for _ in gen_all_pyth_trips(10**7)) обходит это ограничение и возвращает через 30 секунд 23471475 элементов.

Для получения дополнительной информации об используемом алгоритме, ознакомьтесь со статьями по Wolfram и Wikipedia .

13 голосов
/ 22 февраля 2009

Вы должны определить x

for x in range (1, 1000):
    for y in range (x + 1, 1000):
            for z in range(y + 1, 1000):

Еще одна хорошая оптимизация - использовать только x и y и вычислять zsqr = x * x + y * y. Если zsqr - квадратное число (или z = sqrt (zsqr) - целое число), это триплет, иначе нет. Таким образом, вам нужно всего два цикла вместо трех (для вашего примера это примерно в 1000 раз быстрее).

11 голосов
/ 23 февраля 2009

Ранее перечисленные алгоритмы генерации пифагорейских триплетов представляют собой все модификации наивного подхода, основанного на базовом отношении a^2 + b^2 = c^2, где (a, b, c) - триплет натуральных чисел. Оказывается, пифагорейские триплеты удовлетворяют некоторым довольно замечательным отношениям, которые можно использовать для генерации всех пифагорейских триплетов.

Евклид обнаружил первое такое отношение. Он определил, что для каждой пифагорейской тройки (a, b, c), возможно, после переупорядочения a и b существуют относительно простые натуральные числа m и n с m > n, по крайней мере одно из которых является четным положительное целое число k такое, что

a = k (2mn)
b = k (m^2 - n^2)
c = k (m^2 + n^2)

Затем, чтобы сгенерировать пифагорейские триплеты, сгенерируйте относительно простые натуральные числа m и n разной четности, а также натуральное число k и примените приведенную выше формулу.

struct PythagoreanTriple {
    public int a { get; private set; }
    public int b { get; private set; }
    public int c { get; private set; }

    public PythagoreanTriple(int a, int b, int c) : this() {
        this.a = a < b ? a : b;
        this.b = b < a ? a : b;
        this.c = c;
    }

    public override string ToString() {
        return String.Format("a = {0}, b = {1}, c = {2}", a, b, c);
    }

    public static IEnumerable<PythagoreanTriple> GenerateTriples(int max) {
        var triples = new List<PythagoreanTriple>();
        for (int m = 1; m <= max / 2; m++) {
            for (int n = 1 + (m % 2); n < m; n += 2) {
                if (m.IsRelativelyPrimeTo(n)) {
                    for (int k = 1; k <= max / (m * m + n * n); k++) {
                        triples.Add(EuclidTriple(m, n, k));
                    }
                }
            }
        }

        return triples;
    }

    private static PythagoreanTriple EuclidTriple(int m, int n, int k) {
        int msquared = m * m;
        int nsquared = n * n;
        return new PythagoreanTriple(k * 2 * m * n, k * (msquared - nsquared), k * (msquared + nsquared));
    }
}

public static class IntegerExtensions {
    private static int GreatestCommonDivisor(int m, int n) {
        return (n == 0 ? m : GreatestCommonDivisor(n, m % n));
    }

    public static bool IsRelativelyPrimeTo(this int m, int n) {
        return GreatestCommonDivisor(m, n) == 1;
    }
}

class Program {
    static void Main(string[] args) {
        PythagoreanTriple.GenerateTriples(1000).ToList().ForEach(t => Console.WriteLine(t));            
    }
}

Статья в Википедии о Формулы для генерации пифагорейских троек содержит другие подобные формулы.

8 голосов
/ 23 февраля 2009

Алгоритмы могут быть настроены на скорость, использование памяти, простоту и другие вещи.

Вот алгоритм pythagore_triplets, настроенный на скорость, за счет использования памяти и простоты. Если все, что вам нужно, это скорость, это может быть путь.

Расчет list(pythagore_triplets(10000)) занимает на моем компьютере 40 секунд, против 63 секунд для алгоритма ΤΖΩ seconds и, возможно, дней вычисления для алгоритма Тафкаса (и всех других алгоритмов, которые используют 3 встроенных цикла вместо 2).

def pythagore_triplets(n=1000):
   maxn=int(n*(2**0.5))+1 # max int whose square may be the sum of two squares
   squares=[x*x for x in xrange(maxn+1)] # calculate all the squares once
   reverse_squares=dict([(squares[i],i) for i in xrange(maxn+1)]) # x*x=>x
   for x in xrange(1,n):
     x2 = squares[x]
     for y in xrange(x,n+1):
       y2 = squares[y]
       z = reverse_squares.get(x2+y2)
       if z != None:
         yield x,y,z

>>> print list(pythagore_triplets(20))
[(3, 4, 5), (5, 12, 13), (6, 8, 10), (8, 15, 17), (9, 12, 15), (12, 16, 20)]

Обратите внимание, что если вы собираетесь вычислить первый миллиард триплетов, то этот алгоритм завершится сбоем еще до его запуска из-за ошибки нехватки памяти. Таким образом, алгоритм ΤΖΙΟΥΩ probably, вероятно, является более безопасным выбором для больших значений n.

Кстати, вот алгоритм Тафкаса, переведенный на python для целей моих тестов производительности. Его недостаток - 3 петли вместо 2.

def gcd(a, b):
  while b != 0:
    t = b
    b = a%b
    a = t
  return a

def find_triple(upper_boundary=1000):
  for c in xrange(5,upper_boundary+1):
    for b in xrange(4,c):
      for a in xrange(3,b):
        if (a*a + b*b == c*c and gcd(a,b) == 1):
          yield a,b,c
6 голосов
/ 22 февраля 2009
def pyth_triplets(n=1000):
    "Version 1"
    for x in xrange(1, n):
        x2= x*x # time saver
        for y in xrange(x+1, n): # y > x
            z2= x2 + y*y
            zs= int(z2**.5)
            if zs*zs == z2:
                yield x, y, zs

>>> print list(pyth_triplets(20))
[(3, 4, 5), (5, 12, 13), (6, 8, 10), (8, 15, 17), (9, 12, 15), (12, 16, 20)]
Алгоритм

V.1 имеет монотонно увеличивающиеся значения x.

EDIT

Кажется, этот вопрос все еще жив:)
Так как я вернулся и повторно посетил код, я попробовал второй подход, который почти в 4 раза быстрее (около 26% процессорного времени для N = 10000), чем мое предыдущее предложение, поскольку он избегает большого количества ненужных вычислений:

def pyth_triplets(n=1000):
    "Version 2"
    for z in xrange(5, n+1):
        z2= z*z # time saver
        x= x2= 1
        y= z - 1; y2= y*y
        while x < y:
            x2_y2= x2 + y2
            if x2_y2 == z2:
                yield x, y, z
                x+= 1; x2= x*x
                y-= 1; y2= y*y
            elif x2_y2 < z2:
                x+= 1; x2= x*x
            else:
                y-= 1; y2= y*y

>>> print list(pyth_triplets(20))
[(3, 4, 5), (6, 8, 10), (5, 12, 13), (9, 12, 15), (8, 15, 17), (12, 16, 20)]

Обратите внимание, что этот алгоритм имеет увеличивающиеся значения z.

Если алгоритм был преобразован в C - где, будучи ближе к металлу, умножения занимают больше времени, чем сложения - можно было бы минимизировать необходимые умножения, учитывая тот факт, что шаг между последовательными квадратами:

(x + 1) ² - x² = (x + 1) (x + 1) - x² = x² + 2x + 1 - x² = 2x + 1

так что все внутренние x2= x*x и y2= y*y будут преобразованы в сложения и вычитания, подобные этому:

def pyth_triplets(n=1000):
    "Version 3"
    for z in xrange(5, n+1):
        z2= z*z # time saver
        x= x2= 1; xstep= 3
        y= z - 1; y2= y*y; ystep= 2*y - 1
        while x < y:
            x2_y2= x2 + y2
            if x2_y2 == z2:
                yield x, y, z
                x+= 1; x2+= xstep; xstep+= 2
                y-= 1; y2-= ystep; ystep-= 2
            elif x2_y2 < z2:
                x+= 1; x2+= xstep; xstep+= 2
            else:
                y-= 1; y2-= ystep; ystep-= 2

Конечно, в Python дополнительный байт-код на самом деле замедляет алгоритм по сравнению с версией 2, но я бы поспорил (без проверки :), что V.3 быстрее в C.

Приветствует всех:)

4 голосов
/ 14 декабря 2016

Я оправдываю расширенный ответ Кайла Гуллиона, чтобы тройки сортировались по гипотенузе, а затем по самой длинной стороне.

Он не использует numpy, но требует SortedCollection (или SortedList), такой как this

def primitive_triples():
""" generates primitive Pythagorean triplets x<y<z
sorted by hypotenuse z, then longest side y
through Berggren's matrices and breadth first traversal of ternary tree
:see: https://en.wikipedia.org/wiki/Tree_of_primitive_Pythagorean_triples
"""
key=lambda x:(x[2],x[1])
triples=SortedCollection(key=key)
triples.insert([3,4,5])
A = [[ 1,-2, 2], [ 2,-1, 2], [ 2,-2, 3]]
B = [[ 1, 2, 2], [ 2, 1, 2], [ 2, 2, 3]]
C = [[-1, 2, 2], [-2, 1, 2], [-2, 2, 3]]

while triples:
    (a,b,c) = triples.pop(0)
    yield (a,b,c)

    # expand this triple to 3 new triples using Berggren's matrices
    for X in [A,B,C]:
        triple=[sum(x*y for (x,y) in zip([a,b,c],X[i])) for i in range(3)]
        if triple[0]>triple[1]: # ensure x<y<z
            triple[0],triple[1]=triple[1],triple[0]
        triples.insert(triple)

def triples():
""" generates all Pythagorean triplets triplets x<y<z 
sorted by hypotenuse z, then longest side y
"""
prim=[] #list of primitive triples up to now
key=lambda x:(x[2],x[1])
samez=SortedCollection(key=key) # temp triplets with same z
buffer=SortedCollection(key=key) # temp for triplets with smaller z
for pt in primitive_triples():
    z=pt[2]
    if samez and z!=samez[0][2]: #flush samez
        while samez:
            yield samez.pop(0)
    samez.insert(pt)
    #build buffer of smaller multiples of the primitives already found
    for i,pm in enumerate(prim):
        p,m=pm[0:2]
        while True:
            mz=m*p[2]
            if mz < z:
                buffer.insert(tuple(m*x for x in p))
            elif mz == z: 
                # we need another buffer because next pt might have
                # the same z as the previous one, but a smaller y than
                # a multiple of a previous pt ...
                samez.insert(tuple(m*x for x in p))
            else:
                break
            m+=1
        prim[i][1]=m #update multiplier for next loops
    while buffer: #flush buffer
        yield buffer.pop(0)
    prim.append([pt,2]) #add primitive to the list

код доступен в модуле math2 из моей библиотеки Python . Он протестирован с некоторыми сериями OEIS (код здесь внизу), который только что позволил мне найти ошибку в A121727 : -)

2 голосов
/ 22 февраля 2009

Я написал эту программу на Ruby, и она похожа на реализацию на python. Важная строка:

if x*x == y*y + z*z && gcd(y,z) == 1:

Затем вам нужно реализовать метод, который возвращает наибольший общий делитель (gcd) из двух данных чисел. Очень простой пример в Ruby:

def gcd(a, b)
    while b != 0
      t = b
      b = a%b
      a = t
    end
    return a
end

Полный рубиновый метон для нахождения триплетов будет:

def find_triple(upper_boundary)

  (5..upper_boundary).each {|c|
    (4..c-1).each {|b|
      (3..b-1).each {|a|
        if (a*a + b*b == c*c && gcd(a,b) == 1)
          puts "#{a} \t #{b} \t #{c}"
        end
      }
    }
  }
end
1 голос
/ 07 июня 2017
from  math import sqrt
from itertools import combinations

#Pythagorean triplet - a^2 + b^2 = c^2 for (a,b) <= (1999,1999)
def gen_pyth(n):
if n >= 2000 :
  return
ELEM =   [  [ i,j,i*i + j*j ] for i , j in list(combinations(range(1, n +   1 ), 2)) if sqrt(i*i + j*j).is_integer() ]
print (*ELEM , sep = "\n")


gen_pyth(200)
1 голос
/ 21 декабря 2016

Старый вопрос, но я все еще буду вводить свои вещи. Существует два основных способа генерирования уникальных пифагорейских троек. Один по шкале, а другой по этой архаичной формуле.

Для какого масштабирования в основном требуется постоянная n, а затем умножить базовую тройку, скажем, 3,4,5 на n. Таким образом, принимая n равным 2, мы получаем 6,8,10 нашей следующей тройки.

Пересчет

def pythagoreanScaled(n):
    triplelist = []
    for x in range(n):
        one = 3*x
        two = 4*x
        three = 5*x
        triple = (one,two,three)
        triplelist.append(triple)
return triplelist

Метод формулы использует тот факт, что если мы возьмем число x, вычислим 2m, m ^ 2 + 1 и m ^ 2-1, эти три всегда будут трифлетом Пифагора.

Формула

def pythagoreantriple(n):
    triplelist = []
    for x in range(2,n):
        double = x*2
        minus = x**2-1
        plus = x**2+1
        triple = (double,minus,plus)
        triplelist.append(triple)
    return triplelist
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...