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 16000000uint08 sieve[RANGE]    = {0};uint32 prime[RANGE/16] = {0}; // only around RANGE/log(RANGE) are primesuint32 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.
Calculate the square root of the sum of the squares of the prime factors of an input number
optimization;c
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 zerovoid 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#endifSieve 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;}