1/* 2 * constants - implementation of different constants to arbitrary precision 3 * 4 * Copyright (C) 2013 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 29static __CZ__euler_mascheroni = 0; 30static __CZ__euler_mascheroni_prec = 0; 31 32 33define e(){ 34 local k temp1 temp2 ret eps factor upperlimit prec; 35 36 prec = digits(1/epsilon()); 37 if(__CZ__euler_mascheroni != 0 && __CZ__euler_mascheroni_prec >= prec) 38 return __CZ__euler_mascheroni; 39 if(prec<=20) return 2.718281828459045235360287471; 40 if(prec<=1800){ 41 __CZ__euler_mascheroni = exp(1); 42 __CZ__euler_mascheroni_prec = prec; 43 } 44 45 eps=epsilon(1e-20); 46 factor = 1; 47 k = 0; 48 upperlimit = prec * ln(10); 49 while(k<upperlimit){ 50 k += ln(factor); 51 factor++; 52 } 53 epsilon(eps); 54 temp1 = 0; 55 ret = 1; 56 for(k=3;k<=factor;k++){ 57 temp2 = temp1; 58 temp1 = ret; 59 ret = (k-1) *(temp1 + temp2); 60 } 61 62 ret = inverse( ret * inverse(factorial(factor) ) ) ; 63 __CZ__euler_mascheroni = ret; 64 __CZ__euler_mascheroni_prec = prec; 65 return ret; 66} 67 68 69/* Lupas' series */ 70static __CZ__catalan = 0; 71static __CZ__catalan_prec = 0; 72define G(){ 73 local eps a s t n; 74 eps = epsilon(epsilon()*1e-10); 75 if(__CZ__catalan != 0 && __CZ__catalan >= log(1/eps)) 76 return __CZ__catalan; 77 a = 1; 78 s = 0; 79 t = 1; 80 n = 1; 81 while(abs(t)> eps){ 82 a *= 32 * n^3 * (2*n-1); 83 a /=((3-16*n+16*n^2)^2); 84 t = a * (-1)^(n-1) * (40*n^2-24*n+3) / (n^3 * (2*n-1)); 85 s += t; 86 n += 1; 87 } 88 s = s/64; 89 __CZ__catalan = s; 90 __CZ__catalan_prec = log(1/eps); 91 epsilon(eps); 92 return s; 93} 94 95 96config("resource_debug", resource_debug_level),; 97if (config("resource_debug") & 3) { 98 print "e()"; 99 print "G()"; 100} 101