1/*
2 * lnseries - special functions (e.g.: gamma, zeta, psi)
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
25/*
26 * hide internal function from resource debugging
27 */
28static resource_debug_level;
29resource_debug_level = config("resource_debug", 0);
30
31
32static __CZ__int_logs;
33static __CZ__int_logs_limit;
34static __CZ__int_logs_prec;
35
36
37define deletelnseries(){
38  free(__CZ__int_logs,__CZ__int_logs_limit,__CZ__int_logs_prec);
39}
40
41define lnfromseries(n){
42  if(    isnull(__CZ__int_logs)
43      || __CZ__int_logs_limit < n
44      || __CZ__int_logs_prec < log(1/epsilon())){
45
46    lnseries(n+1);
47  }
48  return __CZ__int_logs[n,0];
49}
50
51define lnseries(limit){
52  local k j eps ;
53  if(    isnull(__CZ__int_logs)
54      || __CZ__int_logs_limit < limit
55      || __CZ__int_logs_prec < log(1/epsilon())){
56    __CZ__int_logs = mat[limit+1,2];
57    __CZ__int_logs_limit = limit;
58    __CZ__int_logs_prec = log(1/epsilon());
59
60    /* probably still too much */
61    eps = epsilon(epsilon()*10^(-(5+log(limit))));
62    k =2;
63    while(1){
64      /* the prime itself, compute logarithm */
65      __CZ__int_logs[k,0] = ln(k);
66      __CZ__int_logs[k,1] = k;
67
68      for(j = 2*k;j<=limit;j+=k){
69       /* multiples of prime k, add logarithm of k computed earlier */
70       __CZ__int_logs[j,0] += __CZ__int_logs[k,0];
71       /* First hit, set counter to number */
72       if(__CZ__int_logs[j,1] ==0)
73         __CZ__int_logs[j,1]=j;
74       /* reduce counter by prime added */
75       __CZ__int_logs[j,1] //= __CZ__int_logs[k,1];
76      }
77
78      k++;
79      if(k>=limit) break;
80      /* Erastothenes-sieve: look for next prime. */
81      while(__CZ__int_logs[k,0]!=0){
82        k++;
83        if(k>=limit) break;
84      }
85    }
86    /* Second run to include the last factor */
87    for(k=1;k<=limit;k++){
88      if(__CZ__int_logs[k,1] != k){
89        __CZ__int_logs[k,0] +=__CZ__int_logs[ __CZ__int_logs[k,1],0];
90        __CZ__int_logs[k,1] = 0;
91      }
92    }
93
94    epsilon(eps);
95  }
96  return 1;
97}
98
99
100/*
101 * restore internal function from resource debugging
102 */
103config("resource_debug", resource_debug_level),;
104if (config("resource_debug") & 3) {
105    print "lnseries(limit)";
106    print "lnfromseries(n)";
107    print "deletelnseries()";
108}
109