Compact Simpson's Rule

[Dieter]

(12-13-2015 02:06 AM)Namir Wrote:  This listing calculates the integral for f(X)=1/X (see label E) using Simpson's rule. The code is based on the compact pseudo-code that I posted in a thread in the General area.

Here's my version with yet another method of toggling between 2* and 4*f(x): use a flag. This also frees up one register, and even one more can be saved as once h is calculated the value of b is no longer required. So that's two registers less and some bytes of code saved. And one function call as f(b) is not calculated twice.

Code:

01 LBL "SIMP"
02 LBL A
03 " A^B=?"
04 PROMPT
05 STO 02   ' store upper limit for integral
06 X<>Y
07 STO 01   ' store lower limit for integral
08 -
09 STO 03   ' prestore b-a
10 RCL 01
11 XEQ E    ' calculate f(a)
12 STO 00   ' sum = f(a)
13 RCL 02
14 XEQ E    ' calculate f(b)
15 ST+ 00   ' sum = f(a) + f(b) ... no trick required ;-)
16 " N=?"
17 PROMPT
18 STO 02   ' store number of divisions. Must be even.
19 ST/ 03   ' calculate and store h = (b-a)/n
20 SF 02
21 DSE 02   ' n = n - 1
22 LBL 00   ' start loop
23 RCL 03
24 ST+ 01   ' x := x + h
25 RCL 01
26 XEQ E    ' calculate f(x)
27 ST+ X
28 FS? 02
29 ST+ X
30 ST+ 00   ' sum = sum + [2|4]*f(x)
31 FC?C 02  ' toggle flag 2
32 SF 02
33 DSE 02   ' n := n - 1
34 GTO 00   ' if n > 0 do another loop
35 RCL 00
36 RCL 03
37 *
38 3
39 /
40 RTN

41 LBL E    ' Calculate f(x).
42 1/x      ' Insert your f(x) code here.
43 RTN
 
Memory Map
----------

R00 = sum
R01 = a, x
R02 = b, n, counter
R03 = h

Maybe t's also fun to watch the flashing "2" annunciator. ;-)

Dieter