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 %
Content Disclaimer
Informasi ini disarikan dari Wikipedia dan disajikan kembali untuk tujuan edukasi. Konten tersedia di bawah lisensi CC BY-SA 3.0. Kami tidak bertanggung jawab atas ketidakakuratan data yang bersumber dari kontribusi publik tersebut.
- The information displayed on this website is sourced in part or in whole from Wikipedia and has been adapted for the purpose of restating it. We strive to provide accurate and relevant information, however:
- There is no guarantee of absolute accuracy. Wikipedia is an open, collaborative project that can be edited by anyone, so information is subject to change.
- It is not intended to constitute professional advice. The content displayed is for informational and educational purposes only. For important decisions (e.g., medical, legal, or financial), please consult a professional.
- Content copyright. Wikipedia is licensed under the Creative Commons Attribution-ShareAlike License (CC BY-SA). This means that content may be reused with appropriate attribution and shared under a similar license.
- Responsible use. Any risk arising from the use of information from this website is entirely the responsibility of the user.