Find the a 4 digit number whose square is 8 digits AND last 4 digits are the original number

Question

From the comments on my answer here, the question was asked (paraphrase):

Write a Python program to find a 4 digit whole number, that when multiplied to itself, you get an 8 digit whole number whose last 4 digits are equal to the original number.

I will post my answer, but am interested in a more elegant solutions concise but easily readable solution! (Would someone new-ish to python be able to understand it?)

Answer 1

Here is a 1-liner solution without any modules:

>>> next((x for x in range(1000, 10000) if str(x*x)[-4:] == str(x)), None)
9376

If you consider numbers from 1000 to 3162, their square gives you a 7 digit number. So iterating from 3163 would be a more optimized because the square should be a 8 digit one. Thanks to @adrin for such a good point.

>>> next((x for x in range(3163, 10000) if str(x*x)[-4:] == str(x)), None)
9376

Answer 2

If you are happy with using a 3rd party library, you can use numpy. This version combines with numba for optimization.

import numpy as np
from numba import jit

@jit(nopython=True)
def find_result():
    for x in range(1e7**0.5, 1e9**0.5):  
        i = x**2
        if i % 1e4 == x:
            return (x, i)

print(find_result())
# (9376, 87909376)

Answer 3

[Almost] 1-liner:

from math import sqrt, ceil, floor
print(next(x for x in range(ceil(sqrt(10 ** 7)), floor(sqrt(10 ** 8 - 1))) if x == (x * x) % 10000))

printing:

9376

Timing:

%timeit next(x for x in range(ceil(sqrt(10 ** 7)), floor(sqrt(10 ** 8 - 1))) if x == (x * x) % 10000)
546 µs ± 32.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

@theausome 's answer (the shortest (character-wise)):

%timeit next((x for x in range(3163, 10000) if str(x*x)[-4:] == str(x)), None)
3.09 ms ± 119 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

@jpp 's answer (the fastest):

import numpy as np
from numba import jit

@jit(nopython=True)
def find_result():
    for x in range(1e7**0.5, 1e9**0.5):  
        i = x**2
        if i % 1e4 == x:
            return (x, i)
%timeit find_result()
61.8 µs ± 1.46 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Answer 4

The following solution is not as readable as other answers. But what it lacks in readability, it gains in efficiency.

The brute force approach checks every number in a given range, making it O(10^n) where n is the number of desired digits (worst if we account for multiplication and conversions).

Instead, we can build the desired numbers from right to left, adding digits as long as the generated number forms the trailing digits of its square. This provides a O(n³) algorithm (see the time complexity section at the bottom).

def get_all_numbers(n, _digits=None):
    # We are provided digits that form a number with desired properties
    digits = [] if _digits is None else _digits

    # Base case: if we attained the number of digits, return those digits
    if len(digits) >= n:
        return [int(''.join(digits))]

    solutions = []

    # Add digits from 0 to 9 and check if the new number satisfies our property
    for x in range(10):
        next_digits = [str(x)] + digits if x else ['0'] + digits
        next_number = int(''.join(next_digits))

        # If it does try adding yet another digit
        if int(str(next_number ** 2)[-len(next_digits):]) == next_number:
            solutions.extend(get_all_numbers(n, next_digits))

    # Filter out solutions with too few digits
    # Necessary as we could have prepended only zeros to an existing solution
    return [sol for sol in solutions if sol >= 10 ** (n - 1)]

def get_numbers(n, sqr_length=None):
    if sqr_length:
        return [x for x in get_all_numbers(n) if len(str(x ** 2)) == sqr_length]
    else:
        return get_all_numbers(n)

get_numbers(4, 8) # [9376]

This is not necessary for small numbers of digits, but allows to solve the problem for bigger inputs, where the brute force solution takes forever.

get_numbers(100) # Outputs in a few milliseconds

Time complexity

For a given number of digit n, there can only exist at most two solutions other than 0 and 1. And any solution is built from a solution for a smaller number of digits.

From that we conclude that despite the recursion, the algorithm takes O(n) steps to find a solution.

Now, each step has to execute some multiplication and conversions. Integer conversions are O(n²) and multiplication in Python uses Karatsuba's algorithm which is less than conversion.

Overall this yields a O(n³) time complexity.

This could be improved by using a linear integer conversion algorithm and would then provide a O(n^(1 + log(3))) complexity.

Answer 5

Here's the one-line implementation, and it excludes 97.24% of candidates:

import itertools
step = itertools.cycle((24, 1, 24, 51))

[x for x in range(3176, 10000, next(step)) if str(x*x)[-4:] == str(x) ]

Call the number abcd. You can optimize by restricting the last two digits cd, there are only 4 legal possibilities, excluding 96% of candidates for cd. Similarly we only need to test 31 <= ab < 100, excluding 31% of candidates for ab. Thus we excluded 97.24%

cd_cands = set((n**2) % 100 for n in range(0,99+1) if ((n**2 % 100) == n))
cd_cands = sorted(cd_cands)
[0, 1, 25, 76]

for ab in range(31,99+1):
    for cd in cd_cands:
        n = ab*100 + cd
        if n**2 % 10**4 == n :
            print("Solution: {}".format(n))
            #break if you only want the lowest/unique solution
... 
Solution: 9376

(Sure you could squeeze that into a one-line list comprehension, but it would be ugly)

Now we can separate the multiple for-loops with the following observations: strictly we only need to start testing at the first legal candidate above 3162, i.e. 3176. Then, we increment by successively adding the steps (100-76, 1-0, 25-1, 76-25) = (24, 1, 24, 51)

import itertools
step = itertools.cycle((24, 1, 24, 51))

abcd = 3176
while abcd < 10**4:
    if abcd**2 % 10000 == abcd:
        print("Solution: {}".format(abcd))
    abcd += next(step)

