## 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.

Pages: 1 2

### 27 Responses to “Torn Numbers”

1. Paul said

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 t
```
2. You 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
```
3. (10000 complete checks, checking for Torn number using `isTorn` is only taken 1290 times)
(cannot edit posts? :/ )

4. (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…

5. Andras said

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)) {
}
}
}

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));
}
}
```
6. 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;
}
}
```
7. ```{--
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))
--}
```
8. programmingpraxis said

@JoseJuan: Your solution misses some torn numbers.

9. ```{--
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))
...
--}
```
10. Daniel Näslund said
```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))
```
11. Jamie Hope said

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.

12. ```{--
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

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
```
13. svenningsson said

Computing all Torn numbers in Haskell.

```torn = [ a | n <- [4..], let a = n*n, (l,r) <- splitNum a, n == l + r]

splitNum a = takeWhile ((>0).fst) [ divMod a (10^n) | n <- [2..]]
```
14. matthew said

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))))
```
```# \$ time ./torn.py
# [81, 100, 2025, 3025, 9801, 10000, 88209, 494209, 998001, 1000000, 4941729, 7441984, 23804641, 24502500, 25502500, 28005264, 52881984, 60481729, 99980001, 100000000, 300814336, 493817284, 1518037444, 6049417284, 6832014336, 9048004641, 9999800001, 10000000000, 20408122449, 21948126201, 33058148761, 35010152100, 43470165025, 101558217124, 108878221089, 123448227904, 127194229449, 152344237969, 213018248521, 217930248900, 249500250000, 250500250000, 284270248900, 289940248521, 371718237969, 393900588225, 413908229449, 420744227904, 448944221089, 464194217124, 626480165025, 660790152100, 669420148761, 725650126201, 734694122449, 923594037444, 989444005264, 999998000001, 1000000000000, 19753082469136, 24284602499481, 25725782499481, 30024405179209, 30864202469136, 66616807495744, 87841600588225, 99999980000001, 100000000000000, 123456809876544, 186086811780496, 275246813838096, 371449415558529, 390974415863329, 612685018625625, 637690018875625, 953832821345856]

# real	0m22.424s
# user	0m22.360s
# sys	0m0.051s
```

Nice to see 1729 in there. We could do some strength reduction etc. but doesn’t seem to gain much in Python.

15. 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`.

16. matthew said

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.

17. jose.juan said

(I cannot post my code here)
Using previous checker for squares code here

18. matthew said

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.

19. Filip said

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```

20. Mark said

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]
```
21. Mike said

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 ff
```

Without 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.

22. matthew said

@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.

23. Mike said

@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 - inc
```
24. Mike said

Turns out (pow(10,19)-1)//9 is prime, so there aren’t any torn numbers for i = 19.

25. matthew said

@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 n
```

The outputs don’t come out in order, but sorting them we get:

```4941729
7441984
24502500
25502500
52881984
60481729
99980001
100000000
```

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).

26. matthew said

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.

27. Christophe said
```(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))
```