Approximate Squaring
June 8, 2021
Lagarias and Sloane study the “approximate squaring” map f(x) = x⌈x⌉ and its behavior when iterated in this paper.
Consider the fraction x = n / d when n > d > 1; let’s take 8/7 as an example. In the first step, the smallest integer greater than 8/7 (the “ceiling”) is 2, and 8/7 × 2 = 16/7. In the second step, the ceiling of 16/7 is 3, so we have 16/7 × 3 = 48/7. And in the third step, the ceiling of 48/7 is 7, so we have 48/7 × 7 = 48. Now the denominator is 1 and the result is an integer, so iteration stops, and we say that 8/7 goes to 48 in 3 steps.
Study shows the iteration is chaotic; sometimes the iteration stops in just a few steps, as in 8/7, sometimes it takes longer (6/5 goes to a number with 57735 digits in 18 steps), and sometimes it’s just ridiculous (200/199 goes to a number with 10435 digits). It is conjectured but not proven that iterated approximate squaring always terminates in an integer.
Your task is to write a program that iterates approximate squaring. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
Nifty little drill. I suppose some people have way too much time on their hands, writing papers like this!
Here is my take on this using Julia: https://pastebin.com/TtiyB0Wu
Here’s a Haskell version.
import Data.Ratio (denominator, numerator) import System.Environment (getArgs) approxSq :: Rational -> Rational approxSq x = x * toRational (ceiling x :: Integer) iterApproxSq :: Rational -> Integer iterApproxSq = numerator . head . iterateUntil isInteger approxSq where isInteger = (== 1) . denominator iterateUntil p f = dropWhile (not . p) . drop 1 . iterate f main :: IO () main = do rs <- getArgs mapM_ (print . iterApproxSq . read) rsLisp:
(defun approximate-square (number)
(loop
for n = number then (* n (ceiling n))
when (integerp n)
do (return n)))
Here’s a solution in Python.
def approx_squaring(n, d): while n % d: n *= (n + d - 1) // d return n // d print(approx_squaring(8, 7)) print(approx_squaring(10, 6)) print(approx_squaring(6, 5))Output:
Here’s a solution in C using GMP.
// $ gcc -o approx_squaring approx_squaring.c -lgmp #include <stdlib.h> #include <gmp.h> int main(int argc, char* argv[]) { if (argc != 3) return EXIT_FAILURE; mpz_t n; mpz_init(n); mpz_set_str(n, argv[1], 10); mpz_t d; mpz_init(d); mpz_set_str(d, argv[2], 10); mpz_t q; mpz_init(q); while (!mpz_divisible_p(n, d)) { mpz_cdiv_q(q, n, d); mpz_mul(n, n, q); } mpz_divexact(n, n, d); gmp_printf("%Zd\n", n); mpz_clear(n); mpz_clear(d); mpz_clear(q); return EXIT_SUCCESS; }Example usage: