Jump to content

User:Linas/Mangoldt

From Wikipedia, the free encyclopedia

This note intended for User:Dantheox. The code below can compute the von Mangoldt function up to n=10 million in under a minute of cpu time, and higher with some patience (and/or a faster machine). I wrote it for my personal use only, its under-documented and not really meant for general consumption. Should compile cleanly.

I no longer receive personal e-mail due to spam problems. I get more then 10K pieces of spam a day, and the best spam filters are not good enough to whittle this down to less than 10.

Ask if you have questions. --linas

linas@backlot: /home/linas/linas/fractal/tools/src/libs/farey %cat moebius.h

/*

* moebius.h
*
* Return the moebius function of an integer.
* not intended for large integers, works only for small integers
* due to poor-mans factorization algo.
*
* Linas Vepstas Jaunuary 2005
*/
  1. ifdef __cplusplus

extern "C" {

  1. endif

/** classic Moebius mu function */ int moebius_mu (int n);

/** Mertens function, summatory function of mu */ int mertens_m (int n);

/** The number of prime factors of a number */ int liouville_omega (int n);

/** The Liouville lambda function */ int liouville_lambda (int n);

/** The von Mangoldt Lambda function

*  Returns von Mangoldt Lambda for n, which is
*  log(p) if n is a power of prime p, otherwise
*  returns zero. */

long double mangoldt_lambda (int n); long double mangoldt_lambda_cached (int n);

/** The indexed von Mangoldt Lambda function

*  Returns the n'th non-zero von Mangoldt value
*  */

long double mangoldt_lambda_indexed (int n); unsigned int mangoldt_lambda_index_point (int n);

  1. ifdef __cplusplus

};

  1. endif

linas@backlot: /home/linas/linas/fractal/tools/src/libs/farey %cat moebius.c

/*

* moebius.c
*
* Return the moebius function of an integer.
* not intended for large integers, works only for small integers
* due to poor-mans factorization algo.
*
* Linas Vepstas Jaunuary 2005
*/
  1. include <math.h>
  2. include <malloc.h>
  1. include "moebius.h"
  2. include "cache.h"

static int *sieve = NULL; static int sieve_size = 0; static int sieve_max = 0;

  1. define INIT_PRIME_SIEVE(N) \
       if (!sieve || sieve[sieve_max]*sieve[sieve_max] <(N)) {\
               init_prime_sieve(N); \
       }

/* Initialize and fill in a prime-number sieve */ static void init_prime_sieve (int prod) {

       int n, j;
       int nstart;
       int pos;
       if (sieve && sieve[sieve_max]*sieve[sieve_max] > prod) return;
       int max = 1000.0+sqrt (prod);
       if (!sieve)
       {
               sieve = (int *) malloc (8192*sizeof (int));
               sieve_size = 8192;
               sieve_max = 2;
               sieve[0] = 2;
               sieve[1] = 3;
               sieve[2] = 5;
       }
       pos = sieve_max+1;
       nstart = sieve[sieve_max] + 2;
       /* dumb algo, brute-force test all odd numbers against known primes */
       for (n=nstart; n<=max; n+=2)
       {
               for (j=1; ; j++)
               {
                       int p = sieve[j];
                       if (0 == n%p)
                       {
                               break;
                       }
                       if (p*p > n)
                       {
                               /* If we got to here, n must be prime; save it, move on. */
                               sieve[pos] = n;
                               pos ++;
                               if (pos >= sieve_size)
                               {
                                       sieve_size += 8192;
                                       sieve = (int *)realloc (sieve, sieve_size * sizeof (int));
                               }
                               break;
                       }
               }
       }
       sieve_max = pos-1;
  1. if 0
       for (j=0; j<pos; j++)
       {
               printf ("its %d %d\n", j, sieve[j]);
       }
  1. endif

}

/* ====================================================== */

int moebius_mu (int n) {

       if (1 >= n) return 1;
       if (3 >= n) return -1;
       INIT_PRIME_SIEVE(n);
       /* Implement the dumb/simple moebius algo */
       int cnt = 0;
       int i;
       for (i=0; ; i++)
       {
               int k = sieve[i];
               if (0 == n%k)
               {
                       cnt ++;
                       n /= k;
                       /* If still divisible, its a square */
                       if (0 == n%k) return 0;
               }
               if (1 == n) break;
               if (k*k > n)
               {
                       cnt ++;
                       break;
               }
       }
       if (0 == cnt%2) return 1;
       return -1;

}

