= register int a, aa, k; /* primary indices in loops */ int b, bb, c, cc, d, dd; /* secondary indices */ int n; /* the number of vertices */ int n_factor; /* either ${1\over2}(q-1)$ (type~3) or $q-1$ (type 4) */ register Vertex *v; /* the current vertex of interest */ @ Here we implicitly test that |q| is prime, by finding a primitive root whose powers generate everything. If |q| is not prime, its smallest divisor will cause the inner loop in this step to terminate with |k>=q|, because no power of that divisor will be congruent to~1. @= for (a=2; ; a++) if (q_inv[a]==0) { for (b=a,k=1; b!=1&&k =q) dead_panic(bad_specs+1); /* |q| is not prime */ if (k==q-1) break; /* good, |a| is the primitive root we seek */ } @ As soon as we have discovered a primitive root, it is easy to generate all the inverses. (We could also generate the discrete logarithms if we had a need for them.) We set |q_inv[0]=q|; this will be our internal representation of $\infty$. @= for (b=a,bb=aa; b!=bb; b=(a*b)%q,bb=(aa*bb)%q) q_inv[b]=bb,q_inv[bb]=b; q_inv[1]=1; q_inv[b]=b; /* at this point |b| must equal |q-1| */ q_inv[0]=q; @ The conditions we stated for validity of |q| when |p=2| are equivalent to the existence of $\sqrt{-2}$ and $\sqrt{13}$ modulo~|q|, according to the law of quadratic reciprocity (see, for example, {\sl Fundamental Algorithms}, exercise 1.2.4--47). @ = if (p==2) { if (q_sqrt[13%q]<0 || q_sqrt[q-2]<0) dead_panic(bad_specs+2); /* improper prime to go with |p=2| */ } if ((a=p%q)==0) dead_panic(bad_specs+3); /* |p| divisible by |q| */ if (type==0) type=(q_sqrt[a]>0? 3: 4); n_factor=(type==3? (q-1)/2: q-1); switch (type) { case 1: n=q+1;@+break; case 2: n=q*(q+1)/2;@+break; default: if ((q_sqrt[a]>0 && type!=3) || (q_sqrt[a]<0 && type!=4)) dead_panic(bad_specs+4); /* wrong type for |p| modulo |q| */ if (q>1289) dead_panic(bad_specs+5); /* way too big for types 3, 4 */ n=n_factor*q*(q+1); break; } if (p>=(long)(0x3fffffff/n)) dead_panic(bad_specs+6); /* $(p+1)n\ge2^{30}$ */ @* The vertices. Graphs of type 1 will have vertices from the set $\{0,1,\ldots,q-1,\infty\}$, namely the integers modulo~|q| with an additional ``infinite'' element thrown in. The idea will be to operate on these quantities by adding constants, and/or multiplying by constants, and/or taking reciprocals, modulo~|q|. Graphs of type 2 will have vertices that are unordered pairs of distinct elements from that same set. Graphs of types 3 and 4 will have vertices that are $2\times2$ matrices having nonzero determinants modulo~|q|. The determinants of type~3 matrices will, in fact, be nonzero quadratic residues. We consider two matrices to be equivalent if one is obtained from the other by multiplying all entries by a constant (modulo~|q|); therefore we will normalize all the matrices so that the second row is either $(0,1)$ or has the form $(1,x)$ for some~$x$. The total number of equivalence classes of type~4 matrices obtainable in this way is $(q+1)q(q-1)$, because we can choose the second row in $q+1$ ways, after which there are two cases: Either the second row is $(0,1)$, and we can select the upper right corner element arbitrarily and choose the upper left corner element nonzero; or the second row is $(1,x)$, and we can select the upper left corner element arbitrarily and then choose an upper right corner element to make the determinant nonzero. For type~3 the counting is similar, except that ``nonzero'' becomes ``nonzero quadratic residue,'' hence there are exactly half as many choices. It is easy to verify that the equivalence classes of matrices that correspond to vertices in these graphs of types 3 and~4 are closed under matrix multiplication. Therefore the vertices may be regarded as the elements of finite groups. The type~3 group for a given |q| is often called the linear fractional group $LF(2,{\bf F}_q)$, or the projective special linear group $PSL(2,{\bf F}_q)$, or the linear simple group $L_2(q)$; it can also be regarded as the group of $2\times2$ matrices with determinant~1 (mod~$q$), when the matrix $A$ is considered equivalent to $-A$. (This group is a simple group for all primes |q>2|.) The type~4 group is officially known as the projective general linear group of degree~2 over the field of |q|~elements, $PGL(2,{\bf F}_q)$. @ = new_graph=gb_new_graph(n); if (new_graph==NULL) dead_panic(no_room); /* out of memory before we try to add edges */ sprintf(new_graph->id,"raman(%d,%d,%u,%u)",p,q,type,reduce); strcpy(new_graph->format,"ZZZIIZIZZZZZZZ"); v=new_graph->vertices; switch(type) { case 1: @ ;@+break; case 2: @ ;@+break; default: @ ;@+break; } @ Type 1 graphs are the easiest to label. We store a serial number in utility field |x.i|, using $q$ to represent $\infty$. @ = new_graph->format[4]='Z'; for (a=0;a name=gb_save_string(name_buf); v->x.i=a; v++; } v->name=gb_save_string("INF"); v->x.i=q; v++; @ @= static char name_buf[]="(1111,1111;1,1111)"; /* place to form vertex names */ @ The type 2 labels run from $\{0,1\}$ to $\{q-1,\infty\}$; we put the coefficients into |x.i| and |y.i|, where they might prove useful in some applications. @ = for (a=0;a name=gb_save_string(name_buf); v->x.i=a;@+v->y.i=aa; v++; } @ For graphs of types 3 and 4, we set the |x.i| and |y.i| fields to the elements of the first row of the matrix, and we set the |z.i| field equal to the ratio of the elements of the second row (again with $q$ representing~$\infty$). The vertices in this case will consist of |q(q+1)| blocks of vertices having a given second row and a given element in the upper left or upper right position. Within each block of vertices, the determinants will be respectively congruent modulo~|q| to $1^2$, $2^2$, \dots,~$({q-1\over2})^2$ in the case of type~3 graphs, or to 1,~2, \dots,~$q-1$ in the case of type~4. @= new_graph->format[5]='I'; for (c=0;c<=q;c++) for (b=0;b z.i=c; if (c==q) { /* second row of matrix is $(0,1)$ */ v->y.i=b; v->x.i=(type==3? q_sqr[a]: a); /* determinant is $a^2$ or $a$ */ sprintf(name_buf,"(%d,%d;0,1)",v->x.i,b); } else { /* second row of matrix is $(1,c)$ */ v->x.i=b; v->y.i=(b*c+q-(type==3? q_sqr[a]: a))%q; sprintf(name_buf,"(%d,%d;1,%d)",b,v->y.i,c); } /* determinant is $a^2$ or $a$ */ v->name=gb_save_string(name_buf); v++; } @* Group generators. We will define a set of |p+1| permutations $\{\pi_0, \pi_1,\ldots,\pi_p\}$ of the vertices, such that the arcs of our graph will go from $v$ to $v\pi_k$ for |0<=k<=p|. Thus, each path in the graph will be defined by a product of permutations; the cycles of the graph will correspond to vertices that are left fixed by a product of permutations. The graph will be undirected, because the inverse of each $\pi_k$ will also be one of the permutations of the generating set. In fact, each permutation $\pi_k$ will be defined by a $2\times2$ matrix; for graphs of types 3 and~4, the permutations will therefore correspond to certain vertices, and the vertex $v\pi_k$ will simply be the product of matrix $v$ by matrix $\pi_k$. For graphs of type 1, the permutations will be defined by linear fractional transformations, which are mappings of the form $$v\;\longmapsto\; {av+b\over cv+d}\bmod q\,.$$ This transformation applies to all $v\in\{0,1,\ldots,q-1,\infty\}$, under the usual conventions that $x/0=\infty$ when $x\ne0$ and $(x\infty+x')/(y\infty+y')=x/y$. The composition of two such transformations is again a linear fractional transformation, corresponding to the product of the two associated matrices $\bigl({a\,b\atop c\,d}\bigr)$. Graphs of type 2 will be handled just like graphs of type 1, except that we will compute the images of two distinct points $v=\{v_1,v_2\}$ under the linear fractional transformation. The two images will be distinct, because the transformation is invertible. When |p=2|, a special set of three generating matrices $\pi_0$, $\pi_1$, $\pi_2$ can be shown to define Ramanujan graphs; these matrices are described below. Otherwise |p| is odd, and the generators are based on the theory of integral quaternions. Integral quaternions are quadruples of the form $\alpha=a_0+a_1i+a_2j+a_3k$, where $a_0$, $a_1$, $a_2$, and~$a_3$ are integers; we multiply them by using the associative but noncommutative multiplication rules $i^2=j^2=k^2=ijk=-1$. If we write $\alpha=a+A$, where $a$ is the ``scalar'' $a_0$ and $A$ is the ``vector'' $a_1i+a_2j+a_3k$, the product of quaternions $\alpha=a+A$ and $\beta=b+B$ can be expressed as $$(a+A)(b+B)=ab-A\cdot B+aB+bA+A\times B\,,$$ where $A\cdot B$ and $A\times B$ are the usual dot product and cross product of vectors. The conjugate of $\alpha=a+A$ is $\overline\alpha=a-A$, and we have $\alpha\overline\alpha=a_0^2+a_1^2+a_2^2+a_3^2$. This important quantity is called $N(\alpha)$, the norm of $\alpha$. It is not difficult to verify that $N(\alpha\beta)=N(\alpha)N(\beta)$, because we have $\overline{\mathstrut\alpha\beta}=\overline{\mathstrut\beta} \,\overline{\mathstrut\alpha}$ and $\alpha x=x\alpha$ when $x$ is scalar. Integral quaternions have a beautiful theory; for example, there is a nice variant of Euclid's algorithm by which we can compute the greatest common left divisor of any two integral quaternions, and this makes it possible to prove that integral quaternions whose coefficients are relatively prime can be uniquely factored into quaternions whose norm is prime. However, the details of that theory are beyond the scope of this documentation. It will suffice for our purposes to observe that we can use quaternions to define the finite groups $PSL(2,{\bf F}_q)$ and $PGL(2,{\bf F}_q)$ in a different way from the definitions given earlier: Suppose we consider two quaternions to be equivalent if their coefficients are equal modulo~|q|, or if one is a nonzero scalar multiple of the other (modulo~|q|). Thus, for example, if $q=3$ we consider $1+4i-j$ to be equivalent to $1+i+2j$, and also equivalent to $2+2i+j$. It turns out that there are exactly $(q+1)q(q-1)$ such equivalence classes, and they form a group under quaternion multiplication that is the same as the projective group of $2\times2$ matrices under matrix multiplication, modulo~|q|. One way to prove this is by means of the one-to-one correspondence $$a_0+a_1i+a_2j+a_3k\;\longleftrightarrow\; \left(\matrix{a_0+a_1g+a_3h&a_2+a_3g-a_1h\cr -a_2+a_3g-a_1h&a_0-a_1g-a_3h\cr}\right)\,,$$ where $g$ and $h$ are integers with $g^2+h^2\=-1$ (mod~|q|). Jacobi proved that the number of ways to represent any odd number |p| as a sum of four squares $a_0^2+a_1^2+a_2^2+a_3^2$ is 8 times the sum of divisors of~|p|. [This fact appears in the concluding sentence of his monumental work {\sl Fundamenta Nova Theori\ae\ Functionum Ellipticorum}, K\"onigsberg, 1829.] In particular, when |p| is prime, the number of such representations is $8(p+1)$; in other words, there are exactly $8(p+1)$ quaternions $\alpha=a_0+a_1i+a_2j+a_3k$ with $N(\alpha)=p$. These quaternions form |p+1| equivalence classes under multiplication by the eight ``unit quaternions'' $\{\pm1,\pm i,\pm j,\pm k\}$; we will select one element from each equivalence class, and the resulting |p+1| quaternions will correspond to |p+1| matrices, which will generate the |p+1| arcs leading from each vertex in the graphs to be constructed. @= typedef struct { long a0,a1,a2,a3; /* coefficients of a quaternion */ unsigned bar; /* the index of the inverse (conjugate) quaternion */ } quaternion; @ A global variable |gen_count| will be declared below, indicating the number of generators found so far. When |p| isn't prime, we will find more than |p+1| solutions; we allocate one extra slot in the |gen| table to hold a possible overflow entry. @ = gen=gb_alloc_type(p+2,@[quaternion@],working_storage); if (gen==NULL) late_panic(no_room+2); /* not enough memory */ gen_count=0;@+max_gen_count=p+1; if (p==2) @ @; else @ ; if (gen_count!=max_gen_count) late_panic(bad_specs+7); /* |p| is not prime */ @ @ = static quaternion *gen; /* table of the |p+1| generators */ @ As mentioned above, quaternions of norm |p| come in sets of 8, differing from each other only by unit multiples; we need to choose one of the~8. Suppose $a_0^2+a_1^2+a_2^2+a_3^2=p$. If $p\bmod4=1$, exactly one of the $a$'s will be odd; so we call it $a_0$ and assign it a positive sign. When $p\bmod4=3$, exactly one of the $a$'s will be even; we call it $a_0$, and if it is nonzero we make it positive. If $a_0=0$, we make sure that one of the others---say the rightmost appearance of the largest one---is positive. In this way we obtain a unique representative from each set of 8 equivalent quaternions. For example, the four quaternions of norm 3 are $\pm i\pm j+k$; the six of norm~5 are $1\pm2i$, $1\pm2j$, $1\pm2k$. In the program here we generate solutions to $a^2+b^2+c^2+d^2=p$ when $a\not\=b\=c\=d$ (mod~2) and $b\le c\le d$. The variables |aa|, |bb|, and |cc| hold the respective values $p-a^2-b^2-c^2-d^2$, $p-a^2-3b^2$, and $p-a^2-2c^2$. The |for| statements use the fact that $a^2$ increases by $4(a+1)$ when $a$ increases by~2. @ = {@+long sa,sb,sc; /* $p-a^2$, $p-a^2-b^2$, $p-a^2-b^2-c^2$ */ int pp=(p>>1)&1; /* 0 if $p\bmod4=1$, \ 1 if $p\bmod4=3$ */ for (a=1-pp,sa=p-a;sa>0;sa-=(a+1)<<2,a+=2) for (b=pp,sb=sa-b,bb=sb-b-b;bb>=0;bb-=12*(b+1),sb-=(b+1)<<2,b+=2) for (c=b,cc=bb,sc=(sb+cc)>>1;cc>=0;cc-=(c+1)<<3,sc-=(c+1)<<2,c+=2) for (d=c,aa=cc;aa>=0;aa-=(d+1)<<2,d+=2) if (aa==0) @ ; @ ; } @ If |a>0| and |0 = static unsigned gen_count; /* the next available quaternion slot */ static unsigned max_gen_count; /* $p+1$, stored as a global variable */ static void deposit(a,b,c,d) long a,b,c,d; /* a solution to $a^2+b^2+c^2+d^2=p$ */ { if (gen_count>=max_gen_count) /* oops, we already found |p+1| solutions */ gen_count=max_gen_count+1; /* this will happen only if |p| isn't prime */ else { gen[gen_count].a0=gen[gen_count+1].a0=a; gen[gen_count].a1=b;@+gen[gen_count+1].a1=-b; gen[gen_count].a2=c;@+gen[gen_count+1].a2=-c; gen[gen_count].a3=d;@+gen[gen_count+1].a3=-d; if (a) { gen[gen_count].bar=gen_count+1; gen[gen_count+1].bar=gen_count; gen_count+=2; } else { gen[gen_count].bar=gen_count; gen_count++; } } } @ @= { deposit(a,b,c,d); if (b) { deposit(a,-b,c,d);@+deposit(a,-b,-c,d); } if (c) deposit(a,b,-c,d); if (b = {@+register int g,h; int a00,a01,a10,a11; /* entries of $2\times2$ matrix */ for (k=q-1;q_sqrt[k]<0;k--) ; /* find the largest quadratic residue, |k| */ g=q_sqrt[k];@+h=q_sqrt[q-1-k]; for (k=p;k>=0;k--) { a00=(gen[k].a0+g*gen[k].a1+h*gen[k].a3)%q; if (a00<0) a00+=q; a11=(gen[k].a0-g*gen[k].a1-h*gen[k].a3)%q; if (a11<0) a11+=q; a01=(gen[k].a2+g*gen[k].a3-h*gen[k].a1)%q; if (a01<0) a01+=q; a10=(-gen[k].a2+g*gen[k].a3-h*gen[k].a1)%q; if (a10<0) a10+=q; gen[k].a0=a00;@+gen[k].a1=a01;@+gen[k].a2=a10;@+gen[k].a3=a11; } } @ When |p=2|, the following three appropriate generating matrices have been found by P.~Chiu: $$\left(\matrix{1&0\cr 0&-1\cr}\right)\,,\qquad \left(\matrix{2+s&t\cr t&2-s\cr}\right)\,,\qquad\hbox{and}\qquad \left(\matrix{2-s&-t\cr-t&2+s\cr}\right)\,,$$ where $s^2\=-2$ and $t^2\=-26$ (mod~$q$). The determinants of these matrices are respectively $-1$, $32$, and~$32$; the product of the second and third matrices is 32 times the identity matrix. Notice that when 2 is a quadratic residue (this happens when $q=8k+1$), the determinants are all quadratic residues, so we get a graph of type~3; when 2 is a quadratic nonresidue (which happens when $q=8k+3$), the determinants are all nonresidues, so we get a graph of type~4. @ = {@+int s=q_sqrt[q-2], t=(q_sqrt[13%q]*s)%q; gen[0].a0=1;@+gen[0].a1=gen[0].a2=0;@+gen[0].a3=q-1;@+gen[0].bar=0; gen[1].a0=gen[2].a3=(2+s)%q; gen[1].a1=gen[1].a2=t; gen[2].a1=gen[2].a2=q-t; gen[1].a3=gen[2].a0=(q+2-s)%q; gen[1].bar=2;@+gen[2].bar=1; gen_count=3; } @* Constructing the edges. The remaining task is to use the permutations defined by the |gen| table to create the arcs of the graph and their inverses. The |ref| fields in each arc will refer to the permutation leading to the arc. In most cases each vertex |v| will have degree exactly |p+1|, and the edges emanating from it will appear in a linked list having the respective |ref| fields 0,~1, \dots,~|p| in order. However, if |reduce| is nonzero, self-loops and multiple edges will be eliminated, so the degree may be less than |p+1|; in this case the |ref| fields will still be in ascending order, but some generators won't be referenced. There is also a subtle case where |reduce=0| but the degree of a vertex might actually be greater than |p+1|. We want the graph |g| generated by |raman| to satisfy the conventions for undirected graphs stated in |gb_graph|; therefore, if any of the generating permutations has a fixed point, we will create two arcs for that fixed point, and the corresponding vertex |v| will have an edge running to itself. Since each edge consists of two arcs, such an edge will produce two consecutive entries in the list |v->arcs|. If the generating permutation happens to be its own inverse, there will be two consecutive entries with the same |ref| field; this means there will be more than |p+1| entries in |v->arcs|, and the total number of arcs |g->m| will exceed |(p+1)n|. Self-inverse generating permutations arise only when |p=2| or when $p$ is expressible as a sum of three odd squares (hence $p\bmod8=3$); and such permutations will have fixed points only when |type<3|. Therefore this anomaly does not arise often. But it does occur, for example, in the smallest graph generated by |raman|, namely when |p=2|, |q=3|, and |type=1|, when there are 4~vertices and 14 (not~12) arcs. @d ref a.i /* the |ref| field of an arc refers to its permutation number */ @ = for (k=p;k>=0;k--) {@+int kk; if ((kk=gen[k].bar)<=k) /* we assume that |kk=k| or |kk=k-1| */ for (v=new_graph->vertices;v vertices+n;v++) { register Vertex* u; @ ; if (u==v) { if (!reduce) { gb_new_edge(v,v,1); v->arcs->ref=kk;@+(v->arcs+1)->ref=k; /* see the remarks above regarding the case |kk=k| */ } } else {@+register Arc* ap; if (u->arcs && u->arcs->ref==kk) continue; /* |kk=k| and we've already done this two-cycle */ else if (reduce) for (ap=v->arcs;ap;ap=ap->next) if (ap->tip==u) goto done; /* there's already an edge between |u| and |v| */ gb_new_edge(v,u,1); v->arcs->ref=k;@+u->arcs->ref=kk; if ((ap=v->arcs->next)!=NULL && ap->ref==kk) { v->arcs->next=ap->next;@+ap->next=v->arcs;@+v->arcs=ap; } /* now the |v->arcs| list has |ref| fields in order again */ done:; } } } @ For graphs of types 3 and 4, our job is to compute a $2\times2$ matrix product, reduce it modulo~|q|, and find the appropriate equivalence class~|u|. @ = if (type<3) @ @; else {@+long a0=gen[k].a0,a1=gen[k].a1,a2=gen[k].a2,a3=gen[k].a3; a=v->x.i;@+b=v->y.i; if (v->z.i==q) c=0,d=1; else c=1,d=v->z.i; @ ; a=(cc? q_inv[cc]: q_inv[dd]); /* now |a| is a normalization factor */ d=(a*dd)%q;@+c=(a*cc)%q;@+b=(a*bb)%q;@+a=(a*aa)%q; @ ; } @ @ = aa=(a*a0+b*a2)%q; bb=(a*a1+b*a3)%q; cc=(c*a0+d*a2)%q; dd=(c*a1+d*a3)%q; @ @ = if (c==0) d=q,aa=a; else { aa=(a*d-b)%q; if (aa<0) aa+=q; b=a; } /* now |aa| is the determinant of the matrix */ u=new_graph->vertices+((d*q+b)*n_factor+(type==3? q_sqrt[aa]: aa)-1); @* Linear fractional transformations. Given a nonsingular $2\times2$ matrix $\bigl({a\,b\atop c\,d}\bigr)$, the linear fractional transformation $z\mapsto(az+b)/(cz+d)$ is defined modulo~$q$ by the following subroutine. We assume that the matrix $\bigl({a\,b\atop c\,d}\bigr)$ appears in row |k| of the |gen| table. @ = static long lin_frac(a,k) long a; /* the number being transformed; $q$ represents $\infty$ */ unsigned k; /* index into |gen| table */ {@+register long q=q_inv[0]; /* the modulus */ long a00=gen[k].a0, a01=gen[k].a1, a10=gen[k].a2, a11=gen[k].a3; /* the coefficients */ register num, den; /* numerator and denominator */ if (a==q) num=a00, den=a10; else num=(a00*a+a01)%q, den=(a10*a+a11)%q; if (den==0) return q; else return (num*q_inv[den])%q; } @ We are computing the same values of |lin_frac| over and over again in type~2 graphs, but the author was too lazy to optimize this. @ = if (type==1) u=new_graph->vertices+lin_frac(v->x.i,k); else { a=lin_frac(v->x.i,k);@+aa=lin_frac(v->y.i,k); u=new_graph->vertices+(a