Calculate the square root of the sum of the squares of the prime factors of an input number

Question

This question is related to another one that I have asked here.

I need to find ways to improve the runtime performance of the following piece of C-code.

I want to increase RANGE to around 1G, and still complete the execution within a reasonable time.

---

Auxiliary data structures:

#define RANGE 16000000

uint08 sieve[RANGE]    = {0};
uint32 prime[RANGE/16] = {0}; // only around RANGE/log(RANGE) are primes
uint32 numOfPrimes     =  0 ;

Initialize the auxiliary data structures:

void CalcAuxiliaryData()
{
    uint32 i,j;

    uint32 root = (uint32)sqrt((double)RANGE);

    for (i=2; i<=root; i++)
    {
        if (sieve[i] == 0)
            for (j=i+i; j<RANGE; j+=i)
                sieve[j] = 1;
    }

    for (i=2; i<RANGE; i++)
    {
        if (sieve[i] == 0)
            prime[numOfPrimes++] = i;
    }
}

Calculate the square root of the sum of the squares of the prime factors of an input number:

uint32 CalcDiagonalLen(uint32 n)
{
    uint32 i;

    uint64 square;
    uint32 length;

    if (sieve[n] == 0) // quickly resolve the case of a prime number
        return n;

    square = 0;
    for (i=0; i<numOfPrimes && n>1; i++)
    {
        uint32 p = prime[i];
        uint64 pp = (uint64)p*p;
        while (n%p == 0)
        {
            n /= p;
            square += pp;
        }
    }

    length = (uint32)sqrt((double)square);
    if ((uint64)length*length == square)
        return length;

    return 0; // indicate that the result is not integer
}

Calculate the square root of the sum of the squares of the prime factors of each number:

int main()
{
    uint32 i;

    uint32 diagonal_len;

    CalcAuxiliaryData();

    for (i=2; i<RANGE; i++)
    {
        diagonal_len = CalcDiagonalLen(i);
        if (diagonal_len != 0)
            printf("%u %u\n",i,diagonal_len);
    }

    return 0;
}

---

Function CalcDiagonalLen holds the bottleneck, but any suggestions will be greatly appreciated.

Answer 1

You can use the sieve approach to create an array of square sums. For any new prime that you sieve, you have:

sqsum[m * p] = sqsum[m] + p*p;

You can build that array as you go. The primes need a square sum value, too.

uint8_t sieve[RANGE] = {0};
uint64_t sqsum[RANGE] = {0};

void CalcAuxiliaryData()
{
    uint32_t i, j, k;

    for (i = 2; i <= RANGE; i++) {
        if (sieve[i] == 0) {
            uint64_t ii = (uint64_t) i * i;

            sqsum[i] = ii;

            j = 2; 
            k = i + i;
            while (k < RANGE) {
                sieve[k] = 1;
                sqsum[k] = sqsum[j++] + ii;
                k += i;
            }
        }
    }
}

uint32_t CalcDiagonalLen(uint32_t n)
{
    uint64_t square;
    uint32_t length;

    if (sieve[n] == 0) return n;

    square = sqsum[n];
    length = sqrt(square);
    if ((uint64_t) length * length == square) return length;

    return 0;
}

This approach adds more time to the initialisation, because you can't cut sieving short at sqrt(RANGE), but the sum lookup is fast.

You could also create a lookup array that tells you whether a number is a perfect square or not, but the sqrt function isn't the bottleneck here. And you could make the array itself a direct lookup table for the diagonal length.

Answer 2

