Torn Numbers
September 16, 2014
In 1917, Henry Ernest Dudeney published a book Amusements in Mathematics of arithmetic puzzles. Today’s exercise solves puzzle 113 from that book:
A number n is a torn number if it can be chopped into two parts which, when added together and squared, are equal to the original number. For instance, 88209 is a torn number because (88 + 209)2 = 2972 = 88209.
Your task is to write a program to find torn numbers. 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 
In Python.
from itertools import count, islice def torn_bf(): for m in count(1): n = m ** 2 s = str(n) for p in range(1, len(s)): if (int(s[:p]) + int(s[p:])) ** 2 == n: yield n for t in islice(torn_bf(), 20): print tYou can reduce brute force x10 but must solve one quadratic equation for each digit partition (but, probably it can reduce too).
For only AABBB partition only 10000 checks (instead 100000 checks).
a5 :: Double -> Double -> Double -> Double -> Double a5 a1 a2 a3 a4 = 50 - 0.1 * a1 - a2 - 10 * a3 - 0.1 * a4 - 0.1 * sqrt (250000-9990*a2-99900*a3-999*a1) main = do let ds = [0..9] ss = [round (a + 10 * b + 100 * c + 1000 * d) + 10000 * e | a <- ds, b <- ds, c <- ds, d <- ds , let e = round $ a5 a b c d , e `elem` [1..9]] mapM_ print $ filter isTorn ss(10000 complete checks, checking for Torn number using `isTorn` is only taken 1290 times)
(cannot edit posts? :/ )
(arghhh!)
Obviously, computing better the quadratic solution (e.g. I’m using `round` instead one error bound) it can be reduced; also, quadratic domain could be inspected to inspect only regions with possible solutions. But not like trivial…
Java
public class TornNumbers { public static List<Integer> tornNumbers(int max) { List<Integer> tornNumbers = newArrayList(); for (int i = 10; i <= max; i++) { if (isTornNumber(i)) { tornNumbers.add(i); } } return tornNumbers; } public static boolean isTornNumber(Integer number) { boolean isTorn = false; String numberString = String.valueOf(number); for (int i = 1; i < numberString.length(); i++) { Integer a = Integer.valueOf(numberString.substring(0, i)); Integer b = Integer.valueOf(numberString.substring(i)); if ((a + b) * (a + b) == number) { System.out.println(number + ": " + a + " + " + b + " (" + (a + b) + ")"); isTorn = true; } } return isTorn; } @Test public void testTornNumbers() throws Exception { System.out.println(tornNumbers(88209)); } }Perl version – note that I’m avoiding multiplication except in the output!
use strict; my $z = 0; foreach my $n (1..1e6) { $_ = $z += $n+$n - 1; my $c = q(); while( s{^(\d)}{} ) { $c .= $1; next unless $c+$_ eq $n; printf "%10d %20d (%10d,%10d)\n", $n, $n*$n, $c, $_+0; last; } }{-- Computing big Torn numbers (e.g. 393.900.588.225 ) (but not all, e.g. only one for each suffix). Can be improved. --} import Data.Maybe import System.Environment import Data.Real.Constructible p10 :: Construct -> Construct p10 0 = 1 p10 x = 10 * p10 (x - 1) b' :: Construct -> Construct -> Maybe Construct b' a x = if q < 0 then Nothing else Just $ -a + p10 x / 2 - sqrt q / 2 where q = 4*a + p10 (2*x) - 4 * p10 x * a b a x = case b' a x of Nothing -> Nothing Just y' -> if fromIntegral y == y' then Just y else Nothing where y = round y' torns :: Integer -> [(Integer,(Integer,Integer))] torns s = filter ((>1).fst) $ concatMap (\x -> catMaybes [b a (fromIntegral x) >>= \i->return (i * 10^x + round a, (i, round a)) | a <- [0..10^x-1]]) [0..s] main = getArgs >>= mapM_ print . torns . read . head {-- josejuansquare$ time -f "%E" ./s 7 (2025,(20,25)) (88209,(88,209)) (4941729,(494,1729)) (7441984,(744,1984)) (24502500,(2450,2500)) (23804641,(238,4641)) (300814336,(3008,14336)) (493817284,(4938,17284)) (28005264,(28,5264)) (1518037444,(1518,37444)) (20408122449,(20408,122449)) (21948126201,(21948,126201)) (33058148761,(33058,148761)) (35010152100,(35010,152100)) (43470165025,(43470,165025)) (101558217124,(101558,217124)) (108878221089,(108878,221089)) (123448227904,(123448,227904)) (127194229449,(127194,229449)) (152344237969,(152344,237969)) (213018248521,(213018,248521)) (217930248900,(217930,248900)) (249500250000,(249500,250000)) (393900588225,(39390,588225)) --}@JoseJuan: Your solution misses some torn numbers.
{-- All Torn numbers Using relation (a+b)^2 = a + 10^x b --} import Data.Maybe import System.Environment import Data.Real.Constructible import Control.Parallel import Control.Parallel.Strategies p10 :: Construct -> Construct p10 0 = 1 p10 x = 10 * p10 (x - 1) b' :: Construct -> Construct -> Maybe (Construct, Construct) b' a x = if q < 0 then Nothing else Just (w-z,w+z) where q = 4*a + p10 (2*x) - 4 * p10 x * a z = sqrt q / 2 w = -a + p10 x / 2 b x a = case b' a x of Nothing -> [] Just (a, b) -> ww a ++ ww b ww y' = if fromIntegral y == y' then [y] else [] where y = round y' torns s = filter ((>1).fst) $ concat $ (parMap rseq) (\x -> concatMap (\a -> map (\n ->(n*10^x+a,(n,a))) $ b (fromIntegral x) (fromIntegral a)) [0..10^x-1]) [0..s] main = getArgs >>= mapM_ print . torns . read . head {-- josejuansquare$ time -f "%E" ./s 7 +RTS -N2 (100,(10,0)) (81,(8,1)) (10000,(100,0)) (9801,(98,1)) (2025,(20,25)) (3025,(30,25)) (1000000,(1000,0)) (998001,(998,1)) (88209,(88,209)) (494209,(494,209)) (100000000,(10000,0)) (99980001,(9998,1)) (4941729,(494,1729)) (60481729,(6048,1729)) (7441984,(744,1984)) (52881984,(5288,1984)) (24502500,(2450,2500)) (25502500,(2550,2500)) (10000000000,(100000,0)) (9999800001,(99998,1)) (23804641,(238,4641)) (9048004641,(90480,4641)) (300814336,(3008,14336)) (6832014336,(68320,14336)) (493817284,(4938,17284)) (6049417284,(60494,17284)) (1000000000000,(1000000,0)) (999998000001,(999998,1)) (28005264,(28,5264)) (989444005264,(989444,5264)) (1518037444,(1518,37444)) (923594037444,(923594,37444)) (20408122449,(20408,122449)) (734694122449,(734694,122449)) (21948126201,(21948,126201)) (725650126201,(725650,126201)) (33058148761,(33058,148761)) (669420148761,(669420,148761)) (35010152100,(35010,152100)) (660790152100,(660790,152100)) (43470165025,(43470,165025)) (626480165025,(626480,165025)) (101558217124,(101558,217124)) (464194217124,(464194,217124)) (108878221089,(108878,221089)) (448944221089,(448944,221089)) (123448227904,(123448,227904)) (420744227904,(420744,227904)) (127194229449,(127194,229449)) (413908229449,(413908,229449)) (152344237969,(152344,237969)) (371718237969,(371718,237969)) (213018248521,(213018,248521)) (289940248521,(289940,248521)) (217930248900,(217930,248900)) (284270248900,(284270,248900)) (249500250000,(249500,250000)) (250500250000,(250500,250000)) (100000000000000,(10000000,0)) (99999980000001,(9999998,1)) ... --}def toint(x): return 0 if x == '' else int(x) def tornnumber(x): s = str(x) for mid in range(len(s)): a = toint(s[:mid]) b = toint(s[mid:]) if (a + b)**2 == x: print x map(tornnumber, range(1000000))A simple optimization over the suggested solution is to immediately filter out numbers which aren’t perfect squares:
(define (torn? n) (let ((r (sqrt n))) (if (not (integer? r)) #f (let loop ((ten-power 10)) (if (< n ten-power) #f (let ((a (quotient n ten-power)) (b (modulo n ten-power))) (if (= (+ a b) r) #t (loop (* ten-power 10)))))))))On this machine using Kawa, the original solution takes about 4 minutes 45 seconds to check all numbers from 10 to 100 million; substituting this one instead finishes in 45 seconds.
{-- Improved (and parallelized) version of "scanning suffixes" ... (8934133805179209,(89341338,5179209)) (66616807495744,(666168,7495744)) (8434234407495744,(84342344,7495744)) (123456809876544,(1234568,9876544)) (7901234409876544,(79012344,9876544)) ... Find all Torn's numbers with suffixes with certains digits (e.g. from 1 to 7) --} import Data.Maybe import System.Environment import System.IO import Data.Real.Constructible import Control.Concurrent.Async import Control.Concurrent import Control.Monad digits x = mapM_ (\n -> forkIO (b $ fromIntegral n)) [0..x10-1] where x10 = 10^x :: Integer x20 = 10^(2*x) :: Integer cx = fromIntegral x cx10 = fromIntegral x10 cx20 = fromIntegral x20 cx10h = cx10 / 2 b' :: Construct -> Maybe (Construct, Construct) b' a = if q < 0 then Nothing else Just (w-z, w+z) where q = 4 * a * (1 - cx10) + cx20 z = sqrt q / 2 w = cx10h - a b n = case b' (fromIntegral n) of Nothing -> return () Just (a, b) -> ww a >> ww b where ww :: Construct -> IO () ww y' = when (y' > 0 && fromIntegral y == y') $ print (y*x10+n,(y,n)) >> hFlush stdout where y = round y' :: Integer torns a b = mapConcurrently digits [a..b] main = do (a:b:_) <- getArgs torns (read a) (read b)Computing all Torn numbers in Haskell.
Another nice problem. We want a(k+1) + b = (a+b)^2, with k = 9,99,999 etc. Put c = a+b, and substitute c-a for b, rearrange to ak = c(c-1) ie. c(c-1) mod k == 0. Maximum possible value for c is k+1 (when a == k+1 and b == 0). Also, k is a multiple of 9, therefore either c or c-1 is also (since they are coprime), so we don’t have to try every value for c. Results are produced according to tear position, so sort the output:
def torn (n): lim = pow(10,n); for i in range(1,n+1): k = pow(10,i)-1 # 9,99,999,... c = 9; inc = 1 # c or c-1 is a multiple of 9 while c <= k+1: # Using range blows up memory d = c*c if d > lim: break if (d-c)%k == 0: yield(d) c += inc; inc = 9-inc print(sorted(list(torn(15))))Nice to see 1729 in there. We could do some strength reduction etc. but doesn’t seem to gain much in Python.
Yes @matthew I see it also using `n = q^2 = (a + b)^2 = a + 10^x b` then substracting equations `n = a + 10^x b` and `q = a + b` is `n – q = 999..999 b`.
Indeed, and having another look this morning: for a given k, ak = c(c-1) is a Diophantine equation of the simple form Ax^2 + Bx + Cy = 0, which can then be solved using general methods (eg. see http://www.alpertron.com.ar/QUAD.HTM, which also has an handy online solver – for k = 999, we get (with a couple of simple cases) x = 999 u + 297, y = 999 u^2 + 593 u + 88 or x = 999 u + 703, y = 999 u^2 + 1405 u + 494 which gives our 88209 and 494209 solutions respectively for u = 0.
(I cannot post my code here)
Using previous checker for squares code here
Actually, I’m not sure if there is a general method – that solver seems to be doing essentially what my Python program does (it has handy “show derivation steps” feature), though we can generalize the solution found.
C# Solution:
public partial class Form1 : Form
{
public Form1()
{
InitializeComponent();
}
Thread myThread;
private void button1_Click(object sender, EventArgs e)
{
lstTornNumbers.Items.Clear();
label1.Text = "Current number: " + 0;
myThread = new Thread(CalculateTornNumbers);
myThread.Start();
button1.Enabled = false;
button2.Enabled = true;
}
void CalculateTornNumbers()
{
int i = 0;
while(myThread.ThreadState != ThreadState.AbortRequested)
{
label1.Invoke((MethodInvoker)delegate
{
label1.Text = "Current number: " + i;
});
int locResult;
if (int.TryParse(Math.Sqrt(i).ToString(), out locResult))
{
for(int x =1;x<i.ToString().Length;x++)
{
int a = Convert.ToInt32(i.ToString().Substring(0, x));
int b = Convert.ToInt32(i.ToString().Substring(x, i.ToString().Length - x));
if ((a + b) == locResult)
{
lstTornNumbers.Invoke((MethodInvoker) delegate
{
ListViewItem myItem = new ListViewItem(new string[] { i.ToString(), a.ToString(), b.ToString() });
lstTornNumbers.Items.Add(myItem);
lstTornNumbers.Refresh();
});
}
}
}
i++;
}
}
private void button2_Click(object sender, EventArgs e)
{
myThread.Abort();
button1.Enabled = true;
button2.Enabled = false;
}
}
</code
Not efficient but works
print [n for n in range(10**5) for i in range(1, len(str(n))) if (int(str(n)[:i]) + int(str(n)[i:]))**2 == n]python version
uses a couple optimizations:
1. only checks perfect squares
2. uses (n + 1)**2 = n**2 + (2*n + 1) instead of squaring root_n
in the loop
3. if n has d digits, the left part, a, has at most (d+1)//2 digits.
Or, conversely, the right hand part, b, has at least d//2 digits.
4. both parts have to be less than root_n
5. root_n % 9 must be 0 or 1 (see matthew’s solution above)
def torn(end=1e8, parts=False): root_n = 9 n = 81 ff = True pstart = 10 shiftpoint = 1000 while not end or n < end: # keep track of min digits for right part while n > shiftpoint: pstart *= 10 shiftpoint *= 100 p = pstart b = n % p while b < root_n: a = n // p if a + b == root_n: yield (n, a, b, root_n) if parts else n break p *= 10 b = n % p # calculate the next square and it's root # root_n % 9 must be 0 or 1 (see mathew's solution above) if ff: n += 2*root_n + 1 root_n += 1 else: n += 16*root_n + 64 root_n += 8 ff = not ffWithout matthew’s optimization, it would calculate torn(1e15) in about 55 seconds. The optimization dropped that to about 15 seconds. Changing a,b = divmod(n, p) to b = n % p and a = n //p shaved another 2 seconds off.
Profiling showed that the inner “while b < root_n" loop runs an average of 1.1 times per n.
@MIke: nice. Your solution seems to run about the same speed as mine when generating all torn numbers, though it’s faster for generating up to a limit (7 secs vs. 20 secs up to 10^15) – my program generates the actual solutions fairly quickly, but then spends most of its time checking there aren’t any more. Maybe I can improve my bounds checking, but it’s not obvious to me just how.
@Matthew: Looking at your code again, I tried to find better ways to initialize and increment the value of c. For example, if c*c has i+1 digits, then c has at least (i+1)//2 digits. Or c could be initialized to the first value > sqrt(k+1), such that c%9 ==0.
Then it occurred to me that if c(c-1) % k == 0 then c(c-1) % bf(k) == 0, where bf(k) is the biggest factor of k. So, c can be initialized to and incremented by bf(k). This speeds up the search considerably. It now calculates torn(15) in 0.5 seconds, and torn(18) in 2 seconds. (At torn(19), my table of primes runs out so I didn’t go any further.)
Here’s the code:
from utils import primes def bf(n): f = 1 for p in primes: if p*p > n: return n if n != 1 else f while n%p == 0: f = p n //= p raise ValueError('overflowed primes with n = {}'.format(n)) def torn (n): lim = pow(10,n) for i in range(1,n+1): k = pow(10,i)-1 # 9, 99, 999, ... f = bf(k) # f is the largest factor of k c = f # c or c-1 is a multiple of 'f' inc = 1 while c <= k+1: d = c*c if d > lim: break if (d-c)%k == 0: yield(d) c += inc inc = f - incTurns out (pow(10,19)-1)//9 is prime, so there aren’t any torn numbers for i = 19.
@Mike: That’s brilliant – and it makes me realize we can go further still: we want c(c-1) = ak. Suppose d is a divisor of k and suppose it’s coprime to e = k//d, then if we can find x and y such that xe – yd = 1, then we have a = xy and c = xe. But we can solve xe -yd = 1 with the extended Euclidean algorithm (and we can find a solution just when d and e are indeed coprime). So the problem comes down to finding the divisors of k – and there are some tricks we can use when k is of the form 9….9). Here’s a program that doesn’t use any tricks:
def egcd (a,b): x,y,z,w = 1,0,0,1 while b != 0: q,r = divmod(a,b) a,b = b,r x,y,z,w = z,w,x-q*z,y-q*w return a,x,y def torn(k): i = 2 while i <= k: if k%i == 0: (a,x,y) = egcd(i,k//i) if a == 1: while x<=0: x += k//i yield i*i*x*x i += 1 yield (k+1)*(k+1) for n in sorted(list(torn(9999))): print nThe outputs don’t come out in order, but sorting them we get:
The numbers come in pairs corresponding to the 2 solutions to Jose’s quadratic equation – and also to divisor pair d and k//d, so there are some more simplifications that can be done there (and it would be nice to get rid of that special case for 10…0).
OK, so here’s another version:
def egcd (a,b): x,y,z,w = 1,0,0,1 while b != 0: q,r = divmod(a,b) a,b = b,r x,y,z,w = z,w,x-q*z,y-q*w return a,x,y def torn(n): k = pow(10,n)-1 i = 3 while i*i <= k: if k%i == 0: j = k//i (d,x,y) = egcd(i,j) if d == 1: while x<=0: x += j yield i*i*x*x while y<=0: y += i yield j*j*y*y i += 2 yield k*k yield (k+1)*(k+1) for n in range(1,100): print(sorted(list(torn(n))))Still nothing fancy for finding divisors, though now we only check odd numbers (since k is odd) and we only iterate up to sqrt(k), but we get two solutions each time around. Looks like divisors are relevant to todays problem too.
(use 'clojure.math.numeric-tower) (defn tear [n] (let [digits (map #(Integer. (str %)) (seq (str n)))] digits)) (defn untear [ns] (let [total (count ns)] (loop [factor (expt 10 (dec total)) xs ns number 0] (if (empty? xs) number (recur (quot factor 10) (rest xs) (+ (* factor (first xs)) number)))))) (defn sublists [xs] (let [length (count xs)] (loop [xs xs ys '() l (dec length)] (if (== l 0) ys (let [first (take l xs) secnd (drop l xs) yys (cons first (list secnd)) yss (concat (list yys) ys)] (recur xs yss (dec l))))))) (defn torn? [n] (let [subnums (sublists (tear n)) torns (map (fn [[ part1 part2]] (let [untorn1 (untear part1) untorn2 (untear part2) sum (+ untorn1 untorn2) exp (expt sum 2)] (if (== n exp) true false))) subnums)] (some true? torns))) (filter torn? (range 1 1000))