# Find 7 digit prime numbers with this property;

When you subtract the sum of the squares of the digits of the number from the original number it gives you another prime number squared.

## migrated from meta.mathoverflow.netJul 30 '13 at 1:26

This question came from our discussion, support, and feature requests site for professional mathematicians.

• No offense, but I think it's bad policy to change the question (from 5 to 7-digit numbers) once answers have been attempted. If you wish to know the answer for 7-digit numbers, then post a separate question, or run the code with modified parameters. – Siva Jul 30 '13 at 2:00

The primes from $1000$ to $10^8$ with that property (and the corresponding root) are:

(2903,53)
(3533,59)
(3803,61)
(5197,71)
(9533,97)
(18973,137)
(24763,157)
(37321,193)
(73561,271)
(96953,311)
(113621,337)
(124777,353)
(129097,359)
(134837,367)
(139241,373)
(398341,631)
(830003,911)
(1100509,1049)
(1585201,1259)
(1661789,1289)
(2211257,1487)
(4541309,2131)
(4871077,2207)
(4897709,2213)
(5340949,2311)
(5958751,2441)
(7393123,2719)
(8185501,2861)
(8744003,2957)
(11485559,3389)
(15343039,3917)
(15343079,3917)
(16008193,4001)
(16459441,4057)
(16736519,4091)
(17648693,4201)
(17901563,4231)
(21077401,4591)
(22915511,4787)
(22915591,4787)
(25877861,5087)
(26553629,5153)
(28058447,5297)
(28058467,5297)
(30261149,5501)
(32137747,5669)
(33443269,5783)
(34304561,5857)
(36954461,6079)
(40360753,6353)
(43073159,6563)
(43309751,6581)
(47375963,6883)
(47596441,6899)
(48762557,6983)
(53743763,7331)
(56505481,7517)
(62141909,7883)
(65496941,8093)
(67782521,8233)
(76510171,8747)
(80982263,8999)
(82864927,9103)
(83302253,9127)
(84401197,9187)
(85027057,9221)
(86360003,9293)
(90992939,9539)
(95394619,9767)
(96373819,9817)

module Main (main) where

import Math.NumberTheory.Primes
import Math.NumberTheory.Powers

import System.Environment (getArgs)

