c_cpp/etc/calc/cal/brentsolve.cal

   1 /*
   2  *  brentsolve- Root finding with the Brent-Dekker trick.
   3  *
   4  * Copyright (C) 2013 Christoph Zurnieden
   5  *
   6  * brentsolve is open software; you can redistribute it and/or modify it under
   7  * the terms of the version 2.1 of the GNU Lesser General Public License
   8  * as published by the Free Software Foundation.
   9  *
  10  * brentsolve is distributed in the hope that it will be useful, but WITHOUT
  11  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  12  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
  13  * Public License for more details.
  14  *
  15  * A copy of version 2.1 of the GNU Lesser General Public License is
  16  * distributed with calc under the filename COPYING-LGPL.  You should have
  17  * received a copy with calc; if not, write to Free Software Foundation, Inc.
  18  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  19  *
  20  * @(#) $Revision: 30.3 $
  21  * @(#) $Id: brentsolve.cal,v 30.3 2013/08/11 08:41:38 chongo Exp $
  22  * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/brentsolve.cal,v $
  23  *
  24  * Under source code control:   2013/08/11 01:31:28
  25  * File existed as early as:    2013
  26  */
  27
  28 static resource_debug_level;
  29 resource_debug_level = config("resource_debug", 0);
  30
  31
  32 /*
  33   A short explanation is at http://en.wikipedia.org/wiki/Brent%27s_method
  34   I tried to follow the description at wikipedia as much as possible to make
  35   the slight changes I did more visible.
  36   You may give http://people.sc.fsu.edu/~jburkardt/cpp_src/brent/brent.html a
  37   short glimpse (Brent's originl Fortran77 versions and some translations of
  38   it).
  39 */
  40
  41 static true = 1;
  42 static false = 0;
  43 define brentsolve(low, high,eps){
  44   local a b c d fa fb fc fa2 fb2 fc2 s fs  tmp tmp2  mflag i places;
  45   a = low;
  46   b = high;
  47   c = 0;
  48
  49   if(isnull(eps))
  50     eps = epsilon(epsilon()*1e-3);
  51   places = highbit(1 + int( 1/epsilon() ) ) + 1;
  52
  53   d = 1/eps;
  54
  55   fa = f(a);
  56   fb = f(b);
  57
  58   fc = 0;
  59   s = 0;
  60   fs = 0;
  61
  62   if(fa * fb >= 0){
  63     if(fa < fb){
  64       epsilon(eps);
  65       return a;
  66     }
  67     else{
  68       epsilon(eps);
  69       return b;
  70     }
  71   }
  72
  73   if(abs(fa) < abs(fb)){
  74     tmp = a; a = b; b = tmp;
  75     tmp = fa; fa = fb; fb = tmp;
  76   }
  77
  78   c = a;
  79   fc = fa;
  80   mflag = 1;
  81   i = 0;
  82
  83   while(!(fb==0) && (abs(a-b) > eps)){
  84     if((fa != fc) && (fb != fc)){
  85       /* Inverse quadratic interpolation*/
  86       fc2 = fc^2;
  87       fa2 = fa^2;
  88       s = bround(((fb^2*((fc*a)-(c*fa)))+(fb*((c*fa2)-(fc2*a)))+(b*((fc2*fa)
  89                     -(fc*fa2))))/((fc - fb)*(fa - fb)*(fc - fa)),places++);
  90     }
  91     else{
  92       /* Secant Rule*/
  93       s =bround( b - fb * (b - a) / (fb - fa),places++);
  94     }
  95     tmp2 = (3 * a + b) / 4;
  96     if( (!( ((s > tmp2) && (s < b))||((s < tmp2) && (s > b))))
  97         || (mflag && (abs(s - b) >= (abs(b - c) / 2)))
  98         || (!mflag && (abs(s - b) >= (abs(c - d) / 2)))) {
  99       s = (a + b) / 2;
 100       mflag = true;
 101     }
 102     else{
 103       if( (mflag && (abs(b - c) < eps))
 104           || (!mflag && (abs(c - d) < eps))) {
 105         s = (a + b) / 2;
 106         mflag = true;
 107       }
 108       else
 109         mflag = false;
 110     }
 111     fs = f(s);
 112     c = b;
 113     fc = fb;
 114     if (fa * fs < 0){
 115       b = s;
 116       fb = fs;
 117     }
 118     else {
 119       a = s;
 120       fa = fs;
 121     }
 122
 123     if (abs(fa) < abs(fb)){
 124       tmp = a; a = b; b = tmp;
 125       tmp = fa; fa = fb; fb = tmp;
 126     }
 127     i++;
 128     if (i > 1000){
 129       epsilon(eps);
 130       return newerror("brentsolve: does not converge");
 131     }
 132   }
 133   epsilon(eps);
 134   return b;
 135 }
 136
 137 /*
 138   A variation of the solver to accept functions named differently from "f". The
 139   code should explain it.
 140 */
 141 define brentsolve2(low, high,which,eps){
 142   local a b c d fa fb fc fa2 fb2 fc2 s fs  tmp tmp2  mflag i places;
 143   a = low;
 144   b = high;
 145   c = 0;
 146
 147   switch(param(0)){
 148     case 0:
 149     case 1: return newerror("brentsolve2: not enough argments");
 150     case 2: eps = epsilon(epsilon()*1e-2);
 151             which = 0;break;
 152     case 3: eps = epsilon(epsilon()*1e-2);break;
 153     default: break;
 154   };
 155   places = highbit(1 + int(1/epsilon())) + 1;
 156
 157   d = 1/eps;
 158
 159   switch(which){
 160     case 1:  fa = __CZ__invbeta(a);
 161              fb = __CZ__invbeta(b); break;
 162     case 2:  fa = __CZ__invincgamma(a);
 163              fb = __CZ__invincgamma(b); break;
 164
 165     default: fa = f(a);fb = f(b);   break;
 166   };
 167
 168   fc = 0;
 169   s = 0;
 170   fs = 0;
 171
 172   if(fa * fb >= 0){
 173     if(fa < fb)
 174       return a;
 175     else
 176       return b;
 177   }
 178
 179   if(abs(fa) < abs(fb)){
 180     tmp = a; a = b; b = tmp;
 181     tmp = fa; fa = fb; fb = tmp;
 182   }
 183
 184   c = a;
 185   fc = fa;
 186   mflag = 1;
 187   i = 0;
 188
 189   while(!(fb==0) && (abs(a-b) > eps)){
 190
 191     if((fa != fc) && (fb != fc)){
 192       /* Inverse quadratic interpolation*/
 193       fc2 = fc^2;
 194       fa2 = fa^2;
 195       s = bround(((fb^2*((fc*a)-(c*fa)))+(fb*((c*fa2)-(fc2*a)))+(b*((fc2*fa)
 196                     -(fc*fa2))))/((fc - fb)*(fa - fb)*(fc - fa)),places);
 197       places++;
 198     }
 199     else{
 200       /* Secant Rule*/
 201       s =bround( b - fb * (b - a) / (fb - fa),places);
 202       places++;
 203     }
 204     tmp2 = (3 * a + b) / 4;
 205     if( (!( ((s > tmp2) && (s < b))||((s < tmp2) && (s > b))))
 206         || (mflag && (abs(s - b) >= (abs(b - c) / 2)))
 207         || (!mflag && (abs(s - b) >= (abs(c - d) / 2)))) {
 208       s = (a + b) / 2;
 209       mflag = true;
 210     }
 211     else{
 212       if( (mflag && (abs(b - c) < eps))
 213           || (!mflag && (abs(c - d) < eps))) {
 214         s = (a + b) / 2;
 215         mflag = true;
 216       }
 217       else
 218         mflag = false;
 219     }
 220     switch(which){
 221       case 1:  fs = __CZ__invbeta(s); break;
 222       case 2:  fs = __CZ__invincgamma(s); break;
 223
 224       default: fs = f(s);             break;
 225     };
 226     c = b;
 227     fc = fb;
 228     if (fa * fs < 0){
 229       b = s;
 230       fb = fs;
 231     }
 232     else {
 233       a = s;
 234       fa = fs;
 235     }
 236
 237     if (abs(fa) < abs(fb)){
 238       tmp = a; a = b; b = tmp;
 239       tmp = fa; fa = fb; fb = tmp;
 240     }
 241     i++;
 242     if (i > 1000){
 243       return newerror("brentsolve2: does not converge");
 244     }
 245   }
 246   return b;
 247 }
 248
 249
 250 /*
 251  * restore internal function from resource debugging
 252  */
 253 config("resource_debug", resource_debug_level),;
 254 if (config("resource_debug") & 3) {
 255     print "brentsolve(low, high,eps)";
 256     print "brentsolve2(low, high,which,eps)";
 257 }