User:Linas/Mangoldt

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 */

extern "C" {
 * 1) ifdef  __cplusplus
 * 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);

}; linas@backlot: /home/linas/linas/fractal/tools/src/libs/farey %cat moebius.c
 * 1) ifdef  __cplusplus
 * 1) endif

/* * 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 
 * 2) include 


 * 1) include "moebius.h"
 * 2) include "cache.h"

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

if (!sieve || sieve[sieve_max]*sieve[sieve_max] <(N)) {\ init_prime_sieve(N); \ }
 * 1) define 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;

for (j=0; j= 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 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 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 ; } linas@backlot: /home/linas/linas/fractal/tools/src/libs/farey %
 * 1) endif

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;

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

/** 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;

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

/** 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 


 * 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) */ 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 cache[en] = 0.0L;    \ c->ticky[en] = 0;      \ }      \                c->nmax = newsize-1;    \ return 0;      \ }      \        \        return (c->ticky[n]);   \ }
 * 1) define CACHE_CHECK(TYPE_NAME,TYPE) \

/** * TYPE_NAME##_d_cache_fetch - fetch value from cache */ TYPE TYPE_NAME##_one_d_cache_fetch (TYPE_NAME##_cache *c, unsigned int n)     \{       \ if (c->disabled) return 0.0L;  \ return c->cache[n];    \ }
 * 1) define CACHE_FETCH(TYPE_NAME,TYPE) \

/** * TYPE_NAME##_d_cache_store - store value in cache */ 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 CACHE_STORE(TYPE_NAME,TYPE) \

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

CACHE_CHECK(TYPE_NAME,TYPE) \ CACHE_FETCH(TYPE_NAME,TYPE) \ CACHE_STORE(TYPE_NAME,TYPE)
 * 1) define DEFINE_CACHE(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 
 * 2) include 
 * 3) include 


 * 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; i0; i--) {               long double x = ((double) i)/((double) nmax);

// #define 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);
 * 1) ifdef TOTIENT_SERIES

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 x = 1.0L-tp; long double y = divisor_series (x); // y *= (1.0L-x)*(1.0L-x); y *= tp;
 * 1) ifdef DIVISOR_SERIES

printf ("%d    %Lg     %26.18Lg\n", i, x, y);
 * 1) endif

// #define C_DIVISOR_SERIES long double y = c_divisor_series (x); long double z = c_erdos_series (x);
 * 1) ifdef C_DIVISOR_SERIES

printf ("%d    %Lg     %26.18Lg        %26.18Lg        %26.18Lg\n", i, x, y, z, y-z);
 * 1) endif

// #define ERDOS_SERIES long double y = z_erdos_series (x); y *= (1.0L-x);
 * 1) ifdef ERDOS_SERIES

printf ("%d    %Lg     %26.18Lg\n", i, x, y);
 * 1) endif

// #define MOEBIUS_SERIES long double y = moebius_series (x);
 * 1) ifdef MOEBIUS_SERIES

printf ("%d    %Lg     %26.18Lg\n", i, x, y); fflush (stdout);
 * 1) endif

x *= 0.000001; // long double y = mangoldt_series (x); long double y = mangoldt_idx_series (x); y -= 0.337877; y += 0.898*x;
 * 1) define MANGOLDT_SERIES
 * 2) ifdef MANGOLDT_SERIES

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 %