生成唯一的有序勾股三元组

毕达哥拉斯三胞胎就是一个很好的例子。 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个版本中,这比版本2快5倍,但是 仍然 立方对 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 只有一次。它不仅比 n 在1000中,它是二次开的 n 因此,加速速度增加为 n 生长。

我经常遇到这种改进,除了最简单的用法(例如遍历数组),我不相信“计数循环”。

更新: 显然,我应该指出一些关于v4的东西,这些东西很容易被忽略。

1. 两个 的 while 循环由以下值控制: Z (一个直接,另一个间接通过 Z )内部 虽然 实际上是在加速外层 虽然 而不是与之正交。 查看循环在做什么是很重要的,而不仅仅是计算有多少个循环。

2. v4中的所有计算都是严格的整数算术。相比之下,浮点转换和浮点计算的成本很高。

3. v4在常量内存中运行,只需要三个整型变量。没有要分配和初始化的数组或哈希表(并且可能会导致内存不足错误)。

4. 最初的问题允许 X , Y 和 X 在同一范围内变化。v1..v4遵循这个模式。

下面是一组不太科学的时序(使用在我的旧笔记本电脑上的Eclipse和其他运行的东西……),其中“使用x,y,z”是通过实例化一个三值的三重对象来实现的,并把它放在一个数组中。(对于这些跑步, n 设为10000,每种情况下产生12471个三倍。)

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源代码的请求,为其他变化,我的时间,如果我误执行任何东西。

比目前任何解决方案都快得多。通过三元树找到三连体。

Wolfram 说:

霍尔(1970)和罗伯茨(1977)证明了这是一个原始的毕达哥拉斯三重奏,如果且仅如果

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

其中m是矩阵u,a,d的有限积。

在这里我们有一个公式来产生每一个原始三重态。

在上面的公式中,斜边一直在增长,所以很容易检查最大长度。

在蟒蛇中:

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个元素

作为参考, accepted answer (in python) 12471个元素用了38秒。

只是为了好玩,把上限定为一百万 list(gen_all_pyth_trips(10**6)) 在2.66秒内返回1980642个元素(在3秒内几乎是200万个三倍)。 list(gen_all_pyth_trips(10**7)) 当列表变大时,我的电脑就要崩溃了,它会消耗掉内存的最后一点。做一些类似的事情 sum(1 for _ in gen_all_pyth_trips(10**7)) 绕过这个限制,在30秒内返回23471475个元素。

有关所用算法的更多信息,请参阅 沃尔夫拉姆 和 Wikipedia .

您应该定义x<y<z。

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倍)。

前面列出的生成算法 Pythagorean triplets 是否所有对幼稚方法的修改都源于基本关系? a^2 + b^2 = c^2 哪里 (a, b, c) 是正整数的三联体。结果表明,勾股儿三连体满足一些相当显著的关系,这些关系可用于生成所有勾股儿三连体。

Euclid 发现了第一个这样的关系。他决定每三个毕达哥拉斯 (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)

然后生成勾股三元组,生成相对素数正整数。 米 和 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));            
    }
}

维基百科关于 Formulas for generating Pythagorean triples 包含其他此类公式。

算法可以针对速度、内存使用、简单性和其他方面进行调整。

这里是一个 pythagore_triplets 算法以节省内存和简单性为代价进行了速度调整。如果你想要的只是速度,这可能是你要走的路。

计算 list(pythagore_triplets(10000)) 我的电脑需要40秒,而Tafkas算法需要63秒,Tafkas算法可能需要几天的计算时间(以及所有其他使用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)]

请注意,如果要计算前十亿个三元组,那么这个算法甚至会在启动之前崩溃,因为内存不足。因此,对于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

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 价值观。

编辑

这个问题似乎还活着:)
由于我返回并重新访问了代码,我尝试了第二种方法,它的速度几乎是之前建议的4倍(n=10000的CPU时间的26%),因为它避免了许多不必要的计算:

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___,在C____

(x+1)_2-x_2=(x+1)(x+1)-x_2=x_2+2x+1-x_2=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中更快。

为大家干杯:)

我只是把凯尔·古利昂的答案扩展了一下,这样三个字母就按押金排列,然后是最长的一边。

它不使用numpy,但需要SortedCollection(或SortedList),例如 this one

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 module 属于 my Python library . 它是根据一系列的OEIS(代码 here 在底部),这让我 find a mistake in A121727 -)

我用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

是的,有。

好吧,现在你想知道为什么。为什么不限制它使z>y?尝试

for z in range (y+1, 1000)

旧问题,但我还是会输入我的资料。 有两种生成独特的毕达哥拉斯三倍体的一般方法。一种是按比例缩放,另一种是使用这个古老的公式。

