/*
 * The Npic library
 *
 * Copyright (C) 2003 Edouard Thiel <Edouard.Thiel@lif.univ-mrs.fr>
 *
 * This library is free software under the terms of the 
 * GNU Lesser General Public License (LGPL) version 2.1.
*/

/*
 * mask_gsym.c - 22/01/2004
 *
 * Computation of all G-symmetries for generating a half or full mask.
 *
 * References:
 *
 * [1] S.M. Johnson, 1963. "Generation of Permutations by Adjacent Transpositions".
 *     Math. Comput. 17, pp. 282-285.
 *
 * [2] R. Seroul, 1995. "Informatique pour mathematiciens". InterEditions,
 *     ISSN 978-2-7296-0568-1.
 *
 * [3] R. Seroul, 2000. "Permutations: Johnson's Algorithm". Programming 
 *     for Mathematicians, Berlin, Springer-Verlag, pp. 213-218, section 8.15,
 *     ISSN 978-3-540-66422-2.
*/


#include <npic.h>


/*----------------------------------------------------------------------------*/

/* 
 * Johnson's algorithm to compute all permutations in 1..dim
 * see [1][2][3].
 * This function also check if (x[1], .., x[dim]) is unique.
*/

void NpicGsymPermuInit (Npic_gsym *gs)
{
    int i;
    
    /* Initialize Johnson's algorithm */
    gs->pos[0] = gs->pos[gs->dim+1] = gs->dim+1;
    for (i = 1; i <= gs->dim; i++) { 
        gs->dir[i] = -1; gs->pos[i] = i; 
    }
}

int NpicGsymPermuNext (Npic_gsym *gs)
{
    int i, j, m, t1, t2, tmp, ok;
    
    do {  
        m = 0; t1 = 0;
    
        /* Search m */
        for (i = 1; i <= gs->dim; i++)
            if (gs->pos[i] > m &&
                gs->pos[i] > gs->pos[i + gs->dir[i]])
        {
            m = gs->pos[i];
            t1 = i;
        }
    
        /* No more permutation */
        if (m == 0) return 0;

        /* do 1 permutation */
        t2 = t1 + gs->dir[t1];
        tmp = gs->pos[t1]; gs->pos[t1] = gs->pos[t2]; gs->pos[t2] = tmp;
        tmp = gs->dir[t1]; gs->dir[t1] = gs->dir[t2]; gs->dir[t2] = tmp;
        tmp = gs->x[t1]; gs->x[t1] = gs->x[t2]; gs->x[t2] = tmp;
    
        /* change directions */
        if (m < gs->dim)
        for (i = 1; i <= gs->dim; i++)
            if (gs->pos[i] > m) gs->dir[i] = -gs->dir[i];
        
        /* Test if a distinct (x[1], .., x[dim]) was found.
         * The idea is to verify ascending order on pos[] 
         * for the x[] which are equal, because the solution is
         * always unique.
        */
        ok = 1;
        for (i = 1; i <= gs->dim; i++)
        for (j = i+1; j <= gs->dim; j++)
            if (gs->x[i] == gs->x[j] && gs->pos[i] > gs->pos[j]) 
            ok = 0;
        
    } while (ok == 0);
    
    return 1;
}


/* 
 * Compute all the distinct sign combinations of (x[1], .., x[dim])
 * which are in the half-mask in scan line order (when half==1),
 * or in the full-mask (when half==0).
 *
 * The scan line order is : for t[dim], for t[dim-1], .., for t[1].
 * A point (x[1], .., x[j], 0 ... 0) belongs to the half mask iff the
 * last non-null coordinate is positive, i.e x[j] > 0.
*/

void NpicGsymCombiInit (Npic_gsym *gs)
{
    int i, j;
    
    /* Initialize sign combinations */
    gs->k = gs->kmax = 1 << gs->dim;
    
    /* Compute binary mask to forbid sign combinations 
     * - for x[i] == 0 (so as to have distinct points)
     * - for x[j] where j is such that x[j] is the last non null x[]
     *   (separates half masks if wanted) in scan line order.
    */
    gs->bin_mask = 0; j = 0;
    for (i = 1; i <= gs->dim; i++)
        if (gs->x[i] == 0) 
             gs->bin_mask |= 1 << (i-1);
        else j = i;
    if (gs->half && j > 0) gs->bin_mask |= 1 << (j-1);
}

int NpicGsymCombiNext (Npic_gsym *gs)
{
    int i;
    while (gs->k > 0)
    {
        gs->k--; 
        
        /* Check if gs->k is not forbidden 
         * Remark: 0 is always allowed
        */
        if ((gs->k & gs->bin_mask) == 0)
        {
            /* Do 1 sign combination */
            for (i = 1; i <= gs->dim; i++)
                if (gs->k & (1 << (i-1)))
                     gs->x[i] = -abs(gs->x[i]);
                else gs->x[i] =  abs(gs->x[i]);
            return 1;
        }
    }
    return 0;
}


/*
 * Compute all the distinct G-symmetries of (x[1], .., x[dim]) 
 * in the half (or full) mask.
 * The G-symmetries are all the permutations and sign combinations.
 *
 * USAGE:
 *    Npic_gsym gs;
 *    gs.x[1] = 1; gs.x[2] = 3; gs.x[3] = 0; gs.dim = 3; gs.half = 1;
 *    NpicGsymInit (&gs);
 *    while (NpicGsymNext (&gs))
 *        NpicGsymPrint (&gs);
*/

void NpicGsymInit (Npic_gsym *gs)
{
    NpicGsymCombiInit (gs);
    NpicGsymPermuInit (gs);
}

int NpicGsymNext (Npic_gsym *gs)
{
    if (gs->k == 0)
    {
        if (NpicGsymPermuNext (gs) == 0) return 0;
        NpicGsymCombiInit (gs);
    }   
    
    return NpicGsymCombiNext (gs);
}    
    


/* 
 * Print positions pos[], directions dir[] and coordinates x[]
*/

void NpicGsymPrint (Npic_gsym *gs)
{
    int i;
    for (i = 1; i <= gs->dim; i++)
        printf ("%3d%c ", gs->pos[i], gs->dir[i] == 1 ? '>':'<');
    printf ("   ");
    for (i = 1; i <= gs->dim; i++)
        printf ("%3d ", gs->x[i]);
    printf ("\n");
}


/*----------------------------------------------------------------------------*/