gb_raman.w

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

  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 the formulas
$$g=\sqrt{\mathstrut k}\qquad\hbox{and}\qquad 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 long g,h;
  long 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 Patrick~Chiu:
@^Chiu, Patrick@>
$$\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@>=
{@+long 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 might 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 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 {\sc 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--) {@+long 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,1L);
          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,1L);
        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 a00=gen[k].a0,a01=gen[k].a1,a10=gen[k].a2,a11=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)*(a00,a01;a10,a11)|@>;
  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*a00+b*a10)%q;
bb=(a*a01+b*a11)%q;
cc=(c*a00+d*a10)%q;
dd=(c*a01+d*a11)%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$ */
  long 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 long 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.