HP42s first major program (Double Integral) Best way to approach?

[Albert Chan]

This version uses my implementation of EKmc(m) = E(m) - c, K(m) - c, where c = Pi/2

I added HV(d,D) on top of EKmc(m) to test its accuracy.

HV(d, D) = D/3 (d^2 (E+K) + D^2 (E-K))
              = D/3 (d^2 ((E-c)+(K-c) + Pi) + D^2 ((E-c)-(K-c)))

Note: this is a recursive algorithm, requiring Free42 LSTO command.

Code:

00 { 167-Byte Prgm }
01▸LBL "HV"
02 LSTO "c"
03 X<>Y
04 LSTO "b"
05 R↓
06 XEQ 03       ; HV(c,a)
07 X<> "b"
08 RCL "c"
09 XEQ 03       ; HV(c,b)
10 RCL "b"
11 RTN
12▸LBL 03       ; HV(d,D)
13 X>Y?         ; = D/3*(dd*(e+k) + DD*(e-k))
14 X<>Y         ; = D/3*((dd+DD)*(e-c) + (dd-DD)*(k-c) + dd*(2c))
15 X↑2
16 LSTO "d"     ; d := d^2
17 X<>Y
18 LSTO "D"
19 X=Y?
20 GTO 01       ; HV(d,d) = 2/3*d^3
21 X↑2
22 ÷            ; m = (dd)/(DD)
23 XEQ "EKmc"   ; e-c    k-c
24 RCL "d"
25 RCL "D"
26 X↑2
27 -
28 STO× ST Z    ; dd-DD  e-c    (dd-DD)*(k-c)
29 X<> ST L
30 RCL+ "d"     ; dd+DD  e-c    (dd-DD)*(k-c)
31 ×
32 +
33 PI
34 RCL× "d"
35▸LBL 00
36 +
37 RCL× "D"
38 3
39 ÷
40 RTN
41▸LBL 01       ; (D) -> 2/3*D^3
42 X↑2
43 ENTER
44 GTO 00
45▸LBL "EKmc"   ; E(m)-c, K(m)-c, where c=PI/2
46 ENTER
47 ABS
48 1
49 +
50 LASTX        ; 1  1+|m|  m
51 X=Y?
52 GTO 02       ; EKmc base case
53 STO ST Y
54 RCL- ST Z    ; 1-m  1    m
55 SQRT
56 LSTO "b"     ; b    1   m
57 +
58 X↑2
59 ÷
60 LSTO "h"     ; h = (1-b)/(1+b) = m/(1+b)^2
61 X↑2
62 XEQ "EKmc"   ; EKmc(h*h) = e, k
63 RCL× "b"     ; EKmc(m) = (hk-k) + (e+be) - bhc, (hk+k) + hc
64 LASTX
65 +
66 X<>Y
67 RCL× "h"
68 LASTX
69 RCL+ ST Y
70 X<>Y
71 LASTX
72 -
73 R↑
74 +
75 PI
76 2
77 ÷
78 RCL× "h"
79 STO+ ST Z
80 RCL× "b"
81 -
82 RTN
83▸LBL 02
84 PI           ; PI    1   1+|m|   m
85 R↑
86 ×
87 8
88 ÷
89 ENTER
90 +/-          ; -pi/8*m   pi/8*m
91 END

Let's try HV(1,1000):

HV(1,1000) = 785.3980652226656130844841050849260, error = -2 ULP
BORE       = 785.3980652226656130844841051786603, error = -937345 ULP
Mathemtaica: 785.3980652226656130844841050849257781988463583655713086050...

* Mathematica: HV = Pi/4 d^2 D Hypergeometric2F1[-1/2,1/2,2,(d/D)^2] /. {d->1, D->1000}

Update: HV now has same functionality of BORE, and better accuracy.
Formula changed to avoid subtraction cancellation, when m = (d/D)² ≈ 1

HV(d, D) = D/3*((dd+DD)*(e-c) + (dd-DD)*(k-c) + dd*(2c))
→ HV1(m) = 1/3*((m+1)*(e-c) + (m-1)*(k-c) + m*pi)

First 2 terms both negative, with -sum less than 1/2 of the third.
We know this because 2/3 ≤ HV1/m ≤ pi/4, and (third term)/m = pi/3

→ (third term)/HV1 = [4/3, pi/2] ≈ [1.3333, 1.5708]
→ -(sum of first 2 terms)/(third term) ≈ [0.25000, 0.36338]