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