c_cpp/etc/calc/cal/lnseries.cal

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