/* ====================================================== */

long double mangoldt_lambda (int n) {

       if (1 >= n) return 0.0L;
       INIT_PRIME_SIEVE(n);
       /* Implement the dumb/simple factorization algo */
       int i;
       for (i=0; ; i++)
       {
               int k = sieve[i];
               if (0 == n%k)
               {
                       n /= k;
                       while (0 == n%k) n /= k;
                       if (1 == n) return logl ((long double)k);
                       return 0.0L;
               }
               if (k*k > n) return logl ((long double) n);
       }
       return 0.0L;

}

/* ====================================================== */

long double mangoldt_lambda_cached (int n) {

       DECLARE_LD_CACHE (mangoldt_cache);
       if(ld_one_d_cache_check (&mangoldt_cache, n))
       {
               return ld_one_d_cache_fetch(&mangoldt_cache, n);
       }
       else
       {
               long double val = mangoldt_lambda(n);
               ld_one_d_cache_store (&mangoldt_cache, val, n);
               return val;
       }

}

/* ====================================================== */

DECLARE_LD_CACHE (mangoldt_idx_cache); DECLARE_UI_CACHE (mangoldt_pow_cache); static int man_last_val =1; static int man_last_idx =0;

long double mangoldt_lambda_indexed (int n) {

       if(ld_one_d_cache_check (&mangoldt_idx_cache, n))
       {
               return ld_one_d_cache_fetch(&mangoldt_idx_cache, n);
       }
       else
       {
               ui_one_d_cache_check (&mangoldt_pow_cache, n);
               while (1)
               {
                       man_last_val++;
                       long double val = mangoldt_lambda(man_last_val);
                       if (val != 0.0L)
                       {
                               man_last_idx++;
                               ld_one_d_cache_store (&mangoldt_idx_cache, val, man_last_idx);
                               ui_one_d_cache_store (&mangoldt_pow_cache, man_last_val, man_last_idx);
                               if (n == man_last_idx)
                                       return val;
                       }
               }
       }

}

unsigned int mangoldt_lambda_index_point (int n) {

       if(ui_one_d_cache_check (&mangoldt_pow_cache, n))
       {
               return ui_one_d_cache_fetch(&mangoldt_pow_cache, n);
       }
       else
       {
               ld_one_d_cache_check (&mangoldt_idx_cache, n);
               while (1)
               {
                       man_last_val++;
                       long double val = mangoldt_lambda(man_last_val);
                       if (val != 0.0L)
                       {
                               man_last_idx++;
                               ld_one_d_cache_store (&mangoldt_idx_cache, val, man_last_idx);
                               ui_one_d_cache_store (&mangoldt_pow_cache, man_last_val, man_last_idx);
                               if (n == man_last_idx)
                                       return man_last_val;
                       }
               }
       }

}

/* ====================================================== */

int mertens_m (int n) {

       int i;
       int acc = 0;
       for (i=1; i<=n; i++)
       {
               acc += moebius_mu (i);
       }
       return acc;

}

/* ====================================================== */ /* count the number of prime factors of n */

int liouville_omega (int n) {

       int i;
       int acc = 2;
       if (1 >= n) return 1;
       if (2 >= n) return 2;
       INIT_PRIME_SIEVE(n);
       i=0;
       while (1)
       {
               int d = sieve[i];
               if (0 == n%d)
               {
                       acc ++;
                       n /= d;
               }
               else
               {
                       i++;
               }
               if (d*d > n) break;
       }
       return acc;

}

int liouville_lambda (int n) {

       int omega = liouville_omega (n);
       if (0 == omega%2) return 1;
       return -1;

}

/* ====================================================== */

// #define TEST 1

  1. ifdef TEST

int test_moebius(void) {

       int n;
       int have_error = 0;
       int nmax = 40000;
       for (n=1; n<nmax; n++)
       {
               /* Perform a Dirichlet sum */
               int sum = 0;
               int d;
               for (d=1; ; d++)
               {
                       if (2*d > n) break;
                       if (n%d) continue;
                       sum += moebius_mu (d);
                       // printf ("%d divides %d and sum=%d\n", d, n, sum);
               }
               if (1 != n) sum += moebius_mu (n);
               if (0 != sum)
               {
                       printf ("ERROR for moebius mu at n=%d sum=%d \n", n, sum);
                       have_error ++;
               }
       }
       if (0 == have_error)
       {
               printf ("PASS: tested moebius function w/ dirichlet up to %d\n", nmax);
       }
       return have_error;

}

