Special pythagorean triplet
- Task
- The following problem is taken from Project Euler
- Related task
ALGOL 68
Using Euclid's formula, as in the XPL0 sample. <lang algol68>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 sqrt 1000 = ENTIER sqrt( 1000 );
FOR n TO sqrt 1000 DO # m and must have different parity, i.e. one even, one odd #
FOR m FROM n + 1 BY 2 TO sqrt 1000 DO
# a = m^2 - n^2, b = 2mn, c = m^2 + n^2 ( Euclid's formula ), so #
# a + b + c = m^2 - n^2 + 2mn + m^2 + n^2 = 2( m^2 + mn ) = 2m( m + n ) #
IF m * ( m + n ) = 500 THEN
INT m2 = m * m, n2 = n * n;
INT a = m2 - n2;
INT b = 2 * m * n;
INT c = m2 + n2;
print( ( "a = ", whole( a, 0 ), ", b = ", whole( b, 0 ), ", c = ", whole( c, 0 ), newline ) );
print( ( "a * b * c = ", whole( a * b * c, 0 ), newline ) )
FI
OD
OD
END</lang>
- Output:
a = 375, b = 200, c = 425 a * b * c = 31875000
ALGOL W
...but doesn't stop on the first solution (thus verifying there is only one). <lang algolw>% 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 .</lang>
- Output:
a = 200, b = 375, c = 425 a + b + c = 1000 a * b * c = 31875000
F#
<lang fsharp> // Special pythagorean triplet. Nigel Galloway: August 31st., 2021 let rec fN i g e l=if g=i then None else match compare(e+g*g)(l*l) with 0->Some(i,g,l) |1->fN i (g-1) e (l+1) |_->None seq{1..332}|>Seq.choose(fun n->let g=(999-n)/2 in fN n g (n*n) (1000-n-g))|>Seq.iter(fun(n,g,l)->printfn $"%d{n*n}(%d{n})+%d{g*g}(%d{g})=%d{l*l}(%d{l})") </lang>
- Output:
40000(200)+140625(375)=180625(425)
Go
<lang go>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
}
}
}
}</lang>
- Output:
a = 200, b = 375, c = 425 a + b + c = 1000 a * b * c = 31875000 Took 77.664µs
Julia
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)
Nim
My solution from Project Euler:
<lang Nim>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}"</lang>
- Output:
Product: 31875000 a: 200 b: 375 c: 425
Phix
brute force (83000 iterations)
Not that this is in any way slow (0.1s, or 0s with the displays removed), or 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)
It would of course be 100 iterations if we quit/exit once found.
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
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...
<lang pli>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;
DO M = N + 1 TO SQRT$1000;
/* 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</lang>
- Output:
31875000
Python
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
<lang perl6>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²
}
}</lang>
- Output:
200² + 375² = 425² 200 + 375 + 425 = 1000 200 × 375 × 425 = 31875000
REXX
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.
<lang rexx>/*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 s hi . /*obtain optional argument from the CL.*/
if s== | s=="," then s= 1000 /*Not specified? Then use the default.*/
if hi== | hi=="," then hi= 1000 /* " " " " " " */
hi2= hi-2
do j=1 for hi; @.j= j*j /*precompute squares up to HI. */
end /*j*/
- = 0 /*#: the number of solutions found. */
do a=1 for hi2; aa= @.a /*go hunting for solutions to equations*/
do b=a+1 for hi2-a; ab= a + b /*calculate sum of two (A,B) squares.*/
if ab>hi 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 for hi2-b /*test integers that satisfy equations.*/
if @.c > aabb then iterate b /*Square> A^2+B^2? Then stop with C's.*/
if @.c \== aabb then iterate /*Square\=A^2+B^2? Then keep searching*/
abc= ab + c /*compute the sum of A + B + C */
if abc == s then call show /*Does A+B+C = S? Then solution found*/
if abc > s then iterate b /* " " > S? Then stop with C's.*/
end /*c*/
end /*b*/
end /*a*/
say # ' solutions found.' exit 0 /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ show: #= # + 1; say ' a=' a " b=" b ' c=' c; exit 0
/*replace EXIT 0 with RETURN to find more possible solutions──►──┘ */</lang>
- output when using the default inputs:
a= 200 b= 375 c= 425
Ring
<lang ring> load "stdlib.ring" see "working..." + nl
time1 = clock() for a = 1 to 1000
for b = a+1 to 1000
for c = b+1 to 1000
if (pow(a,2)+pow(b,2)=pow(c,2))
if a+b+c = 1000
see "a = " + a + " b = " + b + " c = " + c + nl
exit 3
ok
ok
next
next
next
time2 = clock() time3 = time2/1000 - time1/1000 see "Elapsed time = " + time3 + " s" + nl see "done..." + nl </lang>
- Output:
working... a = 200 b = 375 c = 425 Elapsed time = 497.61 s done...
Wren
Very simple approach, only takes 0.013 seconds even in Wren. <lang ecmascript>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
}</lang>
- 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:
<lang ecmascript>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
}
}</lang>
XPL0
<lang XPL0>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);
]</lang>
- Output:
31875000