/**************************************************************
          *                                                            *
          *     cycles.cpp:  contains procedures which find a set      *
          *     of generating cycles for the homology of the Seifert   *
          *     surface obtained by Seifert's algorithm, and then      *
          *     compute the linking form                               *
          *                                                            *
          **************************************************************/



#include "link_diagram.h"
#include <algorithm>



void Link::tree()
{
    int i, j;

    /*
     *  compute homology cycles for Seifert surface; begin by constructing maximal tree.
     *  Much use is made of a structure QCQ-triple (quadrant-circuit-quadrant):
     *  these provide the edges of the graph whose vertices are Seifert circuits
     *  and whose edges are crossings of the diagram.
     *
     *  We deal with circuit #0 (the root circuit) first.
     */
    
    Seifert_circuit *s0 = seifert_circuit_list[0];
    QCQ_triple *qc = new QCQ_triple;
    qc->circuit = s0;
    qc->forward_quadrant = NULL;
    qc->backward_quadrant = NULL;
    qc->next = NULL;
    qc->prev = NULL;
    s0->path_from_root.push_back( qc );
    
    /*
     *  now find standard path from root circuit to each other circuit
     */    
    
    for ( int ccirc=1; ccirc<num_seifert_circuits; ++ccirc )
    {

        Seifert_circuit *s0 = seifert_circuit_list[ccirc];
        Seifert_circuit *s = s0;
        
        while ( s->index != 0 )
        {

            QCQ_triple *qc = new QCQ_triple;
            qc->forward_quadrant = s->previous_quadrant;
            qc->circuit = s->previous_circuit;
            
            s0->path_from_root.push_back( qc );
           
            s = s->previous_circuit;

        }
        
        /*
         *  We want to reverse the path_from_root vector so that it starts at circuit #0.
         *  Also we need to add an entry for the (target) circuit s0, and enter the
         *  backward_quadrants, so as to get a doubly linked list.  In a sense,
         *  the union of the histories is a maximal tree for the graph whose 
         *  vertices are circuits and whose edges are crossings of the diagram.
         */
        
        reverse( s0->path_from_root.begin(), s0->path_from_root.end() );
        
        s0->path_from_root[0]->backward_quadrant = NULL;  //  can't go back from the beginning
        
        for ( i=1; i<s0->path_from_root.size(); ++i )
            s0->path_from_root[i]->backward_quadrant = s0->path_from_root[i-1]->forward_quadrant->clockwise->clockwise;
        
        QCQ_triple *qc = new QCQ_triple;
        
        qc->circuit = s0;
        qc->forward_quadrant = NULL;      //  can't go forward from the end
        qc->backward_quadrant = s0->path_from_root.back()->forward_quadrant->clockwise->clockwise;
        
        s0->path_from_root.push_back( qc );
        
        /*
         *  make doubly linked list of QCQ_triples
         */

        s0->path_from_root.front()->prev = NULL;
        s0->path_from_root.back()->next = NULL;
        
        for ( i=0; i<s0->path_from_root.size()-1; ++i )
            s0->path_from_root[i]->next = s0->path_from_root[i+1];
        for ( i=1; i<s0->path_from_root.size(); ++i )
            s0->path_from_root[i]->prev = s0->path_from_root[i-1];
                    
    }

/*
 * (diagnostic)
    for ( i=1; i<num_seifert_circuits; ++i )
    {
        Seifert_circuit *s = seifert_circuit_list[i];
        
        for ( j=0; j<s->path_from_root.size(); ++j )
        {
            printf( "circ=%3d  ", s->path_from_root[j]->circuit->index );
            if ( s->path_from_root[j]->forward_quadrant != NULL )
                printf( "forward=%3d", s->path_from_root[j]->forward_quadrant->crossing->index );
            if ( s->path_from_root[j]->backward_quadrant != NULL )
                printf( "backward=%3d", s->path_from_root[j]->backward_quadrant->crossing->index );
            printf( "\n" );
        }
    }
*/    


}


