// code associated with the paper https://arxiv.org/abs/2006.14232

#include <boost/numeric/interval.hpp>

using namespace boost::numeric;
using namespace interval_lib;

// Fix rounding policies for the transcendental functions
typedef interval<double, policies<save_state<rounded_transc_std<double>>, checking_base<double>>> I;

const I r=sqrt(I(2))-I(1); // radius

// global variables
I x; // proportion of large discs
I delta; // optimal density
I V111, V11r, V1r1, V1rr, Vr1r, Vrrr; // base vertex potentials
I a1, ar;
double m1, mr; // coefficient of angle deviation from base vertex potential
double Z1, Zr; // ceiling of vertex potential
double eps;
int nt=0; // number of tight triangles
int nb=0; // number of bad (not feasible) triangles
int ng=0; // number of good (validated) triangles


//***************************************************//
// useful routines

I val(char q)
{
    return q=='1'?I(1):r;
}

I Vijk(char i,char j, char k)
{
    std::string t;
    t+=i;
    t+=j;
    t+=k;
    if (t=="111") return V111;
    if (t=="11r") return V11r;
    if (t=="r11") return V11r;
    if (t=="r1r") return Vr1r;
    if (t=="1r1") return V1r1;
    if (t=="1rr") return V1rr;
    if (t=="rr1") return V1rr;
    if (t=="rrr") return Vrrr;
}

// angle opposite to side of length a
I angle(I a, I b, I c)
{
    return acos((b/c+c/b-square(a)/(b*c))/I(2));
}

// area of the triangle
I area(I a, I b, I c)
{
    return sqrt((a + b + c)*(a + b - c)*(a - b + c)*(- a + b + c))/I(4);
}

// area covered by disc sectors. ri is the radius in vertex opposed to edge i
I covered(I a, I b, I c, I ra, I rb, I rc)
{
    return (square(ra)*angle(a,b,c)+square(rb)*angle(b,c,a)+square(rc)*angle(c,a,b))/I(2);
}

// excess of the triangle
I excess(I a, I b, I c, I ra, I rb, I rc)
{
    return delta*area(a,b,c)-covered(a,b,c,ra,rb,rc);
}

//***************************************************//
// compute the radius of the support circle
// Input: side length a, b and c ot a triangle, radii ra, rb and rc of discs centered on vertices
I radius(I a, I b, I c, I ra, I rb, I rc)
{
    I S=(a+b+c)*(a+b-c)*(a-b+c)*(-a+b+c);
    I a2=square(a);
    I b2=square(b);
    I c2=square(c);
    I ra2=square(ra);
    I rb2=square(rb);
    I rc2=square(rc);
    I A=I(4)*((ra2-ra*rb-ra*rc+rb*rc)*a2 + (rb2-ra*rb+ra*rc-rb*rc)*b2 + (rc2+ra*rb-ra*rc-rb*rc)*c2)-S;
    I B=I(2)*(a2*(a2*ra-b2*(ra+rb)+((I)2*ra2*ra-ra2*rb-ra*rb2-ra2*rc+rb2*rc-ra*rc2+rb*rc2)) + b2*(b2*rb-c2*(rb+rc)+(I(2)*rb2*rb-rb2*rc-rb*rc2-rb2*ra+rc2*ra-rb*ra2+rc*ra2)) + c2*(c2*rc-a2*(rc+ra)+((I)2*rc2*rc-rc2*ra-rc*ra2-rc2*rb+ra2*rb-rc*rb2+ra*rb2)));
    I C=a2*b2*c2-(ra2+rb2)*a2*b2-(ra2+rc2)*a2*c2-(rb2+rc2)*b2*c2+(ra2*a2+square(ra2)-ra2*rb2-ra2*rc2+rb2*rc2)*a2+(rb2*b2-ra2*rb2+square(rb2)+ra2*rc2-rb2*rc2)*b2+(rc2*c2+ra2*rb2-ra2*rc2-rb2*rc2+square(rc2))*c2;
    if (zero_in(A) || zero_in(C))
    {
        I x=I(4)*A*C/square(B);
        return (upper(x)<0.78 ? -C/B*(I(1)+x) : I(-INFINITY,INFINITY));
    }
    else
    {
        I D=I(4)*S*(a+rb-rc)*(a-rb+rc)*(b+ra-rc)*(b-ra+rc)*(c+ra-rb)*(c-ra+rb);
        return (upper(C)<0 ? (-B+sqrt(D))/(I(2)*A) : (-B-sqrt(D))/(I(2)*A));
    }
}