int test_omega(void) {

       int n;
       int have_error = 0;
       int nmax = 40000;
       for (n=1; n<nmax; n++)
       {
               /* Perform a Dirichlet sum */
               int sum = 0;
               int d;
               for (d=1; ; d++)
               {
                       if (2*d > n) break;
                       if (n%d) continue;
                       sum += liouville_lambda (d);
                       // printf ("%d divides %d and sum=%d\n", d, n, sum);
               }
               if (1 != n) sum += liouville_lambda (n);
               if (0 == sum) continue;
               int ns = sqrt (n);
               if (ns*ns != n)
               {
                       printf ("ERROR at liouville lambda at n=%d sum=%d \n", n, sum);
                       have_error ++;
               }
       }
       if (0 == have_error)
       {
               printf ("PASS: tested liouville omega function w/ dirichlet up to %d\n", nmax);
       }
       return have_error;

}

int main() {

       test_omega ();
       test_moebius ();

}

  1. endif

linas@backlot: /home/linas/linas/fractal/tools/src/libs/farey %

linas@backlot: /home/linas/linas/fractal/tools/src/libs/farey %cat cache.h

/* cache.h

* Generic cache management for commonly computed numbers
*
* Linas Vepstas 2005,2006
*/

/* ======================================================================= */ /* Cache management */

typedef struct {

       unsigned int nmax;
       long double *cache;
       char *ticky;
       short disabled;

} ld_cache;


  1. define DECLARE_LD_CACHE(name) \
       static ld_cache name = {.nmax=0, .cache=NULL, .ticky=NULL, .disabled = 0}

/** ld_one_d_cache_check() -- check if long double value is in the cache

*  Returns true if the value is in the cache, else returns false.
*  This assumes a 1-dimensional cache layout (simple aray)
*/

int ld_one_d_cache_check (ld_cache *c, unsigned int n);

/**

* ld_one_d_cache_fetch - fetch value from cache
*/

long double ld_one_d_cache_fetch (ld_cache *c, unsigned int n);

/**

* ld_one_d_cache_store - store value in cache
*/

void ld_one_d_cache_store (ld_cache *c, long double val, unsigned int n);

/* ======================================================================= */ /* Cache management */

typedef struct {

       unsigned int nmax;
       unsigned int *cache;
       char *ticky;
       short disabled;

} ui_cache;

  1. define DECLARE_UI_CACHE(name) \
       static ui_cache name = {.nmax=0, .cache=NULL, .ticky=NULL, .disabled = 0}

/** ui_one_d_cache_check() -- check if long double value is in the cache

*  Returns true if the value is in the cache, else returns false.
*  This assumes a 1-dimensional cache layout (simple aray)
*/

int ui_one_d_cache_check (ui_cache *c, unsigned int n);

/**

* ui_one_d_cache_fetch - fetch value from cache
*/

unsigned int ui_one_d_cache_fetch (ui_cache *c, unsigned int n);

/**

* ui_one_d_cache_store - store value in cache
*/

void ui_one_d_cache_store (ui_cache *c, unsigned int val, unsigned int n);

linas@backlot: /home/linas/linas/fractal/tools/src/libs/farey %cat cache.c

/* cache.c

* Generic cache management for commonly computed numbers
*
* Linas Vepstas 2005,2006
*/
  1. include <stdlib.h>
  1. include "cache.h"

/* =============================================================== */ /** TYPE_NAME##_d_cache_check() -- check if long double value is in the cache

*  Returns true if the value is in the cache, else returns false.
*  This assumes a 1-dimensional cache layout (simple aray)
*/
  1. define CACHE_CHECK(TYPE_NAME,TYPE) \

int TYPE_NAME##_one_d_cache_check (TYPE_NAME##_cache *c, unsigned int n) \{ \

       if (c->disabled) return 0;      \
       if ((n > c->nmax) || 0==n )     \
       {       \
               unsigned int newsize = 1.5*n+1; \
               c->cache = (TYPE *) realloc (c->cache, newsize * sizeof (TYPE))\                c->ticky = (char *) realloc (c->ticky, newsize * sizeof (char))\        \
               unsigned int en;        \
               unsigned int nstart = c->nmax+1;        \
               if (0 == c->nmax) nstart = 0;   \
               for (en=nstart; en <newsize; en++)      \
               {       \
                       c->cache[en] = 0.0L;    \
                       c->ticky[en] = 0;       \
               }       \
               c->nmax = newsize-1;    \
               return 0;       \
       }       \
       \
       return (c->ticky[n]);   \

}

