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 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        COMMENT
 101  
 102        From a letter written by De Morgan to Whewell in 1861:
 103        “The reason I call x^3-2x-5=0 a celebrated equation is because it was the one on which Wallis chanced 
 104        to exhibit Newton’s method when he first published it, in consequence of which every numerical solver 
 105        has felt bound in duty to make it one of his examples. Invent a numerical method, neglect to show how 
 106        it works on this equation, and you are a pilgrim who does not come in at the little wicket”
 107  
 108        From “Computer Methods for Mathematical Computations”, George E. Forsythe, Michael A. Malcolm, Cleve B. Moler.
 109  
 110        COMMENT
 111  
 112        solve (2, 3, (NUM x) NUM: x * (x * x - 2) - 5, "x^3 - 2x - 5")
 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  Expression: x^3 - 2x - 5
 133  [2, 3]
 134  x: 2.09455148154232659148238654057930
 135  y: +2.890624152673e -51
 136  
 137  COMMENT
Deze website is gearchiveerd door de KB, nationale bibliotheek.

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