//***************************************************//
// Input: idem as radius + rs radius support circle
// Output: signed distance to edge a of the center of the support circle
I dist(I a, I b, I c, I ra, I rb, I rc, I rs, I S)
{
    I a2=square(a);
    I b2=square(b);
    I c2=square(c);
    I ra2=square(ra);
    I rb2=square(rb);
    I rc2=square(rc);
    return ((b2+c2-a2 + I(2)*rs*(rb+rc-I(2)*ra)+rb2+rc2-I(2)*ra2)*a2 + (I(2)*rs+rb+rc)*(rb-rc)*(b2-c2))/(I(2)*a*sqrt(S));
}

//***************************************************//
// determination of eps

bool check(I ra,I rb,I rc, I eps)
{
    I a=rb+rc+I(0,upper(eps));
    I b=rc+ra+I(0,upper(eps));
    I c=ra+rb+I(0,upper(eps));
    I de=((square(a)+square(b)-square(c))*square(rb)+(square(a)+square(c)-square(b))*square(rc))/(I(2)*a)+a*((square(b)+square(c)-square(a))*delta-I(2)*square(ra))/I(2);
    I du=I(2)*a*(ra==I(1)? m1 : mr)+abs(square(a)+square(b)-square(c))*(rb==I(1) ? m1 : mr)+abs(square(a)+square(c)-square(b))*(rc==I(1) ? m1 : mr);
    return lower(de)>=upper(du);
}

bool test_eps(I eps)
{
    return (check(I(1),I(1),I(1),eps) && check(I(1),I(1),r,eps) && check(I(1),r,r,eps) && check(r,I(1),I(1),eps) && check(r,I(1),r,eps) && check(r,r,r,eps));
}

double find_eps()
{
    I eps=I(0,1);
    std::pair<I,I> e;
    while (width(eps)>0.000001)
    {
//        printf("%f,%f\n",lower(eps),upper(eps));
        e=bisect(eps);
        if (test_eps(e.first)) eps=e.second;
        else eps=e.first;
    }
    return upper(eps);
}


//***************************************************//

