Special pythagorean triplet
- Task
- The following problem is taken from Project Euler
- Related task
ALGOL 68
...but doesn't stop on the first solution (thus verifying there is only one). <lang algol68># find the Pythagorian triplet a, b, c where a + b + c = 1000, a < b < c # FOR a TO 1000 OVER 3 DO
INT a2 = a * a;
FOR b FROM a + 1 TO 1000 WHILE INT a plus b = a + b;
INT c = 1000 - a plus b;
c > b
DO
IF a2 + b*b = c*c THEN
print( ( "a = ", whole( a, 0 ), ", b = ", whole( b, 0 ), ", c = ", whole( c, 0 ), newline ) );
print( ( "a + b + c = ", whole( a plus b + c, 0 ), newline ) );
print( ( "a * b * c = ", whole( a * b * c, 0 ), newline ) )
FI
OD
OD</lang>
- Output:
a = 200, b = 375, c = 425 a + b + c = 1000 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)
PL/M
As the original 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;
/* PL/M ONLY DOES UNSIGNED ARITHMETIC SO WE USE THE EQUATIONS FOR */
/* A, B AND C: A = M2 - N2, B = 2MN, C = M2 + N2 TO CALCULATE */
/* A + B + C = M2 - N2 + 2MN + M2 + N2 = 2( M2 + MN ) = 2M( M + N )*/
IF ( 2 * M * ( M + N ) ) = 1000 AND M > N THEN DO;
M2 = M * M;
N2 = N * N;
A = M2 - N2;
B = 2 * M * N;
C = M2 + N2;
CALL SET$LONG$INTEGER( .LA, A );
CALL SET$LONG$INTEGER( .LB, B );
CALL SET$LONG$INTEGER( .LC, C );
CALL LONG$MULTIPLY( .LA, .LB, .ABC );
CALL LONG$MULTIPLY( .ABC, .LC, .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 = 1000" and exit if $a2 + $b² == $c²
}
}</lang>
- Output:
200² + 375² = 425² 200 + 375 + 425 = 1000
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 /* " " " " " " */
do j=1 for hi; @.j= j*j /*precompute squares up to HI. */
end /*j*/
hi2= hi-2
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 leave a /*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 ' a=' a " b=" b ' c=' c exit 0 /*stick a fork in it, we're all done. */</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