Calculate the square root of the sum of the squares of the prime factors of an input number 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&lt;=root; i++)
    {
        if (sieve[i] == 0)
            for (j=i+i; j&lt;RANGE; j+=i)
                sieve[j] = 1;
    }

    for (i=2; i&lt;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&lt;numOfPrimes &amp;&amp; n&gt;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&lt;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.
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 &lt;= RANGE; i++) {
        if (sieve[i] == 0) {
            uint64_t ii = (uint64_t) i * i;

            sqsum[i] = ii;

            j = 2; 
            k = i + i;
            while (k &lt; 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.