main :: IO ()
main = do
args <- getArgs
let (lo,hi) = case args of
_     -> (1000, 10^7)
mapM_ print $funPrimes lo hi digitSquareSum :: Integer -> Integer digitSquareSum = go 0 where go acc 0 = acc go acc m = case m quotRem 10 of (q,r) -> go (acc + r*r) q rootPrime :: Integer -> Maybe Integer rootPrime n = do let m = n - digitSquareSum n r <- exactSquareRoot m guard (isPrime r) return r funPrimes :: Integer -> Integer -> [(Integer, Integer)] funPrimes low high = takeWhile (<= high) (sieveFrom low) >>= \p -> case rootPrime p of Just r -> [(p,r)] Nothing -> [] Between$10^8$and$10^9$, we have (105411479,10267) (110271107,10501) (111492631,10559) (111492671,10559) (116014531,10771) (116014571,10771) (117787913,10853) (117787993,10853) (117961573,10861) (120846227,10993) (134258869,11587) (148767061,12197) (149793487,12239) (156325123,12503) (161010841,12689) (161366489,12703) (165199973,12853) (166900837,12919) (167780441,12953) (169078339,13003) (169234399,13009) (171846139,13109) (184878781,13597) (187334143,13687) (188485741,13729) (199007723,14107) (207043433,14389) (207677111,14411) (208023083,14423) (209641673,14479) (214300417,14639) (218182633,14771) (221384893,14879) (235039859,15331) (235899239,15359) (238486597,15443) (242082611,15559) (247338829,15727) (258020197,16063) (258791927,16087) (258791987,16087) (262083953,16189) (267944477,16369) (270701371,16453) (275792723,16607) (276191479,16619) (286320427,16921) (286320487,16921) (288354637,16981) (288762337,16993) (312193733,17669) (312688781,17683) (316519933,17791) (319730353,17881) (328298413,18119) (343694849,18539) (357928927,18919) (359443943,18959) (365612881,19121) (368909131,19207) (371217443,19267) (371217463,19267) (393744983,19843) (399720367,19993) (417916507,20443) (424319009,20599) (441798719,21019) (454414673,21317) (462981581,21517) (472497401,21737) (475371007,21803) (486776309,22063) (486952877,22067) (487482577,22079) (490932931,22157) (500282917,22367) (503598839,22441) (508187159,22543) (508457681,22549) (511619333,22619) (511709977,22621) (520615621,22817) (520615681,22817) (524730821,22907) (532271219,23071) (562591219,23719) (565155733,23773) (575088641,23981) (575088799,23981) (622153489,24943) (648364667,25463) (658281941,25657) (671794927,25919) (671794987,25919) (702939491,26513) (726464467,26953) (739350793,27191) (744362309,27283) (785737283,28031) (808435721,28433) (811851211,28493) (863360941,29383) (864419033,29401) (878589319,29641) (940097269,30661) (949194871,30809) (959946707,30983) (968641477,31123) (974876041,31223) (978251059,31277) (990801919,31477)$111$more. So these beasts are rare, but there are enough. • Enough for what Daniel ? – Kelly Henrehan Jul 30 '13 at 10:41 • Enough for finding a bunch without starving while one waits. – Daniel Fischer Jul 30 '13 at 10:46 • Are you going to write a paper about them ? – Kelly Henrehan Jul 30 '13 at 10:48 • I don't think so. They're nice curiosities, but I don't see what one could learn from them. – Daniel Fischer Jul 30 '13 at 10:54 • Have you heard of the book Prime curios by G.L. Hanaker . – Kelly Henrehan Jul 30 '13 at 10:56 EDIT: Okay, here's my solution: for i in primes 10000 100000;do x=$(echo $i|python -c 'import sys, math;s = sys.stdin.read()[:-1]; p = math.sqrt(int(s) - sum(map(lambda x : int(x) * int(x), list(s)))); if p % 1 == 0: print str(int(p))') if [ "$x" != "" ];then
if  [ $(factor$x|awk '{ print NF }') -eq "2" ];then
echo $i fi fi done Ugly code, but I get: 18973 24763 37321 73561 96953 Let's check 37321. We get 37321 - 3² - 7² - 3² - 2² - 1², whose square root is 193, which is prime. For 7 digits, here's what I have so far: 1100509 … Let's check 1100509. We get sqrt(1100509 - 1 - 1 - 25 - 81) = 1049, which is prime. Feel free to expend more CPU cycles to get more. • Now these primes should be classified under your and my name , – Kelly Henrehan Jul 30 '13 at 1:45 • 124 is not a prime and may numbers are missing. – Vedran Šego Jul 30 '13 at 1:46 • Many numbers are missing true – Kelly Henrehan Jul 30 '13 at 1:49 • Ah, I made a mistake. Let me fix it. – kumanna Jul 30 '13 at 1:50 • I think I'll stick to mine. And you really shouldn't switch it to 7 digits all of a sudden. Not fair! – kumanna Jul 30 '13 at 1:58 My previous solution was not checking that the numbers are prime, just that$n - \!\!\!\!\!\!\!\!\sum\limits_{\text{digits $d$ of $n$}} d^2\$ are.

So, in Mathematica:

n = 7;
Select[Range[10^(n - 1) + 1, 10^n - 1, 2],
And[
PrimeQ[#1],
PrimeQ[Sqrt[#1 - Apply[Plus, IntegerDigits[#1]^2]]]
] &
]

Output:

{1100509, 1585201, 1661789, 2211257, 4541309, 4871077, 4897709,
5340949, 5958751, 7393123, 8185501, 8744003}
• Now that many examples of this type of prime number have been fleshed out , now what .... somebody writes a paper about these animals ? – Kelly Henrehan Jul 30 '13 at 9:38
• For a paper to be accepted in any decent journal, there must be a bit more. It is very welcome that there is some use for them, or at least something interesting to learn. One cannot just make up a formula or a rule, find a bunch of numbers that fit it, and expect to publish a paper. A few minutes of programming is usually far from enough for a paper. – Vedran Šego Jul 30 '13 at 11:18
• But it's a start , and you have to start somewhere Vedran. – Kelly Henrehan Aug 2 '13 at 12:56