%I #37 Mar 10 2021 17:26:50
%S 1,2,6,22,342,1444,33184,399235,12502550,117906198,7740054144,
%T 74673569118,4724493959332,121637216075836,4503600768557056,
%U 89450720590507768,10119960926575526448,152232968281237988010,16384000020089600480000,552020693349464673399080,35271474934322858202723576
%N Number of fundamentally different graceful labelings of digraphs with n arcs.
%C Consider vertices numbered 0 to n. For 1 <= k <= n, add an arc from v_k to (v_k+k) mod q, where q = n+1. (Thus (n+1)^n possibilities.) Two such digraphs are considered equivalent under the following operations: (i) rename each v to (v+1) mod q; (ii) rename each v to (av) mod q, where a is relatively prime to q; (iii) form the converse digraph by reversing the direction of every arc. The number of equivalence classes is a(n).
%D G. S. Bloom and D. F. Hsu, On graceful digraphs and a problem in network addressing, Congr. Numer., 35 (1982) 91-103.
%D D. E. Knuth, The Art of Computer Programming, forthcoming exercise in Section 7.2.2.3
%H Peter Luschny, <a href="/A341884/b341884.txt">Table of n, a(n) for n = 1..100</a>
%e Each equivalence class has exactly one digraph with v_1=0.
%e For n=3 the six classes of digraphs 0v_2v_3 are: {000,032}, {001,011,021,031}, {002,010,022,030}, {003,033}, {012,020}, {013,023}.
%p # This is a port from the Mathematica program below.
%p sols := proc(alf, bet, q)
%p igcd(alf, q); `if`(modp(bet, %) <> 0, 0, %) end:
%p F := proc(l, a, b, q) local s, r, ll;
%p s := 1; r := 1; ll := modp(l*a, q);
%p while ll > l do
%p s := s + 1;
%p ll := modp(ll*a, q);
%p r := modp(r*a+1, q); od;
%p `if`(ll <> l, 1, sols(a^s - 1, -r*b, q)) end:
%p G := proc(l, a, b, q) local s, r, ll;
%p s := 1; r := 1; ll := modp(-l*a, q);
%p while ll > l do
%p s := s + 1;
%p ll := modp(-ll*a, q);
%p r := modp(r*a+1, q); od;
%p `if`(ll <> l, 1, sols(a^s - 1, -r*b - l*irem(s, 2), q)) end:
%p f := (a, b, q) -> mul(F(l, a, b, q), l = 1..q-1):
%p g := (a, b, q) -> mul(G(l, a, b, q), l = 1..q-1):
%p a := proc(n) local q; q := n + 1;
%p add(`if`(igcd(a, q) > 1, 0, add(f(a, b, q) + g(a, b, q),
%p b = 0..q-1)), a = 1..q-1);
%p % / (2*q*numtheory:-phi(q)) end:
%p seq(a(n), n=1..21); # _Peter Luschny_, Mar 10 2021
%t sols[alf_,bet_,q_]:=Block[{d=GCD[alf,q]},If[Mod[bet,d]!=0,0,d]]
%t (* that many solutions to alf x == bet (modulo q) for 0<=x<q *)
%t f[l_,a_,b_,q_]:=Block[{r,s,ll,atos},
%t s=1;ll=Mod[l*a,q];r=1;
%t While[ll>l, s++;ll=Mod[ll*a,q];r=Mod[r*a+1,q]];
%t If[ll==l,sols[a^s-1,-r b,q],1]]
%t g[l_,a_,b_,q_]:=Block[{r,s,ll,atos},
%t s=1;ll=Mod[-l*a,q];r=1;
%t While[ll>l, s++;ll=Mod[-ll*a,q];r=Mod[r*a+1,q]];
%t If[ll==l,sols[a^s-1,-r b-l Mod[s,2],q],1]]
%t f[a_,b_,q_]:=Product[f[l,a,b,q],{l,q-1}]
%t g[a_,b_,q_]:=Product[g[l,a,b,q],{l,q-1}]
%t x[q_]:=Sum[If[GCD[a,q]>1,0,Sum[f[a,b,q]+g[a,b,q],{b,0,q-1}]],{a,q-1}]/
%t (2q EulerPhi[q])
%t a[n_]:=x[n+1]
%Y Cf. A342357, A071118.
%K nonn
%O 1,2
%A _Don Knuth_, Feb 22 2021
%E a(10) from _Don Knuth_, Mar 10 2021
%E More terms from _Peter Luschny_, Mar 10 2021