Special pythagorean triplet: Difference between revisions
(→{{header|Common Lisp}}: Added Common Lisp version) |
(→{{header|Common Lisp}}: Added more conventional implementation.) |
||
| Line 76: | Line 76: | ||
=={{header|Common Lisp}}== |
=={{header|Common Lisp}}== |
||
| ⚫ | |||
A conventional solution: |
|||
<lang lisp> |
|||
(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))))) |
|||
</lang> |
|||
| ⚫ | |||
<lang lisp> |
<lang lisp> |
||
Revision as of 05:06, 30 November 2021
- Task
- The following problem is taken from Project Euler problem 9
- Related task
ALGOL 68
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.
<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 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</lang>
- 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
...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
Common Lisp
A conventional solution:
<lang lisp> (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)))))
</lang>
This version utilizes SCREAMER which provides constraint solving:
<lang lisp> (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)
</lang>
F#
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. <lang fsharp> // 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) </lang>
- 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
<lang freebasic>' 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</lang>
- 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
<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
jq
Works with gojq, the Go implementation of jq <lang 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)}</lang>
- Output:
{"a":200,"b":375,"c":425,"product":31875000}
Or, with a tiny bit of thought: <lang jq>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)}}</lang>
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
Perl
<lang perl>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
}
}</lang>
- Output:
200² + 375² = 425² 200 + 375 + 425 = 1000
Phix
strictly adhering to the task description
brute force (83000 iterations)
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)
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
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; /* 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</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 (for the previous DO loop).
<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 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</lang>
- output when using the default inputs:
a= 200 b= 375 c= 425
1 solution found.
Ring
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. <lang ring>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..."</lang>
- 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
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