(41C) Pythagorean Triples

[SlideRule]

An extract from PROGRAMMABLE CALCULATORS Implications for the Mathematics Curriculum, Clearinghouse for Science, Mathematics and Environmental Education, Ohio State University, December 1980

"Independent Study with a Programmable Calculator … pg. 49
… The following material is an excerpt from Mathematical Recreations for the Programmable Calculator (Mohler and Hoffman, 1981) a collection of programming problems designed to teach the standard techniques of programming …
Pythagorean Triples … pg. 50
… Write a program which searches out and finds all Pythagorean triples …
PYTR Program for the HP-41C … pg. 77
…

Parts of the extract are difficult to read BUT the HP-41C program listing is only slightly fuzzy near the bottom of the listing.

BEST!
SlideRule

[Thomas Klemm]

On page 26 we find:

Quote:Our calculators used Reverse Polish Notation, which is a great convenience and posed no difficulty to the children.

No one in this forum is surprised by that.

---

I picked the problem to print primitive Pythagorean triples.

Here's a program for the HP-42S:

Code:

00 { 99-Byte Prgm } #
01▸LBL "PPT"        # p
02 1ᴇ3              # 1000 p
03 ÷                # 0.ppp
04 2                # 2 0.ppp
05 +                # i=2.ppp
06 STO 00           # i -> R00
07▸LBL 00           # for i in 2.ppp
08 RCL 00           # i
09 IP               # m
10 STO 02           # m -> R02
11 2                # 2 m
12 MOD              # m%2
13 1                # 1 m%2
14 +                # b=m%2+1
15 RCL 02           # m b
16 1                # 1 m b
17 -                # e=m-1 b
18 .02              # 0.02 e b
19 +                # e.02 b
20 1ᴇ3              # 1000 e.02 b
21 ÷                # 0.eee02
22 +                # j=b.eee02
23 STO 01           # j -> R01
24▸LBL 01           # for j in b.eee02
25 RCL 02           # m
26 RCL 01           # j m
27 IP               # n m
28 STO 03           # n -> R03
29▸LBL 02           # ( m n -- gdc(m, n) )
30 X<>Y             # m n
31 RCL ST Y         # n m n
32 MOD              # m%n n
33 X≠0?             #
34 GTO 02           # n' m'
35 R↓               # gcd(m, n)
36 1                # 1 gcd(m, n)
37 X≠Y?             # gcd(m, n) ≠ 1 ?
38 GTO 03           # skip triple
39 RCL 03           # n
40 X↑2              # n^2
41 RCL 02           # m n^2
42 RCL× 03          # m*n n^2
43 STO+ ST X        # 2*m*n n^2
44 RCL 02           # m 2*m*n n^2
45 X↑2              # m^2 2*m*n n^2
46 ENTER            # m^2 m^2 2*m*n n^2
47 R↑               # n^2 m^2 m^2 2*m*n
48 STO- ST Z        # n^2 m^2 m^2-n^2 2*m*n
49 +                # m^2+n^2 m^2-n^2 2*m*n
50 CLA              # ""
51 ARCL ST Y        # "Y"
52 ├" "             # "Y "
53 ARCL ST Z        # "Y Z"
54 ├" "             # "Y Z "
55 ARCL ST X        # "Y Z X"
56 AVIEW            #                                 
57▸LBL 03           # resume
58 ISG 01           # j -> j+1
59 GTO 01           #
60 ISG 00           # i -> i+1
61 GTO 00           #
62 END              #

But it should also work with the HP-41C after the obvious transformations.

Example

12 XEQ "PPT"

3 4 5
5 12 13
15 8 17
7 24 25
21 20 29
9 40 41
35 12 37
11 60 61
45 28 53
33 56 65
13 84 85
63 16 65
55 48 73
39 80 89
15 112 113
77 36 85
65 72 97
17 144 145
99 20 101
91 60 109
51 140 149
19 180 181
117 44 125
105 88 137
85 132 157
57 176 185
21 220 221
143 24 145
119 120 169
95 168 193
23 264 265

Make sure to SF 21 and hit R/S to get to the next triple.

---

The program is based on the following Python program:

Code:

def primitive(m):
    for n in range(2 if m % 2 else 1, m, 2):
        if gcd(m, n) == 1:
            triple(m, n)

def gcd(a, b):
    while b > 0:
        a, b = b, a % b
    return a

def triple(m, n):
    a = m*m - n*n
    b = 2*m*n
    c = m*m + n*n
    print(f"{a} {b} {c}")

for m in range(2, 13):
    primitive(m)

---

References

• Pythagorean triple
  

[C.Ret]

(06-25-2022 03:40 PM)Thomas Klemm Wrote:  But it should also work with the HP-41C after the obvious transformations.

As usual, Thomas Klemm produces an excellent algorithm. I particularly like how the GCD determination is efficiently incorporated into the code.

So I grabbed my HP-41C and its trusty 82240A to make the obvious transformations needed to aim and print:

01 LBL"PYTH"
02 STO 04  XEQ 04
04 LBL 01
05  RCL 04  STO 05
07  LBL 00
08   DSE 05  XEQ 02  DSE 05  GTO 00
12  DSE 04  GTO 01
14 GTO 04

