Special pythagorean triplet
(Redirected from Special Pythagorean Triplet)
Special pythagorean triplet is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
- Task
Find the Pythagorean triplet for which a + b + c = 1000 and print the product of a, b, and c. The problem was taken from Project Euler problem 9
- Related task
11l[edit]
V PERIMETER = 1000
L(a) 1 .. PERIMETER
L(b) a + 1 .. PERIMETER
V c = PERIMETER - a - b
I a * a + b * b == c * c
print((a, b, c))- Output:
(200, 375, 425)
ALGOL 68[edit]
Uses Euclid's formula, as in the XPL0 sample but also uses the fact that M and N must be factors of half the triangle's perimeter to reduce the number of candidate M's to check. A loop is not needed to find N once a candidate M has been found.
Does not stop after the solution has been found, thus verifying there is only one solution.
BEGIN # find the product of the of the Pythagorian triplet a, b, c where: #
# a + b + c = 1000, a2 + b2 = c2, a < b < c #
INT perimeter = 1000;
INT half perimeter = perimeter OVER 2;
INT max factor := half perimeter;
INT count := 0;
FOR m WHILE m < max factor DO
count +:= 1;
# using Euclid's formula: #
# a = m^2 - n^2, b = 2mn, c = m^2 + n^2 for some integer m, n, so #
# a + b + c = m^2 - n^2 + 2mn + m^2 + n^2 = 2m( m + n ) #
# so m and ( m + n ) are factors of half the perimeter #
IF half perimeter MOD m = 0 THEN
# have a factor of half the perimiter #
INT other factor = half perimeter OVER m;
INT n = other factor - m;
INT m2 = m * m, n2 = n * n;
INT a := IF m > n THEN m2 - n2 ELSE n2 - m2 FI;
INT b := 2 * m * n;
INT c = m2 + n2;
IF ( a + b + c ) = perimeter THEN
# have found the required triple #
IF b < a THEN INT t = a; a := b; b := t FI;
print( ( "a = ", whole( a, 0 ), ", b = ", whole( b, 0 ), ", c = ", whole( c, 0 ) ) );
print( ( "; a * b * c = ", whole( a * b * c, 0 ), newline ) )
FI;
max factor := other factor
FI
OD;
print( ( whole( count, 0 ), " iterations", newline ) )
END- Output:
a = 200, b = 375, c = 425; a * b * c = 31875000 24 iterations
Note that if we stopped after finding the solution, it would be 20 iterations.
ALGOL W[edit]
...but doesn't stop on the first solution (thus verifying there is only one).
% find the Pythagorian triplet a, b, c where a + b + c = 1000 %
for a := 1 until 1000 div 3 do begin
integer a2, b;
a2 := a * a;
for b := a + 1 until 1000 do begin
integer c;
c := 1000 - ( a + b );
if c <= b then goto endB;
if a2 + b*b = c*c then begin
write( i_w := 1, s_w := 0, "a = ", a, ", b = ", b, ", c = ", c );
write( i_w := 1, s_w := 0, "a + b + c = ", a + b + c );
write( i_w := 1, s_w := 0, "a * b * c = ", a * b * c )
end if_a2_plus_b2_e_c2 ;
b := b + 1
end for_b ;
endB:
end for_a .- Output:
a = 200, b = 375, c = 425 a + b + c = 1000 a * b * c = 31875000
AWK[edit]
# syntax: GAWK -f SPECIAL_PYTHAGOREAN_TRIPLET.AWK
# converted from FreeBASIC
BEGIN {
main()
exit(0)
}
function main(a,b,c, limit) {
limit = 1000
for (a=1; a<=limit; a++) {
for (b=a+1; b<=limit; b++) {
for (c=b+1; c<=limit; c++) {
if (a*a + b*b == c*c) {
if (a+b+c == limit) {
printf("%d+%d+%d=%d\n",a,b,c,a+b+c)
printf("%d*%d*%d=%d\n",a,b,c,a*b*c)
return
}
}
}
}
}
}
- Output:
200+375+425=1000 200*375*425=31875000
C++[edit]
#include <cmath>
#include <concepts>
#include <iostream>
#include <numeric>
#include <optional>
#include <tuple>
using namespace std;
optional<tuple<int, int ,int>> FindPerimeterTriplet(int perimeter)
{
unsigned long long perimeterULL = perimeter;
auto max_M = (unsigned long long)sqrt(perimeter/2) + 1;
for(unsigned long long m = 2; m < max_M; ++m)
{
for(unsigned long long n = 1 + m % 2; n < m; n+=2)
{
if(gcd(m,n) != 1)
{
continue;
}
// The formulas below will generate primitive triples if:
// 0 < n < m
// m and n are relatively prime (gcd == 1)
// m + n is odd
auto a = m * m - n * n;
auto b = 2 * m * n;
auto c = m * m + n * n;
auto primitive = a + b + c;
// check all multiples of the primitive at once
auto factor = perimeterULL / primitive;
if(primitive * factor == perimeterULL)
{
// the triplet has been found
if(b<a) swap(a, b);
return tuple{a * factor, b * factor, c * factor};
}
}
}
// the triplet was not found
return nullopt;
}
int main()
{
auto t1 = FindPerimeterTriplet(1000);
if(t1)
{
auto [a, b, c] = *t1;
cout << "[" << a << ", " << b << ", " << c << "]\n";
cout << "a * b * c = " << a * b * c << "\n";
}
else
{
cout << "Perimeter not found\n";
}
}
- Output:
[200, 375, 425] a * b * c = 31875000
Common Lisp[edit]
A conventional solution:
(defun special-triple (sum)
(loop
for a from 1
do (loop
for b from (1+ a)
for c = (- sum a b)
when (< c b) do (return)
when (= (* c c) (+ (* a a) (* b b)))
do (return-from conventional-triple-search (list a b c)))))
This version utilizes SCREAMER which provides constraint solving:
(ql:quickload "screamer")
(in-package :screamer-user)
(defun special-pythagorean-triple (sum)
(let* ((a (an-integer-between 1 (floor sum 3)))
(b (an-integer-between (1+ a) (floor (- sum a) 2)))
(c (- sum a b)))
(when (< c b) (fail))
(if (= (* c c) (+ (* a a) (* b b)))
(list a b c (* a b c))
(fail))))
(print (one-value (special-pythagorean-triple 1000)))
;; => (200 375 425 31875000)
Delphi[edit]
{Define structure to contain triple}
type TTriple = record
A, B, C: integer;
end;
{Make dynamic array of triples}
type TTripleArray = array of TTriple;
procedure PythagorianTriples(Limit: integer; var TA: TTripleArray);
{Find pythagorian Triple up to limit - Return result in list TA}
var Limit2: integer;
var I,J,K: integer;
begin
SetLength(TA,0);
Limit2:=Limit div 2;
for I:=1 to Limit2 do
for J:=I to Limit2 do
for K:=J to Limit do
if ((I+J+K)<Limit) and ((I*I+J*J)=(K*K)) then
begin
SetLength(TA,Length(TA)+1);
TA[High(TA)].A:=I;
TA[High(TA)].B:=J;
TA[High(TA)].C:=K;
end;
end;
procedure ShowPythagoreanTriplet(Memo: TMemo);
var TA: TTripleArray;
var I, Sum, Prod: integer;
begin
{Find triples up to 1100}
PythagorianTriples(1100,TA);
for I:=0 to High(TA) do
begin
{Look for sum of 1000}
Sum:=TA[I].A + TA[I].B + TA[I].C;
if Sum<>1000 then continue;
{Display result}
Prod:=TA[I].A * TA[I].B * TA[I].C;
Memo.Lines.Add(Format('%d + %d + %d = %10.0n',[TA[I].A, TA[I].B, TA[I].C, Sum+0.0]));
Memo.Lines.Add(Format('%d * %d * %d = %10.0n',[TA[I].A, TA[I].B, TA[I].C, Prod+0.0]));
end;
end;
- Output:
200 + 375 + 425 = 1,000 200 * 375 * 425 = 31,875,000 Elapsed Time: 210.001 ms.
F#[edit]
Here I present a solution based on the ideas on the discussion page. It finds all Pythagorean triplets whose elements sum to a given value. It runs in O[n] time. Normally I would exclude triplets with a common factor but for this demonstration I prefer to leave them.
// Special pythagorean triplet. Nigel Galloway: August 31st., 2021
let fG n g=let i=(n-g)/2L in match (n+g)%2L with 0L->if (g*g)%(4L*i)=0L then Some(g,i-(g*g)/(4L*i),i+(g*g)/(4L*i)) else None
|_->if (g*g-2L*i-1L)%(4L*i+2L)=0L then Some(g,i-(g*g)/(4L*i+2L),i+1L+(g*g)/(4L*i+2L)) else None
let E9 n=let fN=fG n in seq{1L..(n-2L)/3L}|>Seq.choose fN|>Seq.iter(fun(n,g,l)->printfn $"%d{n*n}(%d{n})+%d{g*g}(%d{g})=%d{l*l}(%d{l})")
[1L..260L]|>List.iter(fun n->printfn "Sum = %d" n; E9 n)
- Output:
Sum = 1 Sum = 2 Sum = 3 Sum = 4 Sum = 5 Sum = 6 Sum = 7 Sum = 8 Sum = 9 Sum = 10 Sum = 11 Sum = 12 9(3)+16(4)=25(5) Sum = 13 Sum = 14 Sum = 15 Sum = 16 Sum = 17 Sum = 18 Sum = 19 Sum = 20 Sum = 21 Sum = 22 Sum = 23 Sum = 24 36(6)+64(8)=100(10) Sum = 25 Sum = 26 Sum = 27 Sum = 28 Sum = 29 Sum = 30 25(5)+144(12)=169(13) Sum = 31 Sum = 32 Sum = 33 Sum = 34 Sum = 35 Sum = 36 81(9)+144(12)=225(15) Sum = 37 Sum = 38 Sum = 39 Sum = 40 64(8)+225(15)=289(17) Sum = 41 Sum = 42 Sum = 43 Sum = 44 Sum = 45 Sum = 46 Sum = 47 Sum = 48 144(12)+256(16)=400(20) Sum = 49 Sum = 50 Sum = 51 Sum = 52 Sum = 53 Sum = 54 Sum = 55 Sum = 56 49(7)+576(24)=625(25) Sum = 57 Sum = 58 Sum = 59 Sum = 60 100(10)+576(24)=676(26) 225(15)+400(20)=625(25) Sum = 61 Sum = 62 Sum = 63 Sum = 64 Sum = 65 Sum = 66 Sum = 67 Sum = 68 Sum = 69 Sum = 70 400(20)+441(21)=841(29) 441(21)+400(20)=841(29) Sum = 71 Sum = 72 324(18)+576(24)=900(30) Sum = 73 Sum = 74 Sum = 75 Sum = 76 Sum = 77 Sum = 78 Sum = 79 Sum = 80 256(16)+900(30)=1156(34) Sum = 81 Sum = 82 Sum = 83 Sum = 84 144(12)+1225(35)=1369(37) 441(21)+784(28)=1225(35) Sum = 85 Sum = 86 Sum = 87 Sum = 88 Sum = 89 Sum = 90 81(9)+1600(40)=1681(41) 225(15)+1296(36)=1521(39) Sum = 91 Sum = 92 Sum = 93 Sum = 94 Sum = 95 Sum = 96 576(24)+1024(32)=1600(40) Sum = 97 Sum = 98 Sum = 99 Sum = 100 Sum = 101 Sum = 102 Sum = 103 Sum = 104 Sum = 105 Sum = 106 Sum = 107 Sum = 108 729(27)+1296(36)=2025(45) Sum = 109 Sum = 110 Sum = 111 Sum = 112 196(14)+2304(48)=2500(50) Sum = 113 Sum = 114 Sum = 115 Sum = 116 Sum = 117 Sum = 118 Sum = 119 Sum = 120 400(20)+2304(48)=2704(52) 576(24)+2025(45)=2601(51) 900(30)+1600(40)=2500(50) Sum = 121 Sum = 122 Sum = 123 Sum = 124 Sum = 125 Sum = 126 784(28)+2025(45)=2809(53) Sum = 127 Sum = 128 Sum = 129 Sum = 130 Sum = 131 Sum = 132 121(11)+3600(60)=3721(61) 1089(33)+1936(44)=3025(55) Sum = 133 Sum = 134 Sum = 135 Sum = 136 Sum = 137 Sum = 138 Sum = 139 Sum = 140 1600(40)+1764(42)=3364(58) 1764(42)+1600(40)=3364(58) Sum = 141 Sum = 142 Sum = 143 Sum = 144 256(16)+3969(63)=4225(65) 1296(36)+2304(48)=3600(60) Sum = 145 Sum = 146 Sum = 147 Sum = 148 Sum = 149 Sum = 150 625(25)+3600(60)=4225(65) Sum = 151 Sum = 152 Sum = 153 Sum = 154 1089(33)+3136(56)=4225(65) Sum = 155 Sum = 156 1521(39)+2704(52)=4225(65) Sum = 157 Sum = 158 Sum = 159 Sum = 160 1024(32)+3600(60)=4624(68) Sum = 161 Sum = 162 Sum = 163 Sum = 164 Sum = 165 Sum = 166 Sum = 167 Sum = 168 441(21)+5184(72)=5625(75) 576(24)+4900(70)=5476(74) 1764(42)+3136(56)=4900(70) Sum = 169 Sum = 170 Sum = 171 Sum = 172 Sum = 173 Sum = 174 Sum = 175 Sum = 176 2304(48)+3025(55)=5329(73) 3025(55)+2304(48)=5329(73) Sum = 177 Sum = 178 Sum = 179 Sum = 180 324(18)+6400(80)=6724(82) 900(30)+5184(72)=6084(78) 2025(45)+3600(60)=5625(75) Sum = 181 Sum = 182 169(13)+7056(84)=7225(85) Sum = 183 Sum = 184 Sum = 185 Sum = 186 Sum = 187 Sum = 188 Sum = 189 Sum = 190 Sum = 191 Sum = 192 2304(48)+4096(64)=6400(80) Sum = 193 Sum = 194 Sum = 195 Sum = 196 Sum = 197 Sum = 198 1296(36)+5929(77)=7225(85) Sum = 199 Sum = 200 1600(40)+5625(75)=7225(85) Sum = 201 Sum = 202 Sum = 203 Sum = 204 2601(51)+4624(68)=7225(85) Sum = 205 Sum = 206 Sum = 207 Sum = 208 1521(39)+6400(80)=7921(89) Sum = 209 Sum = 210 1225(35)+7056(84)=8281(91) 3600(60)+3969(63)=7569(87) 3969(63)+3600(60)=7569(87) Sum = 211 Sum = 212 Sum = 213 Sum = 214 Sum = 215 Sum = 216 2916(54)+5184(72)=8100(90) Sum = 217 Sum = 218 Sum = 219 Sum = 220 400(20)+9801(99)=10201(101) Sum = 221 Sum = 222 Sum = 223 Sum = 224 784(28)+9216(96)=10000(100) Sum = 225 Sum = 226 Sum = 227 Sum = 228 3249(57)+5776(76)=9025(95) Sum = 229 Sum = 230 Sum = 231 Sum = 232 Sum = 233 Sum = 234 4225(65)+5184(72)=9409(97) 5184(72)+4225(65)=9409(97) Sum = 235 Sum = 236 Sum = 237 Sum = 238 Sum = 239 Sum = 240 225(15)+12544(112)=12769(113) 1600(40)+9216(96)=10816(104) 2304(48)+8100(90)=10404(102) 3600(60)+6400(80)=10000(100) Sum = 241 Sum = 242 Sum = 243 Sum = 244 Sum = 245 Sum = 246 Sum = 247 Sum = 248 Sum = 249 Sum = 250 Sum = 251 Sum = 252 1296(36)+11025(105)=12321(111) 3136(56)+8100(90)=11236(106) 3969(63)+7056(84)=11025(105) Sum = 253 Sum = 254 Sum = 255 Sum = 256 Sum = 257 Sum = 258 Sum = 259 Sum = 260 3600(60)+8281(91)=11881(109)
I present results with timing for increasing powers of 10 to demonstrate its O[n] running.
E9 1000L 40000(200)+140625(375)=180625(425) Real: 00:00:00.001 E9 10000L 4000000(2000)+14062500(3750)=18062500(4250) Real: 00:00:00.000 E9 100000L 400000000(20000)+1406250000(37500)=1806250000(42500) 478515625(21875)+1296000000(36000)=1774515625(42125) Real: 00:00:00.001 E9 1000000L 40000000000(200000)+140625000000(375000)=180625000000(425000) 47851562500(218750)+129600000000(360000)=177451562500(421250) Real: 00:00:00.005 E9 10000000L 54931640625(234375)+23814400000000(4880000)=23869331640625(4885625) 4000000000000(2000000)+14062500000000(3750000)=18062500000000(4250000) 4785156250000(2187500)+12960000000000(3600000)=17745156250000(4212500) Real: 00:00:00.040 E9 100000000L 5493164062500(2343750)+2381440000000000(48800000)=2386933164062500(48856250) 400000000000000(20000000)+1406250000000000(37500000)=1806250000000000(42500000) 478515625000000(21875000)+1296000000000000(36000000)=1774515625000000(42125000) Real: 00:00:00.382 E9 1000000000L 549316406250000(23437500)+238144000000000000(488000000)=238693316406250000(488562500) 40000000000000000(200000000)+140625000000000000(375000000)=180625000000000000(425000000) 47851562500000000(218750000)+129600000000000000(360000000)=177451562500000000(421250000) Real: 00:00:03.704
FreeBASIC[edit]
' version 06-10-2021
' compile with: fbc -s console
' -------------------------------------------------------------
' Brute force
Dim As UInteger a, b, c
Dim As Double t1 = Timer
Print "Brute force" : Print
For a = 1 To 1000
For b = a+1 To 1000
For c = b+1 To 1000
If (a*a+b*b) = (c*c) Then
If a+b+c = 1000 Then
Print a,b,c
End If
End If
Next
Next
Next
Print
Print Using "runtime: ###.## msec.";(Timer-t1)*1000
Print String(60,"-") : Print
' -------------------------------------------------------------
' limit for a = 1000\3 and b = 1000\2, c = 1000 - a - c
Print "Set limits for a and b, c = 1000 - a - b" : Print
t1 = Timer
For a = 3 To 333
For b = a+1 To 500
For c= b+1 To (1000-a-b)
If a+b+c = 1000 Then
If (a*a+b*b) = (c*c) Then
Print a,b,c
Exit For,For,For
End If
End If
Next
Next
Next
Print
Print Using "runtime: ###.## msec.";(Timer-t1)*1000
Print String(60,"-") : Print
' -------------------------------------------------------------
' primative pythagoras triples
' m and n are positive integer and odd
' m > n, m start at 3
' if n > 1 then GCD(m,n) must be 1
' a = m * n
' b = (m ^ 2 - n ^ 2) \ 2
' c = (m ^ 2 + n ^ 2) \ 2
' swap a and b if a > b
Function gcd(x As UInteger, y As UInteger) As UInteger
While y
Dim As UInteger t = y
y = x Mod y
x = t
Wend
Return x
End Function
Function check(n As UInteger) As Boolean
If n And 1 = 0 Then Return FALSE
For i As UInteger = 3 To Sqr(n)
If n Mod i = 0 Then Return TRUE
Next
Return FALSE
End Function
Dim As UInteger m, n, temp
Print "Using primative pythagoras triples" : Print
t1 = Timer
For m = 3 To 201 Step 2
For n = 1 To m Step 2
If n > 1 And (gcd(m, n) <> 1) Then
Continue For
End If
a = m * n
b = (m * m - n * n) \ 2
c = b + n * n
If a > b Then Swap a, b
temp = a + b + c
If 1000 Mod temp = 0 Then
' temp must be odd and split in two number (no prime)
If check(temp) Then
temp = 1000 \ temp
a = a * temp
b = b * temp
c = c * temp
Print a,b,c
Exit For, For
End If
End If
Next
Next
Print
Print Using "runtime: ###.## msec.";(Timer-t1)*1000
Print String(60,"-") : Print
' empty keyboard buffer
While InKey <> "" : Wend
Print : Print "hit any key to end program"
Sleep
End
- Output:
Brute force 200 375 425 runtime: 555.19 msec. ------------------------------------------------------------ Set limits for a and b, c = 1000 - a - b 200 375 425 runtime: 68.91 msec. ------------------------------------------------------------ Using primative pythagoras triples 200 375 425 runtime: 2.49 msec. ------------------------------------------------------------
Go[edit]
package main
import (
"fmt"
"time"
)
func main() {
start := time.Now()
for a := 3; ; a++ {
for b := a + 1; ; b++ {
c := 1000 - a - b
if c <= b {
break
}
if a*a+b*b == c*c {
fmt.Printf("a = %d, b = %d, c = %d\n", a, b, c)
fmt.Println("a + b + c =", a+b+c)
fmt.Println("a * b * c =", a*b*c)
fmt.Println("\nTook", time.Since(start))
return
}
}
}
}
- Output:
a = 200, b = 375, c = 425 a + b + c = 1000 a * b * c = 31875000 Took 77.664µs
jq[edit]
Works with gojq, the Go implementation of jq
range(1;1000) as $a
| range($a+1;1000) as $b
| (1000 - $a - $b) as $c
| select($a*$a + $b*$b == $c*$c)
| {$a, $b, $c, product: ($a*$b*$c)}- Output:
{"a":200,"b":375,"c":425,"product":31875000}
Or, with a tiny bit of thought:
range(1;1000/3) as $a
| range($a+1;1000/2) as $b
| (1000 - $a - $b) as $c
| select($a*$a + $b*$b == $c*$c)
| {$a, $b, $c, product: ($a*$b*$c)}}Julia[edit]
julia> [(a, b, c) for a in 1:1000, b in 1:1000, c in 1:1000 if a < b < c && a + b + c == 1000 && a^2 + b^2 == c^2]
1-element Vector{Tuple{Int64, Int64, Int64}}:
(200, 375, 425)
or, with attention to timing:
julia> @time for a in 1:1000
for b in a+1:1000
c = 1000 - a - b; a^2 + b^2 == c^2 && @show a, b, c
end
end
(a, b, c) = (200, 375, 425)
0.001073 seconds (20 allocations: 752 bytes)
Mathematica / Wolfram Language[edit]
sol = First@
FindInstance[
a + b + c == 1000 && a > 0 && b > 0 && c > 0 &&
a^2 + b^2 == c^2, {a, b, c}, Integers];
Join[List @@@ sol, {{"Product:" , a b c /. sol}}] // TableForm
- Output:
a 200 b 375 c 425 Product: 31875000
Nim[edit]
My solution from Project Euler:
import strformat
from math import floor, sqrt
var
p, s, c : int
r: float
for i in countdown(499, 1):
s = 1000 - i
p = 1000 * (500 - i)
let delta = float(s * s - 4 * p)
r = sqrt(delta)
if floor(r) == r:
c = i
break
echo fmt"Product: {p * c}"
echo fmt"a: {(s - int(r)) div 2}"
echo fmt"b: {(s + int(r)) div 2}"
echo fmt"c: {c}"
- Output:
Product: 31875000 a: 200 b: 375 c: 425
Perl[edit]
use strict;
use warnings;
for my $a (1 .. 998) {
my $a2 = $a**2;
for my $b ($a+1 .. 999) {
my $c = 1000 - $a - $b;
last if $c < $b;
print "$a² + $b² = $c²\n$a + $b + $c = 1000\n" and last if $a2 + $b**2 == $c**2
}
}
- Output:
200² + 375² = 425² 200 + 375 + 425 = 1000
Phix[edit]
strictly adhering to the task description[edit]
brute force (83000 iterations)[edit]
Not that this is in any way slow (0.1s, or 0s with the displays removed), and not that it deliberately avoids using sensible loop limits, you understand.
with javascript_semantics constant n = 1000 integer count = 0 for a=1 to floor(n/3) do for b=a+1 to floor((n-a)/2) do count += 1 integer c = n-(a+b) if a*a+b*b=c*c then printf(1,"a=%d, b=%d, c=%d, a*b*c=%d\n",{a,b,c,a*b*c}) end if end for end for printf(1,"%d iterations\n",count)
- Output:
a=200, b=375, c=425, a*b*c=31875000 83000 iterations
smarter (166 iterations)[edit]
It would of course be 100 iterations if we quit once found (whereas the above would be 69775).
with javascript_semantics constant n = 1000 integer count = 0 for a=2 to floor(n/3) by 2 do count += 1 integer nn2a = n*(n/2-a), na = n-a if remainder(nn2a,na)=0 then integer b = nn2a/na, c = n-(a+b) printf(1,"a=%d, b=%d, c=%d, a*b*c=%d\n",{a,b,c,a*b*c}) end if end for printf(1,"%d iterations\n",count)
- Output:
a=200, b=375, c=425, a*b*c=31875000 166 iterations
PL/M[edit]
Based on the XPL0 solution.
As the original 8080 PL/M compiler only has unsigned 8 and 16 bit integer arithmetic, the PL/M long multiplication routines and also a square root routine based that in the PL/M sample for the Frobenius Numbers task are used - which makes this somewhat longer than it would otherwose be...
100H: /* FIND THE PYTHAGOREAN TRIPLET A, B, C WHERE A + B + C = 1000 */
/* CP/M BDOS SYSTEM CALL */
BDOS: PROCEDURE( FN, ARG ); DECLARE FN BYTE, ARG ADDRESS; GOTO 5; END;
/* I/O ROUTINES */
PRINT$CHAR: PROCEDURE( C ); DECLARE C BYTE; CALL BDOS( 2, C ); END;
PRINT$STRING: PROCEDURE( S ); DECLARE S ADDRESS; CALL BDOS( 9, S ); END;
PRINT$NL: PROCEDURE; CALL PRINT$STRING( .( 0DH, 0AH, '$' ) ); END;
/* LONG MULTIPLICATION */
/* LARGE INTEGERS ARE REPRESENTED BY ARRAYS OF BYTES WHOSE VALUES ARE */
/* A SINGLE DECIMAL DIGIT OF THE NUMBER */
/* THE LEAST SIGNIFICANT DIGIT OF THE LARGE INTEGER IS IN ELEMENT 1 */
/* ELEMENT 0 CONTAINS THE NUMBER OF DIGITS THE NUMBER HAS */
DECLARE LONG$INTEGER LITERALLY '(21)BYTE';
DECLARE DIGIT$BASE LITERALLY '10';
/* PRINTS A LONG INTEGER */
PRINT$LONG$INTEGER: PROCEDURE( N$PTR );
DECLARE N$PTR ADDRESS;
DECLARE N BASED N$PTR LONG$INTEGER;
DECLARE ( D, F ) BYTE;
F = N( 0 );
DO D = 1 TO N( 0 );
CALL PRINT$CHAR( N( F ) + '0' );
F = F - 1;
END;
END PRINT$LONG$INTEGER;
/* SETS A LONG$INTEGER TO A 16-BIT VALUE */
SET$LONG$INTEGER: PROCEDURE( LN, N );
DECLARE ( LN, N ) ADDRESS;
DECLARE V ADDRESS;
DECLARE LN$PTR ADDRESS, LN$BYTE BASED LN$PTR BYTE;
DECLARE LN$0 ADDRESS, LN$0BYTE BASED LN$0 BYTE;
LN$0, LN$PTR = LN;
LN$0BYTE = 1;
LN$PTR = LN$PTR + 1;
LN$BYTE = ( V := N ) MOD DIGIT$BASE;
DO WHILE( ( V := V / DIGIT$BASE ) > 0 );
LN$PTR = LN$PTR + 1;
LN$BYTE = V MOD DIGIT$BASE;
LN$0BYTE = LN$0BYTE + 1;
END;
END SET$LONG$INTEGER;
/* IMPLEMENTS LONG MULTIPLICATION, C IS SET TO A * B */
/* C CAN BE THE SAME LONG$INTEGER AS A OR B */
LONG$MULTIPLY: PROCEDURE( A$PTR, B$PTR, C$PTR );
DECLARE ( A$PTR, B$PTR, C$PTR ) ADDRESS;
DECLARE ( A BASED A$PTR, B BASED B$PTR, C BASED C$PTR ) LONG$INTEGER;
DECLARE MRESULT LONG$INTEGER;
DECLARE RPOS BYTE;
/* MULTIPLIES THE LONG INTEGER IN B BY THE INTEGER A, THE RESULT */
/* IS ADDED TO C, STARTING FROM DIGIT START */
/* OVERFLOW IS IGNORED */
MULTIPLY$ELEMENT: PROCEDURE( A, B$PTR, C$PTR, START );
DECLARE ( B$PTR, C$PTR ) ADDRESS;
DECLARE ( A, START ) BYTE;
DECLARE ( B BASED B$PTR, C BASED C$PTR ) LONG$INTEGER;
DECLARE ( CDIGIT, D$CARRY, BPOS, CPOS ) BYTE;
D$CARRY = 0;
CPOS = START;
DO BPOS = 1 TO B( 0 );
CDIGIT = C( CPOS ) + ( A * B( BPOS ) ) + D$CARRY;
IF CDIGIT < DIGIT$BASE THEN D$CARRY = 0;
ELSE DO;
/* HAVE DIGITS TO CARRY */
D$CARRY = CDIGIT / DIGIT$BASE;
CDIGIT = CDIGIT MOD DIGIT$BASE;
END;
C( CPOS ) = CDIGIT;
CPOS = CPOS + 1;
END;
C( CPOS ) = D$CARRY;
/* REMOVE LEADING ZEROS BUT IF THE NUMBER IS 0, KEEP THE FINAL 0 */
DO WHILE( CPOS > 1 AND C( CPOS ) = 0 );
CPOS = CPOS - 1;
END;
C( 0 ) = CPOS;
END MULTIPLY$ELEMENT ;
/* THE RESULT WILL BE COMPUTED IN MRESULT, ALLOWING A OR B TO BE C */
DO RPOS = 1 TO LAST( MRESULT ); MRESULT( RPOS ) = 0; END;
/* MULTIPLY BY EACH DIGIT AND ADD TO THE RESULT */
DO RPOS = 1 TO A( 0 );
IF A( RPOS ) <> 0 THEN DO;
CALL MULTIPLY$ELEMENT( A( RPOS ), B$PTR, .MRESULT, RPOS );
END;
END;
/* RETURN THE RESULT IN C */
DO RPOS = 0 TO MRESULT( 0 ); C( RPOS ) = MRESULT( RPOS ); END;
END;
/* INTEGER SUARE ROOT: BASED ON THE ONE IN THE PL/M FOR FROBENIUS NUMBERS */
SQRT: PROCEDURE( N )ADDRESS;
DECLARE ( N, X0, X1 ) ADDRESS;
IF N <= 3 THEN DO;
IF N = 0 THEN X0 = 0; ELSE X0 = 1;
END;
ELSE DO;
X0 = SHR( N, 1 );
DO WHILE( ( X1 := SHR( X0 + ( N / X0 ), 1 ) ) < X0 );
X0 = X1;
END;
END;
RETURN X0;
END SQRT;
/* FIND THE PYTHAGORIAN TRIPLET */
DECLARE ( A, B, C, M, N, M2, N2, SQRT$1000 ) ADDRESS;
DECLARE ( LA, LB, LC, ABC ) LONG$INTEGER;
SQRT$1000 = SQRT( 1000 );
DO N = 1 TO SQRT$1000; /* M AND N MUST HAD DIFFERENT PARITY, */
DO M = N + 1 TO SQRT$1000 BY 2; /* I.E. ONE ODD, ONE EVEN */
/* NOTE: A = M2 - N2, B = 2MN, C = M2 + N2 */
/* A + B + C = M2 - N2 + 2MN + M2 + N2 = 2( M2 + MN ) = 2M( M + N )*/
IF ( M * ( M + N ) ) = 500 THEN DO;
M2 = M * M;
N2 = N * N;
CALL SET$LONG$INTEGER( .A, M2 - N2 );
CALL SET$LONG$INTEGER( .B, 2 * M * N );
CALL SET$LONG$INTEGER( .C, M2 + N2 );
CALL LONG$MULTIPLY( .A, .B, .ABC );
CALL LONG$MULTIPLY( .ABC, .C, .ABC );
CALL PRINT$LONG$INTEGER( .ABC );
CALL PRINT$NL;
END;
END;
END;
EOF- Output:
31875000
Python[edit]
Python 3.8.8 (default, Apr 13 2021, 15:08:03) Type "help", "copyright", "credits" or "license" for more information. >>> [(a, b, c) for a in range(1, 1000) for b in range(a, 1000) for c in range(1000) if a + b + c == 1000 and a*a + b*b == c*c] [(200, 375, 425)]
Raku[edit]
hyper for 1..998 -> $a {
my $a2 = $a²;
for $a + 1 .. 999 -> $b {
my $c = 1000 - $a - $b;
last if $c < $b;
say "$a² + $b² = $c²\n$a + $b + $c = {$a+$b+$c}\n$a × $b × $c = {$a×$b×$c}"
and exit if $a2 + $b² == $c²
}
}
- Output:
200² + 375² = 425² 200 + 375 + 425 = 1000 200 × 375 × 425 = 31875000
REXX[edit]
Some optimizations were done such as pre-computing the squares of all possible A's, B's, and C's.
Also, there were multiple shortcuts to limit an otherwise exhaustive search; Once a sum or a square was too big,
the next integer was used (for the previous DO loop).
/*REXX pgm computes integers A, B, C that solve: 0<A<B<C; A+B+C = 1000; A^2+B^2 = C^2 */
parse arg sum hi n . /*obtain optional argument from the CL.*/
if sum=='' | sum=="," then sum= 1000 /*Not specified? Then use the default.*/
if hi=='' | hi=="," then hi= 1000 /* " " " " " " */
if n=='' | n=="," then n= 1 /* " " " " " " */
hh= hi - 2 /*N: number of solutions to find/show.*/
do j=1 for hi; @.j= j*j /*pre─compute squares ──► HI, inclusive*/
end /*j*/
#= 0; pad= left('', 9) /*#: the number of solutions found. */
do a=2 for hh%2 by 2; aa= @.a /*search for solutions to the equations*/
do b=a+1; ab= a + b /*compute the sum of 2 numbers (A & B).*/
if ab>hh then iterate a /*Sum of A+B>HI? Then stop with B's */
aabb= aa + @.b /*compute the sum of: A^2 + B^2 */
do c=b+1 while @.c <= aabb /*test integers that satisfy equations.*/
if @.c\==aabb then iterate /* " \=A^2+B^2? Then keep searching*/
abc= ab + c /*compute the sum of: A + B + C */
if abc > sum then iterate b /*Is A+B+C > SUM? Then stop with C's.*/
if abc == sum then call show /*Does " = SUM? Then solution found*/
end /*c*/
end /*b*/
end /*a*/
done: if #==0 then #= 'no'; say pad pad pad # ' solution's(#) "found."
exit 0 /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
s: if arg(1)==1 then return arg(3); return word(arg(2) 's', 1) /*simple pluralizer*/
show: #= #+1; say pad 'a=' a pad "b=" b pad 'c=' c; if #>=n then signal done; return
- output when using the default inputs:
a= 200 b= 375 c= 425
1 solution found.
Ring[edit]
Various algorithms presented, some are quite fast. Timings are from Tio.run. On the desktop of a core i7-7700 @ 3.60Ghz, it goes about 6.5 times faster.
tf = 1000 # time factor adjustment for different Ring versions
? "working..."
? "turtle method:" # 3 nested loops is not a good approach,
# even when cheating by cherry-picking the loop start/end points...
st = clock()
for a = 100 to 400
for b = 200 to 800
for c = b to 1000 - a - b
if a + b + c = 1000
if a * a + b * b = c * c exit 3 ok
ok
next
next
next
et = clock() - st
? "a = " + a + " b = " + b + " c = " + c
? "Elapsed time = " + et / (tf * 1000) + " s" + nl
? "brutally forced method:" # eliminating the "c" loop speeds it up a bit
st = clock()
for a = 1 to 1000
for b = a to 1000
c = 1000 - a - b
if a * a + b * b = c * c exit 2 ok
next
next
et = clock() - st
? "a = " + a + " b = " + b + " c = " + c
? "Elapsed time = " + et / tf + " ms" + nl
# some basic info about this task:
p = 1000 pp = p * p >> 1 # perimeter, half of perimeter^2
maxc = ceil(sqrt(pp) * 2) - p # minimum c = 415 ceiling of 414.2135623730950488016887242097
maxa = (p - maxc) >> 1 # maximum a = 292, shorter leg will shrink
minb = p - maxc - maxa # minimum b = 293, longer leg will lengthen
? "polished brute force method:" # calculated realistic limits for the loops,
# cached some vars that didn't need recalcs over and over
st = clock()
minb = maxa + 1 maxb = p - maxc pma = p - 1
for a = 1 to maxa
aa = a * a
c = pma - minb
for b = minb to maxb
if aa + b * b = c * c exit 2 ok
c--
next
pma--
next
et = clock() - st
? "a = " + a + " b = " + b + " c = " + c
? "Elapsed time = " + et / tf + " ms" + nl
? "quick method:" # down to one loop, using some math insight
st = clock()
n2 = p >> 1
for a = 1 to n2
b = p * (n2 - a)
if b % (p - a) = 0 exit ok
next
b /= (p - a)
? "a = " + a + " b = " + b + " c = " + (p - a - b)
et = clock() - st
? "Elapsed time = " + et / tf + " ms" + nl
? "even quicker method:" # generate primitive Pythagorean triples,
# then scale to fit the actual perimeter
st = clock()
md = 1 ms = 1
for m = 1 to 4
nd = md + 2 ns = ms + nd
for n = m + 1 to 5
if p % (((n * m) + ns) << 1) = 0 exit 2 ok
nd += 2 ns += nd
next
md += 2 ms += md
next
et = clock() - st
a = n * m << 1 b = ns - ms c = ns + ms d = p / (((n * m) + ns) << 1)
? "a = " + a * d + " b = " + b * d + " c = " + c * d
? "Elapsed time = " + et / tf + " ms" + nl
? "alternate method:" # only uses addition / subtraction inside the loop.
# makes a guess, then tweaks the guess until correct
st = clock()
a = maxa b = minb
g = p * (a + b) - a * b # guess
while g != pp
if pp > g
b++ g += p - a # step "b" only when the "a" step went too far
ok
a-- g -= p - b # step "a" on every iteration
end
et = clock() - st
? "a = " + a + " b = " + b + " c = " + (p - a - b)
? "Elapsed time = " + et / tf + " ms" + nl
see "done..."- Output:
working... turtle method: a = 200 b = 375 c = 425 Elapsed time = 13.36 s brutally forced method: a = 200 b = 375 c = 425 Elapsed time = 983.94 ms polished brute force method: a = 200 b = 375 c = 425 Elapsed time = 216.66 ms quick method: a = 200 b = 375 c = 425 Elapsed time = 0.97 ms even quicker method: a = 200 b = 375 c = 425 Elapsed time = 0.18 ms alternate method: a = 200 b = 375 c = 425 Elapsed time = 0.44 ms done...
Wren[edit]
Very simple approach, only takes 0.013 seconds even in Wren.
var a = 3
while (true) {
var b = a + 1
while (true) {
var c = 1000 - a - b
if (c <= b) break
if (a*a + b*b == c*c) {
System.print("a = %(a), b = %(b), c = %(c)")
System.print("a + b + c = %(a + b + c)")
System.print("a * b * c = %(a * b * c)")
return
}
b = b + 1
}
a = a + 1
}
- Output:
a = 200, b = 375, c = 425 a + b + c = 1000 a * b * c = 31875000
Incidentally, even though we are told there is only one solution, it is almost as quick to verify this by observing that, since a < b < c, the maximum value of a must be such that 3a + 2 = 1000 or max(a) = 332. The following version ran in 0.015 seconds and, of course, produced the same output:
for (a in 3..332) {
var b = a + 1
while (true) {
var c = 1000 - a - b
if (c <= b) break
if (a*a + b*b == c*c) {
System.print("a = %(a), b = %(b), c = %(c)")
System.print("a + b + c = %(a + b + c)")
System.print("a * b * c = %(a * b * c)")
}
b = b + 1
}
}
XPL0[edit]
int N, M, A, B, C;
for N:= 1 to sqrt(1000) do
for M:= N+1 to sqrt(1000) do
[A:= M*M - N*N; \Euclid's formula
B:= 2*M*N;
C:= M*M + N*N;
if A+B+C = 1000 then
IntOut(0, A*B*C);
]- Output:
31875000