【问题标题】:How do I solve this billion iteration sum? Project Euler problem 94, Python我如何解决这个十亿迭代和?欧拉计划问题 94,Python
【发布时间】:2020-11-20 16:57:56
【问题描述】:

问题 -- 问题 94,Project Euler -- Python -- 几乎等边三角形

*很容易证明,不存在边长为整数、面积为整数的等边三角形。然而,几乎等边三角形 5-5-6 的面积为 12 个正方形单位。

我们将定义一个几乎等边三角形为两个边相等且第三条边相差不超过一个单位的三角形。

求所有边长和面积为整数且周长不超过十亿 (1,000,000,000) 的几乎等边三角形的周长之和。*

输出:518408346

我有一个非常低效的解决方案,如下:

  1. 我尝试暴力破解(花费了相当长的时间)(使用 for 循环)
  2. 然后我在条件中添加了检查它是否是三角形(两条边之和大于第三条)
  3. 接下来,我在周长变量中定义了值,并为面积变量使用了简单的几何图形
  4. 接下来,为了检查它们是否有效,对于周界,我使用了一种类型方法来检查它是否是 int 并且因为我不能对 area 做同样的事情(因为即使它是一个整数,它也是一部分float 类),我使用模 1 == 0 来检查它是否是整数
  5. 如果它满足条件,我将它添加到我的最终“周长总和”变量中。

代码如下:

import time
start = time.time()
plst = 0
o = 100000
for i in range(1,o):
    j = i+1
    if 2*i > j and j > 0:
        p = 2*i + j
        a = 0.5 * j* ((i**2 - (j/2)**2)**0.5)
        # if p > 100:
        #     break
        # else:
        #     pass
        if type(p) == int and a % 1 == 0:
            print(i)
            print('\t', p)
            plst += p
        else:
            pass
    else:
        pass
for k in range(1,o):
    b = k-1
    if 2*k > b and b > 0:
        p = 2*k + b
        a = 0.5 * b* ((k**2 - (b/2)**2)**0.5)
        # if p > 100:
        #     break
        # else:
        #     pass
        if type(p) == int and a % 1 == 0:
            print(k)
            print('\t', p)
            plst += p
        else:
            pass
    else:
        pass
print(plst)
print(f'{time.time() - start} seconds')

注意:

  1. 是的,我知道我的 o 值为 100000。

  2. 是的,我注释掉了检查周长变量是否超过十亿的条件,因为我在达到十亿之前遇到了问题。

  3. 不要介意注释代码。 我有两个问题需要帮助 按优先顺序排列

  4. 当我运行代码时,直到 10000,它运行良好且快速,但是当我跳到 100000 时,我意识到输出中弹出了一些不必要的值。(如 93686、93686、 93687)

    A.我称之为不必要的,因为当我搜索时,一个网站显示了一个最终会产生答案的值表,而 93686 不是其中的一部分。(这是网站,https://blog.dreamshire.com/project-euler-94-solution/

    B.但是我在我的程序中找不到错误,我尝试了这些“不必要”数字的面积,令人惊讶的是,面积竟然是一个整数,这让我感到惊讶,这导致我得到了错误的答案

那么有人能解释一下错误是什么吗?

2。有人可以给我一个有效的解决方案,并用 Python 进行适当的解释吗?

我使用 python 3.8

【问题讨论】:

标签: python python-3.x


【解决方案1】:

您可以完全不使用平方根来获得结果。事实上,它可以只使用加法和乘法来获得:

如果我们将一个几乎等边三角形定义为 (a,a,c),那么 c 是 a+1 或 a-1。 (a,a,c) 三角形的高度可以使用由底边 (c/2) 的一半构成的直角三角形来计算,并以 a 为斜边。

                                               Pythagorean
                        /|\                |\  Triangle  
                       / | \               | \   
 c = a +/- 1        a /  |  \ a            |  \ a 
                     /  h|   \           h |   \  
                    /    |    \            |    \
                   /_____|_____\           |_____\
                         c                   c/2
                     c/2   c/2 

要生成整数曲面,hc/2 必须是整数,并且 (h,c/2,a) 必须是毕达哥拉斯三元组。

基于此,我们可以遍历c 的可能(偶数)值,同时依次推进相应的h。随着c 的增加,h 也会增加,这将允许我们简单地将 h^2 与 a^2 - (c/2)^2 匹配以检测毕达哥拉斯三元组。

def almostEqui(maxPerimeter):
    h = 1
    for c in range(2,maxPerimeter//3+2,2):
        cc4 = c*c//4
        for a in (c-1,c+1):
            hh = a*a - cc4
            while h*h < hh:
                h += 1
            if h*h == hh and 2*a+c <= maxPerimeter:
                yield a,a,c
        h -= 1 

def euler94():
    return sum(a+b+c for a,b,c in almostEqui(1000000000))

euler94() 函数在我的笔记本电脑上在 2.21 分钟(133 秒)内返回 518408346。这比计算平方根大约快 2.5 倍。

[编辑] 一种更快的方法是生成毕达哥拉斯原始三元组,并检查它们在加倍时是否形成几乎等边三角形:

这是一个基于欧几里得 Primitive Pythagorean triples 的互质公式的生成器函数:

from math import gcd
def genTriples(k): # primitive pythagorean triangles a,b,c where a<b<c<k
    n,m = 1,2
    while m*m+1<k:                  # while c<k (for largest m producing c)
        if n>=m: n,m = m%2,m+1      # n reached m, advance m, reset n
        c = m*m+n*n                 # compute c 
        if c >= k: n=m;continue     # skip remaining n when c >= k
        if gcd(n,m) == 1:           # trigger on coprimes
            yield m*m-n*n,2*m*n,c   # return a,b,c triple
        n += 2                      # advance n, odd with evens

输出:

from time import time
start=time()

N = 1000000000
result = 0
for a,b,c in genTriples(N//3+2):
    if abs(2*a-c)==1:               # (c,c,2*a) almost equilateral triangle
        result += 2*a+2*c
    if abs(2*b-c)==1:               # (c,c,2*b) almost equilateral triangle
        result += 2*b+2*c

print(time()-start)  # 51 seconds
print(result)        # 518408346

这比以前的方法快 2 倍以上。

【讨论】:

  • 我有一个类似的解决方案可以在一分钟内运行。我注意到我的代码更快的一种方法是它不检查等边三角形。如果等边三角形的边长为整数,则它们的面积永远不会是整数。因此,您需要更改此行“for a in (c-1,c+1):”以跳过 c。
【解决方案2】:

“不必要的”数字是由于浮点精度和舍入误差造成的。例如,对于边三角形(93686、93686、93687),分别计算公式的各个部分,我们得到area = 93687 * 81134.16730176011 / 2。这显然不是整数,但 area % 1 == 0 返回 True

您可以尝试使用Decimal 来提高精度。或者,您可以使用integer is either a perfect square or its square root is an irrational number 的定理。要使面积成为整数,平方根内的项必须是完全平方。 b 也必须是偶数,否则a**2 - (b/2)**2 不能是整数。

使用函数is_square 来检查项是否为正方形,您的代码可以像这样进行更正和稍微优化:

from math import isqrt

def is_square(i):
    return i == isqrt(i) ** 2

for a in range(3, o, 2):  
    if is_square(a ** 2 - ((a+1)//2) ** 2):
        p = 3 * a + 1
        print(a, p)
        plst += p
    if is_square(a ** 2 - ((a-1)//2) ** 2):
        p = 3 * a - 1
        print(a, p)
        plst += p

(条件type(p) == int是多余的,因为p总是一个整数;如果我们在a = 2开始循环,2*a &gt; bb &gt; 0也是多余的)

这可能是您实现蛮力算法的最快速度,我将其留给其他人来演示必须使用数学存在的更有效的算法。

【讨论】:

  • 嘿斯图尔特,非常感谢这个解决方案,它工作正常,感谢您清除我的疑问(浮点精度)。这给我带来了很多困惑。您的解决方案非常详细
【解决方案3】:

你可以决定使用函数:

def y(x):
  s = x**2 - 1
  if s % 3 == 0:
    h = (s / 3) ** 0.5
    if int(h) == h and h > 0:
      return int(h)

def hypoteneus(x):
  s = 2 * x
  uppr, lowr = s + 1, s - 1
  if uppr % 3 == 0:
    return uppr // 3
  if lowr % 3 == 0:
    return lowr // 3

def triangle(x):
  hyp = hypoteneus(x)
  height = y(x)
  if height is not None and hyp is not None:
    base = (hyp ** 2 - height **2) ** 0.5
    if base == int(base) and base > 0:
      return hyp, hyp, int(base) * 2

def  almost_eq_triangles(x):
  per = 0
  p1 = 0
  for i in range(1, int(1e9)):
    b = triangle(i)
    if b:
      p1 = int(sum(b))
      per +=p1
      if p1 > x:
        print('-'*(68 + len(str(p1))))
        print(f'{i:<10}: {str(b):<40}\t Perimeter: {p1}\tSum:{per}')
        print(f'\n\t\t{"Total sum of Perimeters (exclude last)"}: {per - p1}')
        print(f'\t\t{"Total sum of Perimeters (include last)"}: {per}')
        break
      print(f'{i:<10}: {str(b):<40}\t Perimeter: {p1}')

例如,如果我想要周长小于 10000000 的所有几乎等边三角形,我们可以这样做:

     almost_eq_triangles(10000000)
7         : (5, 5, 6)                                    Perimeter: 16
26        : (17, 17, 16)                                 Perimeter: 50
97        : (65, 65, 66)                                 Perimeter: 196
362       : (241, 241, 240)                              Perimeter: 722
1351      : (901, 901, 902)                              Perimeter: 2704
5042      : (3361, 3361, 3360)                           Perimeter: 10082
18817     : (12545, 12545, 12546)                        Perimeter: 37636
70226     : (46817, 46817, 46816)                        Perimeter: 140450
262087    : (174725, 174725, 174726)                     Perimeter: 524176
978122    : (652081, 652081, 652080)                     Perimeter: 1956242
3650401   : (2433601, 2433601, 2433602)                  Perimeter: 7300804
----------------------------------------------------------------------------
13623482  : (9082321, 9082321, 9082320)                  Perimeter: 27246962

                Total sum of Perimeters (exclude last): 9973078
                Total sum of Perimeters (include last): 37220040

【讨论】:

  • 谢谢你,你能告诉我这个程序的运行时间是什么吗?另外我认为你错过了一些三角形,因为输出是 518408346。
  • @TejasA 我实际上对 518408346 很好奇。因为这个数字不会产生一个几乎等边的三角形
  • 518408346 不会产生几乎等边三角形,输出是 518408346,因为问题是要找到所有几乎等边三角形的总和,而 518408346 是总和
  • 如果你只对这个数字感兴趣,你可以通过使用不同的公式来暗示其他人
猜你喜欢
  • 1970-01-01
  • 2011-05-20
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2011-02-08
  • 1970-01-01
相关资源
最近更新 更多