void update_mZ(char q)
{
    std::string t[3]={"11","1r","rr"};
    // u: lower bound on angles 1q1, 1qr and rqr in any FM-triangle
    I* u=(I*)calloc(3,sizeof(I));
    for (int i=0;i<3;i++)
        u[i]=lower(asin(min(val(t[i][0])/(val(q)+val(t[i][0])+I(2)*r),val(t[i][1])/(val(q)+val(t[i][1])+I(2)*r))));
    // a: angles 1q1, 1qr and rqr in a tight triangle
    I* a=(I*)calloc(3,sizeof(I));
    for (int i=0;i<3;i++)
        a[i]=angle(val(t[i][0])+val(t[i][1]),val(q)+val(t[i][0]),val(q)+val(t[i][1]));
    // ranges through coronas k such that sum of angle may be less or equal 2*pi
    int k[3]={0,0,0}; //int* k=(int*)calloc(3,sizeof(int));
    I z=0; // sum of angles in corona k
    int i=0; // index of the "digit" of k to increment
    I mq=I(0);
    while (i<3)
    {
        // increment corona k
        i=0;
        while (i<3 && z+u[i]>pi_twice<I>())
        {
            z-=I(k[i])*u[i];
            k[i]=0;
            i+=1;
        }
        if (i==3) break; // all the coronas have been considered -> finished
        k[i]+=1;
        z+=u[i];
        // skip non-coronable k
        if (k[1]%2!=0 || (k[0]*k[2]!=0 && k[1]==0) || (k[0]+k[1]+k[2]<3)) continue;
        // skip corona k such that the vertex inequality is ensured by definition of potentials
        if (q=='r' && k[0]==4 && k[1]==0 && k[2]==0) continue; // r-corona 1111
        if (q=='r' && k[0]==0 && k[1]==0 && k[2]==6) continue; // r-corona rrrrrr
        if (q=='1' && k[0]==3 && k[1]==4 && k[2]==0) continue; // 1-corona 1111r1r
        if (q=='1' && k[0]==6 && k[1]==0 && k[2]==0) continue; // 1-corona 111111
        if (q=='1' && k[0]==0 && k[1]==8 && k[2]==0) continue; // 1-corona 1r1r1r1r
        // compute the coefficient of mq
        I deno=abs(pi_twice<I>()-I(k[0])*a[0]-I(k[1])*a[1]-I(k[2])*a[2]);
        // should not happen (skipped above)
        if (zero_in(deno)) throw std::string("ERROR!");
        // compute the sum of vertex base potentials
        I s= (q=='1' ? a1-I(k[0])*V111-I(k[1])*V11r-I(k[2])*Vr1r : ar-I(k[0])*V1r1-I(k[1])*V1rr-I(k[2])*Vrrr);
        // updates the lower bound on mq
        mq=max(mq,s/deno);
    }
    // ceiling Zq
    I Zq=I(INFINITY);
    for (int i=0;i<3;i++)
    {
        mq=max(mq,-Vijk(t[i][0],q,t[i][1])/a[i]);
        Zq=min(Zq,Vijk(t[i][0],q,t[i][1])/a[i]);
    }
    Zq=-pi_twice<I>()*Zq;
    // update m1/Z1 or mr/Zr
    if (q=='1') {m1=upper(mq); Z1=upper(Zq);}
    else  {mr=upper(mq); Zr=upper(Zq);}
    free(u);
    free(a);
}

//***************************************************//
// update variables when x changes

void update()
{
    if (lower(x)<0.5) // more small discs
    {
        //update target density (formula modified to have only one occurence of x, hence small interval)
        delta=pi<I>()*(I(1)-square(r))/(I(4)*(I(1)-square(r)*sqrt(I(3))))*(I(1)+(I(4)*square(r)*(I(1)-square(r)*sqrt(I(3)))/(I(1)-square(r))-I(2)*square(r)*sqrt(I(3)))/(I(4)*x*(I(1)-square(r)*sqrt(I(3)))+I(2)*square(r)*sqrt(I(3))));
        // update base vertex potentials which depends on p
        V1r1=excess(I(2)*r,I(2)*r,I(2)*r,r,r,r)/I(2);
        // arbitrary "magic" value of V1rr
        V1rr=((I(86)/I(5269)*x - I(98)/I(5063))*x + I(10)/I(1139))*x - I(4)/I(2831);
        V1rr=I(median(V1rr));
    }
    else // more large discs
    {
        //update target density (formula modified to have only one occurence of x, hence small interval)
        delta=pi<I>()*(I(1)-square(r))/(I(4)*(sqrt(I(3))-I(1)))*(I(1)+(square(r)*I(4)*(sqrt(I(3))-I(1))/(I(1)-square(r))-(I(4)-I(2)*sqrt(I(3))))/(I(4)*x*(sqrt(I(3))-I(1))+I(4)-I(2)*sqrt(I(3))));
        // update base vertex potentials which depends on p
        V1r1=excess(I(2),I(1)+r,I(1)+r,r,I(1),I(1))-excess(I(2),I(2),I(2),I(1),I(1),I(1))/I(2);
        // arbitrary "magic" value of V1rr
        V1rr=-I(6)/I(1000);
    }
    // update remaining base vertex potential (to sum up to excess in each of the 4 tight triangles)
    Vrrr=excess(I(2)*r,I(2)*r,I(2)*r,r,r,r)/I(3);
    V111=excess(I(2),I(2),I(2),I(1),I(1),I(1))/I(3);
    Vr1r=excess(I(2)*r,I(1)+r,I(1)+r,I(1),r,r)-I(2)*V1rr;
    V11r=(excess(I(1)+r,I(1)+r,I(2),I(1),I(1),r)-V1r1)/I(2);
    // update lower bound on potential sums around a vertex
    a1=I(8)*V11r;
    ar=I(4)*V1r1;
    // add eta=0.0001 to a1 and ar in order to bound the proportions of bad neighborhoods (Prop. 2)
    // uncomment the two following lines to check Prop. 2, comment them to check Th. 1
    a1+=I(1)/I(10000);
    ar+=I(1)/I(10000);
    // update m1/Z1 and mr/Zr
    update_mZ('1');
    update_mZ('r');
    // find eps such that local inequality is ensured on eps-tight triangles
    eps=find_eps();
}