标度基本上取一个常数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,计算2米,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

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)

第5版给乔尔·尼利。

因为x可以是'n-2'的最大值,y可以是'n-1'的最大值,范围是1..n。因为z max是n,y max是n-1,x可以是sqrt的最大值(n*n-(n-1)*(n-1))=sqrt(2*n-1),并且可以从3开始。

MaxX = ( 2 * N - 1 ) ** 0.5

for x in 3..MaxX {
  y = x+1
  z = y+1
  m = x*x + y*y
  k = z * z
  while z <= N {
     while k < m {
        z = z + 1
        k = k + (2*z) - 1
    }
    if k == m and z <= N then {
        // use x, y, z
    }
    y = y + 1
    m = m + (2 * y) - 1
  }
 }

只是检查一下,但我一直在使用下面的代码来生成毕达哥拉斯三元组。它非常快(我在这里尝试了一些例子,虽然我有点了解它们,写了自己的,回来检查了这里(2年前))。我认为这段代码正确地找到了所有的毕达哥拉斯三重奏(说出你的极限),而且速度也相当快。我用C++来制作它。

乌龙是无符号长龙,我创建了几个平方和根函数 我的根函数基本上说,如果给定数的平方根(使其成为整数(integral))的平方根不等于给定数,那么返回-1,因为它不可根。 _ Square和_-Root按照上面描述的那样做,我知道优化它的另一种方法,但是我还没有做过也没有测试过。

generate(vector<Triple>& triplist, ullong limit) {
cout<<"Please wait as triples are being generated."<<endl;
register ullong a, b, c;
register Triple trip;
time_t timer = time(0);

for(a = 1; a <= limit; ++a) {
    for(b = a + 1; b <= limit; ++b) {
        c = _root(_square(a) + _square(b));

        if(c != -1 && c <= limit) {
            trip.a = a; trip.b = b; trip.c = c;

            triplist.push_back(trip);

        } else if(c > limit)
            break;
    }
}

timer = time(0) - timer;
cout<<"Generated "<<triplist.size()<<" in "<<timer<<" seconds."<<endl;
cin.get();
cin.get();

}

告诉我你们的想法。根据我交给的老师,它生成所有原始和非原始三元组。(如果我没记错的话,她测试了100次)。

前一个编码器提供的v4结果如下

下面是一组不太科学的时序(使用在我的旧笔记本电脑上的Eclipse和其他运行的东西……),其中“使用x,y,z”是通过实例化一个三值的三重对象来实现的,并把它放在一个数组中。(对于这些运行,n设置为10000,每种情况下产生12471个三倍。)

版本4:46秒。 使用平方根:134秒。 阵列和映射:400秒。

我的结果是 生成多少个三倍:10000

正在生成三元组,请稍候。 在2秒内生成12471。

在我开始通过编译器进行优化之前。(我记得以前我用大量的特殊选项和东西从10000秒降到0秒)。我的代码还生成了所有三个值,其中100000作为在3.2分钟内1,2,hyp可以达到的上限(我认为1000000的上限需要一个小时)。

我对代码做了一点修改,并将10000限制降到1秒(没有优化)。除此之外,仔细考虑,我的可以被分解成块,并按照给定的范围(例如100000分为4个相等的块,每3个CPU(1个额外的,以备万一,希望消耗CPU时间),范围为1到25000(从1开始,限制为25000),25000到50000,50000到75000,和75000到结束。我可以这样做,看看它是否加快了速度(我将预先处理线程,而不将它们包括在执行三重函数的实际时间内)。我需要一个更精确的计时器和连接向量的方法。我认为,如果1个3.4 GHz的CPU,有8 GB的RAM,可以在1秒钟内完成10000个LIM,那么3个CPU应该在1/3秒内完成这项任务(我会像ATM那样调到更高的秒数)。

应该注意,对于A、B和C,不需要一直循环到N。

对于A,您只需要从1循环到 int(sqrt(n**2/2))+1 ,对于B, a+1 到 int(sqrt(n**2-a**2))+1 ,对于来自的C int(sqrt(a**2+b**2) 到 int(sqrt(a**2+b**2)+2 .

# To find all pythagorean triplets in a range
import math
n = int(input('Enter the upper range of limit'))
for i in range(n+1):
    for j in range(1, i):
        k = math.sqrt(i*i + j*j)
        if k % 1 == 0 and k in range(n+1):
            print(i,j,int(k))

你可以试试这个

  triplets=[]
    for a in range(1,100):
        for b in range(1,100):
            for c in range(1,100):
                if  a**2 + b**2==c**2:
                    i=[a,b,c]
                    triplets.append(i)
                    for i in triplets:
                          i.sort()
                          if triplets.count(i)>1:
                              triplets.remove(i)
    print(triplets)