c-sieve/sieves.c

#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <stdbool.h>
#include <math.h>
#include <stdio.h>

/* =========================
   Basic factor structures
   ========================= */

typedef struct {
    uint64_t prime;
    uint64_t power;
} Factor;

typedef struct {
    uint64_t num_primes;
    Factor factors[];
} Factors;

typedef struct {
    uint64_t len;
    uint64_t elems[];
} U64Array;

/* =========================
   Utilities
   ========================= */

static inline void *get(void *buf, size_t pitch, uint64_t i) {
    return (uint8_t*)buf + pitch*i;
}

static inline uint64_t get_factors(uint64_t w){
    return offsetof(U64Array, elems) + w*sizeof(uint64_t);
}

static inline uint64_t get_factorizations(uint64_t w){
    return offsetof(Factors, factors) + w*sizeof(Factor);
}

static inline void factor_append(Factors *f, uint64_t p, uint64_t k){
    f->factors[f->num_primes].prime = p;
    f->factors[f->num_primes].power = k;
    f->num_primes++;
}

/* =========================
   Math helpers
   ========================= */

static uint64_t u64_pow(uint64_t a, uint64_t e){
    uint64_t r=1;
    while(e){
        if(e&1) r*=a;
        a*=a;
        e>>=1;
    }
    return r;
}

static uint64_t u64_nth_root(uint64_t x, uint64_t n){
    return pow((double)x,1.0/n);
}

static uint64_t u64_binom(uint64_t n,uint64_t k){
    if(k>n) return 0;
    uint64_t r=1;
    for(uint64_t i=1;i<=k;i++)
        r=r*(n-k+i)/i;
    return r;
}

static int64_t i64_gcd(int64_t a,int64_t b){
    while(b){
        int64_t t=a%b;
        a=b;
        b=t;
    }
    return a;
}

static int64_t i64_lcm(int64_t a,int64_t b){
    return a/i64_gcd(a,b)*b;
}

static int64_t i64_egcd(int64_t a,int64_t b,int64_t *x,int64_t *y){
    if(!b){
        *x=1; *y=0;
        return a;
    }
    int64_t x1,y1;
    int64_t g=i64_egcd(b,a%b,&x1,&y1);
    *x=y1;
    *y=x1-y1*(a/b);
    return g;
}

uint64_t max_prime_divs(uint64_t max){
    static const uint64_t primorials[]={
        6ULL,30ULL,210ULL,2310ULL,30030ULL,
        510510ULL,9699690ULL,223092870ULL
    };
    for(uint64_t i=0;i<sizeof(primorials)/sizeof(*primorials);i++)
        if(max<primorials[i])
            return i+1;
    return 15;
}

uint64_t max_primes_le(uint64_t max){
    return 1.25506*max/log(max);
}

/* =========================
   Omega sieve
   ========================= */

uint8_t *sieve_omega(uint64_t max){
    uint8_t *buf=calloc(max+1,sizeof(uint8_t));
    if(!buf) return NULL;

    for(uint64_t n=2;n<=max;n++){
        if(buf[n]) continue;
        for(uint64_t m=n;m<=max;m+=n)
            buf[m]++;
    }
    return buf;
}


uint64_t *sieve_largest_factors(uint64_t max){
    uint64_t *buf=calloc(max+1,sizeof(uint64_t));
    if(!buf) return NULL;

    for(uint64_t p=2;p<=max;p++){
        if(buf[p]) continue;
        for(uint64_t m=p;m<=max;m+=p)
            buf[m]=p;
    }
    return buf;
}

uint32_t *sieve_smallest_factors(uint64_t max){
    uint32_t *buf=calloc(max+1,sizeof(uint32_t));
    if(!buf) return NULL;

    uint64_t r=sqrt(max);

    for(uint64_t p=2;p<=r;p++){
        if(buf[p]) continue;
        for(uint64_t m=p*p;m<=max;m+=p)
            if(!buf[m]) buf[m]=p;
    }

    for(uint64_t i=2;i<=max;i++)
        if(!buf[i]) buf[i]=i;

    return buf;
}

void fill_factors_from_smallest(Factors *out,uint64_t n,const uint32_t spf[]){
    out->num_primes=0;

    while(n>1){
        uint64_t p=spf[n];
        uint64_t k=0;
        while(spf[n]==p){
            n/=p;
            k++;
        }
        factor_append(out,p,k);
    }
}

/* =========================
   Totient sieve https://cp-algorithms.com/algebra/phi-function.html
   ========================= */

uint64_t *sieve_phi(uint64_t max){
    uint64_t *phi=malloc((max+1)*sizeof(uint64_t));
    if(!phi) return NULL;

    for(uint64_t i=1;i<=max;i++)
        phi[i]=i;

    for(uint64_t p=2;p<=max;p++){
        if(phi[p]!=p) continue;

        for(uint64_t m=p;m<=max;m+=p)
            phi[m]-=phi[m]/p;
    }

    return phi;
}

/* =========================
   Sigma_0 (divisor count)
   ========================= */

uint64_t *sigma0(uint64_t max){
    uint64_t *d=calloc(max+1,sizeof(uint64_t));
    if(!d) return NULL;

    for(uint64_t i=1;i<=max;i++)
        for(uint64_t j=i;j<=max;j+=i)
            d[j]++;

    return d;
}

/* =========================
   Mobius sieve https://mathoverflow.net/questions/99473/calculating-m%C3%B6bius-function
   ========================= */

int8_t *mobius(uint64_t max){
    int8_t *mu=malloc(max+1);
    if(!mu) return NULL;

    for(uint64_t i=1;i<=max;i++)
        mu[i]=1;

    for(uint64_t p=2;p<=max;p++){
        if(mu[p]==1){
            for(uint64_t m=p;m<=max;m+=p)
                mu[m]*=-1;

            uint64_t p2=p*p;
            for(uint64_t m=p2;m<=max;m+=p2)
                mu[m]=0;
        }
    }

    return mu;
}


uint64_t *sieve_primes(uint64_t max,uint64_t *count){
    uint8_t *is_comp=calloc(max+1,1);
    if(!is_comp) return NULL;

    uint64_t cap=max_primes_le(max);
    uint64_t *primes=malloc(cap*sizeof(uint64_t));

    uint64_t pc=0;

    for(uint64_t p=2;p<=max;p++){
        if(is_comp[p]) continue;

        primes[pc++]=p;

        if(p*p<=max)
            for(uint64_t m=p*p;m<=max;m+=p)
                is_comp[m]=1;
    }

    free(is_comp);

    *count=pc;
    return realloc(primes,pc*sizeof(uint64_t));
}

int main(){

    uint64_t max = 10000000;

    uint64_t pc;
    uint64_t *primes=sieve_primes(max,&pc);

    printf("Primes up to %llu:\n",(unsigned long long)max);
    for(uint64_t i=0;i<pc;i++)
        printf("%llu ",(unsigned long long)primes[i]);
    printf("\n");

    free(primes);

    uint64_t *phi=sieve_phi(max);
    printf("phi(36) = %llu\n",(unsigned long long)phi[36]);
    free(phi);

    int8_t *mu=mobius(max);
    printf("mu(30) = %d\n",mu[30]);
    free(mu);

    return 0;
}