Based on @MOehm's [second answer](http://codereview.stackexchange.com/a/75737/34177), I have managed to increase `RANGE` to 2 billion and even 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); } **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 int x = (unsigned int)sqrt((double)sum); if ((unsigned long long)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 = 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:** 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("%10u%10u\n",hits[i].n,hits[i].x); return 0; }