/**

* TYPE_NAME##_d_cache_fetch - fetch value from cache
*/
  1. define CACHE_FETCH(TYPE_NAME,TYPE) \

TYPE TYPE_NAME##_one_d_cache_fetch (TYPE_NAME##_cache *c, unsigned int n) \{ \

       if (c->disabled) return 0.0L;   \
       return c->cache[n];     \

}

/**

* TYPE_NAME##_d_cache_store - store value in cache
*/
  1. define CACHE_STORE(TYPE_NAME,TYPE) \

void TYPE_NAME##_one_d_cache_store (TYPE_NAME##_cache *c, TYPE val, unsigned int n) \ { \

       if (c->disabled) return;        \
       c->cache[n] = val;      \
       c->ticky[n] = 1;        \

}

/* =============================================================== */

  1. define DEFINE_CACHE(TYPE_NAME,TYPE) \
                 CACHE_CHECK(TYPE_NAME,TYPE) \
                 CACHE_FETCH(TYPE_NAME,TYPE) \
                 CACHE_STORE(TYPE_NAME,TYPE)

DEFINE_CACHE(ld, long double) DEFINE_CACHE(ui, unsigned int)

linas@backlot: /home/linas/linas/fractal/tools/src/libs/farey %

linas@backlot: /home/linas/linas/fractal/misc/num-theory %cat series.c

/*

* series.c
*
* Graphs of maclaurin series of totient and other number-theoretic
* arithmetic series.  These are ordinary x-y plots; output is
* ascii list of x-y values
*
* Linas Vepstas December 2004
*/
  1. include <complex.h>
  2. include <math.h>
  3. include <stdio.h>
  1. include "divisor.h"
  2. include "gcf.h"
  3. include "moebius.h"
  4. include "totient.h"

long double totient_series (long double x) {

       long double acc = 0.0;
       long double xp = 1.0;
       int n=1;
       while (1)
       {
               long double term = xp * totient_phi (n);
               acc += term;
               if (term < 1.0e-20*acc) break;
               xp *= x;
               n++;
       }
       return acc;

}

// return the limit as the totient sum goes to x-> 1 void limit (void) {

       long double p = 0.5L;
       long double prev = 0.0;
       int i=1;
       while (1)
       {
               long double x = 1.0L - p;
               long double y = totient_series (x);
               y *= p*p;
               long double guess = y + (y-prev)/3.0L;
               printf ("%d     %Lg     %26.18Lg        %26.18Lg\n", i, x, y,  guess);
               p *= 0.5L;
               i++;
               prev = y;
       }

}

long double divisor_series (long double x) {

       long double acc = 0.0;
       long double xp = 1.0;
       int n=1;
       while (1)
       {
               long double term = xp * divisor (n);
               acc += term;
               if (term < 1.0e-20*acc) break;
               xp *= x;
               n++;
       }
       return acc;

}

long double c_divisor_series (long double x) {

       long double complex acc = 0.0;
       long double complex xi = x*I;
       long double complex xp = x*I;
       int n=1;
       while (1)
       {
               long double complex term = xp * divisor (n);
               acc += term;
               if (cabsl(term) < 1.0e-20*cabsl(acc)) break;
               xp *= xi;
               n++;
       }
       return cabsl (acc);

}

long double c_erdos_series (long double x) {

       long double complex acc = 0.0;
       long double complex xi = x*I;
       long double complex xp = x*I;
       while (1)
       {
               long double complex term = xp / (1.0L-xp);
               acc += term;
               if (cabsl(term) < 1.0e-20*cabsl(acc)) break;
               xp *= xi;
       }
       return cabsl(acc);

}

long double erdos_series (long double x) {

       long double acc = 0.0;
       long double xp = x;
       while (1)
       {
               // long double term = xp / (1.0L-xp);
               long double term = 1.0L/(xp * (xp-1.0L));
               acc += term;
               if (term < 1.0e-20*acc) break;
               xp *= x;
       }
       return acc;

}

long double z_erdos_series (long double x) {

       long double acc = 0.0;
       long double tk = 0.5L;
       long double xp = x;
       while (1)
       {
               long double term = xp / (1-tk);
               acc += term;
               if (term < 1.0e-20*acc) break;
               xp *= x;
               tk *= 0.5L;
       }
       return acc;

}

