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