Differently Perfect

For any positive integer n let f(n) denote the number of sets of 
four distinct positive integers a,b,c,d less than n such that 
ab-cd, ac-bd, and ad-bc are each multiples of n.  In terms of 
congruences, f(n) is the number of sets of four distinct non-zero 
residues a,b,c,d such that

         ab=cd      ac=bd      ad=bc     (mod n)

Obviously f(n) is a measure of the divisibility of n, somewhat
analogous to the sum proper divisors ("aliquot parts") of n, 
denoted by A(n).  Recall that the ancient Greeks classified all 
the integers as either "deficient", "perfect", or "abundant" 
accordingly as A(n) is less than, equal to, or greater than n.  
Similarly we could define a different kind of perfection (and 
deficiency and abundance) according to whether f(n) is less than, 
equal to, or greater than n.  According to this definition, the
only "perfect" number I can find is 42, although proving that
this is the only possible perfect number in this sense seems
difficult.

It would be helpful to have a simple expression for f(n) in
terms of the factorization of n, but the problem is complicated
by the fact that, unlike A(n)+n, the function f(n) is not 
multiplicative.  (I've tried modifying the definition of f
by allowing a,b,c,d to be 0 and/or non-distinct, but nothing
seems to result in a multiplicative function.)

The three basic congruences imply that, for any prime divisor p 
of n, if any one of the numbers a,b,c,d is divisible by p then at 
least three of them must be divisible by p.  Also, multiplying 
the three congruences in various combinations gives

             (bc)a^2 = (bc)d^2

             (bd)a^2 = (bd)c^2    (mod n)

             (cd)a^2 = (cd)b^2

so in the case where a,b,c,d are all coprime to p, it follows that

           a^2 = b^2 = c^2 = d^2     (mod p)

Therefore, modulo any given prime divisor p of n, the set {a,b,c,d}
must be of one of the three forms

                   {0,0,0,0}
                       or
                   {k,0,0,0}                     (mod p)
                       or
        {sqrt(k),sqrt(k),sqrt(k),sqrt(k)}

Now, since we require a,b,c,d to be positive and distinct, it's
clear that there are no acceptable quadruples if n is a prime or 
twice a prime (noting that a given residue k has at most two 
distinct square roots modulo a prime or twice a prime).  Thus 
we have f(p) = f(2p) = 0 for every prime p.  We can also see by 
inspection that f(9)=f(18)=0, but every other composite number n 
evidently has at least one solution.  Here's a short list of the 
non-zero values of f(n) for all n less than 100:

    n  f(n)       n  f(n)        n  f(n)        n  f(n)
   --- ----      --- ----       --- ----       --- ----
     8    1       36   20        57    9        80  291
    12    2       39    6        60  238        81   97
    15    2       40   65        63   36        84  356
    16    3       42   42        64  133        85   16
    20    4       44   10        65   12        87   14
    21    3       45   24        66   70        88  161
    24   33       48  147        68   16        90  288
    25    1       49   15        69   11        91   18
    27    9       50   26        70   84        92   22
    28    6       51    8        72  311        93   15
    30   28       52   12        75  161        95   18
    32   31       54  102        76   18        96 1071
    33    5       55   10        77   15        98  190
    35    6       56   97        78   84        99   60

Taking the lazy approach of just looking for patterns, these numbers
(and a few beyond this range) suggest the following rules for odd
primes p, q, and r:

      f(p)    =       0
      f(2p)   =       0
      f(4p)   =     1(p-1)  +  0
      f(8p)   =    16(p-1) +  1
      f(16p)  =    72(p-1) +  3
      f(32p)  =   520(p-1) + 31
      f(64p)  =  2128(p-1) + 133  =  133 f(8p)

      f(pq)  =   (1/4)(p-1)(q-1)
      f(2pq) =   (7/2)(p-1)(q-1)
      f(4pq) = (115/4)(p-1)(q-1) + f(4p) + f(4q) + f(pq)
             =   29(p-1)(q-1) + (p-1) + (q-1)

      f(9p) = 6(p-1)   (except for p=3)

      f(pqr) = (31/4)(p-1)(q-1)(r-1) + f(pq) + f(pr) + f(qr)

Checking each of these forms for "perfection", we find that the
only solution is for the case 2pq = f(2pq), which factors as

                  (3p-7)(3q-7) = 28

with the solution p=3, q=7, giving the perfect number 42.  It
appears that the formulas for the remaining forms (e.g., f(27p)
and f(pqrs) and so on) cannot give any more "perfect" numbers,
but it would be nice to have a simple formula giving f(n) for
ALL positive integers n.  Does anyone know of such a formula?
It seems to be some combination of the totients of the divisors
of n, but I can't quite see it.