gb_raman.w

模拟器提供了一个简单易用的平台

@<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(%ld,%ld,%lu,%lu)",p,q,type,reduce);
strcpy(new_graph->util_types,"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->util_types[4]='Z';
for (a=0;a<q;a++) {
  sprintf(name_buf,"%ld",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,"{%ld,INF}",a);
    else sprintf(name_buf,"{%ld,%ld}",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 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 are
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->util_types[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,"(%ld,%ld;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;
            /* determinant is $a^2$ or $a$ */
        sprintf(name_buf,"(%ld,%ld;1,%ld)",b,v->y.I,c);
      }
      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, for our purposes,
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 of the basic identity
$\overline{\mathstrut\alpha\beta}=\overline{\mathstrut\beta}
\,\overline{\mathstrut\alpha}$ and the fact that
$\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 having odd norm. This algorithm
makes it possible to prove that integral quaternions whose coefficients are
relatively prime can be canonically 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
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,
when we omit quaternions whose norm is a multiple of~$q$;
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
@^Jacobi, Carl Gustav Jacob@>
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 long 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, so we allocate an extra slot in
the |gen| table to hold a possible overflow entry.

@<Compute |p+1| generators...@>=
gen=gb_typed_alloc(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 $\null\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; /* $p-a^2$, $p-a^2-b^2$ */
  long 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;cc>=0;cc-=(c+1)<<3,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 long gen_count; /* the next available quaternion slot */
static unsigned long 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$ */
{