Write
((x)) = x − ⌊x⌋ − 1/2 ( = 0 if x ∈ ℤ )
and
s(h , k) = Σ_{r=1}^{k-1} ((r/k)) ((hr/k)) ( Dedekind sum ).

1. A classical congruence (Meyer 1957)
   For coprime positive integers h , k

     12 k · s(h , k) ≡ h + h̄    (mod k),          (1)
   

where h̄ is the inverse of h modulo k.
See J. Meyer, J. Reine Angew. Math. 203 (1960), 143–157;
H. Rademacher – E. Grosswald, “Dedekind sums”, Th. 2.

Equation (1) is the only ingredient that is needed.

2. The case n = 3 of your question

Put a, b, c positive with gcd(a,b,c)=1, write a′=gcd(b,c), b′=gcd(a,c), c′=gcd(a,b) and d=a′b′c′. For simplicity first assume a, b, c are pairwise coprime, i.e. a′=b′=c′=1, d=1; then

      S(a,b,c)=s(bc , a)+s(ac , b)+s(ab , c).

Fix the modulus a.
The second and the third summand contain the factor a and therefore are 0 modulo a. By (1) with (h,k)=(bc , a)

      12 a s(bc , a) ≡ bc + (bc)^{-1}  (mod a).

Multiplying by bc gives

      12 a b c s(bc , a) ≡ (bc)^2 + 1  (mod a).

Hence

      12 a b c S(a,b,c) ≡ (bc)^2 + 1      (mod a).

The same argument with the moduli b and c shows

      12 a b c S(a,b,c) ≡ (ab)^2 + 1  (mod b),
      12 a b c S(a,b,c) ≡ (ca)^2 + 1  (mod c).

Combining the three congruences by the Chinese remainder theorem yields

      12 a b c S(a,b,c) ≡ (ab)^2 + (bc)^2 + (ca)^2 + 1   (mod abc),

exactly your formula in the coprime situation.
Re-introducing the factors a′, b′, c′ and d is routine (one applies (1) to the reduced pairs ( bc/d , a/(b′c′) ) etc.), and produces your more general congruence

      12 abc 𝔖(a,b,c) ≡ (ab)^2 + (bc)^2 + (ca)^2 + d^2   (mod abc).

3. Arbitrary n (pairwise coprime)

Let n ≥ 2, pairwise coprime a₁,…,aₙ and
P = a₁ … aₙ,  S = Σ_{i=1}^{n} s(P/a_i , a_i).

Fix j and set x = P/a_j , y = a_j. By (1)

      12 y s(x , y) ≡ x + x̄   (mod y).

Multiplying by P/y = x we get

      12 P s(x , y) ≡ x² + 1   (mod y).

Modulo a_j every other term of 12 P S contains the factor a_j and vanishes, so

      12 P S ≡ (P/a_j)² + 1   (mod a_j)  (j = 1,…,n).

The n congruences again combine (Chinese remainder) to

      12 P S ≡ Σ_{i=1}^{n} (P/a_i)² + 1   (mod P),

i.e. exactly the formula you conjectured for n = 4, and in fact for every n.

4. Literature

• J. Meyer, “Über einige Anwendungen Dedekindscher Summen”,
J. Reine Angew. Math. 203 (1960), 143–157.
• H. Rademacher & E. Grosswald, “Dedekind sums”, Carus Monograph 16 (MAA, 1972), Theorem 2.
• Further appearances of (1) and its easy corollaries can be found in
R. Dieter, J. Number Theory 3 (1971), 433–454;
D. Hickerson, Pacific J. Math. 90 (1980), 207–208; etc.

Hence both your three–variable congruence and its extension to an arbitrary number of pairwise coprime integers are straightforward corollaries of the classical congruence (1). They have therefore been implicitly known since Meyer’s paper, although they seldom appear in print in exactly the form you give.