题目75：找出只能以一种方式构成直角三角形的电线长度的数目

欧拉计划

Singular integer right triangles

It turns out that 12 cm is the smallest length of wire that can be bent to form an integer sided right angle triangle in exactly one way, but there are many more examples.

12 cm: (3,4,5)
24 cm: (6,8,10)
30 cm: (5,12,13)
36 cm: (9,12,15)
40 cm: (8,15,17)
48 cm: (12,16,20)

In contrast, some lengths of wire, like 20 cm, cannot be bent to form an integer sided right angle triangle, and other lengths allow more than one solution to be found; for example, using 120 cm it is possible to form exactly three different integer sided right angle triangles.

120 cm: (30,40,50), (20,48,52), (24,45,51)

Given that L is the length of the wire, for how many values of L ≤ 1,500,000 can exactly one integer sided right angle triangle be formed?

题目：

事实证明 12 cm 是最短的只能够以一种方式弯曲成整数边直角三角形的电线长度。但是还有很多其他例子：

12 cm: (3,4,5)
24 cm: (6,8,10)
30 cm: (5,12,13)
36 cm: (9,12,15)
40 cm: (8,15,17)
48 cm: (12,16,20)

与之相反，有一些长度，例如 20cm，不能够弯曲形成整数边的直角三角形。而其他的一些长度则可以以多于一种的方式弯曲成整数边直角三角形；例如，120cm 可以构成三种不同的整数边直角三角形。

120 cm: (30,40,50), (20,48,52), (24,45,51)

L 为电线的长度，对于 L ≤ 1,500,000，有多少个只能以一种方式构成整数边直角三角形？

jerryxjr1220

这题运用了欧几里得公式：
以下是出自wikipedia
Euclid's formula[1] is a fundamental formula for generating Pythagorean triples given an arbitrary pair of positive integers m and n with m > n. The formula states that the integers

a = m^2 - n^2 ,\ \, b = 2mn ,\ \, c = m^2 + n^2
form a Pythagorean triple. The triple generated by Euclid's formula is primitive if and only if m and n are coprime and m − n is odd. If both m and n are odd, then a, b, and c will be even, and so the triple will not be primitive; however, dividing a, b, and c by 2 will yield a primitive triple if m and n are coprime.[2]

Every primitive triple arises from a unique pair of coprime numbers m, n, one of which is even. It follows that there are infinitely many primitive Pythagorean triples. This relationship of a,b and c to m and n from Euclid's formula is referenced throughout the rest of this article.

Despite generating all primitive triples, Euclid's formula does not produce all triples—for example, (9, 12, 15) cannot be generated using integer m and n. This can be remedied by inserting an additional parameter k to the formula. The following will generate all Pythagorean triples uniquely:

a = k\cdot(m^2 - n^2)   ,\ \, b = k\cdot(2mn) ,\ \, c = k\cdot(m^2 + n^2)
where m, n, and k are positive integers with m > n, m − n odd, and with m and n coprime.

That these formulas generate Pythagorean triples can be verified by expanding a2 + b2 using elementary algebra and verifying that the result coincides with c2. Since every Pythagorean triple can be divided through by some integer k to obtain a primitive triple, every triple can be generated uniquely by using the formula with m and n to generate its primitive counterpart and then multiplying through by k as in the last equation.

Many formulas for generating triples with particular properties have been developed since the time of Euclid.

jerryxjr1220

ans =  161667

import math  
  
L = 1500000  
maxm = int(math.sqrt(L/2))  
triple = [0 for x in range(0,L+1)]  
ans = 0

for m in range(2,maxm):  
    for n in range (1,m):  
        if ((m+n)%2)==1 and math.gcd(m,n)==1:  
            a=m*m-n*n  
            b=2*m*n  
            c=m*m+n*n  
            p=a+b+c  
            while p<=L:  
                triple[p]=triple[p]+1  
                if triple[p]==1:  
                    ans=ans+1  
                if triple[p]==2:  
                    ans=ans-1  
                p=p+a+b+c  
print('ans = ',ans) 

najin

用的matlab 结果是
ans =

      161667

时间已过 0.105196 秒。
>>

fc1735

from math import gcd

