zeroin.a68
1 COMMENT
2
3 @section Synopsis
4
5 Zeroin is a classic root-finding algorithm.
6
7 After MCA 2310 in 'ALGOL 60 Procedures in Numerical Algebra' by Th. J. Dekker.
8
9 This Algol 68 version originates from the legacy "REVISED MC ALGOL 68 TEST SET":
10
11 Dick Grune, The Revised MC ALGOL 68 Test Set, IW XX/79,
12 Mathematical Centre, Amsterdam.
13
14 The Mathematical Centre ("Stichting Mathematisch Centrum" or SMC) was a Dutch
15 non-profit institution aiming at the promotion of pure mathematics and its
16 applications.
17
18 SMC is now "Stichting Centrum Wiskunde & Informatica" (CWI). The test set
19 is available as an open access publication from the CWI repository:
20
21 https://ir.cwi.nl/pub/
22
23 Selected (modified) "Revised MC ALGOL 68 Test Set" programs are distributed
24 with Algol 68 Genie with kind permission of Dick Grune.
25
26 COMMENT
27
28 BEGIN MODE NUM = LONG LONG REAL;
29
30 PROC zero in = (REF NUM x, y, PROC (NUM) NUM f, tol) BOOL:
31 BEGIN NUM a := x, b := y; NUM fa := f (a), fb := f (b);
32 NUM c := a, fc := fa;
33 WHILE
34 (ABS fc < ABS fb |
35 (a := b, fa := fb); (b := c, fb := fc); (c := a, fc := fa));
36 NUM tolb := tol (b), m := (c + b) * .5;
37 ABS (m - b) > tolb
38 DO NUM p := (b - a) * fb, q := fa - fb;
39 (p < 0 | (p := -p, q := -q));
40 (a := b, fa := fb);
41 fb := f (b :=
42 IF p <= ABS q * tolb
43 THEN (c > b | b + tolb | b - tolb)
44 ELSE (p < (m - b) * q | p / q + b | m)
45 FI);
46 IF ABS (SIGN fb + SIGN fc) = 2
47 THEN (c := a, fc := fa)
48 FI
49 OD;
50 (x := b, y := c); ABS (SIGN fb + SIGN fc) < 2
51 END;
52
53 # Pretty print NUM. #
54
55 OP PRETTY = (NUM z) STRING:
56 IF ABS (z - ROUND z) < 10 * long long small real
57 THEN whole (z, 0)
58 ELIF ABS (z) >= 0.1 THEF ABS (z) <= 10
59 THEN STRING buf;
60 FOR digits TO long real width
61 WHILE puts (buf, fixed (z, 0, digits));
62 NUM y;
63 gets (buf, y);
64 z /= y
65 DO ~ OD;
66 buf
67 ELSE STRING buf, INT expw = 4;
68 FOR digits TO long real width
69 WHILE puts (buf, float (z, 4 + digits + expw, digits, expw));
70 NUM y;
71 gets (buf, y);
72 z /= y
73 DO ~ OD;
74 buf
75 FI;
76
77 NUM eps = long small real;
78
79 PROC solve = (NUM u, v, PROC (NUM) NUM f, STRING s) VOID:
80 IF u >= v
81 THEN error
82 ELSE print (("Expression: ", s, new line));
83 print (("[", PRETTY u, ", ", PRETTY v, "]", new line));
84 IF NUM x, y;
85 zero in (x := u, y := v, f, (NUM p) NUM : eps + eps * ABS p)
86 THEN print (("x: ", PRETTY x, new line));
87 print (("y: ", PRETTY f(x), new line))
88 ELSE print (("zeroin cannot find a root", new line))
89 FI;
90 new line (standout)
91 FI;
92
93 MODE FUN = PROC (NUM) NUM;
94 PROC pow = (FUN f, INT n, NUM x) NUM: f (x) ** n;
95 OP ** = (FUN f, INT n) FUN: pow (f, n, );
96
97 solve (-long long pi, long long pi, (NUM x) NUM : (long long sin **2) (x) - 1, "sin^2 x - 1");
98 solve (-1, 1, (NUM x) NUM: x * x + 1, "x^2 + 1");
99
100 EXIT
101
102 error: print (("Error detected", new line))
103 END
104
105 COMMENT
106
107 Output:
108
109 Expression: sin^2 x -1
110 [-3.14159265358979323846264338327950, 3.14159265358979323846264338327950]
111 x: 1.57079632679489661923132169163975
112 y: 0
113
114 Expression: x^2 + 1
115 [-1, 1]
116 zeroin cannot find a root
117
118 COMMENT