and again that can be reduced to the one-liner(/two-liner) shown at top.

UPDATE: As @mathmandan shows, we can improve further to step = itertools.cycle((1, 624)), eliminating 99.68% of candidates. And start searching at 625 * 6 = 3750

Answer 6

Here is a dynamic programming version:

We build from right to left, using the knowledge that for a number to square to itself, each less significant digit has to as well (On post, this is the same approach taken by @Olivier Melançon):

def dynamic(currents, tens, goal):
    if tens == goal:
        return [i for i in currents if len(str(i)) == tens]
    else:
        out = []
        for x in currents:
            for i in range(0,10):
                val = x + i *10**tens
                if val**2 % 10**(tens+1) == val:
                    out.append(val)
        currents = out
    tens +=1
    return dynamic(currents, tens, goal)

We call it with the 'current goal', the current tens, and the goal tens:

dynamic([0],0,4)
#[9376]

Works nicely on super large numbers, in less than a second:

dynamic([0],0,100)
#[3953007319108169802938509890062166509580863811000557423423230896109004106619977392256259918212890625,6046992680891830197061490109937833490419136188999442576576769103890995893380022607743740081787109376]

Answer 7

Simple one-liner that uses the modulo operator % instead of strings

print [x for x in range(3163, 10000) if x*x % 10000 == x]
# [9376]

The low end of the range 3163 is the smallest four digit number whose square is an eight digit number.

Answer 8

The solution I came up with is:

# Loop over all 4 digit numbers
for x in range(1000, 10000):
  # Multiply x*x
  square = x*x
  # Convert to a string
  square = str(square)
  # Check if the square is 8 digits long
  if len(square) == 8:
    # Check if the last 4 digets match x
    if square.endswith(str(x)):
      # print the number and it's square
      print('x    = {}\nx**2 = {}'.format(str(x), square))

Which outputs:

x    = 9376
x**2 = 87909376

Answer 9

We only need to test 1 in 625 candidate numbers.

Either solution A:

upper_limit = 10**4
lower_limit = int(10**3.5) + 1
rem = lower_limit % 625
if rem > 0:
    lower_limit = lower_limit - rem + 625
for n in xrange(lower_limit, upper_limit, 625):
    if n % 16 in [1, 15]:
        print {1: n, 15: n+1}[n%16]
        break

or Solution B:

print (1 * (-39) * 16 + 0 * 1 * 625) % 10000

Read on for explanation.

Starting from the brute-force list-comprehension that tests all candidates:

from math import ceil

[n for n in xrange(ceil(10**3.5), 10000) if (n*n) % 10000 == n]

(The ceil rounds the square root of 10**7 up to nearest integer).

(Notice that 10000 is the first number whose square has 9 digits, and the loop will not test that number.)

... or if you'd like to early-terminate as soon as you find a solution:

for n in xrange(ceil(10**3.5), 10000):
    if (n*n) % 10000 == n:
        print n
        break

But consider: we are looking for numbers n**2 - n = n*(n-1) which are divisible by 10**4 = 2**4 * 5**4.

• either n or n-1 is odd; so the other one would have to be divisible by the full 2**4 = 16. Similarly, you can't have both n and (n-1) be divisible by 5;
• so we need either n or (n-1) to be divisible by 5**4 = 625.
• if either one (n or (n-1)) is divisible by both 625 and 16, then that number is divisible by 10000. No such number has four digits, so it must be that either n or (n-1) is divisible by 625, and the other one is divisible by 16.
• We can thus restrict our search space to looking only at multiples of 625 which have four digits; we just have to be careful to remember that the multiple of 625 might be either n or (n-1); the other one has to be divisible by 16.

So:

upper_limit = 10**4
lower_limit = ceil(10**3.5)
rem = lower_limit % 625
if rem > 0:
    lower_limit = lower_limit - rem + 625
for n in xrange(lower_limit, upper_limit, 625):
    if (n-1) % 16 == 0:
        print n
        break
    if (n+1) % 16 == 0:
        print n+1
        break

Or if you test n instead of (n-1), and combine both condition branches into n % 16 in [1, 15], and for compactness, you could print {1: n, 15: n+1}[n%16].

This is Solution A. (Also, you can certainly replace n%16 with n & 0xf if you prefer.)

But wait! All of this can in fact be done using the...

Chinese Remainder Theorem

We want to find n such that either: n = 0 (mod 625) and n - 1 = 0 (mod 16), or: n - 1 = 0 (mod 625) and n = 0 (mod 16).

So in each case we have two equations, with coprime moduli, solved for the same number n:

n = 0 (mod 625) and n = 1 (mod 16), or else n = 1 (mod 625) and n = 0 (mod 16).

Now (in both cases) we would use the Extended Euclidean Algorithm to find m1 and m2 such that 16*m1 + 625*m2 = 1. It turns out that -39*16 + 1*625 = 1, which leads to the Solution B above, from the second case. (Note: the first case would yield instead 625, whose square does end in 0625, but doesn't count as a solution.)

For completeness, here is an implementation for the Extended Euclidean Algorithm. The second and third outputs are the m1 and m2; in our case 1 and -39 in some order.

def extended_gcd(a, b):
    last_remainder, remainder = abs(a), abs(b)
    x, last_x, y, last_y = 0, 1, 1, 0
    while remainder:
        last_remainder, (quotient, remainder) = remainder, divmod(last_remainder, remainder)
        x, last_x = last_x - quotient*x, x
        y, last_y = last_y - quotient*y, y
    return last_remainder, last_x * ((-1)**(a < 0)), last_y * ((-1)**(b < 0))

Answer 10

One liner here, just

print(9376)