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