Area of a Polygon for the HP-41
Overview
1°) General Case
2°) Brahmagupta's Formula
a) Focal Programs
b) M-Code Routines ( Heron & Brahmagupta
Formulae )
3°) Convex Cyclic Polygons
a) Area ( 2 programs )
b) Diagonals Lengths
4°) Regular Star Polygons
-> The latest addition is the second program in paragraph 3-a)
1°) General Case
-Here, we assume that the coordinates of all the vertices are given.
-The following program calculates the area of a polygon A1
A2 ...... An where Ai(xi,yi)
are the n vertices.
-The result is positive if the n points are stored counterclockwise.
Data Registers: • R00 = n ( All these registers are to be initialized before executing "PGA" )
• R01 = x1 • R03 = x2
.......... • R2n-1 = xn
• R02 = y1 • R04 = y2
.......... • R2n = yn
Flags: /
Subroutines: /
-Synthetic registers M N O may be replaced by any unused data registers.
01 LBL "PGA" 02 RCL 00 03 ST+ X 04 STO M 05 RCL 02 06 RCL IND Y 07 STO O 08 - 09 DSE M 10 RCL 01 11 RCL IND M 12 STO N 13 + 14 * 15 DSE M 16 LBL 01 17 RCL IND M 18 ENTER^ 19 X<> O 20 - 21 DSE M 22 RCL IND M 23 ENTER^ 24 X<> N 25 + 26 * 27 - 28 DSE M 29 GTO 01 30 2 31 / 32 CLA 33 END |
( 56 bytes / SIZE 2n+1 )
STACK | INPUT | OUTPUT |
X | / | Area |
Example: The 6 vertices of an
hexagon are defined by their coordinates (x,y)
xi | -2 | -1 | 2 | 3 | 2 | -1 |
yi | 0 | -2 | -1 | 1 | 2 | 3 |
-We store these 12 numbers in registers
R01 R03
R05 R07 R09 R11
R02 R04
R06 R08 R10 R12
respectively
-Then, 6 STO 00 XEQ "PGA" >>>> Area = 16.0000 ( in 2.7 seconds )
Note:
-Execution time is of the order of 43 seconds if n = 100
2°) Brahmagupta's Formula
a) Focal Programs
-Brahmagupta's formula generalizes Heron's formula for a convex cyclic quadrilateral ( i-e that may be inscribed in a circle )
-Let a , b , c , d its sides lengths, the area A of the cyclic quadrilateral is
A = ( ( s-a ).( s-b ).( s-c ).( s-d ) )1/2
where s = semi-perimeter = ( a + b + c + d )/2
Data Registers: /
Flags: /
Subroutines: /
01 LBL "BRM" 02 STO M 03 RDN 04 ST+ M 05 RDN 06 ST+ M 07 X<>Y 08 ST+ M 09 SIGN 10 ST+ X 11 ST/ M 12 CLX 13 X<> M 14 ST- Y 15 ST- Z 16 ST- T 17 ST- L 18 X<> L 19 * 20 * 21 * 22 SQRT 23 END |
( 43 bytes / SIZE 000 )
STACK | INPUTS | OUTPUTS |
T | a | / |
Z | b | / |
Y | c | / |
X | d | Area |
Example: a = 4 , b = 5 , c = 6 , d = 7
4 ENTER^
5 ENTER^
6 ENTER^
7 XEQ "BRM" >>>> Area = 28.98275349
Data Registers: /
Flags: /
Subroutines: /
01 LBL "BHM" 02 STO L 03 X<>Y 04 ST- Y 05 ST+ L 06 X<> L 07 X^2 08 R^ 09 STO L 10 R^ 11 ST+ Y 12 ST- L 13 X<> L 14 X^2 15 ST- Z 16 RDN 17 X^2 18 RCL Z 19 X^2 20 - 21 * 22 SQRT 23 4 24 / 25 END |
( 43 bytes / SIZE 000 )
STACK | INPUTS | OUTPUTS |
T | a | / |
Z | b | / |
Y | c | / |
X | d | Area |
Example: a = 4 , b = 5 , c = 6 , d = 7
5 ENTER^
6 ENTER^
7 XEQ "BRM" >>>> Area = 28.98275350
Notes:
-This formula may also be applied for a triangle with d = 0
-If we also need the radius R of the circumcircle, it may be computed by
R = (1/4) [ ( a.b + c.d ) ( a.c + b.d ) ( a.d + b.c ) ]1/2
/ Area
Data Registers: R00 = Area , R01 = a
, R02 = b , R03 = c , R04 = d
Flags: /
Subroutines: /
01 LBL "BRMR" 02 STO 00 03 STO 04 04 RDN 05 ST+ 00 06 STO 03 07 RDN 08 ST+ 00 09 STO 02 10 X<>Y 11 ST+ 00 |
12 STO 01
13 SIGN 14 ST+ X 15 ST/ 00 16 X<> 00 17 ST- Y 18 ST- Z 19 ST- T 20 ST- L 21 X<> L 22 * |
23 * 24 * 25 SQRT 26 STO 00 27 RCL 01 28 RCL 02 29 * 30 RCL 03 31 RCL 04 32 * 33 + |
34 RCL 01 35 RCL 03 36 * 37 RCL 02 38 RCL 04 39 * 40 + 41 * 42 RCL 01 43 RCL 04 44 * |
45 RCL 02 46 RCL 03 47 * 48 + 49 * 50 SQRT 51 4 52 / 53 RCL 00 54 ST/ Y 55 END |
( 76 bytes / SIZE 005 )
STACK | INPUTS | OUTPUTS |
T | a | / |
Z | b | / |
Y | c | R |
X | d | Area |
Example: a = 4 , b = 5 , c = 6 , d = 7
4 ENTER^
5 ENTER^
6 ENTER^
7 XEQ "BRM" >>>> Area = 28.98275349
X<>Y R = 3.973161451
b) M-Code Routines (
Heron & Brahmagupta Formulae )
-"HERON" and "BRHM" use the Héron's & Brahmagupta's formulae.
-They are so similar that it is preferable to combine them.
-There is no check for alpha data and synthetic register Q is used.
-13-digit routines will give a better precision...
08E "N"
00F "O"
012 "R"
005 "E"
008 "H"
104 CLRF 8
033 JNC+06
08D "M"
008 "H"
012 "R"
002 "B"
108 SETF 8
2A0 SETDEC
0F8 C=X
128 L=C
2BE C=-C
10E A=C ALL
0B8 C=Y
01D C=
060 A+C
078 C=Z
025 C=
060 AB+C
04E C=0 ALL
10C ?FSET 8
01B JNC+03
270 RAMSLCT
038 READATA
0E8 X=C
025 C=
060 AB+C
089 AB
064 STO Q+
0B8 C=Y
2BE C=-C
10E A=C ALL
138 C=L
01D C=
060 A+C
078 C=Z
025 C=
060 AB+C
0F8 C=X
025 C=
060 AB+C
0D1 RCL
064 Q+
149 C=
060 AB*CM
089 AB
064 STO Q+
078 C=Z
2BE C=-C
10E A=C ALL
138 C=L
01D C=
060 A+C
0F8 C=X
025 C=
060 AB+C
0B8 C=Y
025 C=
060 AB+C
0D1 RCL
064 Q+
149 C=
060 AB*CM
089 AB
064 STO Q+
0F8 C=X
2BE C=-C
10E A=C ALL
138 C=L
01D C=
060 A+C
0B8 C=Y
025 C=
060 AB+C
078 C=Z
025 C=
060 AB+C
0D1 RCL
064 Q+
149 C=
060 AB*CM
04E C
35C
050 =
190
226 16
269 C=
060 AB/C
305 C=
060 sqrt(AB)
0E8 X=C
3E0 RTN
( 96 words )
STACK | INPUTS | OUTPUTS |
T | a | a |
Z | b | b |
Y | c | c |
X | d | Area |
L | / | d |
Example: a = 4 , b = 5 , c = 6 , d = 7
4 ENTER^
5 ENTER^
6 ENTER^
7 XEQ "BRHM" >>>> Area = 28.98275349
-See "Heron's Formula for the HP-41" for an example of area of a triangle...
3°) Convex Cyclic Polygons
a) Area
-"CPLA" computes the area A of a convex cyclic polygon and the circumradius
R, assuming all the sides lengths are known.
-Moreover, we also assume that the center of the circumcircle is inside
the polygon.
-If a1 , a2 , ......... , an
are the sides lengths and µ1 , µ2
, ......... , µn are the corresponding central angles,
we have to solve the system of (n+1) equations:
2.R sin µ1/2 = a1
2.R sin µ2/2 = a2
.........................
2.R sin µn/2 = an
µ1 + µ2 + ... + µn
= 360°
-Substitutions lead to Arc sin ( a1 / 2.R ) + Arc sin ( a2 / 2.R ) + ..................... + Arc sin ( an / 2.R ) = 180° which has only 1 unknown.
-After finding R , the Area is given by
A = ( R2/2 ) ( Sin a1 +
Sin a2 + ........... + Sin an )
Data Registers: R00 to R03 are used by "SOLVE" ( Registers Rbb thru Ree are to be initialized before executing "CPLA" )
• Rbb = a1 • Rbb+1 = a2 ....................... • Ree = an bb > 05 ( R04 & R05: temp )
Flags: /
Subroutine: "SOLVE" ( cf "Linear
and Non-Linear Systems for the HP-41" or another root-finding program )
01 LBL "CPLA" 02 DEG 03 STO 04 04 "T" 05 ASTO 00 06 0 07 LBL 00 08 RCL IND Y 09 X>Y? 10 X<>Y 11 RDN |
12 ISG Y 13 GTO 00 14 101 15 % 16 XEQ "SOLVE" 17 RCL 04 18 0 19 GTO 02 20 LBL "T" 21 STO 05 22 RCL 04 |
23 0 24 LBL 01 25 RCL IND Y 26 RCL 05 27 / 28 ASIN 29 + 30 ISG Y 31 GTO 01 32 PI 33 R-D |
34 - 35 RTN 36 LBL 02 37 RCL IND Y 38 RCL 01 39 / 40 ASIN 41 ST+ X 42 SIN 43 + 44 ISG Y |
45 GTO 02 46 RCL 01 47 2 48 / 49 STO Z 50 X^2 51 * 52 2 53 / 54 CLD 55 END |
( 90 bytes / SIZE 006+??? )
STACK | INPUTS | OUTPUTS |
Y | / | R |
X | bbb.eee | Area |
where bbb > 005
Example: Find the area of a convex cyclic polygon with sides 4 , 5 , 6 , 7 , 8 , 9 , 10
-Store these 7 numbers into R06 thru R12 , then:
6.012 XEQ "CPLA" >>>>
Area = 174.6757942
---Execution time = 89s---
X<>Y R =
8.143816985
Notes:
-The program doesn't work if the center of the circumcircle is outside the polygon.
-Instead of using the secant method, we can also employ Newton's method:
Data Registers: R00 to R03: temp ( Registers Rbb thru Ree are to be initialized before executing "CPLA" )
• Rbb = a1 • Rbb+1 = a2 ....................... • Ree = an bb > 03
Flags: /
Subroutines: /
01 LBL "CPLA" 02 DEG 03 STO 00 04 X<>Y 05 STO 01 06 LBL 01 07 RCL 01 08 STO 02 09 CLX 10 STO 03 11 LBL 02 12 RCL IND 02 13 RCL 00 |
14 ST+ X 15 / 16 ASIN 17 ST+ 03 18 X<> L 19 ENTER^ 20 X^2 21 CHS 22 1 23 + 24 SQRT 25 / 26 + |
27 ISG 02 28 GTO 02 29 RCL 00 30 / 31 R-D 32 RCL 03 33 PI 34 R-D 35 - 36 X<>Y 37 / 38 ST+ 00 39 VIEW 00 |
40 ABS 41 E-7 42 X<Y? 43 GTO 01 44 RCL 01 45 STO 02 46 CLX 47 LBL 03 48 RCL IND 02 49 RCL 00 50 ST+ X 51 / 52 ASIN |
53 ST+ X 54 SIN 55 + 56 ISG 02 57 GTO 03 58 RCL 00 59 STO Z 60 X^2 61 * 62 2 63 / 64 CLD 65 END |
( 91 bytes / SIZE var. )
STACK | INPUTS | OUTPUTS |
Y | bbb.eee | R |
X | r | Area |
where bbb > 003 and r is an estimation of the radius
Same Example: Find the area of a convex cyclic polygon with sides 4 , 5 , 6 , 7 , 8 , 9 , 10
-Store these 7 numbers into R06 thru R12 , then, if you choose r = 10
6.012 ENTER^
10 XEQ "CPLA"
>>>> Area = 174.6757941
---Execution time = 55s---
X<>Y R =
8.143816984
Notes:
-A bad guess may lead to a DATA ERROR
-The successive R-approximations are displayed ( line 39 )
b) Diagonals Lengths
An * .............
A1 * A4 *
A2 *
A3 *
-Let's call d1 = A1 A3
, d2 = A1 A4 , ........................
, dn-3 = A1 An-1
( n-3 ) diagonals passing through A1
and a1
= A1 A2 , a2 = A2
A3 , ........................ , an = An A1
n side-lengths ( n > 3 )
-Applying the law of cosines to triangles like A1 A2
A3 & A1 A3 A4
....................
and taking into account Cos ( A1 A2
A3 ) = - Cos ( A1 A4 A3 )
since the sum of these angles = 180°
lead to the following system of ( n-3 ) equations in ( n-3 ) unknowns
dj = SQRT [ { aj+3 dj+1
( d2j-1 + a2j+1 ) + dj-1
aj+1 ( a2j+3 + d2j+1)
} / ( aj+3 dj+1 + dj-1 aj+1 )
]
j = 1 , 2 , ..... , n-3
if we define a dummy "d0" = a1 and a dummy
"dn-2" = an
-"CPLD" solves this system by a successive approximation method.
-The convergence is linear only.
Data Registers: • R00 = n > 3 ( Registers R00 thru Rnn are to be initialized before executing "CPLD" )
• R01 = a1 • R02 = a2 ....................... • Rnn = an
Rn+1 = a1 = "d0" , Rn+2 = d1 , Rn+3
= d2 , ................. , R2n-2 = dn-3 , R2n-1 = an
= "dn-2"
Flags: /
Subroutines: /
01 LBL "CPLD" 02 CLA 03 RCL 00 04 1 05 + 06 STO N 07 RCL 00 08 ST+ X 09 DSE X 10 STO Z 11 E3 12 / 13 + 14 CLRGX 15 RCL IND 00 16 STO IND Z 17 RCL 01 18 STO IND Z 19 RCL 00 |
20 2 21 - 22 E3 23 / 24 2 25 + 26 STO M 27 LBL 01 28 RCL N 29 2 30 + 31 RCL IND X 32 STO P 33 RCL M 34 ISG X 35 ABS 36 X<>Y 37 RCL IND Y 38 STO Q |
39 * 40 STO Y 41 RCL IND M 42 X^2 43 RCL IND N 44 X^2 45 + 46 * 47 RCL IND M 48 RCL IND N 49 * 50 ST+ Z 51 X<> P 52 X^2 53 RCL Q 54 X^2 55 + 56 RCL P 57 * |
58 + 59 X<>Y 60 / 61 SQRT 62 ISG N 63 CLX 64 ENTER^ 65 X<> IND N 66 - 67 ABS 68 ST+ O 69 ISG M 70 GTO 01 71 RCL 00 72 3 73 - 74 ST- M 75 ST- N 76 CLX |
77 X<> O 78 VIEW X 79 X#0? 80 GTO 01 81 RCL 00 82 ST+ X 83 2 84 - 85 E3 86 / 87 RCL 00 88 + 89 2 90 + 91 CLA 92 CLD 93 END |
( 141 bytes / SIZE 2n )
STACK | INPUT | OUTPUTS |
X | / | bbb.eee |
where bbb.eee is the control number of the ( n-3 ) diagonals passing through A1
Example: Find the diagonals lengths of a convex cyclic hexagon with sides 4 , 5 , 6 , 7 , 8 , 9
n = 6 STO 00 and store
4 , 5 , 6 , 7 , 8 , 9 into R01 thru R06
XEQ "CPLD" >>>> 8.010 ---Execution time = 82s---
-And we get
R08 = d1 = 8.462784374
R09 = d2 = 12.12358501
R10 = d3 = 12.97690534
-There are in fact n(n-3)/2 diagonals whose lengths may be obtained by "rotating" the sides lengths in registers R01 to R06 and we have similarly:
d4 = 9.998827970
d7 = 11.30861231
d5 = 13.01214483
d8 = 13.06010803
n(n-3)/2 = 9 if n = 6
d6 = 11.49035918
d9 = 12.32872367
Notes:
-The successive sums of the differences between 2 consecutive approximations
( in absolute values ) are displayed.
-They should tend to zero.
-However, the termination criterion, line 79 X#0? may lead to
an infinte loop.
-So, it could be replaced for instance by something like: RCL 00
RCL 01 * E9 / X<Y?
-The iteration starts with all diagonals lengths = 0 ( line 14 ) which
is very uninspired !
-Modify the first lines as you like provided they contain the following instructions
before LBL 01
CLA
RCL 01 RCL IND
00 E3
RCL 00 STO IND Y
STO IND Y /
1
RCL 00 RCL
00
2
+
ST+ X
2
+
STO N DSE X
-
STO M
-Since there are 319 registers at most, "CPLD" can find the diagonals lengths
of a 159-gon
-The execution time will not be small without an emulator.
-We can also write an M-code routine to perform the calculations of lines
39 to 61 without disturbing registers M , N , O.
-"DIAG" computes
SQRT [ { x.y ( z2 + t2
) + z.t ( x2 + y2 ) } / ( x.y + z.t ) ]
assuming x , y , z , t are in registers X ,
Y , Z , T upon entry
087 "G"
001 "A"
009 "I"
004 "D"
2A0 SETDEC
0F8 C=X
10E A=C ALL
0B8 C=Y
135 C=
060 A*C
268 Q=C
046 C
270 =
038 T
128 L=C
10E A=C ALL
135 C=
060 A*C
070 N=C ALL
078 C=Z
10E A=C ALL
135 C=
060 A*C
0B0 C=N ALL
025 C=
060 AB+C
278 C=Q
13D C=
060 AB*C
228 P=C
138 C=L
10E A=C ALL
078 C=Z
135 C=
060 A*C
128 L=C
278 C=Q
025 C=
060 AB+C
268 Q=C
0F8 C=X
10E A=C ALL
135 C=
060 A*C
070 N=C ALL
0B8 C=Y
10E A=C ALL
135 C=
060 A*C
0B0 C=N ALL
025 C=
060 AB+C
138 C=L
13D C=
060 AB*C
238 C=P
025 C=
060 AB+C
278 C=Q
269 C=
060 AB/C
305 C=
060 sqrt(AB)
0E8 X=C
3E0 RTN
( 65 words )
-The listing becomes:
01 LBL "CPLD" 02 CLA 03 RCL 00 04 1 05 + 06 STO N 07 RCL 00 08 ST+ X 09 DSE X 10 STO Z 11 E3 12 / 13 + 14 CLRGX |
15 RCL IND 00 16 STO IND Z 17 RCL 01 18 STO IND Z 19 RCL 00 20 2 21 - 22 E3 23 / 24 2 25 + 26 STO M 27 LBL 01 28 RCL N |
29 2 30 + 31 RCL IND X 32 RCL M 33 1 34 + 35 X<>Y 36 RCL IND Y 37 RCL IND M 38 RCL IND N 39 DIAG 40 ISG N 41 CLX 42 ENTER^ |
43 X<> IND N 44 - 45 ABS 46 ST+ O 47 ISG M 48 GTO 01 49 RCL 00 50 3 51 - 52 ST- M 53 ST- N 54 CLX 55 X<> O 56 VIEW X |
57 X#0? 58 GTO 01 59 RCL 00 60 ENTER^ 61 ST+ X 62 2 63 ST+ Z 64 - 65 E3 66 / 67 + 68 CLA 69 CLD 70 END |
( 110 bytes / SIZE 2n )
STACK | INPUT | OUTPUTS |
X | / | bbb.eee |
where bbb.eee is the control number of the ( n-3 ) diagonals passing through A1
Example: The same one: the diagonals lengths of a convex cyclic hexagon with sides 4 , 5 , 6 , 7 , 8 , 9
n = 6 STO 00 and store
4 , 5 , 6 , 7 , 8 , 9 into R01 thru R06
XEQ "CPLD" >>>> 8.010 ---Execution time = 65s---
-And we get
R08 = d1 = 8.462784376
R09 = d2 = 12.12358502
R10 = d3 = 12.97690535
Notes:
-We have saved 17 seconds, significant but not a fantastic improvement
!
-There are small differences in the last places because "DIAG" employs a
few 13-digit routines.
-The modifications suggested for the 1st version also apply with this one.
-After finding the diagonals lengths, we can use Heron's formula to calculate
the polygon area by adding the areas of several triangles.
-Unlike "CPLA", this method will also work if the center of the circumcircle
is outside the polygon.
4°) Regular Star Polygons
-"STLA" computes the area A, the perimeter P, the inradius r and the circumradius R of a regular star polygon { n / k } from its edge length a
An-1 *
........
A0 *
A1 *
A2 * A3
*
-Assuming A0 A1 ...................... An-1
An = A0 is a regular convex polygon ( which is
not obvious on the drawing above... ),
the regular star polygon { n / k } is obtained by connecting
A0 Ak A2k ..............
-In the formulas below, a = A0 Ak = Ak
A2k = .....
Formulae:
A = n R2 Sin (180°/n) Cos (180°k/n)
/ Cos [180°(k-1)/n]
a = 2.R Sin (180°k/n)
r = R Cos (180°k/n)
P = 2.A / r
Data Registers: R00 to R04: temp
Flags: /
Subroutines: /
01 LBL "STLA" 02 STO 00 03 RDN 04 STO 01 05 X<>Y 06 2 07 / 08 STO 03 09 STO 04 10 X^2 |
11 * 12 1 13 CHS 14 ACOS 15 STO 02 16 ST* 00 17 RCL 01 18 / 19 SIN 20 * |
21 RCL 00 22 RCL 01 23 / 24 1 25 P-R 26 ST* 03 27 ST* Z 28 RDN 29 ST/ 03 30 ST/ 04 |
31 X^2 32 / 33 RCL 00 34 RCL 02 35 - 36 RCL 01 37 / 38 COS 39 / 40 ENTER^ |
41 ST+ X 42 RCL 03 43 ST/ Y 44 RCL 04 45 RDN 46 X<> Z 47 END |
( 64 bytes / SIZE 005 )
STACK | INPUTS | OUTPUTS |
T | / | R |
Z | a | r |
Y | n | P |
X | k < n/2 | A |
---Execution time = 3.5s---
Examples:
• a = 1 , n = 5 , k = 2
1 ENTER^
5 ENTER^
2 XEQ "STLA" >>>>
A = 0.310270701
RDN P = 3.819660113
RDN r = 0.162459848
RDN R = 0.525731112
• a = 1 , n = 10 , k = 3
1 ENTER^
10 ENTER^
3 R/S >>>>
A = 0.857567126
RDN P = 4.721359547
RDN r = 0.363271264
RDN R = 0.618033989
• a = PI , n = 41 , k = 13
PI ENTER^
41 ENTER^
13 R/S >>>>
A = 9.855571194
RDN P = 19.37713086
RDN r = 1.017237409
RDN R = 1.871409374
Notes:
-This program works in all angular modes.
-However, DEG mode should be preferable.
-If k = 1, we get the convex regular n-gon. For instance, with a =
1 , n = 5 , k = 1 , "STLA" returns:
A = 1.720477401
P = 5
r = 0.688190960
R = 0.850650808
which corresponds to the regular pentagon.
Reference: