## 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
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

```