void Link::find_cycles()
{
    int i, j, t;
    Quadrant *q1, *q2;

    /*
     *  get one cycle for each crossing not in tree; a cycle is a doubly linked
     *  circular list of QCQ-triples.
     */
    
    for ( int ccross=1; ccross<=crossing_number; ++ccross )
    {
        Crossing *c = crossingList[ccross];
        
        if ( c->nw->tree || c->ne->tree ) continue;  //  we only want edges not in the tree
        
        
        q1 = ( c->nw->circuit_index == -1 ) ? c->ne : c->nw;
        Seifert_circuit *s1 = seifert_circuit_list[q1->circuit_index];
        q2 = q1->clockwise->clockwise;
        Seifert_circuit *s2 = seifert_circuit_list[q2->circuit_index];
        
        /*
         *  Determine the point in the tree where the standard paths to s1, s2 diverge
         */        
        
        t = 0; 
        bool try_further = true;   
/*
printf( "circuit %3d path:  ", s1->index );
for ( i=0; i<s1->path_from_root.size(); ++i ) printf( "%3d", s1->path_from_root[i]->circuit->index );
printf( "\n" );
printf( "circuit %3d path:  ", s2->index );
for ( i=0; i<s2->path_from_root.size(); ++i ) printf( "%3d", s2->path_from_root[i]->circuit->index );
printf( "\n" );
*/
        while ( try_further )
        {
            if ( s1->path_from_root[t]->next == NULL || s2->path_from_root[t]->next == NULL )
            {
                try_further = false;   //  we've got to the end of a branch
                break;
            }
            
            if ( s1->path_from_root[t]->next->circuit->index != s2->path_from_root[t]->next ->circuit->index )
            {
                try_further = false;   //  parting of the ways
                break;
            }
            
            ++t;
        }
            
        
        
        int branch_point = t;
//printf( "branch: %3d\n", s1->path_from_root[t]->circuit->index );           
        
        QCQ_triple *qcx1 = s1->path_from_root[branch_point];
        QCQ_triple *qcx2 = s2->path_from_root[branch_point];

        /*
         *  cycle is in three parts: (i) go through current crossing from s1 to s2;
         *  (ii) go back down tree from s2 to t; (iii) go forwards along tree from t to s1
         */
        
        std::vector<QCQ_triple *> vqc;
        
        QCQ_triple *qc1 = new QCQ_triple;
        
        if ( s1->path_from_root.back() != qcx1 )
        {
            qc1->backward_quadrant = s1->path_from_root.back()->backward_quadrant;
        }
        else
        {
            qc1->backward_quadrant = qcx2->forward_quadrant;
        }
        qc1->circuit = s1;
        qc1->forward_quadrant = q1;

        vqc.push_back( qc1 );
        
        QCQ_triple *qc2 = new QCQ_triple;
        
        qc2->backward_quadrant =  q2;
        qc2->circuit = s2;
        
        if ( s2->path_from_root.back() != qcx2 )
        {
            qc2->forward_quadrant = s2->path_from_root.back()->backward_quadrant;
        }
        else
        {
            qc2->forward_quadrant = qcx1->forward_quadrant;
        }
        
       vqc.push_back( qc2 );
        
        /*
         *  the next part is the interior of the backwards path from s2 to the
         *  branch point, so skip if this part is empty
         */
        
        if ( branch_point < (int) s2->path_from_root.size() - 2 )
            
        {

            for ( QCQ_triple *qc = s2->path_from_root.back()->prev;
                              qc != s2->path_from_root[branch_point]; qc=qc->prev )
            {
                QCQ_triple *nqc = new QCQ_triple;

                nqc->backward_quadrant = qc->forward_quadrant;
                nqc->circuit = qc->circuit;
                nqc->forward_quadrant = qc->backward_quadrant;
 
                vqc.push_back( nqc );
            }
        }
        
        /*
         *  next deal with branch point
         */
        
        QCQ_triple *bqc1 = s1->path_from_root[branch_point];
        QCQ_triple *bqc2 = s2->path_from_root[branch_point];
        
        if ( bqc1->circuit != s1 && bqc2->circuit != s2 )
        {
            QCQ_triple *nqc = new QCQ_triple;
 
            nqc->backward_quadrant = bqc2->forward_quadrant;
            nqc->circuit = bqc1->circuit;
            nqc->forward_quadrant = bqc1->forward_quadrant;

            vqc.push_back( nqc );
        }
        
        /*
         *  finally interior part of forwards path along tree from branch point to s1
         */


        if ( branch_point < (int) s1->path_from_root.size() - 2 )
            
        {

            for ( QCQ_triple *qc = s1->path_from_root[branch_point]->next;
                              qc != s1->path_from_root.back(); qc=qc->next )
            {
                QCQ_triple *nqc = new QCQ_triple;

                nqc->backward_quadrant = qc->backward_quadrant;
                nqc->circuit = qc->circuit;
                nqc->forward_quadrant = qc->forward_quadrant;
 
                vqc.push_back( nqc );
            }
        }
/*
printf( "\n" );
printf( "cycle: " );
for ( i=0; i<vqc.size(); ++i ) printf( "%3d%3d%3d   ", vqc[i]->backward_quadrant->crossing->index,
                                                       vqc[i]->circuit->index,
                                                       vqc[i]->forward_quadrant->crossing->index );
printf( "\n\n" );
*/    
        
        for ( i=0; i<vqc.size()-1; ++i ) vqc[i]->next = vqc[i+1];
        vqc.back()->next = vqc.front();
        

        seifert_cycle_list.push_back( vqc );
        
        vqc.clear();
    }



}



