# Finding vampire numbers

How to find vampire numbers by using Mathematica?

A number $v=xy$ with an even number $n$ of digits formed by multiplying a pair of $n/2$-digit numbers (where the digits are taken from the original number in any order) $x$ and $y$ together.

For example:

    1.    1260=21*60
2.    1395=15*93
3.    1530=30*51
4.    1827=21*87


I tried to find all the 4 digit vampire numbers:

FindInstance[{1000 a + 100 b + 10 c + d == (10 m + n) (10 l + k),
a > 0, b > 0, c > 0, d > 0, m > 0, n > 0, k > 0, l > 0,
{a, b, c, d, m, n, l, k} ∈ Integers,
{m, n, l, k} ∈ {a, b, c, d}},
{a, b, c, d, m, n, l, k}]


but there is an error message

The Vampire Numbers that the OP refers to are a simple case of a more general problem of expressing numbers in terms of their own digits ... sometimes called Narcissistic Numbers.

Why write 1296 when you can write:

For a classification of these different types (including vampire Numbers), see for instance:

http://www.tri.org.au/numQ/pwn/comparison.html

Some years ago, I played around with these structures in Mathematica, and wrote a fun little paper on it for the Journal of Recreational Mathematics (volume 33(4), 2004-2005, pp.250-254) on generalising these structures to includes radicals and factorials etc, which I called: Pretty Wild Narcissistic Numbers, or Numbers that PWN.

For the Vampire case, for any 4 digit number abcd, the possible permutations into 2 and 2 numbers are:

