polygon

# 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

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
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

( 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:

-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:

[1]  http://mathworld.wolfram.com/StarPolygon.html