L = [0] * 1500001
for n in range(1, int(len(L) ** 0.5) // 2):
  for m in range(n + 1, 999999999999999999999, 2):
    if gcd(m, n) != 1:
      continue
    l = (m + n) * m * 2
    if l > len(L):
      break
    for i in range(l, len(L), l):
      L[i] += 1

print(sum(filter(lambda x: x == 1, L)))

https://github.com/devinizz/project_euler/blob/master/page02/75.py

guosl

这题我们要用本原直角三角形的结论，具体内容参考：百度百科的本原直角三角形。

/*
答案：161667
耗时：2.28576秒（单核）
      0.393828秒（六核）
*/
#include <iostream>
#include <cmath>
#include <omp.h>
using namespace std;

inline int gcd(int x, int y)//本原直角三角形的两个参数必须互素
{
  if (y == 0)
    return x;
  int z = x % y;
  while (z != 0)
  {
    x = y;
    y = z;
    z = x % y;
  }
  return y;
}

int chk(int n)//计算边长为n的本原直角三角形的个数
{
  n /= 2;
  int nCount = 0;//计数
  int d = (int)sqrt((double)n);
  for (int i = 2; i <= d; ++i)//枚举本原直角三角形的大参数
  {
    if (n % i == 0)
    {
      int p = n / i;
      if (gcd(i, p) == 1)
      {
        int x = i, y = p - i; //解出本原直角三角形的两个可能的参数
        if (y > 0 && x > y && ((x - y) & 1) != 0)//检查是否合法
          ++nCount;//存在计数累加
      }
    }
  }
  return nCount;//返回计数
}

int c[1500001];//c[i]表示周长为i的直角三角形的个数

int main(void)
{
  double t = omp_get_wtime();
  int nCount = 0;
#pragma omp parallel shared(c) reduction(+:nCount)
  {
#pragma omp for schedule(dynamic)
    for (int L = 12; L <= 1500000; L += 2)//枚举本原直角三角形的周长（本原直角三角形的周长一定是偶数）
    {
      int p = chk(L);//得到边长为L的本原直角三角形的个数
      int k = 1;
      while (k*L <= 1500000)
      {
#pragma omp atomic        
        c[L*(k++)] += p;//累加边长为本原直角三角形边长倍数的直角三角形个数
      }
    }
#pragma omp for
    for (int L = 12; L <= 1500000; L += 2)
    {
      if (c[L] == 1)
        ++nCount;//统计结果
    }
  }
  t = omp_get_wtime() - t;
  cout << nCount << endl << t << endl;
  return 0;
}

guosl

这个题目如果用GPU加速则时间更短。

/*
答案：161667
耗时：0.119166秒（GTX970）
*/
#include <iostream>
#include <cmath>
#include <openacc.h>
using namespace std;

#pragma acc routine
inline int gcd(int x, int y)//本原直角三角形的两个参数必须互素
{
  if (y == 0)
    return x;
  int z = x % y;
  while (z != 0)
  {
    x = y;
    y = z;
    z = x % y;
  }
  return y;
}

#pragma acc routine
int chk(int n)//计算边长为n的本原直角三角形的个数
{
  n /= 2;
  int nCount = 0;//计数
  int d = (int)sqrt((double)n);
  for (int i = 2; i <= d; ++i)//枚举本原直角三角形的大参数
  {
    if (n % i == 0)
    {
      int p = n / i;
      if (gcd(i, p) == 1)
      {
        int x = i, y = p - i; //解出本原直角三角形的两个可能的参数
        if (y > 0 && x > y && ((x - y) & 1) != 0)//检查是否合法
          ++nCount;//存在计数累加
      }
    }
  }
  return nCount;//返回计数
}

int c[1500001];//c[i]表示周长为i的直角三角形的个数

int main(void)
{
  int nCount = 0;
#pragma acc kernels
  {
#pragma acc loop independent
    for (int L = 12; L <= 1500000; L += 2)//枚举本原直角三角形的周长（本原直角三角形的周长一定是偶数）
    {
      int p = chk(L);//得到边长为L的本原直角三角形的个数
      int k = 1;
      while (k*L <= 1500000)
      {
        int i=L*(k++);
#pragma acc atomic
        c[i] += p;//累加边长为本原直角三角形边长倍数的直角三角形个数
      }
    }
#pragma acc loop independent reduction(+:nCount)
    for (int L = 12; L <= 1500000; L += 2)
    {
      if (c[L] == 1)
        ++nCount;//统计结果
    }
  }
  cout << nCount << endl;
  return 0;
}