C-XSC - A C++ Class Library for Extended Scientific Computing
2.5.4
|
00001 /* 00002 ** CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4) 00003 ** 00004 ** Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik, 00005 ** Universitaet Karlsruhe, Germany 00006 ** (C) 2000-2014 Wiss. Rechnen/Softwaretechnologie 00007 ** Universitaet Wuppertal, Germany 00008 ** 00009 ** This library is free software; you can redistribute it and/or 00010 ** modify it under the terms of the GNU Library General Public 00011 ** License as published by the Free Software Foundation; either 00012 ** version 2 of the License, or (at your option) any later version. 00013 ** 00014 ** This library is distributed in the hope that it will be useful, 00015 ** but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00017 ** Library General Public License for more details. 00018 ** 00019 ** You should have received a copy of the GNU Library General Public 00020 ** License along with this library; if not, write to the Free 00021 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00022 */ 00023 00024 /* CVS $Id: l_rmath.cpp,v 1.28 2014/01/30 17:23:46 cxsc Exp $ */ 00025 00026 #include <l_rmath.hpp> 00027 #include <rmath.hpp> 00028 #include <l_interval.hpp> 00029 00030 namespace cxsc { 00031 00032 l_real sqrt(const l_real &x) throw(ERROR_LREAL_STD_FKT_OUT_OF_DEF) 00033 // Blomquist, additional scaling, 10.12.02 00034 { 00035 int stagsave = stagprec, stagmax = 19, stagcalc; 00036 l_real y; 00037 if (sign(x[1])<0) // Blomquist: is faster than x < 0.0; 00038 cxscthrow(ERROR_LREAL_STD_FKT_OUT_OF_DEF("l_real sqrt(const l_real &x)")); 00039 else if (zero_(x) || x == 1.0 ) // Blomquist, faster than x == 0; 00040 y = x; 00041 else 00042 { 00043 l_real x1 = x; 00044 int ex = expo(x1[1]); 00045 ex = 1021-ex; // ex: optimal scaling factor 00046 if (ex%2) ex--; // ex is now even; 00047 times2pown(x1,ex); // scaling x1 with 2^ex 00048 00049 // Einsatz des Newton-Verfahrens y = y-f(y)/f'(y) 00050 if (stagprec < stagmax) stagcalc = stagprec+1; 00051 else stagcalc = stagmax+1; // Blomquist with +1; 30.11.02; 00052 y = sqrt(_real(x1)); 00053 stagprec = 1; 00054 while (stagprec < stagcalc) 00055 { 00056 stagprec += stagprec; // now calculation in double precision: 00057 if (stagprec > stagcalc) stagprec = stagcalc; 00058 // y = 0.5*((x1/y)+y); 00059 y += x1/y; // Blomquist, 30.11.02; 00060 times2pown(y,-1); // Blomquist, this multiplication is faster!! 00061 } 00062 times2pown(y,-ex/2); // backscaling 00063 stagprec = stagsave; // restore old stagprec 00064 adjust(y); 00065 } 00066 return y; 00067 } 00068 00069 l_real sqrtx2y2(const l_real& x, 00070 const l_real& y) throw() // Blomquist 6.12.02 00071 // Calculation of an approximation of sqrt(x^2+y^2). 00072 // In general the maximum precision is stagprec=19, predifined by the used 00073 // sqrt-function declared in l_rmath.hpp. 00074 // If the difference |exa-exb| of the exponents (base 2) is sufficient high, 00075 // precision and accuracy can be choosen greater 19. 00076 { 00077 int stagsave = stagprec, stagcalc, stagmax = 19; 00078 l_real a,b,r,r1; 00079 int exa,exb,ex; 00080 a = x; b = y; 00081 exa = expo(a[1]); 00082 exb = expo(b[1]); 00083 if (exb > exa) 00084 { // Permutation of a,b: 00085 r = a; a = b; b = r; 00086 ex = exa; exa = exb; exb = ex; 00087 } // |a| >= |b| 00088 if (sign(a[1]) < 0) a = -a; // a == abs(a); 00089 if (sign(b[1]==0)) return a; 00090 00091 // |a| = a >= |b| > 0: 00092 ex = 0; // initialization 00093 if (6*exb < 5*exa-1071) 00094 { // approximation with a + 0.5*b^2/a - b^4/(8*a^3), if the next 00095 // Taylor term b^6/(16*a^5) is less than 2^-1074==minreal; 00096 // minreal is the smallest positive number in IEEE-format. 00097 r1 = (b/a); 00098 r = r1*r1; 00099 times2pown(r,-2); // r = r/4; 00100 r = 1 - r; 00101 r1 *= b; 00102 times2pown(r1,-1); // r = r/2; 00103 r *= r1; 00104 r += a; // r = a + 0.5*b^2/a - b^4/(8*a^3); 00105 } else 00106 { // Scaling ubwards or downwards to achiev optimal accuracy and 00107 // to avoid overflow by accu-reading: 00108 stagcalc = stagprec; 00109 if (stagcalc > stagmax) stagcalc = stagmax; // using sqrt(...) a 00110 stagprec = stagcalc; // higher precision than stagmax makes no sense! 00111 00112 ex = 511 - exa; // scaling factor to avoide overflow by accu-reading 00113 if (ex < 0) 00114 { // sqrt(a^2+b^2) = a*sqrt( 1+(b/a)^2 ) 00115 r = b/a; 00116 r = a*sqrt(1+r*r); // sqrt(...) declared in l_rmath.hpp 00117 } else 00118 { 00119 times2pown(a,ex); // exact scaling with ex >= 0 00120 times2pown(b,ex); // exact scaling eith ex >= 0 00121 dotprecision dot(0.0); 00122 accumulate(dot,a,a); 00123 accumulate(dot,b,b); 00124 r = dot; // r with no overflow! 00125 r = sqrt(r); // sqrt(...) declared in l_rmath.hpp 00126 times2pown(r,-ex); // back-scaling 00127 } 00128 stagprec = stagsave; // restore old stagprec 00129 } 00130 return r; 00131 } // sqrtx2y2(...) 00132 00133 l_real sqrt1px2(const l_real& x) throw() 00134 // Inclusion of sqrt(1+x^2); Blomquist, 13.12.02; 00135 // With stagmax=19 we get about 16*19=304 exact decimal digits. 00136 { 00137 l_real y,t; 00138 int stagsave, stagmax=19; 00139 stagsave = stagprec; 00140 if (stagprec > stagmax) stagprec = stagmax; 00141 if (expo(x[1]) > 260) 00142 { // sqrt(1+x^2) = |x| + 1/(2*|x|) 00143 y = abs(x); 00144 t = 1/y; 00145 times2pown(t,-1); 00146 y += t; 00147 } else y = sqrt(1+x*x); 00148 stagprec = stagsave; 00149 y = adjust(y); 00150 return y; 00151 } 00152 00153 l_real power(const l_real& x, int n) 00154 { 00155 int stagsave = stagprec, stagmax = 19; 00156 long int zhi = 2; 00157 l_real y, neu; 00158 bool neg=false; 00159 00160 if (x == 1.0) 00161 y = x; 00162 else if (n == 0) 00163 y = adjust(_l_real(1.0)); 00164 else 00165 { 00166 if (stagprec < stagmax) 00167 stagprec++; 00168 else 00169 stagprec = stagmax; 00170 00171 if (n == 1) 00172 y = x; 00173 else if (n == 2) 00174 y = x*x; 00175 else 00176 { 00177 if (n < 0) 00178 { 00179 neg = true; 00180 n = -n; 00181 } 00182 00183 // Initialisierung 00184 if (n%2) 00185 y = x; 00186 else 00187 y = 1.0; // Praezision wird bei 1 Mult. auf 00188 // aktuellen Wert gesetzt; 00189 // Berechnung durch binaere Darstellung der n 00190 neu = x*x; 00191 do { 00192 if ((n/zhi)%2) 00193 y *= neu; 00194 zhi += zhi; 00195 if (zhi <= n) // letzte Mult. entfaellt --> schneller 00196 neu *= neu; 00197 } while (zhi <= n); 00198 if (neg) 00199 y = 1/(y); 00200 } 00201 stagprec = stagsave; 00202 y = adjust(y); 00203 } 00204 return y; 00205 } 00206 00207 // real staggered constants (the same as in l_interval.hpp): 00208 l_real Ln2_l_real() throw() // ln(2) 00209 { return mid( Ln2_l_interval() ); } 00210 l_real Ln10_l_real() throw() // ln(10) 00211 { return mid( Ln10_l_interval() ); } 00212 l_real Ln10r_l_real() throw() // 1/ln(10) 00213 { return mid( Ln10r_l_interval()); } 00214 l_real Pid4_l_real() throw() // Pi/4 00215 { return mid( Pid4_l_interval() ); } 00216 l_real Sqrt2_l_real() throw() // sqrt(2) 00217 { return mid( Sqrt2_l_interval() ); } 00218 l_real Sqrt5_l_real() throw() // sqrt(5) 00219 { return mid( Sqrt5_l_interval() ); } 00220 l_real Sqrt7_l_real() throw() // sqrt(7) 00221 { return mid( Sqrt7_l_interval() ); } 00222 l_real Ln2r_l_real() throw() // 1/ln(2) 00223 { return mid( Ln2r_l_interval() ); } 00224 l_real Pi_l_real() throw() // Pi 00225 { return mid( Pi_l_interval() ); } 00226 l_real Pid2_l_real() throw() // Pi/2 00227 { return mid( Pid2_l_interval() ); } 00228 l_real Pi2_l_real() throw() // 2*Pi 00229 { return mid( Pi2_l_interval() ); } 00230 l_real Pid3_l_real() throw() // Pi/3 00231 { return mid( Pid3_l_interval() ); } 00232 l_real Pir_l_real() throw() // 1/Pi 00233 { return mid( Pir_l_interval() ); } 00234 l_real Pi2r_l_real() throw() // 1/(2*Pi) 00235 { return mid( Pi2r_l_interval() ); } 00236 l_real SqrtPi_l_real() throw() // sqrt(Pi) 00237 { return mid( SqrtPi_l_interval() ); } 00238 l_real Sqrt2Pi_l_real() throw() // sqrt(2*Pi) 00239 { return mid( Sqrt2Pi_l_interval() ); } 00240 l_real SqrtPir_l_real() throw() // 1/sqrt(Pi) 00241 { return mid( SqrtPir_l_interval() ); } 00242 l_real Sqrt2Pir_l_real() throw() // 1/sqrt(2*Pi) 00243 { return mid( Sqrt2Pir_l_interval() ); } 00244 l_real Pip2_l_real() throw() // Pi^2 00245 { return mid( Pip2_l_interval() ); } 00246 l_real Sqrt2r_l_real() throw() // 1/sqrt(2) 00247 { return mid( Sqrt2r_l_interval() ); } 00248 l_real Sqrt3_l_real() throw() // sqrt(3) 00249 { return mid( Sqrt3_l_interval() ); } 00250 l_real Sqrt3d2_l_real() throw() // sqrt(3)/2 00251 { return mid( Sqrt3d2_l_interval() ); } 00252 l_real Sqrt3r_l_real() throw() // 1/sqrt(3) 00253 { return mid( Sqrt3r_l_interval() ); } 00254 l_real LnPi_l_real() throw() // ln(Pi) 00255 { return mid( LnPi_l_interval() ); } 00256 l_real Ln2Pi_l_real() throw() // ln(2*Pi) 00257 { return mid( Ln2Pi_l_interval() ); } 00258 l_real E_l_real() throw() // e = exp(1) 00259 { return mid( E_l_interval() ); } 00260 l_real Er_l_real() throw() // 1/e 00261 { return mid( Er_l_interval() ); } 00262 l_real Ep2_l_real() throw() // e^2 00263 { return mid( Ep2_l_interval() ); } 00264 l_real Ep2r_l_real() throw() // 1/e^2 00265 { return mid( Ep2r_l_interval() ); } 00266 l_real EpPi_l_real() throw() // e^Pi 00267 { return mid( EpPi_l_interval() ); } 00268 l_real Ep2Pi_l_real() throw() // e^(2*Pi) 00269 { return mid( Ep2Pi_l_interval() ); } 00270 l_real EpPid2_l_real() throw() // e^(Pi/2) 00271 { return mid( EpPid2_l_interval() ); } 00272 l_real EpPid4_l_real() throw() // e^(Pi/4) 00273 { return mid( EpPid4_l_interval() ); } 00274 l_real EulerGa_l_real() throw() // EulerGamma 00275 { return mid(EulerGa_l_interval() ); } 00276 l_real Catalan_l_real() throw() // Catalan 00277 { return mid( Catalan_l_interval() ); } 00278 00279 } // namespace cxsc 00280