perm2 = Map[ Partition[#,2]&, Permutations[{a,b,c,d}]];


The 24 possible products of those couplings are:

tes[{a_, b_, c_, d_}] = Map[ Times@@(Map[FromDigits, #, 1])&, perm2]


{(10 a + b) (10 c + d), (10 a + b) (c + 10 d), (10 a + c) (10 b + d), (10 a + c) (b + 10 d), (10 b + c) (10 a + d), (b + 10 c) (10 a + d), (a + 10 b) (10 c + d), (a + 10 b) (c + 10 d), (10 b + c) (10 a + d), (10 b + c) (a + 10 d), (10 a + c) (10 b + d), (a + 10 c) (10 b + d), (a + 10 c) (10 b + d), (a + 10 c) (b + 10 d), (b + 10 c) (10 a + d), (b + 10 c) (a + 10 d), (10 a + b) (10 c + d), (a + 10 b) (10 c + d), (10 b + c) (a + 10 d), (b + 10 c) (a + 10 d), (10 a + c) (b + 10 d), (a + 10 c) (b + 10 d), (10 a + b) (c + 10 d), (a + 10 b) (c + 10 d)}

The following takes 0.09 seconds to solve:

Cases[
ParallelTable[If[ MemberQ[ tes[IntegerDigits[i]], i], i],  {i, 1000,9999}], _Integer]
// AbsoluteTiming


{0.090674, {1260, 1395, 1435, 1530, 1827, 2187, 6880}}

Consider the trusty DivisorPair function from MrWizard:

DivisorPairs[n_] := Thread[{#,Reverse[#]}][[;;Ceiling[Length[#]/2]]]&[Divisors[n]]


DivisorPairs returns pairs of divisors multiplying together to form n. From these, select those pairs with both divisors having half the number of digits of n. Then check that the sorted digits of n match the sorted digits of one or more of the remaining divisor pairs.

VampireNumberQ[n_]:=
Block[{m, d, p},
m = IntegerLength[n]/2;
d = Sort[IntegerDigits[n]];
p = Select[DivisorPairs[n], IntegerLength[#] == {m, m} &];
Select[p, Sort[Flatten[IntegerDigits[#]]] == d &] =!= {}
]
SetAttributes[VampireNumberQ,Listable]


VampireNumberQ is about 6 times faster than the OEIS fQ.

AbsoluteTiming[With[{r = Range[1000, 9999]}, Pick[r, VampireNumberQ[r]]]]


{0.227699, {1260, 1395, 1435, 1530, 1827, 2187, 6880}}

• May you include the link to MrWizard's DivisorPairs function in your post so it displays as referenced in the links on the right-hand side. Commented Mar 11, 2017 at 14:46

Since FindInstance willl not give us the whole solution set I prefer to use Solve:

{a, b, c, d} /.
Solve[ Or @@ ({1000, 100, 10, 1}.{a, b, c, d} == (10 #1 + #2)(10 #3 + #4)& @@@
Permutations[{a, b, c, d}]) && And @@ Thread[1 <= # <= 9 &@{a, b, c}]&&
0 <= d <= 9, {a, b, c, d}, Integers]

 {{1, 2, 6, 0}, {1, 3, 9, 5}, {1, 4, 3, 5}, {1, 5, 3, 0},
{1, 8, 2, 7}, {2, 1, 8, 7}, {6, 8, 8, 0}}

• +1 for the concept, but you need additional conditions to ensure that the digits {m, n, l, k} are a permutation of the digits {a, b, c, d}. This is why FindInstance returns 1116 even though it is not a vampire number -- notice how it used 36*31. The additional conditions are something like Or @@ Thread[ FromDigits[{a, b, c, d}] == FromDigits /@ Permutations[{m, n, l, k}]]. It may take longer to solve :) Commented Mar 11, 2017 at 13:23
• but there are only seven 4 digit vampire numbers $1260, 1395, 1435, 1530, 1827, 2187, 6880$
– vito
Commented Mar 11, 2017 at 13:24
• @WReach Thanks for pointing out the problem (I haven't grasped at first sight what the OP expects), in fact, a direct addition of your conditions are not recommended, however a slight modification works fine. Commented Mar 11, 2017 at 14:51
• Nice simplification. I'd +1 again if I could :) Should the digit constraints be 1 <= a <= 9 && 0 <= {b, c, d} <= 9 in the more general case? Commented Mar 11, 2017 at 15:09
• @WReach That's right, however more efficient approach would be Select[1000 a + 100 b + 10 c + d /. Solve[Or @@ ({1000, 100, 10, 1}.{a, b, c, d} == (10 #1 + #2) (10 #3 + #4) & @@@ Permutations[{a, b, c, d}]) && 0 <= {a, b, c, d} <= 9, {a, b, c, d}, Integers], # > 1000 &] The more diverse logical conditions we apply the more time Solve takes. Commented Mar 11, 2017 at 17:36

I found the code at oeis.org

fQ[n_] :=
If[OddQ@IntegerLength@n, False,
MemberQ[Map[Sort@Flatten@IntegerDigits@# &,
Select[Map[{#, n/#} &, TakeWhile[Divisors@n, # <= Sqrt@n &]],
SameQ @@ Map[IntegerLength, #] &]], Sort@IntegerDigits@n]];


and it's quite fast

Select[Range[10^5], fQ] // AbsoluteTiming


$\{1.96896,\{1260,1395,1435,1530,1827,2187,6880\}\}$

You may use Interval to define the region for Element (∈).

From Element's documentation note that you need to specify one of the possible domains or provide a region. A list of integers is neither but you can create a one-dimensional region specification using Interval. For example,

NumberLinePlot[Interval[{1, 1}, {2, 2}, {3, 3}, {4, 4}]]


There Interval defines a one-dimensional region that contains numbers 1, 2, 3, and 4. One point to note is that when specifying the element of a one-dimensional region you must provide a 1D vector, {2}, and not a scalar, 2.

Next note that your constraints are incomplete. You need 0 < a <= 9 and the remainder 0 <= x <= 9.

We also need to ensure that all the numbers {a, b, c, b} are used exactly once. The Sort constraint sees to this.

You may use the fourth parameter of FindInstance to return more than one.

Taking this into account and making use of short-hand equivalents for the constraints we get.

res =
FindInstance[
{
1000 a + 100 b + 10 c + d == (10 m + n) (10 l + k),
0 < a <= 9,
0 <= {b, c, d, m, n, l, k} <= 9,
{a, b, c, d, m, n, l, k} ∈ Integers,
{{m}, {n}, {l}, {k}} ∈ Interval @@ Map[{#, #} &, {a, b, c, d}],
Sort@{a, b, c, d} == Sort@{m, n, l, k}
},
{a, b, c, d, m, n, l, k},
4]

{{a -> 1, b -> 2, c -> 6, d -> 0, m -> 6, n -> 0, l -> 2, k -> 1},
{a -> 1, b -> 4, c -> 3, d -> 5, m -> 3, n -> 5, l -> 4, k -> 1},
{a -> 1, b -> 5, c -> 3, d -> 0, m -> 3, n -> 0, l -> 5, k -> 1},
{a -> 6, b -> 8, c -> 8, d -> 0, m -> 8, n -> 0, l -> 8, k -> 6}}


1000 a + 100 b + 10 c + d == (10 m + n) (10 l + k) /. res

{True, True, True, True}


Hope this helps.

• This suffers from the same problem as I discussed in my comment on @Artes answer -- 9801 is not a vampire number, and the product being used is 99*99. Commented Mar 11, 2017 at 13:33
• @WReach See update. Commented Mar 11, 2017 at 13:37
• That's better. But for reasons that are not immediately apparent to me, it is not finding 1395, 1827 or 2187. Commented Mar 11, 2017 at 13:42
• why there are only four solution Solve[{1000 a + 100 b + 10 c + d == (10 m + n) (10 l + k), 0 < a <= 9, 0 <= {b, c, d, m, n, l, k} <= 9, {a, b, c, d, m, n, l, k} \[Element] Integers, {{m}, {n}, {l}, {k}} \[Element] Interval @@ Map[{#, #} &, {a, b, c, d}], Sort@{a, b, c, d} == Sort@{m, n, l, k}}, {a, b, c, d, m, n, l, k}] ? $\{1260,1435,1530,6880\}$
– vito
Commented Mar 11, 2017 at 13:43
• @WReach Yes. That is odd. I placed 3 in position b in all locations and removed it from the variables and it still does not find it. Has a bug been found? Commented Mar 11, 2017 at 13:55

You could simply try all $x$ and $y$ pairs:

n = 2;

nDigitIntegers = Range[10^(n - 1), 10^n - 1];

vampireNumbersQ = Function[{x, y},
IntegerLength[x*y] == 2 n &&
Sort[IntegerDigits[x + y*10^n]] == Sort[IntegerDigits[x*y]]];

xy = Flatten[Table[
If[vampireNumbersQ[x, y],
{x, y} -> x*y, {}],
{x, nDigitIntegers}, {y, nDigitIntegers[[1]], x}
]]; // AbsoluteTiming


{0.0327479, Null}

xy


{{15, 93} -> 1395, {21, 60} -> 1260, {21, 87} -> 1827, {27, 81} -> 2187, {30, 51} -> 1530, {35, 41} -> 1435, {80, 86} -> 6880}

at least for small numbers, this is very fast. I can find all 2-digit vampire numbers in 0.03 seconds. For n=3 it takes about 3s to find:

{387*351==135837, 401*260==104260, 401*350==140350, 410*323==132430, 410*350==143500, 414*351==145314, 414*396==163944, 422*311==131242, 423*315==133245, 425*317==134725, 431*356==153436, 443*281==124483, 461*317==146137, 470*371==174370, 486*261==126846, 491*395==193945, 500*251==125500, 501*210==105210, 501*300==150300, 510*201==102510, 510*246==125460, 510*300==153000, 516*204==105264, 524*461==241564, 534*231==123354, 542*470==254740, 543*231==125433, 581*269==156289, 581*422==245182, 582*489==284598, 585*261==152685, 588*231==135828, 591*327==193257, 591*534==315594, 600*201==120600, 600*210==126000, 602*437==263074, 608*251==152608, 612*468==286416, 615*204==125460, 623*524==326452, 626*341==213466, 627*201==126027, 630*585==368550, 635*215==136525, 635*530==336550, 641*533==341653, 650*281==182650, 650*641==416650, 651*240==156240, 662*593==392566, 671*323==216733, 678*321==217638, 678*420==284760, 686*533==365638, 692*338==233896, 701*158==110758, 701*167==117067, 704*650==457600, 705*150==105750, 707*431==304717, 719*275==197725, 725*161==116725, 725*179==129775, 725*350==253750, 725*431==312475, 750*231==173250, 759*231==175329, 759*681==516879, 761*152==115672, 780*624==486720, 782*221==172822, 800*473==378400, 801*135==108135, 801*225==180225, 801*252==201852, 806*323==260338, 807*255==205785, 810*225==182250, 810*270==218700, 822*276==226872, 824*152==125248, 831*465==386415, 834*570==475380, 840*141==118440, 840*546==458640, 842*269==226498, 845*491==414895, 845*590==498550, 846*540==456840, 851*296==251896, 855*630==538650, 860*251==215860, 860*800==688000, 863*392==338296, 864*216==186624, 864*657==567648, 870*210==182700, 870*435==378450, 875*650==568750, 876*843==738468, 878*431==378418, 881*248==218488, 890*482==428980, 891*432==384912, 891*468==416988, 891*549==489159, 897*201==180297, 899*545==489955, 900*351==315900, 902*875==789250, 906*210==190260, 906*750==679500, 906*894==809964, 908*446==404968, 909*351==319059, 909*891==809919, 915*210==192150, 926*140==129640, 926*176==162976, 926*248==229648, 926*284==262984, 926*320==296320, 926*356==329656, 926*392==362992, 926*464==429664, 926*572==529672, 926*680==629680, 926*788==729688, 926*896==829696, 927*855==792585, 930*150==139500, 936*360==336960, 936*720==673920, 938*146==136948, 941*476==447916, 942*156==146952, 945*891==841995, 948*366==346968, 951*165==156915, 951*336==319536, 951*588==559188, 951*858==815958, 953*563==536539, 957*825==789525, 963*342==329346, 963*765==736695, 969*381==369189, 971*383==371893, 975*321==312975, 981*216==211896, 981*369==361989, 983*650==638950, 984*807==794088, 986*953==939658, 992*776==769792}

I haven't tried n=4, but it should take about 300s, as it has to check 100 times as many combinations.

Been a while since I participated on this site; don't get to use Mathematica that often either.

Figured I'd post this (not particularly fast) solution anyway:

vQ[n_ /; EvenQ@IntegerLength@n] :=
(Times @@FromDigits /@ {#[[;; (IntegerLength@n)/2]], #[[(IntegerLength@n)/2 + 1 ;;]]} &
/@ Permutations[IntegerDigits[n]])~MemberQ~n

vQ[_]:=False