sum of unit fractions

[hv at crypt.org]

hv at crypt.org wrote:
:The approach I used to test optimality of a(N) for N = 1 .. 5 was:
:- set Limit successively to 1, 2, ...
:- calculate Q_{p^k} as the subsequence of 1..Limit for which gpp(n) = p^k
:- set H(Limit) = sum_{i=1..Limit}{ 1/i }
:- set Kept = 0/1, Discarded = 0/1
:- recursively for each p^k such that Q_{p^k} is non-empty, in decreasing order
:  - calculate Req = the required sum modulo p :
:    - if denominator(Kept) is not divisible by p^k then Req = 0
:      else Req = numerator(Kept) ^ -1 (mod p)
:  - for Q_{p^k}, calculate R, the sequence (q_i/p^k) ^ -1 (mod p)
:  - save the values of Kept and Discarded
:  - for each subsequence S of R such that sum(S) == Req (mod p)
:    - for each q_i in Q, add 1/q_i to Kept or Discarded depending
:      on whether the corresponding r_i was used ('kept') in S
:    - if H(Limit) - Discarded = N then terminate successfully
:    - if H(Limit) - Discarded > N then recurse to the next lower prime power
:    - restore Kept and Discarded to the saved values
:  - derecurse
:- and try the next value for Limit

Oops, the calculation of Req is an error; it should be:
     Req = ( - p^k Kept ) (mod p)
which is to say:
   - if denominator(Kept) is not divisible by p^k then Req = 0
     else Req = (- numerator(Kept) * ( denominator(Kept) / p^k ) ^ -1) (mod p)

The danger of believing one's own code: due to the effect of other
optimisations, this error didn't show up until evaluating a(6).

This improved, and hopefully now correct code gives me a better a(6):
a(6) = 469 : {
  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
  29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
  54 55 56 57 58 60 61 62 63 64 65 66 67 68 69 70 72 74 75 76 77 78 81 82 84
  85 86 87 88 90 91 92 93 94 95 96 98 99 100 102 104 105 106 108 110 111 114
  115 116 117 119 120 121 122 123 124 126 128 129 130 132 133 134 136 138 140
  141 143 145 147 148 150 152 153 155 156 159 161 162 164 168 170 171 174 177
  180 182 183 184 185 186 187 188 190 192 196 201 203 204 205 207 209 210 212
  215 217 221 222 228 230 232 238 242 245 246 247 248 250 252 253 255 258 259
  261 266 268 272 276 280 282 285 287 290 294 295 299 301 304 305 306 310 315
  319 322 323 324 328 329 341 344 345 348 354 360 363 368 370 372 375 376 377
  380 384 387 402 406 407 413 414 424 430 434 435 444 451 460 465 469
}

I'll aim to get a confirmed value for a(7), then submit the actual sequence.

Hugo