long double moebius_series (long double x) {

       long double acc = 0.0;
       long double xp = 1.0;
       int n=1;
       while (1)
       {
               long double term = xp * moebius_mu (n);
               acc += term;
               if (xp < 1.0e-18) break;
               xp *= x;
               n++;
       }
       return acc;

}

long double mangoldt_series (long double y) {

       long double acc = 0.0;
       long double x = expl (-y);
       long double xp = x*x;
       int n=2;
       while (1)
       {
               long double term = xp * (mangoldt_lambda_cached(n) - 1.0);
               acc += term;
               // printf ("pnt=%d term=%Lg acc=%Lg\n", n, term, acc);
               if (fabs(term) < 1.0e-16*fabs(acc)) break;
               xp *= x;
               n++;
       }
       return -acc;

}

long double mangoldt_idx_series (long double y) {

       long double acc = 0.0L;
       long double x = expl (-y);
       long double ox = x/(1.0L-x);
       // printf ("y=%Lg  x=%Lg ox=%Lg\n", y,x,ox);
       long double last_xp = 0.0;
       int n=1;
       unsigned int pnt;
       unsigned int last_pnt=2;
       while (1)
       {
               pnt = mangoldt_lambda_index_point (n);
               long double xp = expl (-y*pnt);
               long double term = mangoldt_lambda_indexed(n);
               term = xp * term;
               // printf ("index=%d pnt=%d term=%Lg acc=%Lg\n", n, pnt, term, acc);
               acc += term;
               if (fabs(term) < 1.0e-16*fabs(acc)) break;
               if (2 == pnt)
               {
                       acc -= xp;
                       last_xp = xp;
               }
               else
               {
                       // term = expl(-y*(pnt-last_pnt));
                       term=1.0L;
                       int i;
                       for (i=0; i<pnt-last_pnt; i++)
                       {
                               term *= x;
                       }
                       term = 1.0L - term;
                       term *= last_xp*ox;
                       acc -= term;
                       //printf ("---dex=%d last_xp=%Lg term=%Lg acc=%Lg\n", n, last_xp, term, acc);
                       last_xp = xp;
               }
               last_pnt = pnt;
               n++;
       }
       fprintf (stderr, "last index=%d pnt=%d\n", n, pnt);
       return -acc;

}

int main () {

       int i;
       int nmax = 410;
       long double tp = 0.5;
       // for (i=1; i<=nmax; i++)
       for (i=nmax; i>0; i--)
       {
               long double x = ((double) i)/((double) nmax);

// #define TOTIENT_SERIES

  1. ifdef TOTIENT_SERIES
               long double y = totient_series (x);
               y *= (1.0L-x)*(1.0L-x);
               long double z = 0.607927101 * sin (0.5*M_PI*x);
               long double r = 2.0L*(y/z - 1.0L);
               printf ("%d     %Lg     %26.18Lg        %26.18Lg        %26.18Lg\n", i, x, y, z,r);
  1. endif


// #define DIVISOR_SERIES

  1. ifdef DIVISOR_SERIES
               x = 1.0L-tp;
               long double y = divisor_series (x);
               // y *= (1.0L-x)*(1.0L-x);
               y *= tp;
               printf ("%d     %Lg     %26.18Lg\n", i, x, y);
  1. endif

// #define C_DIVISOR_SERIES

  1. ifdef C_DIVISOR_SERIES
               long double y = c_divisor_series (x);
               long double z = c_erdos_series (x);
               printf ("%d     %Lg     %26.18Lg        %26.18Lg        %26.18Lg\n", i, x, y, z, y-z);
  1. endif

// #define ERDOS_SERIES

  1. ifdef ERDOS_SERIES
               long double y = z_erdos_series (x);
               y *= (1.0L-x);
               printf ("%d     %Lg     %26.18Lg\n", i, x, y);
  1. endif

// #define MOEBIUS_SERIES

  1. ifdef MOEBIUS_SERIES
               long double y = moebius_series (x);
               printf ("%d     %Lg     %26.18Lg\n", i, x, y);
               fflush (stdout);
  1. endif
  1. define MANGOLDT_SERIES
  2. ifdef MANGOLDT_SERIES
               x *= 0.000001;
               // long double y = mangoldt_series (x);
               long double y = mangoldt_idx_series (x);
               y -= 0.337877;
               y += 0.898*x;
               printf ("%d     %Lg     %26.18Lg\n", i, x, y);
               fflush (stdout);
  1. endif


               tp *= 0.5L;
       }

} linas@backlot: /home/linas/linas/fractal/misc/num-theory %