summaryrefslogtreecommitdiff
path: root/support/graphbase/gb_raman.w
blob: 243796ef1587fbddb39e324b65bb417f26b1b1b1 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
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.