summaryrefslogtreecommitdiff
path: root/support/graphbase/gb_raman.w
diff options
context:
space:
mode:
authorNorbert Preining <norbert@preining.info>2019-09-02 13:46:59 +0900
committerNorbert Preining <norbert@preining.info>2019-09-02 13:46:59 +0900
commite0c6872cf40896c7be36b11dcc744620f10adf1d (patch)
tree60335e10d2f4354b0674ec22d7b53f0f8abee672 /support/graphbase/gb_raman.w
Initial commit
Diffstat (limited to 'support/graphbase/gb_raman.w')
-rw-r--r--support/graphbase/gb_raman.w715
1 files changed, 715 insertions, 0 deletions
diff --git a/support/graphbase/gb_raman.w b/support/graphbase/gb_raman.w
new file mode 100644
index 0000000000..243796ef15
--- /dev/null
+++ b/support/graphbase/gb_raman.w
@@ -0,0 +1,715 @@
+% This file is part of the Stanford GraphBase (c) Stanford University 1992
+\def\title{GB\_\thinspace RAMAN}
+@i boilerplate.w %<< legal stuff: PLEASE READ IT BEFORE MAKING ANY CHANGES!
+\let\==\equiv % congruence sign
+
+\prerequisite{GB\_\thinspace GRAPH}
+@* Introduction. This GraphBase module contains the |raman| subroutine,
+which creates a family of ``Ramanajun graphs'' based on a theory
+developed by Alexander Lubotzky, Ralph Phillips, and Peter Sarnak
+[see {\sl Combinatorica \bf8} (1988), 261--277].
+
+Ramanujan graphs are defined by the following properties:
+They are connected, undirected graphs in which every vertex has
+degree~|k|, and every eigenvalue of the adjacency matrix
+is either $\pm k$ or has absolute value $\le2\sqrt{\mathstrut k-1}$.
+Such graphs are known to have good expansion properties, small diameter,
+and relatively small independent sets; they cannot be colored with
+fewer than $k/\bigl(2\sqrt{\mathstrut k-1}\,\bigr)$ colors unless they are
+bipartite. The particular examples of Ramanujan graphs constructed here
+are based on interesting properties of quaternions with integer coefficients.
+
+An example of the use of this procedure can be found in the demo program
+called |girth|.
+
+@(gb_raman.h@>=
+extern Graph *raman();
+
+@ The subroutine call `|raman(p,q,type,reduce)|'
+constructs an undirected graph in which each vertex has degree~|p+1|.
+The number of vertices is~|q+1| if |type=1|, or~${1\over2}q(q+1)$ if |type=2|,
+or ${1\over2}(q-1)q(q+1)$ if |type=3|, or |(q-1)q(q+1)| if
+|type=4|. The graph will be bipartite if and only if it has type~4.
+Parameters |p| and |q| must be distinct prime numbers,
+and |q|~must be odd. Furthermore there are additional restrictions:
+If |p=2|, the other parameter |q| must satisfy $q\bmod8\in\{1,3\}$
+and $q\bmod13\in{1,3,4,9,10,12}$; this rules out about one fourth of
+all primes. Moreover, if |type=3| the value of |p| must be a
+quadratic residue modulo~$q$; in other words, there must be an
+integer~$x$ such that $x^2\=p$ (mod~$q$). If |type=4|, the value of |p|
+must not be a quadratic residue.
+
+If you specify |type=0|, the procedure
+will choose the largest permissible type (either 3 or~4);
+the value of the type selected will
+appear as part of the string placed in the resulting graph's |id| field.
+For example, if |type=0|, |p=2|, and |q=43|, a type~4 graph will be
+generated, because 2 is not a quadratic residue modulo~43. This
+graph will have $44\times43\times42=79464$ vertices, each of degree~3.
+(Notice that graphs of types 3 and~4 can be quite large even when
+|q| is rather small.)
+
+The largest permissible value of |q| is 46337; this is the largest
+prime whose square is less than $2^{31}$. Of course you would use
+it only for a graph of type~1.
+
+If |reduce| is nonzero, loops and multiple edges will be suppressed.
+In this case the degrees of some vertices may turn out to be less than~|p+1|,
+in spite of what was said above.
+
+Although type 4 graphs are bipartite, the vertices
+are not separated into two blocks as in other bipartite
+graphs produced by GraphBase routines.
+
+All edges of the graphs have length 1.
+
+@ If the |raman| routine encounters a problem, it returns |NULL|
+(\.{NULL}), after putting a code number into the external variable
+|panic_code|. This code number identifies the type of failure.
+Otherwise |raman| returns a pointer to the newly created graph, which
+will be represented with the data structures explained in |gb_graph|.
+(The external variable |@!panic_code| is itself defined in
+|gb_graph|.)
+
+@d panic(c) @+{@+panic_code=c;@+gb_alloc_trouble=0;@+return NULL;@+}
+@d dead_panic(c) {@+gb_free(working_storage);@+panic(c);@+}
+@d late_panic(c) {@+gb_recycle(new_graph);@+dead_panic(c);@+}
+@#
+@f Graph int /* |gb_graph| defines the |Graph| type and a few others */
+@f Vertex int
+@f Arc int
+@f Area int
+
+@ The \Cee\ file \.{gb\_raman.c} has the following general shape:
+
+@p
+#include "gb_graph.h" /* we will use the |gb_graph| data structures */
+@#
+@<Type declarations@>@;
+@<Private variables and routines@>@;
+@#
+Graph *raman(p,q,type,reduce)
+ int p; /* one less than the desired degree; must be prime */
+ int q; /* size parameter; must be prime and properly related to |type| */
+ unsigned type; /* selector between different possible constructions */
+ unsigned reduce; /* if nonzero, multiple edges and self-loops won't occur */
+{@+@<Local variables@>@;
+ @<Prepare tables for doing arithmetic modulo~|q|@>;
+ @<Choose or verify the |type|, and determine the number |n| of vertices@>;
+ @<Set up a graph with |n| vertices, and assign vertex labels@>;
+ @<Compute |p+1| generators that will define the graph's edges@>;
+ @<Append the edges@>;
+ if (gb_alloc_trouble)
+ late_panic(alloc_fault);
+ /* oops, we ran out of memory somewhere back there */
+ gb_free(working_storage);
+ return new_graph;
+}
+
+@ @<Local var...@>=
+Graph *new_graph; /* the graph constructed by |raman| */
+Area working_storage; /* place for auxiliary tables */
+
+@* Brute force number theory. Instead of using routines like Euclid's
+algorithm to compute inverses and square roots modulo~|q|, we have
+plenty of time to build complete tables, since |q| is smaller than
+the number of vertices we will be generating.
+
+We will make three tables: |q_sqr[k]| will contain $k^2$ modulo~|q|;
+|q_sqrt[k]| will contain one of the values of $\sqrt{\mathstrut k}$
+if $k$ is a quadratic residue; and |q_inv[k]| will contain the multiplicative
+inverse of~|k|.
+
+@<Private...@>=
+static int *q_sqr; /* squares */
+static int *q_sqrt; /* square roots (or $-1$ if not a quadratic residue) */
+static int *q_inv; /* reciprocals */
+
+@ @<Prepare tables for doing arithmetic modulo~|q|@>=
+if (q<3 || q>46337) panic(very_bad_specs);
+ /* |q| is way too small or way too big */
+if (p<2) panic(very_bad_specs+1); /* |p| is way too small */
+init_area(working_storage);
+q_sqr=gb_alloc_type(3*q,@[int@],working_storage);
+if (q_sqr==0) panic(no_room+1);
+q_sqrt=q_sqr+q;
+q_inv=q_sqrt+q; /* note that |gb_alloc| has initialized everything to zero */
+@<Compute the |q_sqr| and |q_sqrt| tables@>;
+@<Find a primitive root |a|, modulo |q|, and its inverse |aa|@>;
+@<Compute the |q_inv| table@>;
+
+@ @<Compute the |q_sqr| and |q_sqrt| tables@>=
+for (a=1; a<q; a++) q_sqrt[a]=-1;
+for (a=1,aa=1; a<q; aa=(aa+a+a+1)%q,a++) {
+ q_sqr[a]=aa;
+ q_sqrt[aa]=q-a; /* the smaller square root will survive */
+ q_inv[aa]=-1;
+ /* we make |q_inv[aa]| nonzero when |aa| can't be a primitive root */
+}
+
+@ @<Local v...@>=
+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.
+
+@<Find a primitive root |a|, modulo |q|, and its inverse |aa|@>=
+for (a=2; ; a++)
+ if (q_inv[a]==0) {
+ for (b=a,k=1; b!=1&&k<q; aa=b,b=(a*b)%q,k++) q_inv[b]=-1;
+ if (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$.
+
+@<Compute the |q_inv| table@>=
+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).
+
+@<Choose or verify the |type|...@>=
+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)$.
+
+@<Set up a graph...@>=
+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: @<Assign labels from the set $\{0,1,\ldots,q-1,\infty\}$@>;@+break;
+ case 2: @<Assign labels for pairs of distinct elements@>;@+break;
+ default: @<Assign projective matrix labels@>;@+break;
+}
+
+@ Type 1 graphs are the easiest to label. We store a serial number
+in utility field |x.i|, using $q$ to represent $\infty$.
+
+@<Assign labels from the set $\{0,1,\ldots,q-1,\infty\}$@>=
+new_graph->format[4]='Z';
+for (a=0;a<q;a++) {
+ sprintf(name_buf,"%d",a);
+ v->name=gb_save_string(name_buf);
+ v->x.i=a;
+ v++;
+}
+v->name=gb_save_string("INF");
+v->x.i=q;
+v++;
+
+@ @<Private...@>=
+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.
+
+@<Assign labels for pairs...@>=
+for (a=0;a<q;a++)
+ for (aa=a+1;aa<=q;aa++) {
+ if (aa==q) sprintf(name_buf,"{%d,INF}",a);
+ else sprintf(name_buf,"{%d,%d}",a,aa);
+ v->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.
+
+@<Assign projective matrix labels@>=
+new_graph->format[5]='I';
+for (c=0;c<=q;c++)
+ for (b=0;b<q;b++)
+ for (a=1;a<=n_factor;a++) {
+ v->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.
+
+@<Type de...@>=
+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.
+
+@<Compute |p+1| generators...@>=
+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) @<Fill the |gen| table with special generators@>@;
+else @<Fill the |gen| table with representatives of all quaternions
+ having norm~|p|@>;
+if (gen_count!=max_gen_count) late_panic(bad_specs+7); /* |p| is not prime */
+
+@ @<Private...@>=
+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.
+
+@<Fill the |gen| table with representatives...@>=
+{@+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) @<Deposit the quaternions associated with $a+bi+cj+dk$@>;
+ @<Change the |gen| table to matrix format@>;
+}
+
+@ If |a>0| and |0<b<c<d|, we obtain 48 different classes of quaternions
+having the same norm by permuting $\{b,c,d\}$ in six ways and attaching
+signs to each permutation in eight ways. This happens, for example,
+when $p=71$ and $(a,b,c,d)=(6,1,3,5)$. Fewer quaternions arise when
+|a=0| or |0=b| or |b=c| or |c=d|.
+
+The inverse of the matrix corresponding to a quaternion is the matrix
+corresponding to the conjugate quaternion. Therefore a generating
+matrix $\pi_k$ will be its own inverse if and only if it comes from
+a quaternion with |a=0|.
+
+It is convenient to have a subroutine that deposits a new quaternion
+and its conjugate into the table of generators.
+
+@<Private...@>=
+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...@>=
+{
+ 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<c) {
+ deposit(a,c,b,d);@+deposit(a,-c,b,d);@+deposit(a,c,d,b);@+deposit(a,-c,d,b);
+ if (b) {
+ deposit(a,c,-b,d);@+deposit(a,-c,-b,d);@+deposit(a,c,d,-b);@+
+ deposit(a,-c,d,-b);
+ }
+ }
+ if (c<d) {
+ deposit(a,b,d,c);@+deposit(a,d,b,c);
+ if (b) {
+ deposit(a,-b,d,c);@+deposit(a,-b,d,-c);@+deposit(a,d,-b,c);@+
+ deposit(a,d,-b,-c);
+ }
+ if (c) {
+ deposit(a,b,d,-c);@+deposit(a,d,b,-c);
+ }
+ if (b<c) {
+ deposit(a,d,c,b);@+deposit(a,d,-c,b);
+ if (b) {
+ deposit(a,d,c,-b);@+deposit(a,d,-c,-b);
+ }
+ }
+ }
+}
+
+@ Once we've found the generators in quaternion form, we want to
+convert them to $2\times2$ matrices, using the correspondence mentioned
+earlier:
+$$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|).
+Appropriate values for $g$ and~$h$ can always be found by letting
+$g=\sqrt{\mathstrut k}$ and $h=\sqrt{\mathstrut q-1-k}$, where
+$k$ is the largest quadratic residue modulo~|q|. For if $q-1$ is
+not a quadratic residue, and if $k+1$ isn't a residue either, then
+$q-1-k$ must be a quadratic residue because it is congruent to the
+product $(q-1)(k+1)$ of nonresidues. (We will have |h=0| if and
+only if $q\bmod4=1$; |h=1| if and only if $q\bmod8=3$; $h=\sqrt{\mathstrut2}$
+if and only if $q\bmod24=7$ or 15; etc.)
+
+@<Change the |gen| table to matrix format@>=
+{@+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.
+
+@<Fill the |gen| table with special generators@>=
+{@+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 */
+
+@<Append the edges@>=
+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<new_graph->vertices+n;v++) {
+ register Vertex* u;
+ @<Compute the image, |u|, of |v|
+ under the permutation defined by |gen[k]|@>;
+ 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|.
+
+@<Compute the image, |u|, of |v| under the permutation defined by |gen[k]|@>=
+if (type<3) @<Compute the image, |u|, of |v| under the linear fractional
+ transformation defined by |gen[k]|@>@;
+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;
+ @<Compute the matrix product |(aa,bb;cc,dd)=(a,b;c,d)*(a0,a1;a2,a3)|@>;
+ 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;
+ @<Set |u| to the vertex whose label is |(a,b;c,d)|@>;
+}
+
+@ @<Compute the matrix product...@>=
+aa=(a*a0+b*a2)%q;
+bb=(a*a1+b*a3)%q;
+cc=(c*a0+d*a2)%q;
+dd=(c*a1+d*a3)%q;
+
+@ @<Set |u|...@>=
+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.
+
+@<Private...@>=
+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.
+
+@<Compute the image, |u|, of |v| under the linear fractional
+ transformation defined by |gen[k]|@>=
+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<aa? (a*(2*q-1-a))/2+aa-1:
+ (aa*(2*q-1-aa))/2+a-1);
+}
+
+@* Index. Here is a list that shows where the identifiers of this program are
+defined and used.