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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: