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 */
- ifdef __cplusplus
extern "C" {
- 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);
- ifdef __cplusplus
};
- 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 */
- include <math.h>
- include <malloc.h>
- include "moebius.h"
- include "cache.h"
static int *sieve = NULL; static int sieve_size = 0; static int sieve_max = 0;
- 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;
- if 0
for (j=0; j<pos; j++)
{
printf ("its %d %d\n", j, sieve[j]);
}
- 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
- 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 ();
}
- 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;
- 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;
- 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 */
- include <stdlib.h>
- 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) */
- 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 */
- 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 */
- 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; \
}
/* =============================================================== */
- 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 */
- include <complex.h>
- include <math.h>
- include <stdio.h>
- include "divisor.h"
- include "gcf.h"
- include "moebius.h"
- 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
- 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);
- endif
// #define DIVISOR_SERIES
- 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);
- endif
// #define C_DIVISOR_SERIES
- 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);
- endif
// #define ERDOS_SERIES
- ifdef ERDOS_SERIES
long double y = z_erdos_series (x);
y *= (1.0L-x);
printf ("%d %Lg %26.18Lg\n", i, x, y);
- endif
// #define MOEBIUS_SERIES
- ifdef MOEBIUS_SERIES
long double y = moebius_series (x);
printf ("%d %Lg %26.18Lg\n", i, x, y);
fflush (stdout);
- endif
- define MANGOLDT_SERIES
- 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);
- endif
tp *= 0.5L;
}
} linas@backlot: /home/linas/linas/fractal/misc/num-theory %