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