[seqfan] square unit fractions

[hv at crypt.org]

Apologies for the length of this; there are three parts pertaining to the
same function: one on a conjecture, one on a sequence, and one on calculation.

Consider a set S of positive integers, and define:
  u(S) = sum_{s in S}{1/s}
  p(S) = sum_{s in S}{1/p_s}  the s'th prime
  q(S) = sum_{s in S}{1/s^2}

Conjecture:

It is straightforward to show by a greedy algorithm that for any rational r,
there exists an S with u(S) = r. Conversely, for almost all r there does
not exist an S with p(S) = r: at most one rational with any given
denominator can be reached (and none unless the denominator is squarefree).

Letting \rho = \pi^2/6, q(S) falls in the ranges [0, \rho-1), [1, \rho).
I conjecture that for any rational r within those ranges there exists an S
with q(S) = r. Can anyone prove or disprove this? (This would seem an
interesting result, since squares are notably less dense than primes
within the integers.)

Sequence:

I have calculated as far as I can S: q(S) = r minimizing |S|, for all
r in [0, \rho-1) with denominators up to 50; my results to date can be
seen at:
  https://github.com/hvds/seq/blob/master/unitfraction/res_square

Calculations are complete for all r < 0.92(\rho-1); the incomplete
value of lowest denominator is r = 5/8 ~= 0.969(\rho-1), but some
solution has been found for all of them to give an upper bound.

The cardinality seems like it would be an interesting sequence, but I'm
not sure what the shape of it should be. It seems natural to let the
rationals follow the order of the Farey series, but should it then set
the value to -1 for rationals outside the range?  Or leave those values
out of the sequence?

Calculation:

I have two programs, one that does an exhaustive search to establish
lower bounds on |S|, another heuristic one that searches for upper bounds.

There is a construction similar to that for unit fractions to find
solutions for r = 1/a^2 + 1/b^2; the main trick to speed things up is
to recognise that 1/a^2 + 1/b^2 cannot have numerator or denominator
divisible by an odd power of any prime p == 3 (mod 4). Since the
denominator as constructed is divisible only by small primes, I factorize
that fully; for the numerator I check only whether the binary expansion
ends 110* - for the harder cases, the numerator and denominator can easily
reach hundreds or thousands of digits.

I feel sure that there are further constraints that could speed things
up, but that is all I have found so far.

For the heuristics, I select an initial set that gets close to the
result, then assume any remaining set members will divide the
denominator of what remains. It feels like there should be a way to
use Pythagorean triples to improve that search, but I'm not sure how
to take advantage of that.

Any suggestions for improvement are welcome.

Hugo