//***************************************************//
// potential of a triangle

// parameters for edge potential (depends on p)
std::pair <double,double> LQ(char r1, char r2)
{
    if (r1=='1' && r2=='1') // edge 11
        return (lower(x)<0.5 ? std::pair<double,double>(2.5,0.38) : std::pair<double,double>(2.5,0.02));
    else if (r1!=r2)    // edge 1r
        return (lower(x)<0.5 ? std::pair<double,double>(1.83,0.15) : std::pair<double,double>(1.83,0.05));
    else // edge rr
        return (lower(x)<0.5 ? std::pair<double,double>(1.18,0.15) : std::pair<double,double>(1.18,0.08));
}
    
// potential
I U(I a, I b, I c, char ra, char rb, char rc, I rs, I S)
{
    I Ra=val(ra);
    I Rb=val(rb);
    I Rc=val(rc);
    // u<-vertex potential
    I u=min(I(ra=='1'?Z1:Zr),Vijk(ra,rb,rc)+I(ra=='1'?m1:mr)*abs(angle(a,b,c)-angle(Rb+Rc,Ra+Rc,Ra+Rb)));
     u+=min(I(rb=='1'?Z1:Zr),Vijk(rb,rc,ra)+I(rb=='1'?m1:mr)*abs(angle(b,c,a)-angle(Ra+Rc,Ra+Rb,Rb+Rc)));
     u+=min(I(rc=='1'?Z1:Zr),Vijk(rc,ra,rb)+I(rc=='1'?m1:mr)*abs(angle(c,a,b)-angle(Ra+Rb,Rb+Rc,Ra+Rc)));
    // add edge potential to u
    std::pair <double,double> lq=LQ(rb,rc);
    if (upper(a)>lq.first) u+=lq.second*dist(a,b,c,Ra,Rb,Rc,rs,S);
    lq=LQ(ra,rc);
    if (upper(b)>lq.first) u+=lq.second*dist(b,c,a,Rb,Rc,Ra,rs,S);
    lq=LQ(ra,rb);
    if (upper(c)>lq.first) u+=lq.second*dist(c,a,b,Rc,Ra,Rb,rs,S);
    return u;
}

//***************************************************//

