1#include "kernel.pckh"
2
3Sign predicate(side3)(
4    point(p0), point(p1), point(p2), point(p3),
5    point(q0), point(q1), point(q2)  DIM
6) {
7
8    /* note: 1* for ensuring mcc's degree compatibility */
9    scalar l1 = 1*sq_dist(p1,p0);
10    scalar l2 = 1*sq_dist(p2,p0);
11    scalar l3 = 1*sq_dist(p3,p0);
12
13    scalar a10 = 2*dot_at(p1,q0,p0);
14    scalar a11 = 2*dot_at(p1,q1,p0);
15    scalar a12 = 2*dot_at(p1,q2,p0);
16    scalar a20 = 2*dot_at(p2,q0,p0);
17    scalar a21 = 2*dot_at(p2,q1,p0);
18    scalar a22 = 2*dot_at(p2,q2,p0);
19
20    scalar a30 = 2*dot_at(p3,q0,p0);
21    scalar a31 = 2*dot_at(p3,q1,p0);
22    scalar a32 = 2*dot_at(p3,q2,p0);
23
24    /*
25     * [ b00 b01 b02 ]           [  1   1   1  ]-1
26     * [ b10 b11 b12 ] = Delta * [ a10 a11 a12 ]
27     * [ b20 b21 b22 ]           [ a20 a21 a22 ]
28     */
29
30    scalar b00 = a11*a22-a12*a21;
31    scalar b01 = a21-a22;
32    scalar b02 = a12-a11;
33    scalar b10 = a12*a20-a10*a22;
34    scalar b11 = a22-a20;
35    scalar b12 = a10-a12;
36    scalar b20 = a10*a21-a11*a20;
37    scalar b21 = a20-a21;
38    scalar b22 = a11-a10;
39
40    scalar Delta = b00+b10+b20;
41
42    /*
43     *       [ Lambda0 ]   [ b01 b02 ]   [ l1 ]   [ b00 ]
44     * Delta [ Lambda1 ] = [ b11 b12 ] * [    ] + [ b10 ]
45     *       [ Lambda2 ]   [ b21 b22 ]   [ l2 ]   [ b20 ]
46     */
47
48    scalar DeltaLambda0 = b01*l1+b02*l2+b00 ;
49    scalar DeltaLambda1 = b11*l1+b12*l2+b10 ;
50    scalar DeltaLambda2 = b21*l1+b22*l2+b20 ;
51
52    scalar r = Delta*l3 - (
53           a30 * DeltaLambda0 +
54           a31 * DeltaLambda1 +
55           a32 * DeltaLambda2
56    ) ;
57
58    Sign Delta_sign = sign(Delta) ;
59    Sign r_sign     = sign(r) ;
60
61    generic_predicate_result(Delta_sign*r_sign) ;
62
63    begin_sos4(p0,p1,p2,p3)
64       sos(p0, Sign(Delta_sign*sign(
65                 Delta-((b01+b02)*a30+(b11+b12)*a31+(b21+b22)*a32)
66       )))
67       sos(p1, Sign(Delta_sign*sign((a30*b01)+(a31*b11)+(a32*b21))))
68       sos(p2, Sign(Delta_sign*sign((a30*b02)+(a31*b12)+(a32*b22))))
69       sos(p3, NEGATIVE)
70    end_sos
71}
72