Algol 68 Genie

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
Deze website is gearchiveerd door de KB, nationale bibliotheek.

This website is archived by the National Library of the Netherlands.