Simulated annealing (sa) for baroque numbers--a warm up

The idea of applying the simulated annealing method to baroque numbers was with me for a long time. Finally I had invited others on pl.sci.matematyka (on 2007-05-19) to implement it in parallel with me (independently). I wanted to have some company, it'd be interesting. Somehow, after some positive (and negative) feedback, I ended up alone with a first version 2007-July. It was able to repeat the results obtained in 1990-ies. It couldn't do more because I used a 32-bit C++, and I have represented prime powers directly as actual integers, thus limiting them to the range p^e <>32 (I even stayed lazily withing p^e ≤ 2³¹). Then I barely started to work on a more advanced version 2007-Aug before I already had to stop (for reasons not realated directly to the project). Now, (starting near the end of 2008-August but really this September, I hope) I am trying to psyche myself up for another round. This time I am going to remove the said limitation, and I'd like to actually discover new baroque numbers, not known before.

First let me present here a very naive approach, just for the illustration (I would never actually code it; you may program this way only if you enjoy programming for its own sake). Select a natural number topNum, say topNum := 10000 and an array of primes, say:

P := (2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79}

Your space of "vertices" V consists now of natural numbers n, 1 < 28 =" 2²⋅7 is one of them. You want to find baroque numbers, which belong to V. For instance, 28 is one of them; indeed, the sum of divisors of 28 is:

sd(28) = 1+2+4+7+14+28 = 2⋅28

hence brq(28) = 2, which means that 28 is a perfect number.

Now define the penalty function pen(n) as the number of primes p ε P such that there exists natural exponent e for which p^e is a divisor of n but not of sd(n). Observe thatg n is baroque if and only if pen(n) = 0. Thus now we may say that our goal is to find many (preferably all) n ε V for which pen(n) = 0.

The algorithm may start with vertex v(0) := 43 (or any other prime from P). Let's assume that the algorithm has already reached vertex v(n). Selecting the next vertex, v(n+1), involves two stages: we select a candidate, then the candidate is accepted with probability prob(n), or rejected. If it is rejected, then we select another candidate, which is going to be accepted as v(n+1), with the same probability prob(n), or rejected, etc., until certain v(n+1) gets finally accepted.

In order to select v(n+1), the algorithm selects randomly a prime p ε P. The product v(n)⋅p is the first candidate for v(n+1). Then it checks condition v(n)⋅p ≤ topNum. If it is satisfied, then it checks penalty: if pen(v(n)⋅p) < pen(n); then v(n+1) := v(n)⋅p is accepted as the next vertex. If the penalty didn't decrease then we still allow a prob(n) := 1/(1 + 10⋅log(1+n)) chance that v(n)⋅p is accepted. If it is rejected, for one reason or another, but p is a divisor of n then v(n+1) := v(n)/n is tried, i.e. the algorithm checks pen(v(n+1)). If v(n+1) := v(n)/n t is accepted then the task of finding the next vertex is accomplished (and algorithm will look for v(n+2)). Otherwise a new prime p ε P is randomly selected and the described candidate process is repeated, etc, until finally a new v(n+1) gets accepted.

The describes algorithm can be discussed, modified, refined but my goal was just to provide an idea how a simulated annealing might work for baroque numbers. Among the shortcomings of the described version is the necessity of representing v(n) actually and directly (naively) as an integer. This either limits the program to a small range of numbers or forces us to use a multi-precision software library--indeed, baroque numbers ten to be huge. In the next posts I'll describe a bit less naive approach.

Remark Above, I have proposed

prob(n) := 1/(1 + 10⋅log(1+n))
One may also try one of the many other possibilities, e.g.

prob(n) := 10/(10 + n^1/4)