Finding vampire numbers

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

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

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