C-Code

/********************************************
 * giuga.c    Giuga-Zahlen berechnen
 ********************************************/
#include <stdlib.h>
#include <stdio.h>
 
#define MAX_SIEVE 1000000   // Array Limit Sieb
#define MAX_PRIME 100000    // Array Limit Primzahlen

unsigned char *psieve;
unsigned int *pprime, *ppfact;
int nsieve, nprime, npfact;

int alloc_buffers() {
	size_t bufsiz;
	
	// sieve buffer initialized with 0
	bufsiz = MAX_SIEVE / 8;
	if ((psieve = (unsigned char *)calloc(bufsiz,1)) == NULL) {
	  printf("Error sieve buffer allocation (%d byte)\n", (int)bufsiz);
	  return 1;
    }
	bufsiz = MAX_PRIME * sizeof(int);
	if ((pprime = (unsigned int *)malloc(bufsiz)) == NULL) {
	  printf("Error prime buffer allocation (%d byte)\n", (int)bufsiz);
	  return 1;
    }
	if ((ppfact = (unsigned int *)malloc(bufsiz)) == NULL) {
	  printf("Error factor buffer allocation (%d byte)\n", (int)bufsiz);
	  return 1;
    }
    return 0;
}

void free_buffers() {
	if (psieve != NULL) free(psieve);
	if (pprime != NULL)	free(pprime);
	if (ppfact != NULL)	free(ppfact);
}

void set_sieve(unsigned int number) {
	// sets "NOPRIME = 1"
	unsigned char *p = psieve + number/8;
	*p |= 1 << number%8;
}

int get_sieve(unsigned int number) {
	// if "NOPRIME" return code > 0
	unsigned char *p = psieve + number/8;
	return *p & (1 << number%8);
}

int get_primes(unsigned int max_number) { 
    unsigned int np = 2;
    unsigned int m, p;
 
    set_sieve(0);
    set_sieve(1);
    do {
		m = np + np;
		while(m <= max_number) {
			set_sieve(m);
			m += np;
		}
		while(++np <= max_number && get_sieve(np) > 0)
			;
    } while (np <= max_number);
	// fill primes array
    for (m=p=0; p <= max_number; p++) {
        if (get_sieve(p) == 0) {
			*(pprime + m) = p;
			m++;
		}
	}
	return m;		
}

int factorize(unsigned int number) {
	unsigned int pnum; 
	int np=0, nf=0;

	while (number > 1) {
		pnum = *(pprime + np);
		while (number % pnum == 0) {
			*(ppfact + nf++) = pnum;
			number /= pnum;
		} 
		if (++np >= nprime)
			return 0;
	}
	return nf;
}

int main(int argc, char **argv) {
	unsigned int max_number, number;
	unsigned int pfact, gtest, k;

	if (argc <= 1 || (max_number = atoi(argv[1])) == 0)
		return 1;

	printf("Berechne Primzahlen bis %d...\n", MAX_SIEVE);
	alloc_buffers();
    nprime = get_primes(MAX_SIEVE);
	printf("%d Primzahlen gefunden.\n",nprime);

   	printf("Suche Giuga-Zahlen bis %u\n", max_number);
    for (number=2; number <= max_number; number+=2) {
	    gtest = 1;
	    if ((npfact = factorize(number)) > 0) {
		  gtest = pfact = 0;	
		  for (k=0; k < npfact; k++) {
		    if (pfact < *(ppfact + k)) {
			    pfact = *(ppfact + k);
		        if ((gtest += (number>pfact)?(number/pfact-1)%pfact:1) > 0)
		          break;
  		    }
  		  }  
	    }
        if (gtest == 0) {
			// Treffer
			printf("%u = ", number);
		    for (k=0; k < npfact; k++)
				printf("%u%c", *(ppfact + k), (k<npfact-1)?'*':'\n');
		}	
	}	
	free_buffers();
    return 0;
}