-
Notifications
You must be signed in to change notification settings - Fork 0
/
fdom.c
4851 lines (4547 loc) · 217 KB
/
fdom.c
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
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/*TO DO
0. Look at https://math.dartmouth.edu/~jvoight/articles/belyi-triangle-032717.pdf, and figure out how to compute the fundamental domain for any sub-order from the one for the maximal order. This should be much much faster than the current method! Maybe also include the tiling of the original fundamental domain.
1. algorderdisc can be very slow in some cases. Maybe randomize the choice of i1 -> i4?
2. Do we want the debug level here to be the same as for algebras? Currently it is.
3. In afuch_moreprec_shallow, when alg_hilbert is updated to allow for denominators, this method can be simplified.
4. my_alg_changeorder may have a better successor with the updated quaternion algebra methods.
5. Don't check if in normalizer for primes dividing the discriminant or for unit norms.
6. When initialize by a, b allows for denominators, fix afuchlist to not use this.
7. Geodesics that intersect a vertex. Maybe just forbid this, raising an error?
8. Find elements of norm!=1 faster
POSSIBLE FUTURE ADDITIONS:
1. Parallelization of element enumeration along with partial domain computations.
2. Methods of Imbert (see Voight's original paper) to find a minimal presentation in canonical form.
3. Computation of cohomology groups.
*/
/*INCLUSIONS*/
#include <pari.h>
#include <paripriv.h>
#include "fdom.h"
/*AUREL: Possibly make it's own debugging level? Not sure what is preferred.*/
#define DEBUGLEVEL DEBUGLEVEL_alg
/*STATIC DECLARATIONS*/
/*SECTION 1: GEOMETRIC METHODS*/
/*1: LINES AND ARCS*/
static int angle_onarc(GEN a1, GEN a2, GEN a, GEN tol);
static GEN line_from_crrcrr(GEN z1, GEN z2, GEN tol);
static GEN line_int(GEN l1, GEN l2, GEN tol);
static GEN line_int11(GEN l1, GEN l2, GEN tol);
static GEN line_line_detNULL(GEN l1, GEN l2, GEN tol);
/*1: MATRIX ACTION ON GEOMETRY*/
static GEN defp(long prec);
static GEN gdat_initialize(GEN p, long prec);
static GEN klein_safe(GEN M, long prec);
static GEN uhp_safe(GEN p, long prec);
/*1: TRANSFER BETWEEN MODELS*/
static GEN m2r_to_klein(GEN M, GEN p);
/*1: DISTANCES/AREA*/
static GEN hdiscradius(GEN area, long prec);
static GEN hdiscrandom_arc(GEN R, GEN ang1, GEN ang2, long prec);
/*1: OPERATIONS ON COMPLEX REALS*/
static GEN gtocr(GEN z, long prec);
static GEN addcrcr(GEN z1, GEN z2);
static GEN cr_conj(GEN z);
static GEN divcrcr(GEN z1, GEN z2);
static GEN mulcrcr(GEN z1, GEN z2);
static GEN mulcrcr_conj(GEN z1, GEN z2);
static GEN mulcri(GEN z, GEN n);
static GEN mulcrr(GEN z, GEN r);
static GEN negc(GEN z);
static GEN normcr(GEN z);
static GEN rM_upper_r_mul(GEN M, GEN r);
static GEN rM_upper_add(GEN M1, GEN M2);
static GEN RgM_upper_add(GEN M1, GEN M2);
static GEN subcrcr(GEN z1, GEN z2);
static GEN subrcr(GEN r, GEN z);
/*1: TOLERANCE*/
static GEN deftol(long prec);
static GEN deflowtol(long prec);
static int tolcmp(GEN x, GEN y, GEN tol);
static int toleq(GEN x, GEN y, GEN tol);
static int toleq0(GEN x, GEN tol);
static int tolsigne(GEN x, GEN tol);
/*SECTION 2: FUNDAMENTAL DOMAIN GEOMETRY*/
/*2: ISOMETRIC CIRCLES*/
static GEN icirc_angle(GEN c1, GEN c2, long prec);
static GEN icirc_klein(GEN M, GEN tol);
static GEN icirc_elt(GEN X, GEN g, GEN (*Xtoklein)(GEN, GEN), GEN gdat);
static GEN argmod(GEN x, GEN y, GEN tol);
static GEN argmod_complex(GEN c, GEN tol);
/*2: NORMALIZED BOUNDARY*/
static GEN normbound(GEN X, GEN G, GEN (*Xtoklein)(GEN, GEN), GEN gdat);
static GEN normbound_icircs(GEN C, GEN indtransfer, GEN gdat);
static long normbound_icircs_bigr(GEN C, GEN order);
static void normbound_icircs_insinfinite(GEN elts, GEN vcors, GEN vargs, GEN infinite, GEN curcirc, long *found);
static void normbound_icircs_insclean(GEN elts, GEN vcors, GEN vargs, GEN curcirc, long toins, long *found);
static void normbound_icircs_phase2(GEN elts, GEN vcors, GEN vargs, GEN curcirc, GEN firstcirc, GEN tol, long toins, long *found);
static GEN normbound_area(GEN C, long prec);
static long normbound_outside(GEN U, GEN z, GEN tol);
static int normbound_whichside(GEN side, GEN z, GEN tol);
/*2: NORMALIZED BOUNDARY APPENDING*/
static GEN normbound_append(GEN X, GEN U, GEN G, GEN (*Xtoklein)(GEN, GEN), GEN gdat);
static GEN normbound_append_icircs(GEN Uvcors, GEN Uvargs, GEN C, GEN Ctype, long rbigind, GEN gdat);
/*2: NORMALIZED BOUNDARY ANGLES*/
static int cmp_icircangle(void *tol, GEN c1, GEN c2);
static long args_find_cross(GEN args);
static long args_search(GEN args, long ind, GEN arg, GEN tol);
/*2: NORMALIZED BASIS*/
static GEN edgepairing(GEN U, GEN tol);
static GEN normbasis(GEN X, GEN U, GEN G, GEN (*Xtoklein)(GEN, GEN), GEN (*Xmul)(GEN, GEN, GEN), GEN (*Xinv)(GEN, GEN), int (*Xistriv)(GEN, GEN), GEN gdat);
/*2: NORMALIZED BOUNDARY REDUCTION*/
static GEN red_elt_decomp(GEN X, GEN U, GEN g, GEN z, GEN (*Xtoklein)(GEN, GEN), GEN (*Xmul)(GEN, GEN, GEN), GEN gdat);
static GEN red_elt(GEN X, GEN U, GEN g, GEN z, GEN (*Xtoklein)(GEN, GEN), GEN (*Xmul)(GEN, GEN, GEN), int flag, GEN gdat);
/*2: CYCLES AND SIGNATURE*/
static GEN minimalcycles(GEN pair);
static GEN minimalcycles_bytype(GEN X, GEN U, GEN Xid, GEN (*Xmul)(GEN, GEN, GEN), int (*Xisparabolic)(GEN, GEN), int (*Xistriv)(GEN, GEN));
static GEN signature(GEN X, GEN U, GEN mcyc_bytype);
/*2: GEODESICS*/
static long fdom_intersect_sidesmidpt(long i1, long i2, long n);
static GEN fdom_intersect(GEN U, GEN geod, GEN tol, long s1);
static GEN geodesic_klein(GEN X, GEN g, GEN (*Xtoklein)(GEN, GEN), GEN tol);
static GEN geodesic_fdom(GEN X, GEN U, GEN g, GEN Xid, GEN (*Xtoklein)(GEN, GEN), GEN (*Xmul)(GEN, GEN, GEN), GEN (*Xinv)(GEN, GEN), GEN gdat);
/*2: ELLIPTIC ELEMENTS AND PRESENTATION*/
static GEN elliptic(GEN X, GEN U, GEN mcyc_bytype);
static GEN presentation(GEN X, GEN U, GEN mcyc_bytype);
static GEN word_collapse(GEN word);
static GEN word(GEN X, GEN U, GEN P, GEN g, GEN (*Xtoklein)(GEN, GEN), GEN (*Xmul)(GEN, GEN, GEN), GEN (*Xinv)(GEN, GEN), int (*Xistriv)(GEN, GEN), GEN gdat);
/*SECTION 3: QUATERNION ALGEBRA METHODS*/
/*3: INITIALIZE SYMMETRIC SPACE*/
static GEN afuch_moreprec_shallow(GEN X, long inc);
static GEN afuch_make_kleinmats(GEN A, GEN O, GEN p, long prec);
static GEN afuch_make_m2rmats(GEN A, GEN O, long prec);
static GEN afuch_make_qfmats(GEN kleinmats);
static GEN afuch_make_traceqf(GEN X, GEN nm, GEN Onorm);
static GEN Omultable(GEN A, GEN O, GEN Oinv);
static GEN Onorm_makechol(GEN F, GEN Onorm);
static GEN Onorm_makemat(GEN A, GEN O, GEN AOconj);
static GEN Onorm_toreal(GEN A, GEN Onorm);
/*3: ALGEBRA FUNDAMENTAL DOMAIN CONSTANTS*/
static GEN afuchO1area(GEN A, GEN O, GEN Olevel_fact, long computeprec, long prec);
static GEN afuchbestC(GEN A, GEN O, GEN Olevel_nofact, long prec);
static GEN afuchfdomdat_init(GEN A, GEN O, long prec);
/*3: ALGEBRA FUNDAMENTAL DOMAIN METHODS*/
static GEN afuchfdom_worker(GEN X, GEN *startingset);
static GEN afuchfdom_subgroup(GEN X, GEN M);
static int nextsub(GEN S, long n);
/*3: NON NORM 1 METHODS*/
static GEN afuch_makeunitelts(GEN X);
static GEN afuch_makeALelts(GEN X);
static GEN afuch_makenormelts(GEN X);
static GEN AL_make_norms(GEN B, long split, GEN ideals, long prec);
/*3: ALGEBRA BASIC AUXILLARY METHODS*/
static GEN afuchid(GEN X);
static int afuchinnormalizer(GEN X, GEN g);
static int afuchisparabolic(GEN X, GEN g);
static GEN afuchnorm_fast(GEN X, GEN g);
static GEN afuchnorm_chol(GEN F, GEN chol, GEN g);
static GEN afuchnorm_mat(GEN F, GEN Onorm, GEN g);
static GEN afuchnorm_real(GEN X, GEN g);
static GEN afuchtoklein(GEN X, GEN g);
static GEN afuchtrace(GEN X, GEN g);
/*3: FINDING ELEMENTS*/
static GEN afuch_make_qf(GEN X, GEN nm, GEN z, GEN tracepart, GEN realnm);
static GEN afuchfindelts(GEN X, GEN nm, GEN z, GEN C, long maxelts, GEN tracepart, GEN realnm);
static GEN afuchfindoneelt_i(GEN X, GEN nm, GEN C);
/*3: ALGEBRA HELPER METHODS*/
static GEN algconj(GEN A, GEN x);
static GEN alggeta(GEN A);
static GEN voidalgmul(void *A, GEN x, GEN y);
/*3: ALGEBRA ORDER METHODS*/
static GEN algd(GEN A, GEN a);
static GEN algorderdisc(GEN A, GEN O, int reduced, int factored);
/*SECTION 4: FINCKE POHST FOR FLAG=2 WITH PRUNING*/
/*4: SUPPORTING METHODS TO FINCKE POHST*/
static GEN clonefill(GEN S, long s, long t);
static int mpgreaterthan(GEN x, GEN y);
static GEN norm_aux(GEN xk, GEN yk, GEN zk, GEN vk);
/*4: MAIN FINCKE POHST METHODS*/
static GEN smallvectors_prune(GEN q, GEN C, GEN prune);
/*SECTION 5: METHODS DELETED / NAME CHANGED FROM LIBPARI / DID NOT EXIST IN 2.15*/
static GEN qf_RgM_apply_old(GEN q, GEN M);
static void init_qf_apply(GEN q, GEN M, long *l);
static GEN my_alg_changeorder(GEN al, GEN ord);
static GEN elementabsmultable(GEN mt, GEN x);
static GEN elementabsmultable_Fp(GEN mt, GEN x, GEN p);
static GEN algbasismultable(GEN al, GEN x);
static GEN algtracebasis(GEN al);
static GEN elementabsmultable_Z(GEN mt, GEN x);
static GEN FpM_trace(GEN x, GEN p);
static GEN ZM_trace(GEN x);
/*MAIN BODY*/
/*SECTION 1: GEOMETRIC METHODS*/
/*LINES, SEGMENTS, TOLERANCE, POINTS
Line [a, b, c]
Representing ax+by=c. We will normalize so that c=gen_1 or gen_0. It is assumed that at least one of a, b is non-zero. We also assume that a and b are of type t_REAL.
Point
Stored as t_COMPLEX with t_REAL entries.
Segment [a, b, c, x0, x1]
[a, b, c] gives the line, which has start point x0 and ends at x1, which are complex. We do not allow segments going through oo. We also require x0 and x1 to have type t_COMPLEX with t_REAL components.
tol
The tolerance, which MUST be of type t_REAL. The default choice is tol=deftol(prec). Two points are declared as equal if they are equal up to tolerance.
*/
/* GEOMETRIC DATA
We will need to deal with passing from the upper half plane model -> unit ball model -> Klein model, as well as doing computations with tolerance. For this, we will fix a "geometric data" argument:
gdat = [tol, p]
tol
The tolerance, which should be initially set with tol=deftol(prec).
p
The point in the upper half plane that is mapped to 0 in the unit ball model. This should be chosen to have trivial stabilizer in Gamma, otherwise issues may arise. We convert it to have components that are t_REAL.
precision
Can be retrieved by realprec(tol), so we don't store it.
*/
/*1: LINES AND ARCS*/
/*A is the arc on the unit disc from angle a1 to a2, where a1 and a2 are in [0, 2Pi). For another angle a in [0, 2Pi), this returns 1 if a==a1, 2 if a==a2, 3 if a in in the interior of the counterclockwise arc from a1 to a2, and 0 if it is outside of this arc. All angles must be t_REAL, even if they are 0.*/
static int
angle_onarc(GEN a1, GEN a2, GEN a, GEN tol)
{
int a1cmp = tolcmp(a1, a, tol);
if (!a1cmp) return 1;/*Equal to a1*/
int cross = cmprr(a1, a2);/*No need for tolerance here*/
if (cross > 0){/*The arc crosses the x-axis.*/
if (a1cmp < 0) return 3;/*Between a1 and the x-axis*/
int a2cmp = tolcmp(a, a2, tol);
if (!a2cmp) return 2;/*Equal to a2*/
if (a2cmp < 0) return 3;/*Inside*/
return 0;/*Outside*/
}
if (a1cmp > 0) return 0;/*Between x-axis and a1*/
int a2cmp = tolcmp(a, a2, tol);
if (!a2cmp) return 2;/*Equal to a2*/
if (a2cmp > 0) return 0;/*Outside*/
return 3;/*Inside*/
}
/*Given two complex numbers with real components, this returns the line [a, b, c] through them (signifying ax+by=c), where c=0 or 1.*/
static GEN
line_from_crrcrr(GEN z1, GEN z2, GEN tol)
{
pari_sp av = avma;
GEN det = line_line_detNULL(z1, z2, tol);/*Still works for our inputs!*/
if (!det) {/*det 0, hence 0*/
return gerepilecopy(av, mkvec3(gel(z1, 2), negr(gel(z1, 1)), gen_0));
}
GEN scale = invr(det);
GEN a = mulrr(subrr(gel(z2, 2), gel(z1, 2)), scale);
GEN b = mulrr(subrr(gel(z1, 1), gel(z2, 1)), scale);
return gerepilecopy(av, mkvec3(a, b, gen_1));
}
/*Returns the intersection of two lines. If parallel/concurrent, returns NULL. We will ensure that all output numbers are COMPLEX with REAL components. This means that 0 is stored as a real number with precision.*/
static GEN
line_int(GEN l1, GEN l2, GEN tol)
{
pari_sp av = avma;
GEN c1 = gel(l1, 3), c2 = gel(l2, 3);/*Get c values*/
GEN det = line_line_detNULL(l1, l2, tol);/*ad-bc*/
if (!det) return gc_NULL(av);/*Lines are parallel*/
if (!signe(c1)) {/*ax+by = 0*/
if (!signe(c2)) {/*Both must pass through 0*/
GEN zeror = real_0(realprec(tol));/*prec=realprec(tol)*/
return gerepilecopy(av, mkcomplex(zeror, zeror));
}
return gerepilecopy(av, mkcomplex(divrr(negr(gel(l1, 2)), det), divrr(gel(l1, 1), det)));/*-b/det, a/det*/
}
/*Next, cx+dy = 0*/
if (!signe(c2)) return gerepilecopy(av, mkcomplex(divrr(gel(l2, 2), det), divrr(negr(gel(l2, 1)), det)));/*d/det, -c/det*/
/*Now ax+by = cx+dy = 1*/
GEN x = divrr(subrr(gel(l2, 2), gel(l1, 2)), det);/*(d-b)/det*/
GEN y = divrr(subrr(gel(l1, 1), gel(l2, 1)), det);/*(a-c)/det*/
return gerepilecopy(av, mkcomplex(x, y));
}
/*Line intersection where c1=c2=1, the typical case. l1 and l2 can have length 2 in this case.*/
static GEN
line_int11(GEN l1, GEN l2, GEN tol)
{
pari_sp av = avma;
GEN det = line_line_detNULL(l1, l2, tol);/*ad-bc*/
if(!det) return gc_NULL(av);/*Lines are concurrent*/
GEN x = divrr(subrr(gel(l2, 2), gel(l1, 2)), det);/*(d-b)/det*/
GEN y = divrr(subrr(gel(l1, 1), gel(l2, 1)), det);/*(a-c)/det*/
return gerepilecopy(av, mkcomplex(x, y));
}
/*Given two lines given by a1x+b1y=c1 and a2x+b2y=c2, returns a1b2-a2b1, unless this is within tolerance of 0, when we return NULL. Gerepileupto safe, leaves garbage.*/
static GEN
line_line_detNULL(GEN l1, GEN l2, GEN tol)
{
GEN d = subrr(mulrr(gel(l1, 1), gel(l2, 2)), mulrr(gel(l1, 2), gel(l2, 1)));/*ad-bc*/
if (toleq0(d, tol)) return NULL;/*d=0 up to tolerance*/
return d;
}
/*1: MATRIX ACTION ON GEOMETRY*/
/*EQUATIONS FOR ACTION
UPPER HALF PLANE
M=[a, b;c, d] in (P)GL(2, R)^+ acts on z via (az+b)/(cz+d).
UNIT DISC
M=[A, B;conj(B), conj(A)] in (P)SU(1, 1) acts on z in the same way as PGL.
KLEIN
M=[A, B] corresponding to the same (A, B) as for the unit disc action. The corresponding equation on a point in the Klein model is via
Mz=(A^2z+B^2conj(z)+2AB)/(|A|^2+|B|^2+A*z*conj(B)+conj(A)*conj(z)*B
=(A(Az+B)+B(B*conj(z)+A))/(conj(B)(Az+B)+conj(A)(B*conj(z)+A)).
*/
/*Returns the default value of p, which is Pi/8+0.5*I*/
static GEN
defp(long prec)
{
GEN p = cgetg(3, t_COMPLEX);
gel(p, 1) = Pi2n(-3, prec);
gel(p, 2) = real2n(-1, prec);
return p;
}
/*This gives the action in the unit disc model, as described above. Safe version, usable by a gp user.*/
GEN
disc_act(GEN M, GEN z, long prec)
{
pari_sp av = avma;
GEN Msafe = klein_safe(M, prec);
GEN zsafe = gtocr(z, prec);
GEN A = gel(Msafe, 1), B = gel(Msafe, 2);
GEN AzpB = addcrcr(mulcrcr(A, zsafe), B);/*Az+B*/
GEN BczpAc = addcrcr(mulcrcr_conj(z, B), cr_conj(A));/*conj(B)*z+conj(A)*/
if (gequal0(BczpAc)) pari_err(e_PREC,"Catastrophic precision loss, please recompute to a higher precision level.");
return gerepilecopy(av, divcrcr(AzpB, BczpAc));
}
/*Initializes gdat for a given p and precision. */
static GEN
gdat_initialize(GEN p, long prec)
{
GEN ret = cgetg(3, t_VEC);
gel(ret, 1) = deftol(prec);/*Default tolerance*/
gel(ret, 2) = uhp_safe(p, prec);/*Convert p to have real components, and check that it was a valid input.*/
return ret;
}
/*This gives the action in the Klein model, as described above. Returns NULL if massive precision loss.*/
GEN
klein_act_i(GEN M, GEN z)
{
pari_sp av = avma;
GEN A = gel(M, 1), B = gel(M, 2);
GEN AzpB = addcrcr(mulcrcr(A, z), B);/*Az+B*/
GEN BzcpA = addcrcr(mulcrcr_conj(B, z), A);/*B*conj(z)+A*/
GEN num = addcrcr(mulcrcr(A, AzpB), mulcrcr(B, BzcpA));/*A(Az+B)+B(B*conj(z)+A)*/
GEN denom = addcrcr(mulcrcr_conj(AzpB, B), mulcrcr_conj(BzcpA, A));/*(Az+B)conj(B)+(B*conj(z)+A)conj(A)*/
if (gequal0(denom)) return gc_NULL(av);
return gerepilecopy(av, divcrcr(num, denom));
}
/*The safe version of klein_act_i, where we ensure z and M have the correct format.*/
GEN
klein_act(GEN M, GEN z, long prec)
{
pari_sp av = avma;
GEN Msafe = klein_safe(M, prec);
GEN zsafe = gtocr(z, prec);
GEN zimg = klein_act_i(Msafe, zsafe);
if (!zimg) pari_err(e_PREC,"Catastrophic precision loss, please recompute to a higher precision level.");
return gerepileupto(av, zimg);
}
/*Returns M=[A, B] converted to having complex entries with real components of precision prec.*/
static GEN
klein_safe(GEN M, long prec)
{
if (typ(M) != t_VEC || lg(M) != 3) pari_err_TYPE("M must be a length 2 vector with complex entries.", M);
GEN Msafe = cgetg(3, t_VEC);
gel(Msafe, 1) = gtocr(gel(M, 1), prec);
gel(Msafe, 2) = gtocr(gel(M, 2), prec);
return Msafe;
}
/*Gives the action of a matrix in the upper half plane/unit disc model. We assume that the input/output are not infinity, which could happen with the upper half plane model.*/
GEN
pgl_act(GEN M, GEN z)
{
pari_sp av = avma;
GEN numer = gadd(gmul(gcoeff(M, 1, 1), z), gcoeff(M, 1, 2));
GEN denom = gadd(gmul(gcoeff(M, 2, 1), z), gcoeff(M, 2, 2));
return gerepileupto(av, gdiv(numer, denom));
}
/*p should be a point in the upper half plane, stored as a t_COMPLEX with t_REAL components. This converts p to this format, or raises an error if it is not convertible.*/
static GEN
uhp_safe(GEN p, long prec)
{
if (typ(p) != t_COMPLEX) pari_err_TYPE("The point p must be a complex number with positive imaginary part.", p);
GEN psafe = gtocr(p, prec);
if (signe(gel(psafe, 2)) != 1) pari_err_TYPE("The point p must have positive imaginary part.", p);
return psafe;
}
/*1: TRANSFER BETWEEN MODELS*/
/*Given a point z in the unit disc model, this transfers it to the Klein model.*/
GEN
disc_to_klein(GEN z)
{
pari_sp av = avma;
GEN znorm = normcr(z);
GEN scale = divsr(2, addsr(1, znorm));/*2/(1+|z|^2)*/
return gerepileupto(av, mulcrr(z, scale));/*2z/(1+|z|^2)*/
}
/*Given a point z in the unit disc model, this transfers it to the Klein model. This function is built to be used by the user (safety precautions taken).*/
GEN
disc_to_klein_user(GEN z, long prec)
{
pari_sp av = avma;
GEN zsafe = gtocr(z, prec);
return gerepileupto(av, disc_to_klein(zsafe));
}
/*Given a point z in the unit disc model, this transfers it to the upper half plane model. The formula is conj(p)*z-p/(z-1)*/
GEN
disc_to_plane(GEN z, GEN p)
{
pari_sp av = avma;
GEN num = subcrcr(mulcrcr_conj(z, p), p);/*z*conj(p)-p*/
GEN denom = mkcomplex(subrs(gel(z, 1), 1), gel(z, 2));/*z-1*/
return gerepileupto(av, divcrcr(num, denom));
}
/*Given a point z in the Klein model, this transfers it to the unit disc model. The formula is z/(1+sqrt(1-|z|^2)).*/
GEN
klein_to_disc(GEN z, GEN tol)
{
pari_sp av = avma;
GEN znm1 = subsr(1, normcr(z));/*1-|z|^2*/
if (toleq0(znm1, tol)) return gerepilecopy(av, z);/*z->z, so just copy and return it to avoid precision loss with sqrt(0).*/
GEN rt = sqrtr(znm1);/*sqrt(1-|z|^2)*/
GEN scale = invr(addsr(1, rt));/*1/(1+sqrt(1-|z|^2))*/
return gerepileupto(av, mulcrr(z, scale));/*z/(1+sqrt(1-|z|^2))*/
}
/*Given a point z in the Klein model, this transfers it to the unit disc model. This function takes the tolerance from the precision and is built to be used by the user (safety precautions taken).*/
GEN
klein_to_disc_user(GEN z, long prec)
{
pari_sp av = avma;
GEN tol = deftol(prec);
GEN zsafe = gtocr(z, prec);
return gerepileupto(av, klein_to_disc(zsafe, tol));
}
/*Given a point z in the Klein model, this transfers it to the upper half plane model.*/
GEN
klein_to_plane(GEN z, GEN p, GEN tol)
{
pari_sp av = avma;
GEN zdisc = klein_to_disc(z, tol);
return gerepileupto(av, disc_to_plane(zdisc, p));/*Klein -> disc -> plane*/
}
/*Given a point z in the upper half plane model, this transfers it to the unit disc model. The formula is (z-p)/(z-conj(p))*/
GEN
plane_to_disc(GEN z, GEN p)
{
pari_sp av = avma;
GEN num = subcrcr(z, p);/*z-p*/
GEN denom = subcrcr(z, conj_i(p));/*z-conj(p)*/
return gerepileupto(av, divcrcr(num, denom));
}
/*Given a point z in the upper half plane model, this transfers it to the Klein model.*/
GEN
plane_to_klein(GEN z, GEN p)
{
pari_sp av = avma;
GEN zdisc = plane_to_disc(z, p);
return gerepileupto(av, disc_to_klein(zdisc));/*Plane -> disc -> Klein*/
}
/*Coverts M in M(2, R) to [A, B] which corresponds to upper half plane model -action > unit disc/Klein model action when the matrices have positive determinat (det([A, B])=|A|^2-|B|^2). If M1=1/(p-conj(p))[1,-p;1,-conj(p)] and M2=[conj(p), -p;1, -1], then this is via M1*M*M2=[A, B;conj(B), conj(A)]. If M=[a, b;c, d], the explicit formula are A=(a*conj(p)-|p|^2*c+b-pd)/(p-conj(p)), and B=(-ap+p^2c-b+pd)/(p-conj(p)), with p=x+iy -> 1/(p-conj(p))=-i/2y. We ensure that A and B are of type t_COMPLEX with t_REAL coefficients.*/
static GEN
m2r_to_klein(GEN M, GEN p)
{
pari_sp av = avma;
GEN scale = invr(shiftr(gel(p, 2), 1));/*1/2y; -1 acts trivially*/
GEN ampc = subrcr(gcoeff(M, 1, 1), mulcrr(p, gcoeff(M, 2, 1)));/*a-pc*/
GEN bmpd = subrcr(gcoeff(M, 1, 2), mulcrr(p, gcoeff(M, 2, 2)));/*b-pd*/
GEN ampc_pconj = mulcrcr_conj(ampc, p);/*(a-pc)*conj(p)*/
GEN Apre = mulcrr(addcrcr(ampc_pconj, bmpd), scale);/*((a-pc)*conj(p)+b-pd)/2y*/
GEN ampc_p = mulcrcr(ampc, p);/*(a-pc)*p*/
GEN Bpre = mulcrr(negc(addcrcr(ampc_p, bmpd)), scale);/*((-a+pc)p-b+pd)/2y*/
GEN AB = mkvec2(mkcomplex(negr(gel(Apre, 2)), gel(Apre, 1)), mkcomplex(negr(gel(Bpre, 2)), gel(Bpre, 1)));/*times i*/
return gerepilecopy(av, AB);
}
/*1: DISTANCES/AREAS*/
/*Given the area (t_REAL) of a hyperbolic disc, this returns the radius. The formula is area=4*Pi*sinh(R/2)^2, or R=2arcsinh(sqrt(area/4Pi))*/
static GEN
hdiscradius(GEN area, long prec)
{
pari_sp av = avma;
GEN aover4pi = divrr(area, Pi2n(2, prec));/*area/4Pi*/
GEN rt = sqrtr(aover4pi);
GEN half = gasinh(rt, prec);
return gerepileupto(av, shiftr(half, 1));
}
/*Returns a random point z in the Klein model, uniform inside the ball of radius R. See page 19 of Page (before section 2.5).*/
GEN
hdiscrandom(GEN R, long prec)
{
pari_sp av = avma;
GEN zbound = expIPiR(shiftr(randomr(prec), 1), prec);/*Random boundary point. Now we need to scale by a random hyperbolic distance in [0, R]*/
GEN x = mulrr(gsinh(gshift(R, -1), prec), sqrtr(randomr(prec)));
/*We need to return zbound*tanh(2asinh(x))=zbound*2x*sqrt(x^2+1)/(2*x^2+1).*/
GEN xsqr = sqrr(x);
GEN num = shiftr(mulrr(x, sqrtr(addrs(xsqr, 1))), 1);/*2*x*sqrt(x^2+1)*/
GEN denom = addrs(shiftr(xsqr, 1), 1);/*2x^2+1*/
return gerepileupto(av, gmul(zbound, divrr(num, denom)));
}
/*Returns a random point z in the unit disc, uniform inside the ball of radius R, with argument uniform in [ang1, ang2]. See page 19 of Page (before section 2.5).*/
static GEN
hdiscrandom_arc(GEN R, GEN ang1, GEN ang2, long prec)
{
pari_sp av = avma;
GEN arg = gadd(ang1, gmul(randomr(prec), gsub(ang2, ang1)));/*Random angle in [ang1, ang2]*/
GEN zbound = expIr(arg);/*The boundary point. Now we need to scale by a random hyperbolic distance in [0, R]*/
GEN x = mulrr(gsinh(gshift(R, -1), prec), sqrtr(randomr(prec)));
/*We need to return tanh(2asinh(x))=2x*sqrt(x^2+1)/(2*x^2+1).*/
GEN xsqr = sqrr(x);
GEN num = shiftr(mulrr(x, sqrtr(addrs(xsqr, 1))), 1);/*2*x*sqrt(x^2+1)*/
GEN denom = addrs(shiftr(xsqr, 1), 1);/*2x^2+1*/
return gerepileupto(av, gmul(zbound, divrr(num, denom)));
}
/*1: OPERATIONS ON COMPLEX REALS*/
/*z should be a complex number with real components of type prec. This converts it to one. Clean method.*/
static GEN
gtocr(GEN z, long prec)
{
GEN zsafe = cgetg(3, t_COMPLEX);
if (typ(z) == t_COMPLEX) {
gel(zsafe, 1) = gtofp(gel(z, 1), prec);
gel(zsafe, 2) = gtofp(gel(z, 2), prec);
return zsafe;
}
gel(zsafe, 1) = gtofp(z, prec);
gel(zsafe, 2) = real_0(prec);
return zsafe;
}
/*Adds two complex numbers with real components, giving a complex output. Clean method.*/
static GEN
addcrcr(GEN z1, GEN z2)
{
GEN z = cgetg(3, t_COMPLEX);
gel(z, 1) = addrr(gel(z1, 1), gel(z2, 1));
gel(z, 2) = addrr(gel(z1, 2), gel(z2, 2));
return z;
}
/*conjugate of a complex number with two real components. Not gerepileupto safe, shallow method.*/
static GEN
cr_conj(GEN z)
{
GEN znew = cgetg(3, t_COMPLEX);
gel(znew, 1) = gel(z, 1);
gel(znew, 2) = negr(gel(z, 2));
return znew;
}
/*Divides two complex numbers with real components, giving a complex output. Gerepileupto safe, leaves garbage.*/
static GEN
divcrcr(GEN z1, GEN z2)
{
GEN num = mulcrcr_conj(z1, z2);/*z1/z2=(z1*conj(z2))/norm(z2)*/
GEN den = normcr(z2);
GEN z = cgetg(3, t_COMPLEX);
gel(z, 1) = divrr(gel(num, 1), den);
gel(z, 2) = divrr(gel(num, 2), den);
return z;
}
/*Multiplies two complex numbers with real components, giving a complex output. Gerepileupto safe, leaves garbage.*/
static GEN
mulcrcr(GEN z1, GEN z2)
{
GEN ac = mulrr(gel(z1, 1), gel(z2, 1));
GEN ad = mulrr(gel(z1, 1), gel(z2, 2));
GEN bc = mulrr(gel(z1, 2), gel(z2, 1));
GEN bd = mulrr(gel(z1, 2), gel(z2, 2));
GEN z = cgetg(3, t_COMPLEX);
gel(z, 1) = subrr(ac, bd);
gel(z, 2) = addrr(ad, bc);
return z;
}
/*mulcrcr, except we do z1*conj(z2). Gerepileupto safe, leaves garbage.*/
static GEN
mulcrcr_conj(GEN z1, GEN z2)
{
GEN ac = mulrr(gel(z1, 1), gel(z2, 1));
GEN ad = mulrr(gel(z1, 1), gel(z2, 2));
GEN bc = mulrr(gel(z1, 2), gel(z2, 1));
GEN bd = mulrr(gel(z1, 2), gel(z2, 2));
GEN z = cgetg(3, t_COMPLEX);
gel(z, 1) = addrr(ac, bd);
gel(z, 2) = subrr(bc, ad);
return z;
}
/*Multiplies a complex number with real components by an integer. Clean method.*/
static GEN
mulcri(GEN z, GEN n)
{
GEN zn = cgetg(3, t_COMPLEX);
gel(zn, 1) = mulri(gel(z, 1), n);
gel(zn, 2) = mulri(gel(z, 2), n);
return zn;
}
/*Multiplies a complex number with real components by a real number. Clean method.*/
static GEN
mulcrr(GEN z, GEN r)
{
GEN zp = cgetg(3, t_COMPLEX);
gel(zp, 1) = mulrr(gel(z, 1), r);
gel(zp, 2) = mulrr(gel(z, 2), r);
return zp;
}
/*Returns -z, for a complex number with real components. Clean method.*/
static GEN
negc(GEN z)
{
GEN zp = cgetg(3, t_COMPLEX);
gel(zp, 1) = negr(gel(z, 1));
gel(zp, 2) = negr(gel(z, 2));
return zp;
}
/*Norm of a complex number with real components, giving a real output. Gerepileupto safe, leaves garbage.*/
static GEN
normcr(GEN z)
{
GEN x = sqrr(gel(z, 1));
GEN y = sqrr(gel(z, 2));
return addrr(x, y);
}
/*Multiplies a square upper triangular matrix M with upper triangle having real entries by a real number r. Clean method.*/
static GEN
rM_upper_r_mul(GEN M, GEN r)
{
long n = lg(M), i, j;
GEN P = cgetg(n, t_MAT);
for (i = 1; i < n; i++) {
gel(P, i) = cgetg(n, t_COL);
for (j = 1; j <= i; j++) gmael(P, i, j) = mulrr(gcoeff(M, j, i), r);
for (j = i + 1; j < n; j++) gmael(P, i, j) = gen_0;
}
return P;
}
/*Adds two square upper triangular matrices M1, M2 with upper triangle having real entries. Clean method.*/
static GEN
rM_upper_add(GEN M1, GEN M2)
{
long n = lg(M1), i, j;
GEN P = cgetg(n, t_MAT);
for (i = 1; i < n; i++) {
gel(P, i) = cgetg(n, t_COL);
for (j = 1; j <= i; j++) gmael(P, i, j) = addrr(gcoeff(M1, j, i), gcoeff(M2, j, i));
for (j = i + 1; j < n; j++) gmael(P, i, j) = gen_0;
}
return P;
}
/*Adds two square upper triangular matrices. Clean method.*/
static GEN
RgM_upper_add(GEN M1, GEN M2)
{
long n = lg(M1), i, j;
GEN P = cgetg(n, t_MAT);
for (i = 1; i < n; i++) {
gel(P, i) = cgetg(n, t_COL);
for (j = 1; j <= i; j++) gmael(P, i, j) = gadd(gcoeff(M1, j, i), gcoeff(M2, j, i));
for (j = i + 1; j < n; j++) gmael(P, i, j) = gen_0;
}
return P;
}
/*Subtracts two complex numbers with real components, giving a complex output. Clean method.*/
static GEN
subcrcr(GEN z1, GEN z2)
{
GEN z = cgetg(3, t_COMPLEX);
gel(z, 1) = subrr(gel(z1, 1), gel(z2, 1));
gel(z, 2) = subrr(gel(z1, 2), gel(z2, 2));
return z;
}
/*Does r-z, where r is real and z is complex with real components. Clean method.*/
static GEN
subrcr(GEN r, GEN z)
{
GEN zp = cgetg(3, t_COMPLEX);
gel(zp, 1) = subrr(r, gel(z, 1));
gel(zp, 2) = negr(gel(z, 2));
return zp;
}
/*1: TOLERANCE*/
/*Returns the default tolerance given the precision, which is saying that x==y if they are equal up to half of the precision.*/
static GEN
deftol(long prec)
{
return real2n(-(prec2nbits(prec) >> 1), prec);
}
/*Lower tolerance. Useful if we may have more precision loss and are OK with the larger range (e.g. we have a slow routine to check exact equality, and use this to sift down to a smaller set).*/
static GEN
deflowtol(long prec)
{
return real2n(-16, prec);
}
/*Returns -1 if x<y, 0 if x==y, 1 if x>y (x, y are t_REAL/t_INT/t_FRAC). Accounts for the tolerance, so will deem x==y if they are equal up to tol AND at least one is inexact*/
static int
tolcmp(GEN x, GEN y, GEN tol)
{
pari_sp av = avma;
GEN d = gsub(x, y);
return gc_int(av, tolsigne(d, tol));/*Return sign(x-y)*/
}
/*Returns 1 if x==y up to tolerance tol. If x and y are both t_INT/t_FRAC, will only return 1 if they are exactly equal.*/
static int
toleq(GEN x, GEN y, GEN tol)
{
pari_sp av = avma;
GEN d = gsub(x, y);
return gc_int(av, toleq0(d, tol));/*Just compare d with 0.*/
}
/*If x is of type t_INT or t_FRAC, returns 1 iff x==0. Otherwise, x must be of type t_REAL or t_COMPLEX, and returns 1 iff x=0 up to tolerance tol.*/
static int
toleq0(GEN x, GEN tol)
{
switch (typ(x)) {
case t_REAL:
if (abscmprr(x, tol) < 0) return 1;/*|x|<tol*/
return 0;
case t_COMPLEX:;
long i;
for (i = 1; i <= 2; i++) {
switch (typ(gel(x, i))) {
case t_REAL:
if (abscmprr(gel(x, i), tol) >= 0) return 0;/*Too large*/
break;
case t_INT:
if (signe(gel(x, i))) return 0;
break;
case t_FRAC:/*Fraction component, cannot be 0*/
return 0;
default:/*Illegal input*/
pari_err_TYPE("Tolerance equality only valid for type t_INT, t_FRAC, t_REAL, t_COMPLEX", x);
}
}
return 1;/*We passed*/
case t_INT:/*Given exactly*/
return !signe(x);
case t_FRAC:/*t_FRAC cannot be 0*/
return 0;
}
pari_err_TYPE("Tolerance equality only valid for type t_INT, t_FRAC, t_REAL, t_COMPLEX", x);
return 0;/*So that there is no warning*/
}
/*Returns the sign of x, where it is 0 if is equal to 0 up to tolerance.*/
static int
tolsigne(GEN x, GEN tol)
{
switch (typ(x)) {
case t_REAL:
if (abscmprr(x, tol) < 0) return 0;/*|x|<tol*/
case t_INT:/*Given exactly or real and not 0 up to tolerance.*/
return signe(x);
case t_FRAC:/*t_FRAC cannot be 0*/
return signe(gel(x, 1));
}
pari_err_TYPE("Tolerance comparison only valid for type t_INT, t_FRAC, t_REAL", x);
return 0;/*So that there is no warning*/
}
/*SECTION 2: FUNDAMENTAL DOMAIN GEOMETRY*/
/* DEALING WITH GENERAL STRUCTURES
Let X be a structure, from which Gamma, a discrete subgroup of PSL(2, R), can be extracted. Given a vector of elements, we want to be able to compute the normalized boundary, and the normalized basis. In order to achieve this, we need to pass in various methods to deal with operations in Gamma (that should be memory clean), namely:
i) Multiply elements of Gamma: format as GEN Xmul(GEN X, GEN g1, GEN g2).
ii) Invert elements of Gamma: format as GEN Xinv(GEN X, GEN g). We don't actually need the inverse, as long as g*Xinv(X, g) embeds to a constant multiple of the identity we are fine. In particular, we can use conjugation in a quaternion algebra.
iii) Embed elements of Gamma into PGU(1, 1)^+ that act on the Klein model: format as GEN Xtoklein(GEN X, GEN g). The output should be [A, B], where the corresponding element in PGU(1, 1)^+ is [A, B;conj(B), conj(A)], and |A|^2-|B|^2 is positive. If |A|^2-|B|^2!=1, this is OK, and do not rescale it, as this will be more inefficient and may cause precision loss. If you have a natural embedding into PGL(2, R)^+, then you can use m2r_to_klein, which will ensure that det(M in PGL)=|A|^2-|B|^2. Also ensure that A and B are of type t_COMPLEX with t_REAL components, even if one of them is real or even integral.
iv) Identify if an element is trivial in Gamma: format as int Xistriv(GEN X, GEN g). Since we are working in PSL, we need to be careful that -g==g, since most representations of elements in X would be for SL. Furthermore, if we do not return the true inverse in Xinv, then we have to account for this.
v) Determine if an element of X is parabolic or not: format at int Xisparabolic(GEN X, GEN g) (typically you can compare trd(g)^2 to +/-4*nrd(g)).
vi) Input an element representing the trivial element of X.
We do all of our computations in the Klein model.
*/
/*2: ISOMETRIC CIRCLES*/
/* ISOMETRIC CIRCLE FORMATTING
icirc
[a, b, r, p1, p2, ang1, ang2]
a, b
ax+by=1 is the Klein model equation of the boundary. a and b are t_REAL
(x-a)^2+(y-b)^2=r^2 is the unit disc model equation of the boundary
r
r^2+1=a^2+b^2 as it is orthogonal to the unit disc, and is of type t_REAL
p1/p2
Intersection points with the unit disc; see below for the ordering. Of type t_COMPLEX with t_REAL components.
ang1
Assume that the angle from p1 to p2 is <pi (which uniquely determines them). Then ang1 is the argument of p1, in the range [0, 2*pi). Type t_REAL
ang2
The argument of p2, in the range [0, 2*pi). It may be less than ang1 if the arc crosses 1. Type t_REAL.
*/
/*Given two isometric circles c1 and c2, assumed to intersect, this computes the angle they form (counterclockwise from the tangent to c1 to the tangent to c2 at the intersection point).
If c1 is given by (x-a)^2+(y-b)^2=r^2 and c2 by (x-c)^2+(y-d)^2=s^2 (in the unit ball model) and theta is the angle formed, then a long computation shows that cos(theta)=(1-ac-bd)/rs.
This method is not stack clean.
*/
static GEN
icirc_angle(GEN c1, GEN c2, long prec)
{
GEN ac = mulrr(gel(c1, 1), gel(c2, 1));
GEN bd = mulrr(gel(c1, 2), gel(c2, 2));
GEN omacmbd = subsr(1, addrr(ac, bd));/*1-ac-bd*/
GEN cost = divrr(omacmbd, mulrr(gel(c1, 3), gel(c2, 3)));/*cos(theta)*/
return gacos(cost, prec);/*acos is in the interval [0, Pi], which is what we want.*/
}
/*Given M=[A, B] with |A|^2-|B|^2=w acting on the Klein model, this returns the isometric circle associated to it. This has centre -conj(A/B), radius sqrt(w)/|B|. In the Klein model, the centre being a+bi -> xa+yb=1, and this intersects the unit disc at (a+/-br/(a^2+b^2), (b-/+ar)/(a^2+b^2)). Assumptions:
If we pass in an element giving everything (i.e. B=0), we return 0.
A and B are t_COMPLEX with t_REAL components
*/
static GEN
icirc_klein(GEN M, GEN tol)
{
if (toleq0(gel(M, 2), tol)) return gen_0;/*Isometric circle is everything, don't want to call it here.*/
pari_sp av = avma;
long prec = realprec(tol);
GEN A = gel(M, 1), B = gel(M, 2);
GEN centre_pre = divcrcr(A, B);/*The centre is -conj(centre_pre)*/
GEN Bnorm = normcr(B);
GEN w = subrr(normcr(A), Bnorm);
GEN r = sqrtr(divrr(w, Bnorm));/*sqrt(w)/|B| is the radius*/
GEN a = gel(centre_pre, 1), b = gel(centre_pre, 2);/*The coords of the centre.*/
togglesign(a);/*centre_pre=x+iy, then the centre is -x+iy*/
GEN ar = mulrr(a, r), br = mulrr(b, r);
GEN apbr = addrr(a, br), ambr = subrr(a, br);/*a+/-br*/
GEN bpar = addrr(b, ar), bmar = subrr(b, ar);/*b+/-ar*/
GEN theta1 = argmod(apbr, bmar, tol);/*First point angle in [0, 2pi)*/
GEN theta2 = argmod(ambr, bpar, tol);/*Second point angle in [0, 2pi)*/
GEN asqbsq = addrs(sqrr(r), 1);/*a^2+b^2 = r^2+1*/
GEN p1 = mkcomplex(divrr(apbr, asqbsq), divrr(bmar, asqbsq));/*First intersection point.*/
GEN p2 = mkcomplex(divrr(ambr, asqbsq), divrr(bpar, asqbsq));/*Second intersection point.*/
GEN thetadiff = subrr(theta2, theta1);
if (signe(thetadiff) == 1) {/*If the next if block executes, theta2 is actually the "first" angle, so we swap them.*/
if (cmprr(thetadiff, mppi(prec)) > 0) return gerepilecopy(av, mkvecn(7, a, b, r, p2, p1, theta2, theta1));
}
else {/*Same as above*/
if (cmprr(thetadiff, negr(mppi(prec))) > 0) return gerepilecopy(av, mkvecn(7, a, b, r, p2, p1, theta2, theta1));
}
return gerepilecopy(av, mkvecn(7, a, b, r, p1, p2, theta1, theta2));/*Return the output!*/
}
/*Computes the isometric circle for g in Gamma, returning [g, M, icirc], where M gives the action of g on the Klein model.*/
static GEN
icirc_elt(GEN X, GEN g, GEN (*Xtoklein)(GEN, GEN), GEN gdat)
{
pari_sp av = avma;
GEN tol = gdat_get_tol(gdat);
GEN ret = cgetg(4, t_VEC);
gel(ret, 1) = gcopy(g);
gel(ret, 2) = Xtoklein(X, g);
gel(ret, 3) = icirc_klein(gel(ret, 2), tol);
return gerepileupto(av, ret);
}
/*Returns the argument of x+iy in the range [0, 2*pi). Assumes x and y are not both 0 and are t_REAL. Gerepileupto safe, leaves garbage.*/
static GEN
argmod(GEN x, GEN y, GEN tol)
{
long prec = realprec(tol);
int xsign = tolsigne(x, tol);/*Sign of x up to tolerance.*/
if (xsign == 0) {/*Fixing theta when x == 0, so the line is vertical.*/
GEN theta = Pi2n(-1, prec);/*Pi/2*/
if (signe(y) == -1) theta = addrr(theta, mppi(prec));/*3Pi/2*/
return theta;
}
GEN theta = gatan(gdiv(y, x), prec);
if (xsign == -1) return addrr(theta, mppi(prec));/*atan lies in (-Pi/2, Pi/2), so must add by Pi to get it correct.*/
int ysign = tolsigne(y, tol);
if (ysign == 0) return real_0(prec);/*Convert to real number.*/
if (ysign == -1) theta = addrr(theta, Pi2n(1, prec));/*Add 2Pi to get in the right interval*/
return theta;
}
/*argmod, except we take c=x+iy. Gerepileupto safe, leaves garbage.*/
static GEN
argmod_complex(GEN c, GEN tol) { return argmod(gel(c, 1), gel(c, 2), tol); }
/*2: NORMALIZED BOUNDARY*/
/*NORMALIZED BOUNDARY FORMATTING
A normalized boundary is represented by U, where
U=[elts, sides, vcors, vargs, crossind, kact, area, spair, infinite]
elts
Elements of Gamma whose isometric circles give the sides of the normalied boundary. An infinite side corresponds to the element 0. Elements cannot be stored as type t_INT.
sides
The ith entry is the isometric circle corresponding to elts[i], stored as an icirc. Infinite side -> 0. The first side is the closest to the origin.
vcors
Vertex coordinates, stored as a t_COMPLEX with real components. The side sides[i] has vertices vcor[i-1] and vcor[i], appearing in this order going counterclockwise about the origin.
vargs
The argument of the corresponding vertex, normalized to lie in [0, 2*Pi), stored as t_REAL. These are stored in counterclockwise order, BUT vargs itself is not sorted: it is increasing from 1 to crossind, and from crossind+1 to the end.
crossind
the arc from vertex crossind to crossind+1 contains the positive x-axis. If one vertex IS on this axis, then crossind+1 gives that vertex. Alternatively, crossind is the unique vertex such that vargs[crossind]>vargs[crossind+1], taken cyclically.
kact
The action of the corresponding element on the Klein model, for use in klein_act. Infinite side -> 0
area
The hyperbolic area of U, which will be oo unless we are a finite index subgroup of Gamma.
spair
Stores the side pairing of U, if it exists/has been computed. When computing the normalized boundary, this will store the indices of the sides that got deleted.
infinite
Stores the indices of the infinite sides.
*/
/*Initializes the inputs for normbound_icircs. G is the set of elements we are forming the normalized boundary for. Returns NULL if no elements giving an isometric circle are input. Not gerepileupto safe, and leaves garbage.*/
static GEN
normbound(GEN X, GEN G, GEN (*Xtoklein)(GEN, GEN), GEN gdat)
{
long lG = lg(G), i;
GEN C = vectrunc_init(lG);
GEN indtransfer = vecsmalltrunc_init(lG);/*Transfer indices of C back to G*/
for (i = 1; i < lG; i++) {
GEN circ = icirc_elt(X, gel(G, i), Xtoklein, gdat);
if(gequal0(gel(circ, 3))) continue;/*No isometric circle*/
vectrunc_append(C, circ);
vecsmalltrunc_append(indtransfer, i);
}
if(lg(C) == 1) return NULL;
return normbound_icircs(C, indtransfer, gdat);
}
/*Given C, the output of icirc_elt for each of our elements, this computes and returns the normalized boundary. Assumptions:
C[i]=[g, M, icirc], where g=elt, M=Kleinian action, and icirc=[a, b, r, p1, p2, ang1, ang2].
None of the entries should be infinite sides, sort this out in the method that calls this.
C has at least one element, so the boundary is non-trivial.
indtransfer tracks the indices of G in terms of the indices in C, for tracking the deleted elements.
Not gerepileupto safe, and leaves garbage.
*/
static GEN
normbound_icircs(GEN C, GEN indtransfer, GEN gdat)
{