15 LBL 02
16  RCL 04  RCL 05
18  LBL 03
19   STO Z  MOD  X≠0?  GTO 03
23  10^X  X≠Y?  RTN
26  +  CLA
28  RCL 04  ARCL X  "┝,"  ST* Y  X^2  STO Z
34  RCL 05  ARCL X  ACA   ST* Z  X^2  ST+ T
40  -  FMT  ACX  RDN  ACX  RDN  ACX
47  LBL 04
48   ADV
49 END

   

Shorter than the original code, this one just displays triples in a different order.

[Thomas Klemm]

With obvious transformations, I thought more of something like: STO+ ST X \(\mapsto\) ST+ X.

Thanks for improving my program.

Using DSE instead of ISG avoids the calculation of the index.
Also by using it twice in a row we still use a step size of 2:

Code:

07  LBL 00
08   DSE 05  XEQ 02  DSE 05  GTO 00

I should have remembered that with STO Z we also have tuck ~ swap over:

Code:

18  LBL 03
19   STO Z  MOD  X≠0?  GTO 03

I would probably rather use SIGN instead of 10^X to map \(0 \mapsto 1\) but its nice to add both \(1 + 1 = 2\) and use the result as a factor when calculating \(2 \cdot m \cdot n\) in the following steps:

Code:

23  10^X  X≠Y?  RTN
26  +  CLA

I wasn't aware of FMT or probably just forgot about it.
It makes for a nice listing:

Code:

40  -  FMT  ACX  RDN  ACX  RDN  ACX

Well done!

[Ángel Martin]

very timely: https://www.youtube.com/watch?v=QJYmyhnaaek

Great visual content for math, love that channel/

[John Keith]

(06-29-2022 01:02 PM)Ángel Martin Wrote:  very timely: https://www.youtube.com/watch?v=QJYmyhnaaek

Great visual content for math, love that channel/

That was a great video, thanks for posting!

Additionally, Berggren's method does not require GCD and is pretty fast and simple on any calculator that can handle matrices. It generates all primitive triples but in a different order than the complex squaring method.

[C.Ret]

(07-01-2022 11:55 AM)John Keith Wrote:  Additionally, Berggren's method does not require GCD and is pretty fast and simple on any calculator that can handle matrices. It generates all primitive triples but in a different order than the complex squaring method.

Thanks for pointing out this method which easily produces tons of primary Pythagorean triples using a very simple recursive program on an advanced calculator natively manipulating vectors and matrices:

PYT:
« 1 +
 [[ 1 2 2 ][ 2 1 2 ][ 2 2 3]] → T n M          \(T=\left[a_0,b_0,c_0 \right]\)
 «
   "             " 1 n 2 * SUB T →STR + PR1 STR→
   IF n N <
   THEN
      M T 2 DUP2 GET NEG PUT * n PYT      \(\left[a_1,b_1,c_1\right]=\begin{bmatrix}1&-2&2\\2&-1&2\\2&-2&3\\\end{bmatrix}\times \left [ a_0,b_0,c_0 \right ]=\begin{bmatrix}1&2&2\\2&1&2\\2&2&3\\\end{bmatrix}\times \left [ a_0,-b_0,c_0 \right ]\)
      M T             * n PYT      \(\left[a_2,b_2,c_2\right]=\begin{bmatrix}1&2&2\\2&1&2\\2&2&3\\\end{bmatrix}\times \left [ a_0,b_0,c_0 \right ]\)
      M T 1 DUP2 GET NEG PUT * n PYT      \(\left[a_3,b_3,c_3\right]=\begin{bmatrix}-1&2&2\\-2&1&2\\-2&2&3\\\end{bmatrix}\times \left [ a_0,b_0,c_0 \right ]=\begin{bmatrix}1&2&2\\2&1&2\\2&2&3\\\end{bmatrix}\times \left [ -a_0,b_0,c_0 \right ]\)
   END » »

Store max depth into N register:
4 'N' STO

Set printer online, aim to it and print by typing:
[ 3 4 5 ] 0 PYT

   

[Thomas Klemm]

(07-03-2022 08:18 AM)C.Ret Wrote:  using a very simple recursive program on an advanced calculator natively manipulating vectors and matrices

Do you mind posting the program?

Nice receipt by the way.
I'd rather have it than the ones I usually get.

[C.Ret]

(07-03-2022 10:46 AM)Thomas Klemm Wrote:  

(07-03-2022 08:18 AM)C.Ret Wrote:  using a very simple recursive program on an advanced calculator natively manipulating vectors and matrices

Do you mind posting the program?

Sorry, as I posted this morning, I was so caught up in my attempt to come up with a version for HP-41, that I completely forgot to post the code for HP-28S.

Tonight, after a day full of other Sunday activities, my version for HP-41 is still not finalized. I believe I will revise the way I support recurring calls as well as matrix products in my HP-41C.

[John Keith]

(07-03-2022 08:18 AM)C.Ret Wrote:  PYT:
« 1 +
 [[ 1 2 2 ][ 2 1 2 ][ 2 2 3]] → T n M
 «
   " " 1 n 2 * SUB T →STR + PR1 STR→
   IF n N <
   THEN
      M T 2 DUP2 GET NEG PUT * n PYT
      M T             * n PYT
      M T 1 DUP2 GET NEG PUT * n PYT
   END » »

Store max depth into N register:
4 'N' STO

Set printer online, aim to it and print by typing:
[ 3 4 5 ] 0 PYT

That's a great program! I'm glad you posted it before I posted my sad attempt.