Q function of HP-21S unveiled

[mike-stgt]

With the trace feature of Emu42 I was able to find out the two procedures of the HP-21S to compute Q(z) -- one of them was new to me.
Here an HP-42S program using the same algorithms -- [Edit 2025-04-11] do not use this routine, the results show the same error as the HP-21S. See 3rd append of this thread.

Code: Select all

00 { 131-Byte Prgm }
01▶LBL "Q(z)"
02 X<0?
03 AON
04 ABS
05 ENTER
06 X^2
07 -2
08 ÷
09 E^X
10 PI
11 STO+ ST X
12 SQRT
13 ÷
14 STO 02
15 X<>Y
16 2,32
17 X>Y?
18 GTO 02
19 Rv
20 STO 04
21 STO 06
22 ÷
23 STO 05
24 CLX
25 STO 00
26 STO 01
27 SIGN
28 STO 03
29▶LBL 01
30 ISG 00
31 CLD
32 RCL 01
33 RCL× 00
34 RCL 02
35 STO 01
36 RCL× 06
37 +
38 STO 02
39 RCL 03
40 RCL× 00
41 RCL 04
42 STO 03
43 RCL× 06
44 +
45 STO 04
46 ÷
47 ENTER
48 X<> 05
49 -
50 RCL÷ ST L
51 ABS
52 -9
53 10^X
54 X≤Y?
55 GTO 01
56 SIGN
57 RCL 05
58 FS? 48
59 -
60 GTO 09
61▶LBL 02
62 SIGN
63 STO 00
64 Rv
65 ×
66 STO 01
67 LASTX
68 X^2
69 STO 02
70 X<>Y
71▶LBL 03
72 RCL 01
73 RCL× 02
74 2
75 RCL+ 00
76 STO 00
77 ÷
78 STO 01
79 RCL+ ST Y
80 X≠Y?
81 GTO 03
82 2
83 1/X
84 FC? 48
85 +/-
86 +
87 ABS
88▶LBL 09
89 AOFF
90 END

Note: in contrast to the HP-32E the HP-21S calculates "the upper tail" of the Q function.

The formulas? If you are truly interested take a closer look to this user code routine, it is much simpler than analysing some firmware trace.

M.

Edit: replaced lines 13 1/x and 14 * by 13 / only. It was a holdover from the HP-21S, which uses the constant @03966: #0398942280401432.

[mike-stgt]

[...] I was able to find out the two procedures of the HP-21S to compute Q(z) -- one of them was new to me.

What a shame! Since years I overlooked it, it's A&S 26.2.14 -- no reason to hide any longer my insight from trace logs:

my finding

Q(z) part 1 cfe.png (2.45 KiB) Viewed 2844 times

To compute Q(z) the HP-21S uses "Algorithm 304" (almost) what means iterating either formula A&S 26.2.11 or A&S 26.2.14 depending on input. The switchpoint between the two formulas is x=2.32, this is attributed as "speed is sacrificed to accuracy".

Edit: inserted "almost"

[mike-stgt]

Since a couple of days I know the Q(z) function of the HP-21S uses "Algorithm 304" but somehow unprecise. Now I found a suggestion how modify algorithm 304 so to reduce the iterations computing the CDF's tails. Following this approach I created a new program for the HP-42S, adapted the "upper tail switchpoint" to lessen the relative error, and made it system compatible: input saved in LastX, stack registers Y, Z, T and data registers stay untouched (are restored).
The followning graph shows the result's relative error, which is typically below 2E-11, for input < 1 even less, only in the input range 1..2 relative deviation up to 9E-11 may occure.

Relative error (results of HP-28S taken as "true values")

Note: center part in red from -3.5 to 1.77 instead of 2.32 (nor 1.28).

Next grap shows clearly the improvement compared to HP-21S.

Q(z) relative error (referred to HP-28S)

Note: green = HP-21S, red = HP-42S program of OP

The HP-42S program -- QPL applies -- in short: if you improve it you have to inform me, if you publish your improvement my original must be included.

Code: Select all