void Link::compute_linking()
{

    int i, c1, c2, ccirc, lnk, lnkt, t, dir, pos1, pos2, h1, h2, h;
    Quadrant *q11, *q12, *q21, *q22;
    
    for ( c1=0; c1<seifert_cycle_list.size(); ++c1 ) for ( c2=0; c2<seifert_cycle_list.size(); ++c2 )
    {
        std::vector<QCQ_triple *> qc1 = seifert_cycle_list[c1];
        std::vector<QCQ_triple *> qc2 = seifert_cycle_list[c2];
        
        lnk = 0;
        
        /*
         *  For each circuit, see how the cycles qc1, qc2 enter and leave it.  We're only
         *  interested in circuits visited by both cycles.  There are four cases:
         *
         *  (i)    qc1, qc2  share both entry and exit
         *  (ii)   they share one entry/exit point
         *  (iii)  they don't share, but respective entry/exit points are interlaced
         *  (iv)   ditto, but not interlaced (in which case there is no contribution
         *         to linking).  Note that that the choice of generating cycles
         *         guarantees that in cases (iii), (iv) the cycles can only meet at
         *         one vertex (circuit).
         */
        
        for ( ccirc=0; ccirc<seifert_circuit_list.size(); ++ccirc )
        {
            Seifert_circuit *s = seifert_circuit_list[ccirc];
            int sns = s->sense;
//printf( "crossings for circuit %3d:   ", ccirc );
//for ( i=0; i<s->length; ++i ) printf( "%3d", s->quadrant_list[i]->crossing->index );
//printf( "\n" );            
            /*
             *  find number of shared entry/exit points
             */
            
            Quadrant *q11, *q12;
            
            t = 0;
            
            for ( i=0; i<qc1.size(); ++i ) if ( qc1[i]->circuit == s )
            {
                t = 1;
                pos1 = i;
                
                q11 = qc1[i]->backward_quadrant;
                q12 = qc1[i]->forward_quadrant;
            }
            
            if ( t == 0 ) continue;
            
            t = 0;
            
            for ( i=0; i<qc2.size(); ++i ) if ( qc2[i]->circuit == s )
            {
                t = 1;
                pos2 = i;
                
                q21 = qc2[i]->backward_quadrant;
                q22 = qc2[i]->forward_quadrant;
            }
            
            if ( t == 0 ) continue;
//printf( "we are at:  c1=%3d  c2=%3d  ccirc=%3d  sense=%3d ", c1, c2, ccirc, seifert_circuit_list[ccirc]->sense );
//printf( "crossings: %3d%3d%3d%3d\n", q11->crossing->index, q12->crossing->index, q21->crossing->index, q22->crossing->index );           

//printf( "t11=%3d  t12=%3d  t21=%3d  t22=%3d   sns=%3d\n", q11->twist, q12->twist, q21->twist, q22->twist, sns );


            if ( q11 == q22 || q12 == q21 )
            {
                dir = -1;            //  for convenience, "reverse" one cycle, i.e.
                
                Quadrant *q = q21;   //  swap entry and exit points of second cycle
                q21 = q22;
                q22 = q;
            }
            else
            {
                dir = 1;
            }
                
            if ( q11 == q21 && q12 == q22 )    //  both shared
            {
                lnkt = ( q11->twist + q12->twist ) * dir;                
            }
            
            else if ( q11 == q21 )
            {
                int t =  ( cyclic_order( q12, q11, q22 ) * q11->twist == -1 ) ? 0 : 2*q11->twist;
                lnkt = t*dir;
            }
            
            else if ( q12 == q22 )
            {
                int t =  ( cyclic_order( q11, q12, q21 ) * q12->twist == -1 ) ? 0 : 2*q12->twist;
                lnkt = t*dir;
            }
            
            else
            {
                int it = interlace_test( q11, q12, q21, q22 );
                
                if ( it == 0 )
                {
                    lnkt = 0;
                }
                else
                {
                
                    
                    h1 = h2 = 0;
                    
                    QCQ_triple *qc = qc1[pos1];
 
                    h1 = seifert_contig[qc->circuit->index][qc->next->circuit->index];
                    while ( seifert_contig[qc->circuit->index][qc->next->circuit->index] == same_level
                            && qc->next != qc1[pos1] )
                    {

                        h1 = seifert_contig[qc->circuit->index][qc->next->circuit->index];
                        qc = qc->next;
                    }
 
                    qc = qc2[pos2];
 
                    h2 = seifert_contig[qc->circuit->index][qc->next->circuit->index];
                    while ( seifert_contig[qc->circuit->index][qc->next->circuit->index] == same_level
                            && qc->next != qc2[pos2] )
                    {
                        h2 = seifert_contig[qc->circuit->index][qc->next->circuit->index];
                        qc = qc->next;
                    }
 
                    if ( h1 == 0 && h2 == 0 )
                    {
                        printf( "error in interlace()\n" );
                        exit(0);
                    }
 
 
                    if ( sns == 1 )                              //  cycle qc2 is pushed up
                         h = ( h1 < 0 || h2 > 0 ) ? 0 : 1;
 
                    if ( sns == -1 )                             //  cycle qc2 is pushed down
                         h = ( h1 < 0 || h2 > 0 ) ? -1 : 0;
 
                    lnkt = 4*h*it;
                }

            }

            lnk += lnkt;
            
        }
        
        seifert_matrix[c1][c2] += lnk;
        
        qc1.clear();
        qc2.clear();
    }

    for ( c1=0; c1<seifert_cycle_list.size(); ++c1 ) for ( c2=0; c2<seifert_cycle_list.size(); ++c2 )
        if ( seifert_matrix[c1][c2] % 4 != 0 )
    {
        printf( "not multiple of 4\n" );
        exit(0);
    }
    
    for ( c1=0; c1<seifert_cycle_list.size(); ++c1 ) for ( c2=0; c2<seifert_cycle_list.size(); ++c2 )
        seifert_matrix[c1][c2] /= 4;
    

}
            


int Link::cyclic_order( Quadrant *q1, Quadrant *q2, Quadrant *q3 )
{
    Quadrant *q = q1;
    
    while ( q != q2 && q != q3 ) q = q->seifert_next;
    
    if ( q == q2 ) return 1;
    else return -1;
}



int Link::interlace_test( Quadrant *q1, Quadrant *q2, Quadrant *q3, Quadrant *q4 )
{
    int t1 = cyclic_order( q1, q2, q3 );
    int t2 = cyclic_order( q1, q2, q4 );
    int t3 = cyclic_order( q1, q3, q4 );

    if ( t1 == -1 && t2 == 1 && t3 == 1 ) return 1;
    else if ( t1 == 1 && t2 == -1 && t3 == -1 ) return -1;
    else return 0;
}