-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathPERE12.TXT
101 lines (101 loc) · 3.86 KB
/
PERE12.TXT
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
; calcul ellipse perimeter GaussKummer-like with AGM exact function for
; "convergence boosting"
;
; "HP41 Iterated function GAGM based"
;
; Execution/Inputs
; a ENTER b ENTER XEQ PERE12
; a and b: half parameter of the ellipse
; in case a or b = 0, it would be a flat ellipse. The programm is considering
; this but must not consider this (result will be "perimeter = 4*(a+b)")
; and would be shorter by several steps
;
; Outputs
; results in stack X: perimeter
;
; Modules used
; None
;
; use Register 07..10
; keep free Reg 00-06 if use with SOL from MATH from other programs
; in case this program would be used for example for finding the "a" in case
; the "perimeter" and the "b" are given
;
; create raw file with "rpncomp --raw-output PERE12.TXT", upload
; in V41 emulator and store into virtual drive pyILPER
;
; under CC BY SA CreativeCommons 4.0
; pascaldagornet at yahoo dot de
;
; idea taken from here
; GAGM
; https://semjonadlaj.com/SP/Computer+Algebra+in+Scientific+Computing_37-56.pdf
; Coding steps
; https://www.hpmuseum.org/forum/thread-5820-post-171326.html#pid171326
; P(a,b)=2*( P((a+b)/2,sqrt(a*b)) - pi*a*b/AGM(a,b) )
;
; change log
; date 2023 04 15 creation and release (based on PEREL6 simplified program)
; date 2023 04 16 R11 use deleted
; date 2023 04 17 X#Y?
;
; comments
; all an bn R0xn are from the "n" loop
; all with "N" are updated parameter in the new loop
;
LBL "PERE12"
X=0? ; in case b = 0 this is a flat ellipse
GTO 01 ; ... go out in this special case for calc error avoidance
STO 07 ; b in 07
X<>Y ; ... go out in this special case for calc error avoidance
X=0? ; in case a = 0 this is a flat ellipse
GTO 01 ; ... go out in this special case for calc error avoidance
STO 08 ; a in 08
X^2
X<>Y
X^2
+
STO 09 ; a^2 + b^2 in R09n for temporary storage of this result
1
STO 10 ; 1 in R10n = parameter in this program
;
LBL 00 ; X Y Z T
RCL 08 ; an ? ? ?
RCL X ; an an ? ?
RCL 07 ; bn an an ?
STO T ; bn an an bn
* ; bn*an an bn bn
SQRT ; SQRT(an*bn)=bN an bn bn
X<> 07 ; bn an bn bn
; ... R07 = bN
RCL 07 ; bN bn an bn
RDN ; bn an bn bN
- ; an-bn bn bN bN
2 ; 2 an-bn bn bN
ST* 10 ; 2 an-bn bn bN
; ... R10N= 2*R10n
/ ; (an-bn)/2 bn bN bN
ST- 08 ; (an-bn)/2 bn bN bN
; ... R08N=aN=an-(an-bn)/2=(an+bn)/2
X^2 ; ((an-bn)/2)^2 bn bN bN
RCL 10 ; R10N ((an-bn)/2)^2 bn bN
* ; R10N*((an-bn)/2)^2 bn bN bN
ST- 09 ; R10N*((an-bn)/2)^2 bn bN bN
; ... R09N = R09n- R10N*((an-bn)/2)^2
RDN ; bn bN bN R10N*((an-bn)/2)^2
X#Y? ; ... this comparison of SQRT has a a good convergence so far
; this is not the case if ax is taken
GOTO 00
1/X ; 1/bn bN bN R10N*((an-bn)/2)^2
RCL 09 ; R09N 1/bn bN bN
* ; R09N/bn bN bN bN
PI ; pi R09N/bn bN bN
* ; pi*R09N/bn bN bN bN
RTN
;
LBL 01 ; in case b or a = 0 flat ellipse it goes here out
+ ; a+b
4 ; 4 a+b
* ; 4*(a+b)
RTN
END