00 { 303-Byte Prgm }
01▶LBL "304mx"
02 Rv
03 STO "_Y"
04 Rv
05 STO "_Z"
06 X<>Y
07 STO "_T"
08 RCL "REGS"
09 STO "_R"
10 CLX
11 SIZE 12
12 CLRG
13 R^
14 X≠0?
15 GTO 00
16 2
17 1/X
18 GTO 09
19▶LBL 00
20 X<0?
21 AON
22 ABS
23 STO 00
24 X^2
25 -2
26 ÷
27 E^X
28 PI
29 STO+ ST X
30 SQRT
31 ÷
32 STO 01
33 RCL÷ 00
34 STO 02
35 FC? 48
36 GTO 01
37 1
38 RCL- ST Y
39 1
40 X=Y?
41 GTO 09
42 3,5
43 GTO 02
44▶LBL 01
45 X=0?
46 GTO 09
47 PI
48 SQRT
49▶LBL 02
50 RCL 00
51 X≤Y?
52 GTO 05
53 STO 07
54 X^2
55 3
56 +
57 STO 02
58 2
59 STO 03
60 RCL 01
61 STO 05
62 RCL÷ 00
63 STO 09
64 RCL 02
65 1
66 -
67 RCL× 01
68 STO 06
69 RCL 02
70 RCL× 00
71 STO 08
72 ÷
73 STO 10
74 FC? 48
75 GTO 03
76 1
77 RCL- 09
78 STO 09
79 1
80 RCL- 10
81 STO 10
82▶LBL 03
83 RCL 09
84 RCL 10
85 X=Y?
86 GTO 09
87 RCL 11
88 X<>Y
89 X=Y?
90 GTO 09
91 4
92 STO+ 02
93 RCL 03
94 8
95 -
96 STO 03
97 RCL+ 04
98 STO 04
99 RCL× 05
100 RCL 02
101 RCL× 06
102 +
103 X<> 06
104 STO 05
105 RCL 07
106 RCL× 04
107 RCL 02
108 RCL× 08
109 +
110 X<> 08
111 STO 07
112 RCL 09
113 STO 11
114 RCL 10
115 STO 09
116 RCL 06
117 RCL÷ 08
118 1
119 X<>Y
120 FS? 48
121 -
122 STO 10
123 GTO 03
124▶LBL 05
125 RCL 00
126 RCL× 01
127 STO 03
128 STO 11
129 1
130 STO 02
131▶LBL 06
132 RCL 10
133 RCL 11
134 X=Y?
135 GTO 07
136 STO 10
137 RCL 00
138 X^2
139 2
140 RCL+ 02
141 STO 02
142 ÷
143 RCL× 03
144 STO 03
145 STO+ 11
146 GTO 06
147▶LBL 07
148 2
149 1/X
150 FC? 48
151 +/-
152 +
153 ABS
154▶LBL 09
155 RCL 00
156 FS? 48
157 +/-
158 SIGN
159 CLX
160 RCL "_R"
161 STO "REGS"
162 CLX
163 X<>Y
164 RCL "_T"
165 RCL "_Z"
166 RCL "_Y"
167 R^
168 CLV "_R"
169 CLV "_T"
170 CLV "_Z"
171 CLV "_Y"
172 AOFF
173 .END.

Same as hex dump. Convert it to binary, save with filetype RAW and use menu Edit/Load Object... of Emu42 running an HP-42S.

Code: Select all

CA 00 F6 00 33 30 34 6D 78 75 F3 81 5F 59 75 F3 81 5F 5A 71 F3 81 5F 54 F5 91 52 45 47 53 F3 81 5F 52 77 F3 F7 00 0C 8A 74 63 B1 85 12 00 60 BA 00 01 66 8C 61 30 51 1C 12 00 43 55 72 92 73 52 43 31 F2 D4 00 32 AD 30 B2 90 11 00 F2 D2 72 11 00 78 BA 00 13 1A 15 00 B3 86 02 67 BA 00 72 52 03 20 46 B6 EE 37 51 13 00 40 32 12 00 33 21 35 F2 D4 00 39 22 11 00 41 F2 D3 01 36 22 F2 D3 00 38 43 3A AD 30 B4 8C 11 00 F2 D2 09 39 11 00 F2 D2 0A 3A 04 29 2A 78 BA E5 2B 71 78 BA E0 14 00 92 02 23 18 00 41 33 F2 D1 04 34 F2 D3 05 22 F2 D3 06 40 CE 06 35 27 F2 D3 04 22 F2 D3 08 40 CE 08 37 29 3B 2A 39 26 F2 D4 08 11 00 71 AC 30 41 3A B4 40 06 20 F2 D3 01 33 3B 11 00 32 07 2A 2B 78 B8 92 3A 20 51 12 00 F2 D1 02 32 43 F2 D3 03 33 92 0B B7 18 08 12 00 60 AD 30 54 40 61 0A 20 AC 30 54 7A 77 F3 91 5F 52 F5 81 52 45 47 53 77 71 F3 91 5F 54 F3 91 5F 5A F3 91 5F 59 74 F3 B0 5F 52 F3 B0 5F 54 F3 B0 5F 5A F3 B0 5F 59 8B C0 00 0D

Opinion: Today nobody will use an HP-42S for serious work any more. Thus the routine I present here is pure nostalgia. But it's fun to dig through old papers and apply ancient algorithms on obsolete calculators.
If you wonder about the unrealistic test range -- I also thought for many years the region of interest may be the statistician’s 0 < x < 5, but had to learn now, for solutions of diffusion or heat equations it's rather 10 < x < 14. The used test range from -7 to 23 is just to make sure it works.