Return to Revisions

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

Sieve Interface:

#include <limits.h>

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

// 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/15) /* 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;
}