Sharedwww / tables / fundomain / fundomain.cOpen in CoCalc
Author: William A. Stein
1
/*
2
* This is a program to draw fundamantal domains
3
* this is version 1.0
4
*/
5
6
/*
7
* This program is written in c++, though actually it
8
* is written in c, with just a few c++ commands thrown
9
* in, and could easily be converted back to c, should
10
* you desire.
11
*
12
* This program was written on a linux machine.
13
* It also works on unix systems I've tried it on.
14
* A very slightly modified version works with windows.
15
* You are free to modify it to work on whatever system you're
16
* running.
17
*
18
*/
19
20
/*
21
* Simple instructions to compile and run under unix type systems:
22
*
23
* Give this file the name fundomain.c
24
* Compile with: gcc -o drawfd fundomain.c -lm
25
* This creates the file "drawfd"
26
* then run by typing: drawfd
27
* or, if that doesn't work try: ./drawfd
28
* after pressing return, the program runs, and
29
* you will be prompted for whether you want a picture of
30
* Gamma_0(N) or Gamma_1(N), and for the value of N,
31
* and for the name of the postscript putput file.
32
*
33
*/
34
35
/*
36
* This program was written by Helena A. Verrill, 2000
37
* Date of most recent modification: June 28 2000.
38
*/
39
40
/* #include </usr/include/math.h> */
41
#include <math.h>
42
#include <stdio.h>
43
44
#define PI 3.14159265358979
45
#define infinity 600
46
/* #define scale 150 */
47
48
49
#define max 50
50
51
/* note, in c on windows machine, the following three
52
* matrices were declared to be "const"
53
* but this does not seem to work on unix.
54
*/
55
56
int T[2][2] = {{1,1},{0,1}};
57
58
int Tm[2][2] = {{1,-1},{0,1}};
59
60
int S[2][2] = {{0,1},{-1,0}};
61
62
int proj_st[max][max][3];
63
int scale = 150;
64
65
void findfund(int N, int m[2][2], int done[max][max][2][3], int dist,int sbsp);
66
int find_red_uv(int sign, int m2[2][2], int N, int comp);
67
68
69
FILE *output; /* the output file */
70
FILE *outputtxt; /* the output file */
71
72
void
73
do_point(double x, double y, char *s)
74
{
75
fprintf(output, "%.8lg %.8lg %s ", x, y, s);
76
}
77
78
void
79
my_arc(double thrang, double startang, double radius)
80
{
81
double sign;
82
83
if (radius<0){
84
thrang=thrang*-1;
85
radius=-1*radius;}
86
fprintf(output, "/thr-ang %.8lg def\n",thrang);
87
if (thrang<0) {sign=-1;}
88
if (thrang>0) {sign=1;}
89
fprintf(output, "/sign %.8lg def\n",sign);
90
fprintf(output, "/start-ang %.8lg def\n",startang);
91
fprintf(output, "/radius %.8lg def\n",radius);
92
fprintf(output, "myarc\n");
93
}
94
95
void
96
closed_line(double x1, double y1, double x2, double y2)
97
{
98
fprintf(output, "newpath ");
99
do_point(x1, y1, "moveto");
100
do_point(x2, y2, "lineto");
101
fprintf(output, "stroke\n");
102
}
103
104
void
105
make_title(char thescript, char thetype, int Anumber)
106
{
107
fprintf(output,"/textfontsize %d def\n",4*Anumber);
108
fprintf(output,"/Helvetica-BoldOblique findfont textfontsize scalefont setfont\n");
109
fprintf(output,"100 -100 moveto\n");
110
fprintf(output,"(Gamma%c%c(%d)) show\n\n",thescript,thetype,Anumber);
111
}
112
113
114
115
void
116
open_line(double x1, double y1, double x2, double y2)
117
{
118
do_point(x1, y1, "lineto");
119
do_point(x2, y2, "lineto");
120
}
121
122
double
123
matmulty(double a, double b, double c, double d,double x, double y)
124
{
125
return a*x+b*y;
126
}
127
128
double
129
matmultx(double a, double b, double c, double d,double x, double y)
130
{
131
return c*x+d*y;
132
}
133
134
void
135
triA(double ang1, double ang2, double ang3, double r1, double r2, double r3, double x, double y,double red, double green, double blue)
136
{
137
fprintf(output, "gsave\n");
138
fprintf(output, "newpath\n");
139
fprintf(output, "%.8lg %.8lg moveto\n",x,y);
140
141
fprintf(output, "/r1 %.6lg abs def\n",r1);
142
fprintf(output, "/r2 %.6lg abs def\n",r2);
143
fprintf(output, "/r3 %.6lg abs def\n",r3);
144
fprintf(output, "/a1 %.6lg %d mul def\n",ang1,2*(r1>0)-1);
145
fprintf(output, "/a2 %.6lg %d mul def\n",ang2,2*(r2>0)-1);
146
fprintf(output, "/a3 %.6lg %d mul def\n",ang3,2*(r3>0)-1);
147
fprintf(output, "%.4lg %.4lg %.4lg %.4lg %.4lg mytriA fill\n",red,green,blue,x,y);
148
fprintf(output, "0 0 0 %.4lg %.4lg mytriA stroke grestore\n",x,y);
149
fprintf(output, "\n");
150
}
151
152
void
153
triB
154
(double ang1, double length2, double ang3, double r1, double r3, double x, double y,double red, double green, double blue)
155
{
156
fprintf(output, "gsave\n");
157
fprintf(output, "newpath\n");
158
fprintf(output, "%.8lg %.8lg moveto\n",x,y);
159
fprintf(output, "/r1 %.6lg abs def\n",r1);
160
fprintf(output, "/r3 %.6lg abs def\n",r3);
161
fprintf(output, "/a1 %.6lg %d mul def\n",ang1,2*(r1>0)-1);
162
fprintf(output, "/a3 %.6lg %d mul def\n",ang3,2*(r3>0)-1);
163
fprintf(output, "%.4lg %.4lg %.4lg 0 %.4lg abs %.4lg 0 mytriB fill\n",red,green,blue,length2,x);
164
fprintf(output, "0 0 0 0 %.4lg abs %.4lg 0 mytriB stroke grestore\n",length2,x);
165
fprintf(output, "\n");
166
}
167
168
void
169
triC
170
(double length1, double ang2, double ang3, double r2, double r3, double x, double y,double red, double green, double blue)
171
{
172
fprintf(output, "gsave\n");
173
fprintf(output, "newpath\n");
174
fprintf(output, "%.8lg %.8lg moveto\n",x,y);
175
fprintf(output, "/r2 %.6lg abs def\n",r2);
176
fprintf(output, "/r3 %.6lg abs def\n",r3);
177
fprintf(output, "/a2 %.6lg %d mul def\n",ang2,2*(r2>0)-1);
178
fprintf(output, "/a3 %.6lg %d mul def\n",ang3,2*(r3>0)-1);
179
fprintf(output, "%.4lg %.4lg %.4lg 0 %.4lg abs %.4lg 0 mytriC fill\n",red,green,blue,length1,x);
180
fprintf(output, "0 0 0 0 %.4lg abs %.4lg 0 mytriC stroke grestore\n",length1,x);
181
fprintf(output, "\n");
182
183
}
184
185
void
186
triD
187
(double ang1, double ang2, double r1, double r2, double x, double y,double red, double green, double blue)
188
{
189
190
fprintf(output, "gsave\n");
191
fprintf(output, "newpath\n");
192
fprintf(output, "%.8lg %.8lg moveto\n",x,y);
193
fprintf(output, "/r1 %.6lg abs def\n",r1);
194
fprintf(output, "/r2 %.6lg abs def\n",r2);
195
fprintf(output, "/a1 %.6lg %d mul def\n",ang1,2*(r1>0)-1);
196
fprintf(output, "/a2 %.6lg %d mul def\n",ang2,2*(r2>0)-1);
197
fprintf(output, "%.4lg %.4lg %.4lg %.4lg 0 mytriD fill\n",red,green,blue,x);
198
fprintf(output, "0 0 0 %.4lg 0 mytriD stroke grestore\n",x);
199
fprintf(output, "\n");
200
201
}
202
203
204
void
205
triE
206
(double r2, double ang2, double x, double el,double red, double green, double blue)
207
{
208
fprintf(output, "gsave\n");
209
fprintf(output, "newpath\n");
210
fprintf(output, "%.8lg 0 moveto\n",x);
211
fprintf(output, "/r2 %.6lg abs def\n",r2);
212
fprintf(output, "/a2 %.6lg %d mul def\n",ang2,2*(r2>0)-1);
213
fprintf(output, "%.4lg %.4lg %.4lg 0 %.4lg abs dup %.4lg mytriE fill\n",red,green,blue,el,x);
214
fprintf(output, "0 0 0 0 %.4lg abs dup %.4lg mytriE stroke grestore\n",el,x);
215
fprintf(output, "\n");
216
}
217
218
219
void
220
transtri(double a,double b,double c,double d, double red, double green, double blue)
221
{
222
double xm1,ym1,xp1,yp1,xmh,ymh,xph,yph,xinf,yinf,m1 ,p1 ,mh ,ph, inf;
223
double el1, el2, elI1, elI2, dist1, dist2a, dist2b, dist3;
224
double center1, center2, center3;
225
double ang1, ang2, ang3;
226
double radius1, radius2, radius3;
227
228
if(a*d-c*b>0){
229
xm1 = matmulty(a,b,c,d,-1,1); /* image of -1 */
230
ym1 = matmultx(a,b,c,d,-1,1);
231
xp1 = matmulty(a,b,c,d,1,1); /* image of 1 */
232
yp1 = matmultx(a,b,c,d,1,1);
233
xmh = matmulty(a,b,c,d,-1,2); /* image of -1/2 */
234
ymh = matmultx(a,b,c,d,-1,2);
235
xph = matmulty(a,b,c,d,1,2); /* image of 1/2 */
236
yph = matmultx(a,b,c,d,1,2);
237
xinf = matmulty(a,b,c,d,1,0); /* image of infinity */
238
yinf = matmultx(a,b,c,d,1,0);
239
m1=0; p1=0; mh=0; ph=0; inf=0;
240
if (ym1!=0){
241
m1 = scale * xm1/ym1;}
242
if (yp1!=0){
243
p1 = scale * xp1/yp1;}
244
if (ymh!=0){
245
mh = scale * xmh/ymh;}
246
if (yph!=0){
247
ph = scale * xph/yph;}
248
if (yinf!=0){
249
inf = scale * xinf/yinf;}
250
/* define real and imag parts of elliptic points */
251
el1 = scale * ((a+2*b)*(2*d+c)+3*a*c)/((2*d+c)*(2*d+c)+3*c*c);
252
el2 = scale * ((-a+2*b)*(2*d-c)+3*a*c)/((2*d-c)*(2*d-c)+3*c*c);
253
elI1 = scale * sqrt(3)*2*(a*d-b*c)/((2*d+c)*(2*d+c)+3*c*c);
254
elI2 = scale * sqrt(3)*2*(a*d-b*c)/((2*d-c)*(2*d-c)+3*c*c);
255
256
center1 = (inf+ph)/2;
257
center2 = (inf+mh)/2;
258
center3 = (m1+p1)/2;
259
radius1 = (inf-ph)/2;
260
radius2 = (m1-p1)/2;
261
radius3 = (mh-inf)/2;
262
263
dist1 = el1-center1;
264
dist3 = center2 - el2;
265
dist2a = center3-el1;
266
dist2b = el2-center3;
267
268
ang1 = (acos(dist1/radius1))*180/PI;
269
ang3 = (acos(dist3/radius3))*180/PI;
270
ang2 = -(asin(dist2a/radius2)+asin(dist2b/radius2))*180/PI;
271
272
if (ym1*yp1*ymh*yph*yinf!=0){
273
triA(ang1,ang2,ang3,radius1,radius2,radius3,inf,0,red,green,blue);
274
}
275
if (ym1*yp1==0){
276
triB(ang1,elI2-elI1,ang3,radius1,radius3,inf,0,red,green,blue);
277
}
278
if (yph==0){
279
triC(elI1,ang2,ang3,radius2,radius3,inf,0,red,green,blue);
280
}
281
if (ymh==0){
282
triD(ang1,ang2,radius1,radius2,inf,0,red,green,blue);
283
}
284
if (yinf==0){
285
triE(radius2,ang2,el1,elI1,red,green,blue);
286
}
287
}
288
}
289
290
291
void
292
lots_of_triangles()
293
{
294
double a,b,c,d,e;
295
e=1;
296
for (a = 1; a < 7; a = a + 1) {
297
for (b = -1; b < 6; b = b + 1) {
298
for (c = -1; c < 5; c = c + 1) {
299
for (d = -1; d < 4 ; d = d + 1) {
300
e=e+0.4;
301
transtri(a,b,c,d,1,0.5+cos(e),0.5+sin(e));
302
}}}}
303
}
304
305
void
306
more_trianglesa()
307
{
308
double a,e;
309
e=0;
310
for (a = -5; a < 5; a = a + 1) {
311
e=e+.1;
312
transtri(0,1,-1,a,(1-e)*(1-e)+1-(.5-e)*(.5-e),1-e*e,e*e);
313
transtri(a,1,-1,0,(1-e)*(1-e)+1-(.5-e)*(.5-e),1-e*e,e*e);
314
transtri(1,a,0,1,(1-e)*(1-e)+1-(.5-e)*(.5-e),1-e*e,e*e);
315
transtri(1,0,a,1,(1-e)*(1-e)+1-(.5-e)*(.5-e),1-e*e,e*e);
316
}
317
}
318
319
void
320
more_triangles()
321
{
322
double a,e;
323
e=0;
324
for (a = -17; a < 18; a = a + 1) {
325
e=1-e;
326
transtri(0,1,-1,a,e,0,1);
327
if(a*a!=1){transtri(1,0,a,1,1-e,0,1);}
328
}
329
e=0;
330
for (a = -5; a < 6; a = a + 1) {
331
e=1-e;
332
if(a<-2 || a>1) {transtri(1,1,a,a+1,e,0,1);}
333
if(a<-1 || a>2) {transtri(1,-1,a,-a+1,e,0,1);}
334
}
335
}
336
337
void
338
row()
339
{
340
fprintf(output, "/block{\n");
341
more_triangles();
342
343
fprintf(output, "} def\n");
344
fprintf(output, "\n");
345
fprintf(output, "\n");
346
fprintf(output,"-120 0 translate \n");
347
fprintf(output,"/col 0 def \n");
348
fprintf(output,"4 {\n");
349
fprintf(output,"block\n");
350
fprintf(output,"150 0 translate \n");
351
fprintf(output,"} repeat\n");
352
353
}
354
355
356
void
357
more_trianglesb()
358
{
359
double a,e;
360
e=0;
361
for (a = -5; a < 5; a = a + 1) {
362
e=e+.1;
363
transtri(0,1,-1,a,e,0,1);
364
transtri(a,1,-1,0,1,0,e);
365
transtri(1,a,0,1,1,e,0);
366
transtri(1,0,a,1,0,1,e);
367
}
368
}
369
370
371
int coprime(int a[3])
372
{
373
int b[3];
374
375
/* this function only works if it is fed non negative integers! */
376
if((a[1]*a[1]-1)*(a[0]*a[0]-1)*(a[2]*a[2]-1)==0){return 1;}
377
if(a[0]==0)
378
{
379
if(a[1]==0)
380
{return 0;}
381
else{a[0]=a[1];}
382
}
383
b[0]=a[1]-(a[1]/a[0])*a[0];
384
b[1]=a[2]-(a[2]/a[0])*a[0];
385
b[2]=a[0];
386
return coprime(b);
387
}
388
389
390
void
391
Proj_reps(int N, int mat[10][10])
392
{
393
int j,k,p;
394
395
for(j=1;j<N;++j){
396
for(k=1;k<N;++k){
397
if(mat[j][k]==1){
398
for (p=2;p<N;++p){
399
if(mat[0][p]==1){
400
mat[j*p-((j*p)/N)*N][k*p-((k*p)/N)*N]=0;
401
}
402
}
403
}
404
}
405
}
406
}
407
408
409
410
void
411
matinit(int m1[2][2], int a, int b, int c, int d)
412
{
413
414
m1[0][0]=a;
415
m1[0][1]=b;
416
m1[1][0]=c;
417
m1[1][1]=d;
418
419
}
420
421
422
423
424
425
void
426
matmult(int c[2][2], int a[2][2], int b[2][2])
427
{
428
429
c[0][0]=a[0][0]*b[0][0]+a[0][1]*b[1][0];
430
c[0][1]=a[0][0]*b[0][1]+a[0][1]*b[1][1];
431
c[1][0]=a[1][0]*b[0][0]+a[1][1]*b[1][0];
432
c[1][1]=a[1][0]*b[0][1]+a[1][1]*b[1][1];
433
434
}
435
436
437
void
438
printmat(int m1[2][2])
439
{
440
int i,j;
441
442
for (i=0;i<2;++i)
443
{
444
printf("\n");
445
for (j=0;j<2;++j)
446
{
447
printf(" %d",m1[i][j]);
448
}
449
}
450
printf("\n");
451
}
452
453
454
455
456
457
void
458
bezout(int pair[2])
459
{
460
int u,v,v1,i,j;
461
int newmat[2][2], runningmat[2][2], newrun[2][2];
462
463
464
465
u=pair[0];
466
v=pair[1];
467
468
matinit(newmat,0,1,1,0);
469
printmat(newmat);
470
matinit(runningmat,1,0,0,1);
471
472
printmat(runningmat);
473
474
while(u>0)
475
{
476
477
478
newmat[1][1]=-(v/u);
479
480
matmult(newrun,runningmat,newmat);
481
482
for (i=0;i<2;++i){
483
for (j=0;j<2;++j){
484
runningmat[i][j]=newrun[i][j];
485
486
}
487
}
488
489
printmat(runningmat);
490
491
v1=u;
492
u=v-(v/u)*u;
493
v=v1;
494
495
496
}
497
498
printmat(runningmat);
499
500
501
pair[0]=runningmat[1][0];
502
pair[1]=runningmat[0][0];
503
504
}
505
506
507
508
void
509
listpairs(int N, int pairlist[10][10], int matlist[10][10][2])
510
{
511
int i,j, b[3];
512
int pair[2];
513
514
/* matlist contains the bottom parts of the matrix with top tow [i,j] */
515
516
N=6;
517
518
for (i=0;i<N;++i){
519
for (j=0;j<N;++j){
520
pairlist[i][j]=0;
521
b[0]=i;
522
b[1]=j;
523
b[2]=N;
524
pairlist[i][j]=coprime(b);
525
526
}
527
}
528
529
530
Proj_reps(N,pairlist);
531
532
for (i=0;i<N;++i){
533
for (j=0;j<N;++j){
534
if (i>0 && j>0 && pairlist[i][j]!=0)
535
{ pair[0]=i;
536
pair[1]=j;
537
538
bezout(pair);
539
matlist[i][j][0]=-pair[1];
540
matlist[i][j][1]=pair[0];
541
542
}
543
else {
544
matlist[i][j][0]=0;
545
matlist[i][j][1]=0;
546
}
547
}
548
}
549
550
}
551
552
void
553
proj_eq(int N, int kind)
554
{
555
556
int i,j,k,u,v;
557
558
559
for(i=1;i<N;++i) {
560
for(j=1;j<N;++j) {
561
for(k=0;k<2;++k) {
562
proj_st[i][j][k]=0;
563
}}}
564
565
566
if (kind==0){
567
for(i=1;i<N;++i) {
568
proj_st[0][i][1]=1;
569
proj_st[i][0][0]=1;
570
}
571
572
573
574
for(i=1;i<N;++i) {
575
for(j=1;j<N;++j) {
576
if(proj_st[i][j][0]==0 && proj_st[i][j][1]==0){
577
proj_st[i][j][0]=i;
578
proj_st[i][j][1]=j;
579
for(k=1;k<N;++k){
580
u=k*i-((k*i)/N)*N;
581
v=k*j-((k*j)/N)*N;
582
if(u!=0 && v!=0){
583
if (proj_st[u][v][0]==0&&proj_st[u][v][1]==0){
584
proj_st[u][v][0]=i;
585
proj_st[u][v][1]=j;
586
}}}}}}
587
588
}
589
590
if(kind==1){
591
for(i=1;i<N/2;++i) {
592
proj_st[0][i][1]=i;
593
proj_st[i][0][0]=i;
594
}
595
for(i=N/2;i<N;++i) {
596
proj_st[0][i][1]=N-i;
597
proj_st[i][0][0]=N-i;
598
}
599
600
601
602
for(i=1;i<N;++i) {
603
for(j=1;j<N;++j) {
604
if(j<=N/2) {
605
proj_st[i][j][0]=i;
606
proj_st[i][j][1]=j;}
607
else {
608
proj_st[i][j][0]=N-i;
609
proj_st[i][j][1]=N-j;
610
}
611
}
612
613
}
614
615
if ((N/2)*2==N){
616
for(i=1;i<N/2;++i) {
617
proj_st[N/2][i][0]=N/2;
618
proj_st[N/2][i][1]=i;
619
proj_st[i][N/2][0]=i;
620
proj_st[i][N/2][1]=N/2;
621
}
622
for(i=N/2;i<N;++i) {
623
proj_st[N/2][i][0]=N/2;
624
proj_st[N/2][i][1]=N-i;
625
proj_st[i][N/2][0]=N-i;
626
proj_st[i][N/2][1]=N/2;
627
}
628
}
629
630
631
}
632
633
634
}
635
636
637
638
639
640
641
void
642
shorten(int N, int done[max][max][2][3])
643
{
644
645
int s,i,j, u,v,u1,v1;
646
int m[2][2],m1[2][2],m2[2][2],m3[2][2];
647
648
649
650
651
s=0;
652
653
for(i=0;i<N;++i){
654
for(j=0;j<N;++j){
655
if (done[i][j][0][2] > s) {
656
s=done[i][j][0][2]; u=i; v=j;
657
}
658
}}
659
660
if (s<1) return;
661
for(i=0;i<2;++i){for(j=0;j<2;++j){m[i][j]=done[u][v][i][j];}}
662
663
664
matmult(m1,S,m);
665
666
u1=find_red_uv(1,m1,N,0);
667
v1=find_red_uv(1,m1,N,1);
668
669
if (done[u1][v1][0][0]==m1[0][0] && done[u1][v1][0][1]==m1[0][1] && done[u1][v1][1][0]==m1[1][0] && done[u1][v1][1][1]==m1[1][1]) {
670
if (done[u1][v1][0][2]<s-1) {
671
for(i=0;i<2;++i){for(j=0;j<2;++j){m2[i][j]=done[u1][v1][i][j];}}
672
matmult(m3,S,m2);
673
for(i=0;i<2;++i){for(j=0;j<2;++j){done[u][v][i][j]=m3[i][j];}}
674
done[u][v][0][2]=done[u1][v1][0][2]+1;
675
shorten(N,done);
676
}}
677
678
matmult(m1,T,m);
679
680
u1=find_red_uv(1,m1,N,0);
681
v1=find_red_uv(1,m1,N,1);
682
683
if (done[u1][v1][0][0]==m1[0][0] && done[u1][v1][0][1]==m1[0][1] && done[u1][v1][1][0]==m1[1][0] && done[u1][v1][1][1]==m1[1][1]) {
684
if (done[u1][v1][0][2]<s-1) {
685
for(i=0;i<2;++i){for(j=0;j<2;++j){m2[i][j]=done[u1][v1][i][j];}}
686
matmult(m3,Tm,m2);
687
for(i=0;i<2;++i){for(j=0;j<2;++j){done[u][v][i][j]=m3[i][j];}}
688
done[u][v][0][2]=done[u1][v1][0][2]+1;
689
shorten(N,done);
690
}}
691
692
693
matmult(m1,Tm,m);
694
695
u1=find_red_uv(1,m1,N,0);
696
v1=find_red_uv(1,m1,N,1);
697
698
if (done[u1][v1][0][0]==m1[0][0] && done[u1][v1][0][1]==m1[0][1] && done[u1][v1][1][0]==m1[1][0] && done[u1][v1][1][1]==m1[1][1]) {
699
if (done[u1][v1][0][2]<s-1) {
700
for(i=0;i<2;++i){for(j=0;j<2;++j){m2[i][j]=done[u1][v1][i][j];}}
701
matmult(m3,T,m2);
702
for(i=0;i<2;++i){for(j=0;j<2;++j){done[u][v][i][j]=m3[i][j];}}
703
done[u][v][0][2]=done[u1][v1][0][2]+1;
704
shorten(N,done);
705
}}
706
707
return;
708
709
710
711
712
713
}
714
715
716
717
int
718
find_red_uv(int sign, int m2[2][2], int N, int comp)
719
{
720
721
int uu,vv,u,v;
722
723
u=sign*(m2[0][0]-(m2[0][0]/N)*N);
724
v=sign*(m2[0][1]-(m2[0][1]/N)*N);
725
if(u<0) u=N+u;
726
if(v<0) v=N+v;
727
uu=u;
728
vv=v;
729
return proj_st[uu][vv][comp];
730
731
}
732
733
734
735
void
736
nearest(int m[2][2], int A[2][2], int N, int done[max][max][2][3], int dist,int sbsp)
737
{
738
739
int u,v,i,j;
740
int dist1,dist3,dist4;
741
int m2[2][2],m3[2][2],m4[2][2];
742
int u3,v3,u4,v4;
743
int C[2][2], D[2][2], Cm[2][2], Dm[2][2];
744
745
dist3=max*max;
746
dist4=max*max;
747
748
749
if (A[0][0]==T[0][0] && A[0][1]==T[0][1] && A[0][1]==T[0][1] && A[0][1]==T[0][1]) {
750
for (i=0;i<2;++i){for (j=0;j<2;++j){C[i][j]=S[i][j];}}
751
for (i=0;i<2;++i){for (j=0;j<2;++j){D[i][j]=T[i][j];}}
752
for (i=0;i<2;++i){for (j=0;j<2;++j){Cm[i][j]=S[i][j];}}
753
for (i=0;i<2;++i){for (j=0;j<2;++j){Dm[i][j]=Tm[i][j];}}
754
}
755
756
if (A[0][0]==S[0][0] && A[0][1]==S[0][1] && A[0][1]==S[0][1] && A[0][1]==S[0][1]) {
757
for (i=0;i<2;++i){for (j=0;j<2;++j){C[i][j]=T[i][j];}}
758
for (i=0;i<2;++i){for (j=0;j<2;++j){D[i][j]=Tm[i][j];}}
759
for (i=0;i<2;++i){for (j=0;j<2;++j){Cm[i][j]=Tm[i][j];}}
760
for (i=0;i<2;++i){for (j=0;j<2;++j){Dm[i][j]=T[i][j];}}
761
}
762
763
if (A[0][0]==Tm[0][0] && A[0][1]==Tm[0][1] && A[0][1]==Tm[0][1] && A[0][1]==Tm[0][1]) {
764
for (i=0;i<2;++i){for (j=0;j<2;++j){C[i][j]=S[i][j];}}
765
for (i=0;i<2;++i){for (j=0;j<2;++j){D[i][j]=Tm[i][j];}}
766
for (i=0;i<2;++i){for (j=0;j<2;++j){Cm[i][j]=S[i][j];}}
767
for (i=0;i<2;++i){for (j=0;j<2;++j){Dm[i][j]=T[i][j];}}
768
}
769
770
matmult(m2,m,A);
771
772
773
u=find_red_uv(1,m2,N,0);
774
v=find_red_uv(1,m2,N,1);
775
776
if (done[u][v][0][0]==0 && done[u][v][0][1]==0) {
777
/*&& done[u1][v1][0][0]==0 && done[u1][v1][0][1]==0) {*/
778
779
if (sbsp==0 && m2[1][0]==0) {
780
dist1=0;
781
findfund(N,m2,done,dist1,sbsp);
782
}
783
else {
784
785
786
787
dist1=dist+1;
788
789
/* note, divide by 2 on following row because even though the distance as a graph increases, the size gets smaller at
790
a non linear rate. Possibly should change this to a different function that just divide by 2 */
791
792
if (sbsp==1 && m2[0][0]==1) {
793
dist1 = m2[0][1];
794
if (dist1*dist1 < 5) dist1=dist1/2;
795
if (dist1<0) dist1=-dist1;
796
}
797
798
matmult(m3,m2,C);
799
matmult(m4,m2,D);
800
801
802
u3=find_red_uv(1,m3,N,0);
803
v3=find_red_uv(1,m3,N,1);
804
/* u31=find_red_uv(-1,m3,N,0);
805
v31=find_red_uv(-1,m3,N,1);*/
806
if (done[u3][v3][0][0]!=0 || done[u3][v3][0][1]!=0) {
807
dist3=done[u3][v3][0][2];
808
for (i=0;i<2;++i){for (j=0;j<2;++j){m3[i][j]=done[u3][v3][i][j];}}
809
}
810
/* if (done[u31][v31][0][0]!=0 || done[u31][v31][0][1]!=0) {
811
dist3=done[u31][v31][0][2];
812
for (i=0;i<2;++i){for (j=0;j<2;++j){m3[i][j]=done[u31][v31][i][j];}}
813
}
814
*/
815
816
817
u4=find_red_uv(1,m4,N,0);
818
v4=find_red_uv(1,m4,N,1);
819
/* u41=find_red_uv(-1,m4,N,0);
820
v41=find_red_uv(-1,m4,N,1); */
821
if (done[u4][v4][0][0]!=0 || done[u4][v4][0][1]!=0) {
822
dist4=done[u4][v4][0][2];
823
for (i=0;i<2;++i){for (j=0;j<2;++j){m4[i][j]=done[u4][v4][i][j];}}
824
}
825
/* if (done[u41][v41][0][0]!=0 || done[u41][v41][0][1]!=0) {
826
dist4=done[u41][v41][0][2];
827
for (i=0;i<2;++i){for (j=0;j<2;++j){m4[i][j]=done[u4][v4][i][j];}}
828
}
829
*/
830
831
832
if (dist3 < dist1 && dist3 < dist4) {
833
matmult(m2,m3,Cm);
834
findfund(N, m2, done, dist3+1,sbsp);
835
}
836
837
else if (dist4 <= dist3 && dist4<dist1)
838
{
839
matmult(m2,m4,Dm);
840
findfund(N, m2, done, dist4+1,sbsp);
841
842
}
843
else {
844
findfund(N,m2,done,dist1,sbsp);
845
}
846
847
848
}
849
}
850
851
}
852
853
854
855
void
856
findfund(int N, int m[2][2], int done[max][max][2][3], int dist,int sbsp)
857
{
858
int u,v,Nu,Nv,u1,v1,i,j,ua,va;
859
860
861
u=m[0][0]-(m[0][0]/N)*N;
862
v=m[0][1]-(m[0][1]/N)*N;
863
if(u<0) u=N+u;
864
if(v<0) v=N+v;
865
866
867
868
869
ua=proj_st[u][v][0];
870
va=proj_st[u][v][1];
871
872
u=ua;
873
v=va;
874
875
876
Nu=(N-u)-((N-u)/N)*N;
877
Nv=(N-v)-((N-v)/N)*N;
878
879
u1=proj_st[Nu][Nv][0];
880
v1=proj_st[Nu][Nv][1];
881
882
883
for (i=0;i<2;++i){
884
for (j=0;j<2;++j){
885
886
done[u][v][i][j]=m[i][j];
887
if (done[u1][v1][i][j]!=done[u][v][i][j]);
888
done[u1][v1][i][j]=m[i][j];
889
890
}
891
}
892
done[u][v][0][2]=dist;
893
894
/* fprintf(outputtxt," u1 %d, %d, %d %d %d\n",u,v, done[u][v][0][0], done[u][v][0][1],dist); */
895
896
/* CASE 1 */
897
898
899
nearest(m,T,N,done,dist,sbsp);
900
901
/* CASE 2 */
902
903
nearest(m,S,N,done,dist,sbsp);
904
nearest(m,Tm,N,done,dist,sbsp);
905
906
}
907
908
void
909
topmatter(int N,int sbsp)
910
{
911
double sc;
912
int tra;
913
int scale;
914
915
scale = 150; /* scale tells how accurately the circles are drawn */
916
if (sbsp == 1) scale=400;
917
918
tra=60;
919
sc=(float)3.5/N;
920
if(sbsp==1) {sc = 3; tra=300;}
921
922
923
fprintf(output, "%%!PS \n");
924
if(sbsp==1) {fprintf(output, "0.2 setlinewidth \n");}
925
fprintf(output, "%d 300 translate \n",tra);
926
fprintf(output, "\n");
927
fprintf(output, "%.4lg %.4lg scale\n",sc,sc);
928
fprintf(output, "\n");
929
fprintf(output, "/infinity 200 def\n");
930
fprintf(output, " \n");
931
fprintf(output, "/circle { \n");
932
fprintf(output, " 0 0 rlineto \n");
933
fprintf(output, " 0 chord rlineto \n");
934
fprintf(output, "{\n");
935
fprintf(output, " rotate \n");
936
fprintf(output, " 0 chord rlineto \n");
937
fprintf(output, "} repeat\n");
938
fprintf(output, "}bind def \n");
939
fprintf(output, "\n");
940
fprintf(output, "\n");
941
fprintf(output, "/myarc{\n");
942
fprintf(output, "/N %d def \n",scale);
943
fprintf(output, "/turn thr-ang N div def\n");
944
fprintf(output, "/chord sign 2 radius turn 2 div sin mul mul mul def\n");
945
fprintf(output, "start-ang rotate\n");
946
fprintf(output, "N {turn}repeat N circle\n");
947
fprintf(output, "180 rotate \n");
948
fprintf(output, "} def\n");
949
fprintf(output, "\n");
950
fprintf(output, "\n");
951
fprintf(output,"/mytriA{\n");
952
fprintf(output,"newpath\n");
953
fprintf(output,"moveto\n");
954
fprintf(output,"/thr-ang a1 def\n");
955
fprintf(output,"/sign a1 a1 abs div round def\n");
956
fprintf(output,"/start-ang 0 def\n");
957
fprintf(output,"/radius r1 def\n");
958
fprintf(output,"myarc\n");
959
fprintf(output,"/thr-ang a2 def\n");
960
fprintf(output,"/sign a2 a2 abs div round def\n");
961
fprintf(output,"/start-ang 60 def\n");
962
fprintf(output,"/radius r2 def\n");
963
fprintf(output,"myarc\n");
964
fprintf(output,"/thr-ang a3 def\n");
965
fprintf(output,"/sign a3 a3 abs div round def\n");
966
fprintf(output,"/start-ang 60 def\n");
967
fprintf(output,"/radius r3 def\n");
968
fprintf(output,"myarc\n");
969
fprintf(output,"setrgbcolor\n");
970
fprintf(output,"} def\n");
971
fprintf(output, "\n");
972
fprintf(output, "\n");
973
fprintf(output,"/mytriB{\n");
974
fprintf(output,"newpath\n");
975
fprintf(output,"moveto\n");
976
fprintf(output,"/thr-ang a1 def\n");
977
fprintf(output,"/sign a1 a1 abs div round def\n");
978
fprintf(output,"/start-ang 0 def\n");
979
fprintf(output,"/radius r1 def\n");
980
fprintf(output,"myarc\n");
981
fprintf(output,"60 rotate \n");
982
fprintf(output,"rlineto \n");
983
fprintf(output,"/thr-ang a3 def\n");
984
fprintf(output,"/sign a3 a3 abs div round def\n");
985
fprintf(output,"/start-ang 240 def\n");
986
fprintf(output,"/radius r3 def\n");
987
fprintf(output,"myarc\n");
988
fprintf(output,"setrgbcolor\n");
989
fprintf(output,"} def\n");
990
fprintf(output, "\n");
991
fprintf(output, "\n");
992
fprintf(output,"/mytriC{\n");
993
fprintf(output,"newpath\n");
994
fprintf(output,"moveto\n");
995
fprintf(output,"rlineto \n");
996
fprintf(output,"/thr-ang a2 def\n");
997
fprintf(output,"/sign a2 a2 abs div round def\n");
998
fprintf(output,"/start-ang 240 def\n");
999
fprintf(output,"/radius r2 def\n");
1000
fprintf(output,"myarc\n");
1001
fprintf(output,"/thr-ang a3 def\n");
1002
fprintf(output,"/sign a3 a3 abs div round def\n");
1003
fprintf(output,"/start-ang 60 def\n");
1004
fprintf(output,"/radius r3 def\n");
1005
fprintf(output,"myarc\n");
1006
fprintf(output,"setrgbcolor\n");
1007
fprintf(output,"} def\n");
1008
fprintf(output, "\n");
1009
fprintf(output, "\n");
1010
fprintf(output,"/mytriE{\n");
1011
fprintf(output,"newpath\n");
1012
fprintf(output,"1 copy %d moveto\n",infinity);
1013
fprintf(output,"exch lineto \n");
1014
fprintf(output,"/thr-ang a2 def\n");
1015
fprintf(output,"/sign a2 a2 abs div round def\n");
1016
fprintf(output,"/start-ang 60 def\n");
1017
fprintf(output,"/radius r2 def\n");
1018
fprintf(output,"myarc\n");
1019
fprintf(output,"60 rotate \n");
1020
fprintf(output,"%d exch sub rlineto \n",infinity);
1021
fprintf(output,"setrgbcolor\n");
1022
fprintf(output,"} def \n");
1023
fprintf(output,"\n");
1024
fprintf(output,"\n");
1025
fprintf(output,"/mytriD{\n");
1026
fprintf(output,"newpath\n");
1027
fprintf(output,"moveto\n");
1028
fprintf(output,"/thr-ang a1 def\n");
1029
fprintf(output,"/sign a1 a1 abs div round def\n");
1030
fprintf(output,"/start-ang 0 def\n");
1031
fprintf(output,"/radius r1 def\n");
1032
fprintf(output,"myarc\n");
1033
fprintf(output,"/thr-ang a2 def\n");
1034
fprintf(output,"/sign a2 a2 abs div round def\n");
1035
fprintf(output,"/start-ang 60 def\n");
1036
fprintf(output,"/radius r2 def\n");
1037
fprintf(output,"myarc\n");
1038
fprintf(output,"closepath\n");
1039
fprintf(output,"setrgbcolor\n");
1040
fprintf(output,"240 rotate \n");
1041
fprintf(output,"} def\n");
1042
fprintf(output,"\n");
1043
fprintf(output,"\n");
1044
1045
}
1046
1047
1048
1049
int dopsfile()
1050
{
1051
int i,j,k,l, sbsp, N, gotsbsp, gotN, gotName;
1052
char ch, thescript, thetype;
1053
int kind, gotkind;
1054
int temp;
1055
int done[max][max][2][3],m1[2][2],mt[2][2],mm[2][2];
1056
char GammaNps[256];
1057
1058
gotN = 0;
1059
1060
printf("enter N (0 < N < %d) (or type q at any time to quit)\n",max);
1061
1062
N = 0;
1063
while(gotN==0){
1064
if (scanf("%d",&N)!=EOF)
1065
{if (N>0 && N <50) gotN = 1;}
1066
if (gotN==0){
1067
scanf("%c",&ch);
1068
if (tolower(ch)=='q') return 0;
1069
if (ch!='\n') while (ch!='\n') scanf("%c",&ch);
1070
printf("please enter q or an integer N, 0 < N < 50\n");
1071
}
1072
}
1073
1074
1075
if(N==0) return 0;
1076
if(N>0&&N<=max){
1077
gotsbsp = 0;
1078
printf("upper or lower subscript? (u or l)\n");
1079
scanf("%c",&ch);
1080
1081
while (gotsbsp==0){
1082
if (tolower(ch)=='u') {sbsp=0; thescript='^'; gotsbsp=1;}
1083
else if (tolower(ch)=='l'){sbsp=1; thescript='_'; gotsbsp=1;}
1084
else if (tolower(ch)=='q') return 0;
1085
else if (ch!='\n') printf("please enter u or l\n");
1086
while (ch!='\n') scanf("%c",&ch);
1087
if (gotsbsp==0) scanf("%c",&ch);}
1088
1089
1090
gotkind = 0;
1091
printf("Gamma-0 or Gamma-1? (0 or 1)\n");
1092
scanf("%c",&ch);
1093
1094
while (gotkind==0){
1095
if (tolower(ch)=='0') {kind=0; thetype='0'; gotkind=1;}
1096
else if (tolower(ch)=='1') {kind=1; thetype='1'; gotkind=1;}
1097
else if (tolower(ch)=='q') return 0;
1098
else if (ch!='\n') printf("please enter 0 or 1, or q to quit\n");
1099
while (ch!='\n') scanf("%c",&ch);
1100
if (gotkind==0) scanf("%c",&ch);}
1101
1102
1103
printf("enter ps filename (including .ps part) \n");
1104
scanf("%s",GammaNps);
1105
1106
gotName=0;
1107
while (gotName==0){
1108
output=fopen(GammaNps,"r");
1109
if (fopen(GammaNps,"r")!=NULL){
1110
fclose(output);
1111
printf("%s already exists - ",GammaNps);
1112
printf("are you sure you want to over write it?\n");
1113
printf("enter n for no, and anything else to ");
1114
printf("continue to write the file\n");
1115
scanf("%c",&ch);
1116
while (ch!='\n') scanf("%c",&ch);
1117
scanf("%c",&ch);
1118
if (ch!='n') gotName=1;
1119
else {
1120
printf("enter ps filename (including .ps part) \n");
1121
scanf("%s",GammaNps);
1122
}
1123
}
1124
else gotName=1;
1125
}
1126
1127
printf("\n");
1128
printf("I'm now writing to the file %s\n",GammaNps);
1129
printf("when I've finished I'll ask you to tell me another N.\n");
1130
printf("\n");
1131
1132
1133
/*Output for the start of the PostScript file. */
1134
output = fopen(GammaNps,"w");
1135
if (!output) {
1136
fprintf(stderr, "I couldn't open the output file!\n");
1137
return 1;
1138
}
1139
1140
1141
/* draw base line */
1142
1143
for(i=0;i<max;++i){
1144
for(j=0;j<max;++j){
1145
for(k=0;k<2;++k){
1146
for(l=0;l<2;++l){
1147
done[i][j][k][l]=0;
1148
}
1149
}
1150
}
1151
}
1152
topmatter(N,sbsp);
1153
1154
closed_line(-1*infinity,0,infinity*N,0);
1155
1156
1157
make_title(thescript,thetype,N);
1158
1159
proj_eq(N, kind);
1160
1161
1162
1163
/* for(i=0;i<N;++i){
1164
fprintf(outputtxt,"next \n");
1165
for(j=0;j<N;++j){
1166
fprintf(outputtxt,"[%d, %d]",proj_st[i][j][0],proj_st[i][j][1]);
1167
}
1168
}
1169
1170
*/
1171
1172
1173
if (sbsp==0){
1174
matinit(m1,1,0,0,1);
1175
findfund(N,m1,done,0,sbsp);
1176
}
1177
if (sbsp==1){
1178
matinit(m1,1,-N/3,0,1);
1179
findfund(N,m1,done,N/3,sbsp);
1180
}
1181
1182
temp=0;
1183
1184
1185
1186
shorten(N,done);
1187
1188
1189
1190
for(i=0;i<N;++i){
1191
1192
1193
for(j=0;j<N;++j){
1194
1195
if(done[i][j][0][0]*done[i][j][1][1]-done[i][j][1][0]*done[i][j][0][1]==1){
1196
if(sbsp==0){
1197
transtri(done[i][j][0][0],done[i][j][0][1],done[i][j][1][0],done[i][j][1][1],1,0,0);
1198
1199
temp=temp+1;
1200
}
1201
if(sbsp==1) {
1202
mm[0][0]=done[i][j][0][0];mm[0][1]=done[i][j][0][1];mm[1][0]=done[i][j][1][0];mm[1][1]=done[i][j][1][1];
1203
matmult(mt,S,mm);
1204
transtri(mt[0][0],mt[0][1],mt[1][0],mt[1][1],1,0,0);
1205
}
1206
1207
1208
}
1209
}
1210
}
1211
1212
/* printf("\n index: %d \n",temp); */
1213
1214
fprintf(output, "\n\n showpage\n");
1215
fclose(output);
1216
printf("done!\n");
1217
}
1218
1219
1220
return 1;
1221
}
1222
1223
int
1224
main()
1225
{
1226
1227
int i;
1228
1229
1230
for (i=0;i<5;++i){
1231
if(dopsfile()==0) return 0;
1232
}
1233
1234
return 0;
1235
1236
}
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246