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