Solution to Puzzle Corner problem N/D 2 (Nov/Dec 2018)
contributed by Tony Bielecki (MIT 1981), Nov 14, 2018

Any 3-smooth number can be expressed as 2^p * 3^q. In this problem, we
are looking for a 3-smooth number slightly larger than N = 2^100,000 *
3^100,000. In order for our solution to remain close to N, one of its
exponents, either p or q, must be larger than 100,000 and the other
must be smaller.

Direct evaluation and comparison of such large numbers may be feasible
but is unnecessary since we need only calculate the ratio of the
number to N, which can be expressed as 2^(p-100,000)*3^(q-100,000),
and find the number having the smallest ratio greater than one. These
ratios can be calculated with double-precision floating-point
arithmetic if one orders the operations to avoid overflow

My method begins with a trial of every suitable p and q, beginning
near 100,000 and working outward from there in each direction as we
search for the minimum ratio. First, p is decremented from 99,999 to 0
while q is incremented from 100,001 as needed to keep the ratio
greater than one. The ratio is updated at each step simply by division
by 2 and, when necessary, multiplication by 3. At each step, the ratio
is compared to the minimum thus far discovered, and any improvement to
the ratio is stored along with its parameters p and q.

When search over the range of p has been completed, the roles of p and
q are swapped; q is decremented from 99,999 to 0 while p is
incremented from 100,001 as needed, and the calculations involve
division by 3 and multiplication by 2.

The minimum 3-smooth number greater than N thus obtained is
2^225743*3^20665, which is greater than N by a factor of approximately
1.0000037.

We are not quite finished. We still need to confirm that the minimum
can be accurately determined within the precision of the floating
point arithmetic we employed. The worst-case computation in our search
involved about 360,000 floating point operations, which in
double-precision incurs an estimated cumulative rounding error of
about 1e-13. This is insignificant on the scale of the ratios
calculated. To prove that we have correctly identified the minimum, we
examine the second-best trial, 2^49492*3^131867, with its ratio of
approximately 1.0000073. The roundoff error is much less than the
difference between the best and second-best ratios.

Epilogue: I wonder why the exponent 100,000 was chosen. It is not even
close to the limits of double-precision arithmetic. With the
sequential method which I used, one can solve the problem for N
perhaps as high as 6^1000,000,000. For that, I obtained
2^1630138897*3^602426621 with a ratio of 1.000000000104. At that
point, the problem is skimming just an order of magnitude over the
estimated roundoff error, so I don't have a great deal of confidence
in that result.

With some effort to minimize the floating-point operation count (e.g.,
computing high powers by successive squaring rather than stepwise
multiplication), one could reduce the rounding error and take N even
higher, maybe 6^1000,000,000,000. That's just with regular
double-precision processing, nothing really exotic. However, it would
take almost forever, and I am not about to try it.

The problem may be solved in a similar way with logarithms. Rather
than multiplying and dividing by twos and threes, we add and subtract
log(2) and log(3). The overall program flow is the same, as is the
result. The floating-point error propagation may be different; I have
not looked into that.

# Example Python3 program:
# 3-smooth-ratio.py
# python3 program to solve Puzzle Corner problem N/D 2 (Nov/Dec 2018)
# solution by Tony Bielecki, Nov 14, 2018

# power: the exponent used to define N: N=6**power
power=100000

# the initializing value in these 2 lines must be somewhat greater than 1
best_ratio = 2.0
second_best_ratio = 2.0

# these 4 lines not needed, just a formal declaration
best_p = 0
best_q = 0
second_best_p = 0
second_best_q = 0

# first search: decreasing powers of 2, increasing powers of 3

# initialize ratio and exponents p and q (equivalent to starting at N)
ratio = 1.0
p=power
q=power

# scan all p from power to 0
while p > 0:
    ratio = ratio/2
    p=p-1
    while ratio <= 1:
        # scale up ratio as needed
        ratio = ratio*3
        q=q+1
    if ratio < best_ratio:
        # demote previous best parameters
        second_best_ratio = best_ratio
        second_best_p = best_p
        second_best_q = best_q
        # insert best parameters
        best_ratio = ratio
        best_p = p
        best_q = q
    elif ratio < second_best_ratio:
        # replace second best parameters
        second_best_ratio = ratio
        second_best_p = p
        second_best_q = q

# Note: the best and second-best parameters from the first search
# are retained (not reinitialized) as we proceed...

# second search: increasing powers of 2, decreasing powers of 3

# reinitialize ratio, p, and q (return to N and go the other way)
ratio = 1.0
p=power
q=power

# scan all q from power to 0
while q > 0:
    Ratio = ratio/3
    q=q-1
    while ratio <= 1:
        # scale up ratio as needed
        ratio = ratio*2
        p=p+1
    if ratio < best_ratio:
        # demote previous best parameters
        second_best_ratio = best_ratio
        second_best_p = best_p
        second_best_q = best_q
        # insert best parameters
        best_ratio = ratio
        best_p = p
        best_q = q
    elif ratio < second_best_ratio:
        # replace second best parameters
        second_best_ratio = ratio
        second_best_p = p
        second_best_q = q

# report the results
print("best p =",best_p)
print("best q =",best_q)
print("best ratio =",best_ratio)
print("\nsecond best p =",second_best_p)
print("second best q =",second_best_q)
print("second best ratio =",second_best_ratio)
print("second best minus best ratio =",second_best_ratio-best_ratio)
# end of program