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