%I #25 Jul 10 2024 15:02:55
%S 2,0,0,10,164,2670,56724,1867860,104538928,10461483366,1912179618740,
%T 644464839239880,402785011941549964,468944407349226545614,
%U 1021179521951204217530900,4174755063830188009750183026
%N Number of connected unlabeled symmetric relations (graphs with loops) having n nodes such that complement is also connected.
%H Andrew Howroyd, <a href="/A054922/b054922.txt">Table of n, a(n) for n = 1..50</a>
%H V. A. Liskovets, <a href="http://www.cs.uwaterloo.ca/journals/JIS/VOL3/LISK/Derseq.html">Some easily derivable sequences</a>, J. Integer Sequences, 3 (2000), #00.2.2.
%F a(n) = 2*A054921(n) - A000666(n).
%t A000666 = Cases[Import["https://oeis.org/A000666/b000666.txt", "Table"], {_, _}][[All, 2]];
%t A054921 = Cases[Import["https://oeis.org/A054921/b054921.txt", "Table"], {_, _}][[All, 2]];
%t a[n_] := 2*A054921[[n + 1]] - A000666[[n + 1]];
%t Array[a, 50] (* _Jean-François Alcover_, Aug 31 2019 *)
%o (Python)
%o from functools import lru_cache
%o from itertools import combinations
%o from math import prod, factorial, gcd
%o from fractions import Fraction
%o from sympy.utilities.iterables import partitions
%o from sympy import mobius, divisors
%o def A054922(n):
%o @lru_cache(maxsize=None)
%o def b(n): return int(sum(Fraction(1<<sum(p[r]*p[s]*gcd(r,s) for r,s in combinations(p.keys(),2))+sum(((q>>1)+1)*r+(q*r*(r-1)>>1) for q, r in p.items()),prod(q**r*factorial(r) for q, r in p.items())) for p in partitions(n)))
%o @lru_cache(maxsize=None)
%o def c(n): return n*b(n)-sum(c(k)*b(n-k) for k in range(1,n))
%o return (sum(mobius(d)*c(n//d) for d in divisors(n,generator=True))//n<<1)-b(n) # _Chai Wah Wu_, Jul 10 2024
%Y Cf. A000666, A054291.
%K nonn,easy
%O 1,1
%A _N. J. A. Sloane_, May 24 2000
%E More terms from _Vladeta Jovovic_, Jul 17 2000