-
Notifications
You must be signed in to change notification settings - Fork 1
/
qq_geometry.c
2216 lines (2096 loc) · 104 KB
/
qq_geometry.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
//Geometry methods
#include <pari.h>
#include "qquadraticdecl.h"
//The length (lg, so technically length+1) of a circle/line and arc/segment, and a normalized boundary
#define CIRCLEN 4
#define ARCLEN 9
#define NORMBOUND 9
//STATIC DECLARATIONS
static GEN mobius_arcseg(GEN M, GEN c, int isarc, GEN tol, long prec);
static GEN mobius_circle(GEN M, GEN c, GEN tol, long prec);
static GEN mobius_line(GEN M, GEN l, GEN tol, long prec);
static GEN normalizedbasis_shiftpoint(GEN c, GEN r, int initial, long prec);
//Circle is stored as [centre, radius, 0], where 0 means a bona fide cicle (and not a line)
//Circle arc is [centre, radius, start pt, end pt, start angle, end angle, dir, 0]. It is the arc counterclockwise from startpt to endpt, and dir=1 means oriented counterclockwise, and dir=-1 means oriented clockwise. This can also be left uninitialized if arc is not oriented. The final 0 represents that we have an arc and not a segment.
//Line is stored as [slope, intercept, 1], where the line is y=slope*x+intercept unless slope=oo, where it is x=intercept instead, and the 1 just means a bona fide line (not circle).
//Line segment is stored as [slope, intercept, startpt, endpt, 0, ooendptor, dir, 1] (the extra 0 is to format it the same as a circle arc). The final 1 is to signal a segment. dir=1 means we go from startpt to endpt in the upper half plane, and -1 means we go through infinity (only when neither endpoint is infinite). If one endpoint is oo, ooendptor stores which way we get to it. If ooendptor=1, this means the segment travels either vertically up or right, and -1 means the arc is vertically down or left.
//GEN tol -> The tolerance.
//BASIC LINE, CIRCLE, AND POINT OPERATIONS
//Creates the arc on circle c going from p1 to p2 counterclockwise. if dir=1, the arc is oriented counterclockwise, else it is clockwise (i.e. clockwise from p2 to p1). If dir=0, we take it to be unoriented
GEN arc_init(GEN c, GEN p1, GEN p2, int dir, long prec){
pari_sp top=avma;
GEN ang2=radialangle(c, p2, gen_0, prec);//No tolerance need
GEN arc=cgetg(ARCLEN, t_VEC);
gel(arc, 1)=gcopy(gel(c, 1));
gel(arc, 2)=gcopy(gel(c, 2));
gel(arc, 3)=gcopy(p1);
gel(arc, 4)=gcopy(p2);
gel(arc, 5)=radialangle(c, p1, gen_0, prec);//start angle
gel(arc, 6)=shiftangle(ang2, gel(arc, 5), gen_0, prec);//end angle; no need for tolerance.
if(dir==1) gel(arc, 7)=gen_1;
else if(dir==0) gel(arc, 7)=gen_0;
else gel(arc, 7)=gen_m1;
gel(arc, 8)=gen_0;
return gerepileupto(top, arc);
}
//arc_init with checking of c
GEN arc_init_tc(GEN c, GEN p1, GEN p2, int dir, long prec){
int i1=geom_check(c);
if(i1!=0 && i1!=2) pari_err_TYPE("Please input a circle", c);
return arc_init(c, p1, p2, dir, prec);
}
//Returns the midpoint of the arc between p1 and p2 (counterclockwise) on c.
GEN arc_midpoint(GEN c, GEN p1, GEN p2, GEN tol, long prec){
pari_sp top=avma;
GEN pts=circleline_int(c, perpbis(p1, p2), tol, prec);
GEN a1=radialangle(c, p1, gen_0, prec);
GEN a2=shiftangle(radialangle(c, p2, gen_0, prec), a1, gen_0, prec);//No tolerance concerns
GEN angint1=shiftangle(radialangle(c, gel(pts, 1), gen_0, prec), a1, gen_0, prec);//the angle formed by pts[1] to c with base a1. Again, no tolerance need.
if(gcmp(a2, angint1)==1) return gerepilecopy(top, gel(pts, 1));//No need for tolerance, as if this is an issue our points p1 and p2 would be equal up to tolerance.
return gerepilecopy(top, gel(pts, 2));
}
//arc_midpoint with checking c and presetting tol
GEN arc_midpoint_tc(GEN c, GEN p1, GEN p2, long prec){
int i1=geom_check(c);
if(i1!=0 && i1!=2) pari_err_TYPE("Please input a circle", c);
pari_sp top=avma;
GEN tol=deftol(prec);
return gerepileupto(top, arc_midpoint(c, p1, p2, tol, prec));
}
//Returns the angle of intersection between circles c1 and c2 which intersect at p, with the angle formed by rotating the tangent to c1 at p counterclockwise to the tangent to c2 at p.
GEN circle_angle(GEN c1, GEN c2, GEN p, GEN tol, long prec){
pari_sp top=avma;
GEN s1=circle_tangentslope(c1, p, prec);
GEN s2=circle_tangentslope(c2, p, prec);
GEN ang=anglediff(atanoo(s2, prec), atanoo(s1, prec), tol, prec), pi=mppi(prec);//Difference in angles in [0, 2*Pi]
int topi=tolcmp(ang, pi, tol, prec);//We want to be in the range [0, Pi), so potentially subtract off Pi
if(topi==-1) return gerepileupto(top, ang);
else if(topi==0){avma=top;return gen_0;}//Same angle
return gerepileupto(top, gsub(ang, pi));//Subtract off pi
}
//Circle_angle with typechecking
GEN circle_angle_tc(GEN c1, GEN c2, GEN p, long prec){
pari_sp top=avma;
int i1=geom_check(c1);
if(i1!=0 && i1!=2) pari_err_TYPE("Please input a circle", c1);
int i2=geom_check(c2);
if(i2!=0 && i2!=2) pari_err_TYPE("Please input a circle", c2);
GEN tol=deftol(prec);
return gerepileupto(top, circle_angle(c1, c2, p, tol, prec));
}
//Circle with centre cent passing through p
GEN circle_fromcp(GEN cent, GEN p, long prec){
pari_sp top=avma;
GEN pmcent=gsub(p, cent);
GEN ret=cgetg(4, t_VEC);
gel(ret, 1)=gcopy(cent);
gel(ret, 2)=gabs(pmcent, prec);
gel(ret, 3)=gen_0;
return gerepileupto(top, ret);
}
//Circle through 3 points (with one allowed to be oo, making a line instead)
GEN circle_fromppp(GEN p1, GEN p2, GEN p3, GEN tol, long prec){
if(typ(p1)==t_INFINITY) return line_frompp(p2, p3);
if(typ(p2)==t_INFINITY) return line_frompp(p1, p3);
if(typ(p3)==t_INFINITY) return line_frompp(p1, p2);//Lines
pari_sp top=avma;
GEN l1=perpbis(p1, p2), l2=perpbis(p1, p3);
GEN centre=line_int(l1, l2, tol, prec);//centre is intersection of perp bisectors.
if(typ(centre)==t_INFINITY) return gerepileupto(top, line_frompp(p1, p2));//p1, p2, p3 collinear
return gerepileupto(top, circle_fromcp(centre, p1, prec));//The circle!
}
//circle_fromppp with presetting of tolerance.
GEN circle_fromppp_tc(GEN p1, GEN p2, GEN p3, long prec){
pari_sp top=avma;
GEN tol=deftol(prec);
return gerepileupto(top, circle_fromppp(p1, p2, p3, tol, prec));
}
//Returns the slope of the tangent to c at p
GEN circle_tangentslope(GEN c, GEN p, long prec){
pari_sp top=avma;
GEN c1mp=gsub(gel(c, 1), p);//c[1]-p
GEN c1mpr=real_i(c1mp);
GEN c1mpi=imag_i(c1mp);
return gerepileupto(top, divoo(c1mpr, gneg(c1mpi)));//divoo(real(c[1])-real(p),imag(p)-imag(c[1])));
}
//Checks c is a circle or circle arc and calls circle_tangentslope
GEN circle_tangentslope_tc(GEN c, GEN p, long prec){
int i1=geom_check(c);
if(i1!=0 && i1!=2) pari_err_TYPE("Please input a circle", c);
return circle_tangentslope(c, p, prec);
}
//Crossratio, allows one entry to be infinite.
GEN crossratio(GEN a, GEN b, GEN c, GEN d){
pari_sp top=avma;
if (typ(a)==t_INFINITY) return gerepileupto(top,divoo(gsub(d, b), gsub(c, b)));
if (typ(b)==t_INFINITY) return gerepileupto(top,divoo(gsub(c, a), gsub(d, a)));
if (typ(c)==t_INFINITY) return gerepileupto(top,divoo(gsub(d, b), gsub(d, a)));
if (typ(d)==t_INFINITY) return gerepileupto(top,divoo(gsub(c, a), gsub(c, b)));
return gerepileupto(top,divoo(gmul(gsub(c, a), gsub(d, b)), gmul(gsub(c, b), gsub(d, a))));
}
//Angle between l1, l2, formed by rotating l1 counterclockwise (result is in [0, Pi)
GEN line_angle(GEN l1, GEN l2, long prec){
pari_sp top=avma;
return gerepileupto(top, gmod(gsub(atanoo(gel(l2, 1), prec), atanoo(gel(l1, 1), prec)), mppi(prec)));
}
//The line through p with slope s
GEN line_fromsp(GEN s, GEN p){
if(typ(s)==t_INFINITY){//oo slope
GEN ret=cgetg(4, t_VEC);
gel(ret, 1)=mkoo();
gel(ret, 2)=greal(p);//x-intercept
gel(ret, 3)=gen_1;
return ret;
}
pari_sp top=avma;
GEN strealp=gmul(s, real_i(p)), imagp=imag_i(p);
GEN ret=cgetg(4, t_VEC);
gel(ret, 1)=gcopy(s);
gel(ret, 2)=gsub(imagp, strealp);//y=intercept
gel(ret, 3)=gen_1;
return gerepileupto(top, ret);
}
//Line through two points
GEN line_frompp(GEN p1, GEN p2){
pari_sp top=avma;
return gerepileupto(top, line_fromsp(slope(p1, p2), p1));
}
//Mx, where M is a 2x2 matrix and x is complex or infinite.
GEN mat_eval(GEN M, GEN x){
pari_sp top = avma;
if(typ(x)==t_INFINITY) return gerepileupto(top,divoo(gcoeff(M, 1, 1), gcoeff(M, 2, 1)));
return gerepileupto(top,divoo(gadd(gmul(gcoeff(M, 1, 1), x), gcoeff(M, 1, 2)), gadd(gmul(gcoeff(M, 2, 1), x), gcoeff(M, 2, 2))));
}
//This checks that the inputs are valid.
GEN mat_eval_tc(GEN M, GEN x){
if(typ(M)!=t_MAT) pari_err(e_TYPE,"mat_eval. M should be a 2x2 matrix",M);
if(lg(M)!=3 || lg(gel(M,1))!=3) pari_err(e_DIM,"mat_eval. M should be a 2x2 matrix");
return mat_eval(M,x);
}
//Midpoint of p1 and p2
GEN midpoint(GEN p1, GEN p2){
pari_sp top=avma;
return gerepileupto(top, gdivgs(gadd(p1, p2), 2));
}
//Returns M(c), for c a circle/line/arc/segment
GEN mobius(GEN M, GEN c, GEN tol, long prec){
int i=geom_check(c);
switch(i){
case 0: return mobius_circle(M, c, tol, prec);
case 1: return mobius_line(M, c, tol, prec);
case 2: return mobius_arcseg(M, c, 1, tol, prec);
case 3: return mobius_arcseg(M, c, 0, tol, prec);
}
pari_err_TYPE("Please input a circle/line/arc/segment", c);
return gen_0;//Never reach this far
}
//mobius, where we check the type of M and set the default tolerance
GEN mobius_tc(GEN M, GEN c, long prec){
if(typ(M)!=t_MAT) pari_err_TYPE("M should be a 2x2 matrix", M);
if(lg(M)!=3 || lg(gel(M,1))!=3) pari_err(e_DIM,"M should be a 2x2 matrix");
pari_sp top=avma;
GEN tol=deftol(prec);
return gerepileupto(top, mobius(M, c, tol, prec));
}
//Mobius map acting on an arc or segment
static GEN mobius_arcseg(GEN M, GEN c, int isarc, GEN tol, long prec){
pari_sp top=avma;//We start by finding the new circle/line(ignoring the arc/segment bits)
GEN endpt1, endpt2, extpt;//Endpoints and an extra point (used to have 3 points to translate)
endpt1=mat_eval(M, gel(c, 3));
endpt2=mat_eval(M, gel(c, 4));
//Now we must find an extra point extpt
if(isarc==1) extpt=mat_eval(M, arc_midpoint(c, gel(c, 3), gel(c, 4), tol, prec));//arc
else{//segment
if(typ(gel(c, 3))==t_INFINITY){//Start point infinity
GEN u;
if(gequal(gel(c, 6), gen_1)) u=gen_m2;//segment goes vertically up or right
else u=gen_2;//segment goes vertically down or left
if(typ(gel(c, 1))==t_INFINITY) extpt=mat_eval(M, gadd(gel(c, 4), gmul(u, gen_I())));//Vertical line, M(c[4]+U*I)
else extpt=mat_eval(M, gadd(gel(c, 4), gmul(u, gaddsg(1, gmul(gel(c, 1), gen_I())))));//non-vertical line, M(c[4]+u+u*c[1]*I)
}
else if(typ(gel(c, 4))==t_INFINITY){//End point infinity
GEN u;
if(gequal(gel(c, 6), gen_1)) u=gen_2;//segment goes vertically up or right
else u=gen_m2;//segment goes vertically down or left
if(typ(gel(c, 1))==t_INFINITY) extpt=mat_eval(M, gadd(gel(c, 3), gmul(u, gen_I())));//Vertical line, M(c[3]+U*I)
else extpt=mat_eval(M, gadd(gel(c, 3), gmul(u, gaddsg(1, gmul(gel(c, 1), gen_I())))));//non-vertical line, M(c[3]+u+u*c[1]*I)
}
else{//Start/end points in the plane
if(gequal(gel(c, 7), gen_1)) extpt=mat_eval(M, midpoint(gel(c, 3), gel(c, 4)));//Does not go through oo, can take midpoint
else extpt=mat_eval(M, mkoo());//Use oo, since the line goes through there
}
}
//Now we have the 3 new points used to define the new arc/segment. Let's finish the process.
GEN ret=cgetg(ARCLEN, t_VEC);//The returned arc/segment
GEN newcirc=circle_fromppp(endpt1, endpt2, extpt, tol, prec);//The new circle/line
gel(ret, 1)=gel(newcirc, 1);//Slope/centre
gel(ret, 2)=gel(newcirc, 2);//x or y intercept/radius
gel(ret, 3)=endpt1;//Start point
gel(ret, 4)=endpt2;//end point. These may be in the wrong order, so will fix later if so.
if(gequal0(gel(newcirc, 3))){//Circle
gel(ret, 8)=gen_0;//arc
gel(ret, 5)=radialangle(newcirc, endpt1, tol, prec);//angle 1
gel(ret, 6)=shiftangle(radialangle(newcirc, endpt2, tol, prec), gel(ret, 5), tol, prec);//angle 2
if(isarc) gel(ret, 7)=gel(c, 7);//Temporary; match the old to the new direction
else gel(ret, 7)=gen_1;//Temporary; since for segments we have a bona fide start and end, we start with the direction being 1.
if(!onarc(ret, extpt, tol, prec)){//Must swap start/endpoints, angles and the direction
gel(ret, 3)=endpt2;
gel(ret, 4)=endpt1;
GEN tempang=gel(ret, 5);//angle to endpt1
gel(ret, 5)=shiftangle(gel(ret, 6), gen_0, tol, prec);//The angle to endpt2 shifted to [0, 2*Pi)
gel(ret, 6)=shiftangle(tempang, gel(ret, 5), tol, prec);//Angle to endpt1 shifted with a base of the angle to endpt2
gel(ret, 7)=gneg(gel(ret, 7));//We now run backwards
}
}
else{//Line
if(isarc==1 && gequal(gel(c, 7), gen_m1)){//We need to reverse the order of the points, because the arc ran backwards.
gel(ret, 3)=endpt2;
gel(ret, 4)=endpt1;
}
gel(ret, 8)=gen_1; //segment
gel(ret, 5)=gen_0;//Unused
if(typ(endpt1)==t_INFINITY || typ(endpt2)==t_INFINITY){//oo endpoint
gel(ret, 6)=gen_1;//Temporary, assume we go up/right
gel(ret, 7)=gen_0;//Not used
if(!onseg(ret, extpt, tol, prec)) gel(ret, 6)=gen_m1;//We were wrong, and go down/left.
}
else{//Both endpoints finite
gel(ret, 6)=gen_0;
gel(ret, 7)=gen_1;//Temporary, assume we go through the plane only
if(!onseg(ret, extpt, tol, prec)) gel(ret, 7)=gen_m1;//We were wrong, and go through oo.
}
}
return gerepilecopy(top, ret);
}
//Mobius map acting on circle
static GEN mobius_circle(GEN M, GEN c, GEN tol, long prec){
pari_sp top=avma;
GEN p1=mat_eval(M, gadd(gel(c, 1), gel(c, 2)));//M(c[1]+c[2])
GEN p2=mat_eval(M, gadd(gel(c, 1), gmul(gel(c, 2), gen_I())));//M(c[1]+c[2]*I)
GEN p3=mat_eval(M, gsub(gel(c, 1), gel(c, 2)));//M(c[1]-c[2])
return gerepileupto(top, circle_fromppp(p1, p2, p3, tol, prec));
}
//Mobius map acting on line
static GEN mobius_line(GEN M, GEN l, GEN tol, long prec){
pari_sp top=avma;
GEN p1, p2, p3, I=gen_I();
if(typ(gel(l, 1))==t_INFINITY){//Vertical line
p1=mat_eval(M, gel(l, 2));//M(x-intercept)
p2=mat_eval(M, gadd(gel(l, 2), I));//M(x-intercept+I)
p3=mat_eval(M, gsub(gel(l, 2), I));//M(x-intercept-I)
}
else{//Non-vertical line
GEN slopeIp1=gaddgs(gmul(gel(l, 1), I), 1);//1+Slope*I
GEN p1base=gmul(gel(l, 2), I);//y-intercept
GEN p2base=gadd(p1base, slopeIp1);//y-intercept+1+slope*I
GEN p3base=gadd(p2base, slopeIp1);//y-intercept+2+2*slope*I
p1=mat_eval(M, p1base);p2=mat_eval(M, p2base);p3=mat_eval(M, p3base);
}
return gerepileupto(top, circle_fromppp(p1, p2, p3, tol, prec));
}
//Perpendicular bisector of distinct points
GEN perpbis(GEN p1, GEN p2){
pari_sp top=avma;
return gerepileupto(top, line_fromsp(divoo(gen_m1, slope(p1, p2)), midpoint(p1, p2)));
}
//Angle between p and the centre of c, in the range [0, 2*Pi)
GEN radialangle(GEN c, GEN p, GEN tol, long prec){
pari_sp top=avma;
return gerepileupto(top, shiftangle(garg(gsub(p, gel(c, 1)), prec), gen_0, tol, prec));
}
//Radialangle with default tolerance.
GEN radialangle_tc(GEN c, GEN p, long prec){
int i1=geom_check(c);
if(i1!=0 && i1!=2) pari_err_TYPE("Please input a circle/arc", c);
pari_sp top=avma;
GEN tol=deftol(prec);
return gerepileupto(top, radialangle(c, p, tol, prec));
}
/*
slope, intercept, startpt, endpt, 0, ooendptor, dir, 1
GEN seg_init(GEN p1, GEN p2, int thruoo, long prec){
GEN l=line_frompp(p1, p2);
}*/
//The slope of the line through p1, p2
GEN slope(GEN p1, GEN p2){
pari_sp top=avma;
return gerepileupto(top, divoo(imag_i(gsub(p2, p1)), real_i(gsub(p2, p1))));
}
//INTERSECTION OF LINES/CIRCLES
//Returns the intersection points of two arcs
GEN arc_int(GEN c1, GEN c2, GEN tol, long prec){
pari_sp top=avma;
GEN ipts=circle_int(c1, c2, tol, prec);
if(lg(ipts)==1){avma=top;return cgetg(1, t_VEC);}//No intersection
if(lg(ipts)==2){//One intersection point (tangent circles)
if(!onarc(c1, gel(ipts, 1), tol, prec)){avma=top;return cgetg(1, t_VEC);}//Not on arc 1
if(!onarc(c2, gel(ipts, 1), tol, prec)){avma=top;return cgetg(1, t_VEC);}//Not on arc 2
return gerepilecopy(top, ipts);//On arc
}
//Two intersections
int i1=onarc(c1, gel(ipts, 1), tol, prec);
if(i1==1) i1=onarc(c2, gel(ipts, 1), tol, prec);//Now i1==1 iff the ipts[1] is on both c1 and c2
int i2=onarc(c1, gel(ipts, 2), tol, prec);
if(i2==1) i2=onarc(c2, gel(ipts, 2), tol, prec);//Now i2==1 iff the ipts[2] is on both c1 and c2
if(i1==1){
if(i2==1) return gerepilecopy(top, ipts);//Both pts on the arcs
GEN ret=cgetg(2, t_VEC);//Just point 1
gel(ret, 1)=gcopy(gel(ipts, 1));
return gerepileupto(top, ret);
}
//Now i1=0
if(i2==0){avma=top;return cgetg(1, t_VEC);}//Not on either arc
GEN ret=cgetg(2, t_VEC);//Just point 2
gel(ret, 1)=gcopy(gel(ipts, 2));
return gerepileupto(top, ret);
}
//arc_int with typecheck
GEN arc_int_tc(GEN c1, GEN c2, long prec){
int i1=geom_check(c1), i2=geom_check(c2);
if(i1!=2 || i2!=2) pari_err_TYPE("Please input circle arcs", c1);
pari_sp top=avma;
GEN tol=deftol(prec);
return gerepileupto(top, arc_int(c1, c2, tol, prec));
}
//Returns the intersection points of an arc and a segment
GEN arcseg_int(GEN c, GEN l, GEN tol, long prec){
pari_sp top=avma;
GEN ipts=circleline_int(c, l, tol, prec);
if(lg(ipts)==1){avma=top;return cgetg(1, t_VEC);}//No intersection
if(lg(ipts)==2){//One intersection point (tangent circles)
if(!onarc(c, gel(ipts, 1), tol, prec)){avma=top;return cgetg(1, t_VEC);}//Not on arc
if(!onseg(l, gel(ipts, 1), tol, prec)){avma=top;return cgetg(1, t_VEC);}//Not on segment
return gerepilecopy(top, ipts);//On both
}
//Two intersections
int i1=onarc(c, gel(ipts, 1), tol, prec);
if(i1==1) i1=onseg(l, gel(ipts, 1), tol, prec);//Now i1==1 iff the ipts[1] is on both c and l
int i2=onarc(c, gel(ipts, 2), tol, prec);
if(i2==1) i2=onseg(l, gel(ipts, 2), tol, prec);//Now i2==1 iff the ipts[2] is on both c and l
if(i1==1){
if(i2==1) return gerepilecopy(top, ipts);//Both pts on both
GEN ret=cgetg(2, t_VEC);//Just point 1
gel(ret, 1)=gcopy(gel(ipts, 1));
return gerepileupto(top, ret);
}
//Now i1=0
if(i2==0){avma=top;return cgetg(1, t_VEC);}//Not on either
GEN ret=cgetg(2, t_VEC);//Just point 2
gel(ret, 1)=gcopy(gel(ipts, 2));
return gerepileupto(top, ret);
}
//arc_int with typecheck
GEN arcseg_int_tc(GEN c, GEN l, long prec){
int i1=geom_check(c), i2=geom_check(l);
if(i1!=2 || i2!=3) pari_err_TYPE("Please input a circle arc and line segment", c);
pari_sp top=avma;
GEN tol=deftol(prec);
return gerepileupto(top, arcseg_int(c, l, tol, prec));
}
//Returns the set of points in the intersection of circles c1, c2
GEN circle_int(GEN c1, GEN c2, GEN tol, long prec){
pari_sp top=avma;
GEN a1=real_i(gel(c1, 1)), b1=imag_i(gel(c1, 1)), r1=gel(c1, 2);//x, y coords and radius of c1
GEN a2=real_i(gel(c2, 1)), b2=imag_i(gel(c2, 1)), r2=gel(c2, 2);//x, y coords and radius of c2
GEN a1ma2=gsub(a1, a2), b1mb2=gsub(b1, b2), x1, x2, y1, y2;
int oneint=0;
if(gcmp(gabs(a1ma2, prec), gabs(b1mb2, prec))>=0){//We want to divide by the larger of the two quantities to maximize precision and avoid errors when the centres are on the same line.
if(toleq(a1ma2, gen_0, tol, prec)==1){avma=top;return cgetg(1, t_VEC);}//Same centre, cannot intersect.
//u=(r1^2-r2^2+b2^2-b1^2+a2^2-a1^2)/(2*a2-2*a1)-a1;
GEN u=gsub(gdiv(gsub(gadd(gsub(gadd(gsub(gsqr(r1), gsqr(r2)), gsqr(b2)), gsqr(b1)), gsqr(a2)), gsqr(a1)), gmulgs(a1ma2, -2)), a1);
GEN v=gneg(gdiv(b1mb2, a1ma2));//v=(b1-b2)/(a2-a1), and x=a1+u+vy
GEN uvmb1=gsub(gmul(u, v), b1);//uv-b1
GEN vsqrp1=gaddgs(gsqr(v), 1);//v^2+1
GEN rtpart=gsub(gsqr(uvmb1), gmul(vsqrp1, gadd(gsqr(b1), gsub(gsqr(u), gsqr(r1)))));//(u*v-b1)^2-(v^2+1)*(b1^2+u^2-r1^2)
oneint=tolcmp(rtpart, gen_0, tol, prec);//Comparing rtpart to 0
if(oneint==-1){avma=top;return cgetg(1, t_VEC);}//rtpart must be square rooted, so if it's negative the circles do not intersect
if(oneint==0){//One intersection, so we take rtpart=0. This is CRUCIAL, as taking the square root kills our precision if we don't do this here.
y1=gdiv(gneg(uvmb1), vsqrp1);//y1=(b1-u*v)/(1*v^2+1)
x1=gadd(gadd(a1, u), gmul(v, y1));//x1=a1+u+v*y1
}
else{
y1=gdiv(gadd(gneg(uvmb1), gsqrt(rtpart, prec)), vsqrp1);//y1=(b1-u*v+sqrt((u*v-b1)^2-*(v^2+1)*(b1^2+u^2-r1^2)))/(1*v^2+1)
y2=gadd(gneg(y1), gdiv(gmulgs(uvmb1, -2), vsqrp1));//y2=-y1+(2*b1-2*u*v)/(v^2+1)
GEN a1pu=gadd(a1, u);
x1=gadd(a1pu, gmul(v, y1));//x1=a1+u+v*y1
x2=gadd(a1pu, gmul(v, y2));//x1=a1+u+v*y2
}
}
else{
if(toleq(b1mb2, gen_0, tol, prec)==1){avma=top;return cgetg(1, t_VEC);}//Same centre, cannot intersect.
//u=(r1^2-r2^2+b2^2-b1^2+a2^2-a1^2)/(2*b2-2*b1)-b1;
GEN u=gsub(gdiv(gsub(gadd(gsub(gadd(gsub(gsqr(r1), gsqr(r2)), gsqr(b2)), gsqr(b1)), gsqr(a2)), gsqr(a1)), gmulgs(b1mb2, -2)), b1);
GEN v=gneg(gdiv(a1ma2, b1mb2));//v=(a1-a2)/(b2-b1), and y=b1+u+vx
GEN uvma1=gsub(gmul(u, v), a1);//uv-a1
GEN vsqrp1=gaddgs(gsqr(v), 1);//v^2+1
GEN rtpart=gsub(gsqr(uvma1), gmul(vsqrp1, gadd(gsqr(a1), gsub(gsqr(u), gsqr(r1)))));//(u*v-a1)^2-(v^2+1)*(a1^2+u^2-r1^2))
oneint=tolcmp(rtpart, gen_0, tol, prec);//Comparing rtpart to 0
if(oneint==-1){avma=top;return cgetg(1, t_VEC);}//rtpart must be square rooted, so if it's negative the circles do not intersect
if(oneint==0){//One intersection, so we take rtpart=0. This is CRUCIAL, as taking the square root kills our precision if we don't do this here.
x1=gdiv(gneg(uvma1), vsqrp1);//x1=(a1-u*v)/(v^2+1);
y1=gadd(gadd(b1, u), gmul(v, x1));//y1=b1+u+v*x1;
}
else{
x1=gdiv(gadd(gneg(uvma1), gsqrt(rtpart, prec)), vsqrp1);//x1=(a1-u*v+sqrt((u*v-a1)^2-(v^2+1)*(a1^2+u^2-r1^2)))/(v^2+1);
x2=gadd(gneg(x1), gdiv(gmulgs(uvma1, -2), vsqrp1));//x2=-x1+(2*a1-2*u*v)/(v^2+1)
GEN b1pu=gadd(b1, u);
y1=gadd(b1pu, gmul(v, x1));//y1=b1+u+v*x1;
y2=gadd(b1pu, gmul(v, x2));//y2=b1+u+v*x2;
}
}
if(oneint==0){//One point of intersection (0 pts of intersection was already dealt with and returned)
GEN y1I=gmul(gen_I(), y1);
GEN ret=cgetg(2, t_VEC);
gel(ret, 1)=gadd(x1, y1I);
return gerepileupto(top, ret);
}
GEN y1I=gmul(gen_I(), y1), y2I=gmul(gen_I(), y2);
GEN ret=cgetg(3, t_VEC);
gel(ret, 1)=gadd(x1, y1I);
gel(ret, 2)=gadd(x2, y2I);
return gerepileupto(top, ret);
}
//Checks that c1, c2 are circles and calls circle_int with the default tolerance.
GEN circle_int_tc(GEN c1, GEN c2, long prec){
int i1=geom_check(c1), i2=geom_check(c2);
if(i1!=0 || i2!=0) pari_err_TYPE("Please input circles", c1);
pari_sp top=avma;
GEN tol=deftol(prec);
return gerepileupto(top, circle_int(c1, c2, tol, prec));
}
//Returns the intersection points of c and l
GEN circleline_int(GEN c, GEN l, GEN tol, long prec){
pari_sp top=avma;
if(typ(gel(l, 1))==t_INFINITY){
GEN x1=gel(l, 2);
GEN rtpart=gsub(gsqr(gel(c, 2)), gsqr(gsub(x1, real_i(gel(c, 1)))));//c[2]^2-(x1-real(c[1]))^2
if(gsigne(rtpart)==-1){avma=top;return cgetg(1, t_VEC);}//No intersections.
GEN y1=gadd(imag_i(gel(c, 1)), gsqrt(rtpart, prec));//y1=imag(c[1])+sqrt(c[2]^2-(x1-real(c[1]))^2)
if(toleq(rtpart, gen_0, tol, prec)){//Only one intersection point
GEN ret=cgetg(2, t_VEC);
gel(ret, 1)=cgetg(3, t_COMPLEX);
gel(gel(ret, 1), 1)=gcopy(x1);
gel(gel(ret, 1), 2)=gcopy(y1);
return gerepileupto(top, ret);
}
//Two intersection points
GEN y1py2=gmulgs(imag_i(gel(c, 1)), 2);//2*imag(c[1])
GEN ret=cgetg(3, t_VEC);
gel(ret, 1)=cgetg(3, t_COMPLEX);gel(ret, 2)=cgetg(3, t_COMPLEX);
gel(gel(ret, 1), 1)=gcopy(x1);
gel(gel(ret, 1), 2)=gcopy(y1);
gel(gel(ret, 2), 1)=gcopy(x1);
gel(gel(ret, 2), 2)=gsub(y1py2, y1);
return gerepileupto(top, ret);
}
//Now y=mx+b with m finite
GEN A=gaddgs(gsqr(gel(l, 1)), 1);//l[1]^2+1
GEN l2mic1=gsub(gel(l, 2), imag_i(gel(c, 1)));//l[2]-imag(c[1])
GEN B=gadd(gmulgs(real_i(gel(c, 1)), -2), gmulsg(2, gmul(gel(l, 1), l2mic1)));//-2*real(c[1])+2*l[1]*(l[2]-imag(c[1]))
GEN C=gadd(gsqr(real_i(gel(c, 1))), gsub(gsqr(l2mic1), gsqr(gel(c, 2))));//real(c[1])^2+(l[2]-imag(c[1]))^2-c[2]^2
GEN rtpart=gsub(gsqr(B), gmulsg(4, gmul(A, C)));
int rtpartsig=tolcmp(rtpart, gen_0, tol, prec);
if(rtpartsig==-1){avma=top;return cgetg(1, t_VEC);}//No intersection
if(rtpartsig==0){//One root, and rtpart=0
GEN x1=gdiv(B, gmulgs(A, -2));//-B/(2A)
GEN y1part=gmul(gel(l, 1), x1);//l[1]*x1
GEN ret=cgetg(2, t_VEC);
gel(ret, 1)=cgetg(3, t_COMPLEX);
gel(gel(ret, 1), 1)=gcopy(x1);
gel(gel(ret, 1), 2)=gadd(y1part, gel(l, 2));//y1=l[1]*x1+l[2];
return gerepileupto(top, ret);
}
//Two roots
GEN x1=gdiv(gsub(gsqrt(rtpart, prec), B), gmulgs(A, 2));//x1=(-B+sqrt(B^2-4*A*C))/(2*A);
GEN x2=gsub(gneg(gdiv(B, A)), x1);//-B/A-x1
GEN y1part=gmul(gel(l, 1), x1);//l[1]*x1
GEN y2part=gmul(gel(l, 1), x2);//l[1]*x2
GEN ret=cgetg(3, t_VEC);
gel(ret, 1)=cgetg(3, t_COMPLEX);gel(ret, 2)=cgetg(3, t_COMPLEX);
gel(gel(ret, 1), 1)=gcopy(x1);
gel(gel(ret, 1), 2)=gadd(y1part, gel(l, 2));//l[1]*x1+l[2];
gel(gel(ret, 2), 1)=gcopy(x2);
gel(gel(ret, 2), 2)=gadd(y2part, gel(l, 2));//l[1]*x2+l[2];
return gerepileupto(top, ret);
}
//Checks that c is a circle and l is a line and calls circleline_int with the default tolerance.
GEN circleline_int_tc(GEN c, GEN l, long prec){
int i1=geom_check(c), i2=geom_check(l);
if(i1!=0 || i2!=1) pari_err_TYPE("Please input a line and circle", c);
pari_sp top=avma;
GEN tol=deftol(prec);
return gerepileupto(top, circleline_int(c, l, tol, prec));
}
//Returns the intersection of two circles/lines/arcs/segments (in any combination)
GEN genseg_int(GEN s1, GEN s2, GEN tol, long prec){
int t1=geom_check(s1), t2=geom_check(s2);
if(t1==-1 || t2==-1) pari_err_TYPE("Please input two circles/lines/arcs/segments", s1);
switch(t1){
case 0:
switch(t2){
case 0: return circle_int(s1, s2, tol, prec);//Circle circle
case 1: return circleline_int(s1, s2, tol, prec);//Circle line
case 2: return arc_int(s1, s2, tol, prec);//Circle arc
case 3: return arcseg_int(s1, s2, tol, prec);//Circle segment
}
case 1:
switch(t2){
case 0: return circleline_int(s2, s1, tol, prec);//Line circle
case 1: return line_int(s1, s2, tol, prec);//Line line
case 2: return arcseg_int(s2, s1, tol, prec);//Line arc
case 3: return seg_int(s1, s2, tol, prec);//Line segment
}
case 2:
switch(t2){
case 0: return arc_int(s1, s2, tol, prec);//Arc circle
case 1: return arcseg_int(s1, s2, tol, prec);//Arc line
case 2: return arc_int(s1, s2, tol, prec);//Arc arc
case 3: return arcseg_int(s1, s2, tol, prec);//Arc segment
}
case 3:
switch(t2){
case 0: return arcseg_int(s2, s1, tol, prec);//segment circle
case 1: return seg_int(s1, s2, tol, prec);//segment line
case 2: return arcseg_int(s2, s1, tol, prec);//segment arc
case 3: return seg_int(s1, s2, tol, prec);//segment seg
}
}
pari_err_TYPE("ERROR: UNEXPECTED INPUT TYPES", mkvec2(s1, s2));
return gen_0;
}
//genseg_int with default tolerance.
GEN genseg_int_tc(GEN s1, GEN s2, long prec){
pari_sp top=avma;
GEN tol=deftol(prec);
return gerepileupto(top, genseg_int(s1, s2, tol, prec));
}
//The intersection of two lines
GEN line_int(GEN l1, GEN l2, GEN tol, long prec){
GEN s1=gel(l1, 1), s2=gel(l2, 1);//Slopes
if(toleq(s1, s2, tol, prec)) return mkoo();//Parallel or equal
pari_sp top=avma;
if(typ(s1)==t_INFINITY){//l1 vertical
GEN ypart=gmul(s2, gel(l1, 2));//s2*l1[2]
GEN ipt=cgetg(3, t_COMPLEX);
gel(ipt, 1)=gcopy(gel(l1, 2));
gel(ipt, 2)=gadd(ypart, gel(l2, 2));//s2*l1[2]+l2[2]
return gerepileupto(top, ipt);
}
if(typ(s2)==t_INFINITY){//l2 vertical
GEN ypart=gmul(s1, gel(l2, 2));//s1*l2[2]
GEN ipt=cgetg(3, t_COMPLEX);
gel(ipt, 1)=gcopy(gel(l2, 2));
gel(ipt, 2)=gadd(ypart, gel(l1, 2));//s1*l2[2]+l1[2]
return gerepileupto(top, ipt);
}
GEN x=gdiv(gsub(gel(l2, 2), gel(l1, 2)), gsub(s1, s2));//(l2[2]-l1[2])/(s1-s2)
GEN ypart=gmul(s1, x);//s1*x
GEN ipt=cgetg(3, t_COMPLEX);
gel(ipt, 1)=gcopy(x);
gel(ipt, 2)=gadd(ypart, gel(l1, 2));//s1*x+l1[2]
return gerepileupto(top, ipt);
}
//Checks that l1, l2 are lines and calls line_int with the default tolerance.
GEN line_int_tc(GEN l1, GEN l2, long prec){
int i1=geom_check(l1), i2=geom_check(l2);
if(i1!=1 || i2!=1) pari_err_TYPE("Please input lines", l1);
pari_sp top=avma;
GEN tol=deftol(prec);
return gerepileupto(top, line_int(l1, l2, tol, prec));
}
//p is assumed to be on the circle defined by c; this checks if it is actually on the arc (running counterclockwise from c[3] to c[4]).
int onarc(GEN c, GEN p, GEN tol, long prec){
if(lg(c)==CIRCLEN) return 1;//Allow input of just a circle, so the return is trivially 1
if(toleq(gel(c, 3), p, tol, prec)) return 1;//p=the start point. We have this done seperately in case rounding errors take the angle to <c[5], as this will cause issues with the shifting angle.
pari_sp top=avma;
GEN angle=shiftangle(radialangle(c, p, tol, prec), gel(c, 5), tol, prec);//Getting the angle in the range [c[5], c[5]+2*Pi)
if(tolcmp(angle, gel(c, 6), tol, prec)<=0){avma=top;return 1;}//On the arc
avma=top;
return 0;//Beyond the arc.
}
//onarc with checking c and passing the default tolerance.
int onarc_tc(GEN c, GEN p, long prec){
int i1=geom_check(c);
if(i1!=0 && i1!=2) pari_err_TYPE("Please input a circle arc and point", c);
pari_sp top=avma;
GEN tol=deftol(prec);
int oa=onarc(c, p, tol, prec);
avma=top;
return oa;
}
//p is assumed to be on the line defined by l; this checks if it is actually on the segment l
int onseg(GEN l, GEN p, GEN tol, long prec){
if(lg(l)==CIRCLEN) return 1;//Allow input of a line, so return is trivially 1
if(typ(p)==t_INFINITY){//p is the point at oo
if(!gequal0(gel(l, 6)) || gequal(gel(l, 7), gen_m1)) return 1;//oo is an endpoint OR the seg passes through oo
return 0;//If not, does not pass through oo
}
pari_sp top=avma;
//Okay, now p is not oo and l is a line segment
if(typ(gel(l, 1))==t_INFINITY){//Vertical line
if(typ(gel(l, 3))==t_INFINITY){//Start point in oo
if(equali1(gel(l, 6))){//p must lie BELOW l[4]
if(tolcmp(imag_i(p), imag_i(gel(l, 4)), tol, prec)<=0){avma=top;return 1;}//Lies below l[4]
avma=top;return 0;//lies above l[4]
}
//p must lie ABOVE l[4]
if(tolcmp(imag_i(p), imag_i(gel(l, 4)), tol, prec)>=0){avma=top;return 1;}//Lies above l[4]
avma=top;return 0;//lies below l[4]
}
if(typ(gel(l, 4))==t_INFINITY){//End point is oo
if(equali1(gel(l, 6))){//p must lie ABOVE l[3]
if(tolcmp(imag_i(p), imag_i(gel(l, 3)), tol, prec)>=0){avma=top;return 1;}//Lies above l[3]
avma=top;return 0;//lies below l[3]
}
//p must lie BELOW l[3]
if(tolcmp(imag_i(p), imag_i(gel(l, 3)), tol, prec)<=0){avma=top;return 1;}//Lies below l[3]
avma=top;return 0;//lies above l[3]
}
//Start and end points are finite
int i1=tolcmp(imag_i(gsub(p, gel(l, 3))), gen_0, tol, prec);//sign of imag(p)-imag(l[3])
int i2=tolcmp(imag_i(gsub(p, gel(l, 4))), gen_0, tol, prec);//sign of imag(p)-imag(l[4])
avma=top;
if(i1==0 || i2==0) return 1;//endpoint
if(i1==i2){//p on the same side of l[3] and l[4], so return 1 iff l passes through oo
if(gequal(gel(l, 7), gen_1)) return 0;//Not through oo
return 1;//through oo
}
//p is between l[3] and l[4], so return 1 iff l does not pass through oo
if(gequal(gel(l, 7), gen_1)) return 1;//not through oo
return 0;//through oo
}
//Non-vertical line
if(typ(gel(l, 3))==t_INFINITY){//Start point in oo
if(equali1(gel(l, 6))){//p must lie LEFT OF l[4]
if(tolcmp(real_i(p), real_i(gel(l, 4)), tol, prec)<=0){avma=top;return 1;}//Lies left of l[4]
avma=top;return 0;//lies right of l[4]
}
//p must lie RIGHT OF l[4]
if(tolcmp(real_i(p), real_i(gel(l, 4)), tol, prec)>=0){avma=top;return 1;}//Lies right of l[4]
avma=top;return 0;//lies left of l[4]
}
if(typ(gel(l, 4))==t_INFINITY){//End point is oo
if(equali1(gel(l, 6))){//p must lie RIGHT OF l[3]
if(tolcmp(real_i(p), real_i(gel(l, 3)), tol, prec)>=0){avma=top;return 1;}//Lies right of l[3]
avma=top;return 0;//lies below l[3]
}
//p must lie LEFT OF l[3]
if(tolcmp(real_i(p), real_i(gel(l, 3)), tol, prec)<=0){avma=top;return 1;}//Lies left of l[3]
avma=top;return 0;//lies above l[3]
}
//Start and end points are finite
int i1=tolcmp(real_i(gsub(p, gel(l, 3))), gen_0, tol, prec);//sign of real(p)-real(l[3])
int i2=tolcmp(real_i(gsub(p, gel(l, 4))), gen_0, tol, prec);//sign of real(p)-real(l[4])
avma=top;
if(i1==0 || i2==0) return 1;//endpoint
if(i1==i2){//p on the same side of l[3] and l[4], so return 1 iff l passes through oo
if(gequal(gel(l, 7), gen_1)) return 0;//Not through oo
return 1;//through oo
}
//p is between l[3] and l[4], so return 1 iff l does not pass through oo
if(gequal(gel(l, 7), gen_1)) return 1;//not through oo
return 0;//through oo
}
//onseg with checking l and passing the default tolerance.
int onseg_tc(GEN l, GEN p, long prec){
int i1=geom_check(l);
if(i1!=1 && i1!=3) pari_err_TYPE("Please input a line segment and point", l);
pari_sp top=avma;
GEN tol=deftol(prec);
int os=onseg(l, p, tol, prec);
avma=top;
return os;
}
//Returns the intersection points of l1 and l2 (length 0 or 1 vector)
GEN seg_int(GEN l1, GEN l2, GEN tol, long prec){
pari_sp top=avma;
GEN in=line_int(l1, l2, tol, prec);//Line intersection
if(onseg(l1, in, tol, prec) && onseg(l2, in, tol, prec)){//On the segments
GEN ret=cgetg(2, t_VEC);
gel(ret, 1)=gcopy(in);
return gerepileupto(top, ret);
}
avma=top;
return cgetg(1, t_VEC);//Not on at least one segment.
}
//Checks that l1, l2 are line segments and calls seg_int with the default tolerance.
GEN seg_int_tc(GEN l1, GEN l2, long prec){
int i1=geom_check(l1), i2=geom_check(l2);
if(i1!=3 || i2!=3) pari_err_TYPE("Please input line segments", l1);
pari_sp top=avma;
GEN tol=deftol(prec);
return gerepileupto(top, seg_int(l1, l2, tol, prec));
}
//DISTANCES
//z1 and z2 are complex numbers, this computes the hyperbolic distance between them.
GEN hdist(GEN z1, GEN z2, long prec){
pari_sp top=avma;
GEN x1=gel(z1,1);
GEN y1=gel(z1,2);
GEN x2=gel(z2,1);
GEN y2=gel(z2,2);
GEN x=gaddsg(1,gdiv(gadd(gsqr(gsub(x2,x1)),gsqr(gsub(y2,y1))),gmul(gmulsg(2,y1),y2)));
GEN expd=gadd(x,gsqrt(gsubgs(gsqr(x), 1), prec));
return gerepileupto(top,glog(expd, prec));
}
//hdist with typechecking
GEN hdist_tc(GEN z1, GEN z2, long prec){
if(typ(z1)!=t_COMPLEX || signe(gel(z1,2))!=1) pari_err_TYPE("Please input complex numbers in the upper half plane", z1);
if(typ(z2)!=t_COMPLEX || signe(gel(z2,2))!=1) pari_err_TYPE("Please input complex numbers in the upper half plane", z2);
return hdist(z1, z2, prec);
}
//The hyperbolic distance between z1 and z2 in the unit disc model
GEN hdist_ud(GEN z1, GEN z2, long prec){
pari_sp top=avma;
GEN a = gabs(gsubsg(1, gmul(z1, conj_i(z2))), prec);//|1-z1*conj(z2)|
GEN b = gabs(gsub(z1, z2), prec);//|z1-z2|
GEN num=gadd(a, b);
GEN denom=gsub(a, b);
GEN ret;
pari_CATCH(e_INV){
avma=top;
return mkoo();
}
pari_TRY{
ret=gerepileupto(top, glog(gdiv(num, denom), prec));//log((a+b)/(a-b))
}
pari_ENDCATCH
return ret;
}
//Given the isometric circles (in order) and the vertices (so that circles[i] intersect circles[i+1] is vertices[i]), returns the area of the convex hyperbolic polygon. If one of the sides is oo (an entry of circles is 0), the answer will be oo. If there are n vertices with angles a1,...,an, then the area is is (n-2)Pi-sum(ai)
GEN hpolygon_area(GEN circles, GEN vertices, GEN tol, long prec){
pari_sp top=avma;
long blen=lg(circles);
if(blen==1 || gequal0(gel(circles, 1))) return mkoo();//No cicles or the first is 0, i.e. an infinite side
GEN ang, area=gmulsg(blen-3, mppi(prec));//We subtract off from area.
for(long i=1;i<blen-1;i++){
if(gequal0(gel(circles, i+1))){avma=top;return mkoo();}//The next side is infinite
ang=circle_angle(gel(circles, i+1), gel(circles, i), gel(vertices, i), tol, prec);//Do in opposite order to get correct angle
area=gsub(area, ang);
}
ang=circle_angle(gel(circles, 1), gel(circles, blen-1), gel(vertices, blen-1), tol, prec);//The last side wraps around.
area=gsub(area, ang);
return gerepileupto(top, area);
}
//Normalized boundary area with typechecking
GEN hpolygon_area_tc(GEN circles, GEN vertices, long prec){
pari_sp top=avma;
if(typ(circles)!=t_VEC || typ(vertices)!=t_VEC) pari_err_TYPE("Please enter the edges and vertices as vectors", circles);
if(lg(circles)!=lg(vertices)) pari_err_TYPE("There needs to be the same number of edges as vertices", vertices);
GEN tol=deftol(prec);
return gerepileupto(top, hpolygon_area(circles, vertices, tol, prec));
}
//FUNDAMENTAL DOMAIN COMPUTATION
/*Let Gamma be a Fuschian group; the following methods allow us to compute a fundamental domain for Gamma, following the paper of John Voight. We are assuming that we have an "exact" way to deal with the elements of Gamma. In otherwords, we need separate methods to:
i) Multiply elements of Gamma: format as GEN eltmul(GEN *data, GEN x, GEN y), with data the extra data you need.
ii) Find the image of g in PSL(1, 1): format as GEN gamtopsl(GEN *data, GEN g, long prec), with data the extra data you need (e.g. the quaterniona algebra in the case of Shimura curves).
iii) Identify if g is trivial in Gamma: format as int istriv(GEN *data, GEN g). Need to take care as we care about being trivial in PSL, whereas most representations of Gamma are for SL (so need to check with +/-id).
iv) Pass in the identity element of Gamma and find the area of the fundamental domain. These methods are not passed in; just the values.
We pass references to these methods into the methods here.
We do computations mostly in PSU, and shift from PSL to PSU via phi(z):=(z-p)/(z-conj(p)) for some given upper half plane point p.
*/
//Some of the methods may not be that useful (i.e. assuming we have p only and not the matrices).
//Returns the edge pairing, as VECSMALL v where v[i]=j means i is paired with j. If not all edges can be paired, instead returns [v1, v2, ...] where vi is a VECSMALL that is either [gind, v1ind] or [gind, v1ind, v2ind]. gind is the index of the unpaired side, and the viind are the corresponding unpaired vertices (1 or 2). If rboth=1, returns [paired, unpaired], where this time paired is a vector of vecsmalls [i, j] meaning i is paired to j (differing to the output format when rboth=0 and everything is paired)
GEN edgepairing(GEN U, GEN tol, int rboth, long prec){
pari_sp top=avma;
GEN vangles=gel(U, 4);//Vertex angles
GEN baseangle=gel(vangles, 1);//Base angle
GEN toldata=cgetg(3, t_VEC);//Stores the necessary info for searching with tolerance (tolcmp_sort)
gel(toldata, 1)=tol;
gel(toldata, 2)=cgetg(2, t_VECSMALL);
gel(toldata, 2)[1]=prec;
long lU=lg(gel(U, 1));
GEN unpair=vectrunc_init(lU+1), vim, vimang, pair=vectrunc_init(lU);//Unpair stores the unpaired edges, pair stores the paired edges
long ind1, ind2, i1, i2;
for(long i=1;i<lU;i++){
if(gequal0(gel(gel(U, 5), i))){vectrunc_append(pair, mkvecsmall2(i, i));continue;}//oo side, go next (we say it is paired with itself)
ind1=i;
vim=mat_eval(gel(gel(U, 5), i), gel(gel(U, 3), ind1));//The new vertex
vimang=shiftangle(garg(vim, prec), baseangle, tol, prec);//The new angle
i1=gen_search_old(vangles, vimang, 0, &toldata, &tolcmp_sort);
if(i1!=0) if(!toleq(vim, gel(gel(U, 3), i1), tol, prec)) i1=0;//Just because the angles are equal, the points don't have to be (though this occurence is expected to be extremely rare).
if(i==1) ind2=lU-1;
else ind2=i-1;//The two vertices of the side, this is the second one
vim=mat_eval(gel(gel(U, 5), i), gel(gel(U, 3), ind2));//The second new vertex
if(i1!=0){//If ind2 is paired, it MUST be paired to vertex i1+1
i2=i1+1;
if(i2==lU) i2=1;//Wrap around
if(!toleq(vim, gel(gel(U, 3), i2), tol, prec)) i2=0;
if(i2!=0){
if(i<=i2) vectrunc_append(pair, mkvecsmall2(i, i2));//i<=i2 is put so that we are not pairing things twice.
}
else vectrunc_append(unpair, mkvecsmall2(i, ind2));//Second vertex not paired
}
else{//ind1 not paired
vimang=shiftangle(garg(vim, prec), baseangle, tol, prec);//The second new angle
i2=gen_search_old(vangles, vimang, 0, &toldata, &tolcmp_sort);
if(i2!=0) if(!toleq(vim, gel(gel(U, 3), i2), tol, prec)) i2=0;//Just because the angles are equal, the points don't have to be (though this occurence is expected to be extremely rare).
if(i2!=0) vectrunc_append(unpair, mkvecsmall2(i, ind1));//First vtx not paired
else vectrunc_append(unpair, mkvecsmall3(i, ind1, ind2));//Neither vtx paired
}
}
if(rboth) return gerepilecopy(top, mkvec2(pair, unpair));
if(lg(unpair)==1){
GEN pairvs=cgetg(lU, t_VECSMALL);
for(long i=1;i<lg(pair);i++){
pairvs[gel(pair, i)[1]]=gel(pair, i)[2];
pairvs[gel(pair, i)[2]]=gel(pair, i)[1];
}
return gerepilecopy(top, pairvs);//No unpaired vertices
}
return gerepilecopy(top, unpair);//Unpaired vertices!
}
//Edgepairing with typecheck
GEN edgepairing_tc(GEN U, long prec){
pari_sp top=avma;
if(typ(U)!=t_VEC || lg(U)!=NORMBOUND) pari_err_TYPE("Please enter a normalized boundary", U);
GEN tol=deftol(prec);
return gerepileupto(top, edgepairing(U, tol, 1, prec));
}
//Computes the isometric circle for g, returning [g, image in PSU(1, 1), circle]. Must pass in mats (psltopsu_transmats(p)), and a method that translates g to an element of PSL(2, R).
GEN isometriccircle_mats(GEN g, GEN mats, GEN *data, GEN (*gamtopsl)(GEN *, GEN, long), GEN tol, long prec){
pari_sp top=avma;
GEN ginpsl=gamtopsl(data, g, prec);
GEN ret=cgetg(4, t_VEC);
gel(ret, 1)=gcopy(g);
gel(ret, 2)=psltopsu_mats(ginpsl, mats);
gel(ret, 3)=isometriccircle_psu(gel(ret, 2), tol, prec);
if(!gequal(gel(ret, 3), gen_m1)) return gerepileupto(top, ret);//Everything A-OK
//If we reach here, we did not have enough precision
long newprec=prec;
pari_CATCH(CATCH_ALL){
avma=top;
pari_CATCH_reset();
pari_err(e_MISC,"Could not increase precision enough. Please increase precision/memory");
return gen_0;
}
pari_TRY{
do{
if((prec2nbits(newprec) / prec2nbits(prec)) >=5) pari_err(e_MISC,"Throw");
avma=top;
newprec = nbits2prec(prec2nbits(newprec) + DEFAULTPREC);//Increase precision
tol=deftol(newprec);
if(precision(gel(mats, 3))>0){//p is inexact
GEN p=gtofp(gel(mats, 3), newprec);
mats=psltopsu_transmats(p);//Updating mats to more precision
}
ginpsl=gamtopsl(data, g, newprec);
ret=cgetg(4, t_VEC);
gel(ret, 1)=gcopy(g);
gel(ret, 2)=psltopsu_mats(ginpsl, mats);
gel(ret, 3)=isometriccircle_psu(gel(ret, 2), tol, newprec);
}
while(gequal(gel(ret, 3), gen_m1));
pari_CATCH_reset();
return gerepileupto(top, ret);
}
pari_ENDCATCH
}
//Not sure if useful
//Returns the isometric circle of an element of PSL
GEN isometriccircle_psl(GEN g, GEN p, GEN tol, long prec){
pari_sp top=avma;
GEN mats=psltopsu_transmats(p);
return gerepileupto(top, isometriccircle_psl_mats(g, mats, tol, prec));
}
//Not sure if useful
//Given the mats to translate to PSU, this returns the isometric circle of g in the unit disc model
GEN isometriccircle_psl_mats(GEN g, GEN mats, GEN tol, long prec){
pari_sp top=avma;
GEN gpsu=psltopsu_mats(g, mats);//Image in PSU
return gerepileupto(top, isometriccircle_psu(gpsu, tol, prec));
}
//Given an element g of PSU(1, 1), this returns the isometric circle associated to it.
GEN isometriccircle_psu(GEN g, GEN tol, long prec){
pari_sp top=avma;
if(toleq(gen_0, gcoeff(g, 2, 1), tol, prec)) return gen_0;//Isometric circle is everything, don't want to call it here.
GEN geod=zerovec(8);
gel(geod, 2)=gdivsg(1, gcoeff(g, 2, 1));//Need to take absolute value
gel(geod, 1)=gneg(gmul(gcoeff(g, 2, 2), gel(geod, 2)));//-g[2,2]/g[2,1], the centre of the circle
gel(geod, 2)=gabs(gel(geod, 2), prec);//We do things in this order to save a division.
gel(geod, 7)=gen_1;//Always oriented counterclockwise
pari_CATCH(CATCH_ALL){
avma=top;
pari_CATCH_reset();
return gen_m1;//We increase precision in g and retry.
}
pari_TRY{
GEN ipts=circle_int(geod, mkvec3s(0, 1, 0), tol, prec);//Intersect with x^2+y^2=1
GEN ang=anglediff(garg(gsub(gel(ipts, 2), gel(geod, 1)), prec), garg(gsub(gel(ipts, 1), gel(geod, 1)), prec), tol, prec);
if(gcmp(ang, mppi(prec))==1){//Properly orienting the start and endpoints
gel(geod, 3)=gel(ipts, 2);
gel(geod, 4)=gel(ipts, 1);
}
else{
gel(geod, 3)=gel(ipts, 1);
gel(geod, 4)=gel(ipts, 2);