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 ...but doesn't stop on the first solution (thus verifying there is only one). <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 OVER 2 DO FOR m FROM n + 1 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)
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