As I've seen that my seemingly clever answer runs into serious memory problems sooner or later, I'll post a second answer that promotes a memory efficient, but fast bottom-up approach that basically only requires storage for the sieve of Erathostenes. (I've kept the original representation of one byte per number, but that can be easily compressed by a factor of 16 by using bit arrays and by not storing even numbers.)

Instead of doing the factorisation for each number and enumerating the numbers sequentially, we can vitsit all numbers in another fashion and keep track of our square sums.

We start with a number of 1, a sum of zero and with the first prime, 2. Then we spread out the solution recursively from this point by either going up, which means we consider the number, which is the current number multiplied by the current prime:

n, sum, p -> n*p, sum + p*p, p

Or by going right, which means that we advance the prime number to the next prime (without any evaluation):

n, sum, p -> n, sum, next(p)

The sum of squares is calculated as we go. When n exceeds the ceiling RANGE, recursion stops.

The drawback to this solution is that the numbers are vistied out of order. The first draft code below just printed the hits. The sorting (and assesment of the conjecture) had to be done with an external program.

I've now updated the code so that it keeps a linear, fixed-size array that is later sorted by increasing values of n. That array is checked for overflows. When an overflow occurs, the heuristics RANGE / 1000 must be corrected. The entiries to struct hits are 64-bit ints, which is wasteful for the current RANGE.

As the RANGE increases, it might be worth keeping a red-bleck tree or something similar of the hits, which doesn't require a huge contiguous chunk of memory. A red-black tree would already store the items in a sorted way.

Anyway, here's my new proposal:

#include <stdlib.h>
#include <stdint.h>
#include <stdio.h>
#include <math.h>



#define RANGE 16000000u
#define NHITS (RANGE / 1000)

uint8_t sieve[RANGE] = {0};

struct hit {
    uint64_t n;
    uint64_t l;
};

struct hit hits[NHITS];
size_t nhits = 0;

void init(void)
{
    uint32_t i, k;

    for (i = 2; i*i < RANGE; i++) {
        if (sieve[i] == 0) {
            for (k = i + i; k < RANGE; k += i) sieve[k] = 1;
        }
    }
}

int64_t perfect(uint64_t sq)
{
    uint64_t l = sqrt(sq);
    if (l * l == sq) return l;

    return 0u;
}

int sprawl(uint64_t n, uint64_t sum, uint64_t p)
{
    uint64_t l;

    if (n >= RANGE) return 0;

    l = perfect(sum);
    if (l) {
        if (nhits >= NHITS) {
            fprintf(stderr, "Hits overflow\n");
            exit(1);
        }

        hits[nhits].n = n;
        hits[nhits].l = l;
        nhits++;
    }

    while (p < RANGE) {
        if (sieve[p] == 0) {
            uint64_t nn = n*p;

            if (nn >= RANGE) break;
            sprawl(nn, sum + p*p, p);
        }
        p++;
    }

    return 0;
}

int hitscmp(const void *a, const void *b)
{
    const struct hit *aa = a;
    const struct hit *bb = b;

    return (aa->n > bb->n) - (aa->n < bb->n);
}

int main()
{
    size_t i;

    init();
    sprawl(1, 0, 2);

    qsort(hits, nhits, sizeof(*hits), hitscmp);

    for (i = 0; i < nhits; i++) {
        printf("%16llu%16llu\n", hits[i].n, hits[i].l);
    }

    return 0;
}

That code still processes 16,000,000 numbers in less than 2s and 160,000,000 in less than 20s - sorting doesn't seem to take much of the overall time.

Answer 3

Based on @MOehm's second answer, I have managed to increase RANGE to 4 billion:

Sieve Interface:

#include <limits.h>

// Defined by the user (must be less than 'UINT_MAX')
#define RANGE 4000000000

// The actual length required for the prime-sieve array
#define ARR_LEN (((RANGE-1)/(3*CHAR_BIT)+1))

// Assumes that all entries in 'sieve' are initialized to zero
void Init(char sieve[ARR_LEN]);

// Assumes that 'Init(sieve)' has been called and that '1 < n < RANGE'
int IsPrime(char sieve[ARR_LEN],unsigned int n);

#if RANGE >= UINT_MAX
    #error RANGE exceeds the limit
#endif

Sieve Implementation:

#include <math.h>

#define GET_BIT(sieve,n) ((sieve[(n)/(3*CHAR_BIT)]>>((n)%(3*CHAR_BIT)/3))&1)
#define SET_BIT(sieve,n) sieve[(n)/(3*CHAR_BIT)] |= 1<<((n)%(3*CHAR_BIT)/3)

static void InitOne(char sieve[ARR_LEN],int d)
{
    unsigned int i,j;
    unsigned int root = (unsigned int)sqrt((double)RANGE);

    for (i=6+d; i<=root; i+=6)
    {
        if (GET_BIT(sieve,i) == 0)
        {
            for (j=6*i; j<RANGE; j+=6*i)
            {
                SET_BIT(sieve,j-i);
                SET_BIT(sieve,j+i);
            }
        }
    }
}

void Init(char sieve[ARR_LEN])
{
    InitOne(sieve,-1);
    InitOne(sieve,+1);
}

int IsPrime(char sieve[ARR_LEN],unsigned int n)
{
    return n == 2 || n == 3 || (n%2 != 0 && n%3 != 0 && GET_BIT(sieve,n) == 0);
}

MOehm's Algorithm:

#include <stdio.h>
#include <stdlib.h>

#define MAX_HIT_COUNT (RANGE/20) /* change the division factor if necessary */

typedef struct
{
    unsigned int n;
    unsigned int x;
}
hit;

char sieve[ARR_LEN]     = {0};
hit hits[MAX_HIT_COUNT] = {0};
unsigned int hit_count  =  0;

void Sprawl(unsigned long long n,unsigned long long sum,unsigned long long p)
{
    unsigned long long x = (unsigned long long)sqrt((double)sum);
    if (x*x == sum)
    {
        if (hit_count >= MAX_HIT_COUNT)
        {
            fprintf(stderr,"Hits overflow\n");
            exit(1);
        }
        hits[hit_count].n = (unsigned int)n;
        hits[hit_count].x = (unsigned int)x;
        hit_count++;
    }

    while (p < RANGE)
    {
        if (IsPrime(sieve,(unsigned int)p))
        {
            if (n*p >= RANGE)
                break;
            Sprawl(n*p,sum+p*p,p);
        }
        p++;
    }
}

Test Case:

int hitscmp(const void *a,const void *b)
{
    const hit* aa = (const hit*)a;
    const hit* bb = (const hit*)b;
    return (aa->n > bb->n) - (aa->n < bb->n);
}

int main()
{
    unsigned int i;

    Init(sieve);

    Sprawl(1,0,2);

    qsort(hits,hit_count,sizeof(*hits),hitscmp);

    for (i=0; i<hit_count; i++)
        printf("%12u%12u\n",hits[i].n,hits[i].x);

    return 0;
}