int dicho(I a, I b, I c, char ra, char rb, char rc)
{
    I Ra=val(ra);
    I Rb=val(rb);
    I Rc=val(rc);

    // test if it is epsilon tight
    if (upper(a-Rb-Rc)<eps and upper(b-Ra-Rc)<eps and upper(c-Ra-Rb)<eps) return ++nt;

    // test for feasability in an FM-triangulation
    // the altitude of each vertex must be at least r (no disc crosses the opposite edge)
    if (upper(square(b)-square((square(b)+square(c)-square(a))/(I(2)*c)))<lower(square(r))) return ++nb;
    if (upper(square(c)-square((square(c)+square(a)-square(b))/(I(2)*a)))<lower(square(r))) return ++nb;
    if (upper(square(a)-square((square(a)+square(b)-square(c))/(I(2)*b)))<lower(square(r))) return ++nb;
    // the area S must not be too small
    I S=(a+b+c)*(a+b-c)*(a-b+c)*(-a+b+c);
    if (upper(S)<lower(square(pi_twice<I>()*square(r)))) return ++nb;
    // the radius of support circle must be at most r
    I rs=radius(a,b,c,Ra,Rb,Rc);
    if (lower(rs)>upper(r)) return ++nb;
    // one can assume rs included in [0,r] (otherwise will be only later detected as not feasible)
    rs=max(I(0),min(rs,r));

    // compare excess and potential
    I e=excess(a,b,c,Ra,Rb,Rc);
    I u=U(a,b,c,ra,rb,rc,rs,S);
    if (lower(e)>=upper(u)) return ++ng; // good triangles!
    if (upper(e)<lower(u)) // should never happen // and upper(rs)<=lower(r))
    {
        printf("erreur\n");
        printf("a=RIF([%.20f, %.20f])\n",lower(a),upper(a));
        printf("b=RIF([%.20f, %.20f])\n",lower(b),upper(b));
        printf("c=RIF([%.20f, %.20f])\n",lower(c),upper(c));
        printf("ra=RIF([%.20f, %.20f])\n",lower(Ra),upper(Ra));
        printf("rb=RIF([%.20f, %.20f])\n",lower(Rb),upper(Rb));
        printf("rc=RIF([%.20f, %.20f])\n",lower(Rc),upper(Rc));
        printf("rs=[%.20f, %.20f]\n",lower(rs),upper(rs));
        printf("E=[%.20f, %.20f]\n",lower(e),upper(e));
        printf("U=[%.20f, %.20f]\n",lower(u),upper(u));
        throw std::string("ERROR!");
    }
    // precision does not suffice to compare excess and potential -> refine
    std::pair<I,I> aa=bisect(a);
    std::pair<I,I> bb=bisect(b);
    std::pair<I,I> cc=bisect(c);
    dicho(aa.first ,bb.first ,cc.first ,ra,rb,rc);
    dicho(aa.first ,bb.first ,cc.second,ra,rb,rc);
    dicho(aa.first ,bb.second,cc.first ,ra,rb,rc);
    dicho(aa.first ,bb.second,cc.second,ra,rb,rc);
    dicho(aa.second,bb.first ,cc.first ,ra,rb,rc);
    dicho(aa.second,bb.first ,cc.second,ra,rb,rc);
    dicho(aa.second,bb.second,cc.first ,ra,rb,rc);
    dicho(aa.second,bb.second,cc.second,ra,rb,rc);
}

//***************************************************//
int main()
{
    std::string T[4]={"111","11r","1rr","rrr"};
    I a,b,c;
    char ra,rb,rc;
    int n;

    printf("low(p) triangles\n");
    for(int j=0;j<100;j++)
    {
        x=I(j,j+1)/I(100); // proportion of small discs (interval)
        update();
        printf("%.4f ",lower(x));
//        printf("%.15f %.15f %.15f %.15f\n",median(V1rr),m1,mr,eps);
//        printf("%.15f %.15f\n",median(a1),median(ar));
        fflush(stdout);
        n=0;
        for (int i=0;i<4;i++)
        {
            ra=T[i][0];
            rb=T[i][1];
            rc=T[i][2];
            a=I(lower(val(rb)+val(rc)),upper(val(rb)+val(rc)+I(2)*r));
            b=I(lower(val(ra)+val(rc)),upper(val(ra)+val(rc)+I(2)*r));
            c=I(lower(val(rb)+val(ra)),upper(val(rb)+val(ra)+I(2)*r));
            nt=0;
            nb=0;
            ng=0;
            dicho(a,b,c,ra,rb,rc);
//            printf("%d %d %d ",ng,nt,nb);fflush(stdout);  // nb good/tight/bad triangles of each type
            n+=nt;
            n+=nb;
            n+=ng;
        }
        printf("%d\n",n); // total number of triangles
    }

    return 0;
}