rote_generation.cpp

Go to the documentation of this file.

#include <iostream>
#include <cstdlib>
#include <CGAL/Cartesian.h>
#include <CGAL/Arr_segment_traits_2.h>
#include <CGAL/Arr_extended_dcel.h>
#include <CGAL/Arrangement_2.h>
#include <CGAL/Gmpz.h>
#include <CGAL/Gmpq.h>
#include <boost/math/common_factor_rt.hpp>
#include <stack>

using namespace std;

int m, n;

struct VertexData {};

struct EdgeData {
  int i, j;
  EdgeData() : i(0), j(0) {};
  EdgeData(const int& i_, const int& j_) : i(i_), j(j_) {};
};

struct FaceData {
  bool discovered; 
  char** matrix; 


  FaceData() : discovered(false) {
    typedef char* char_ptr;
    matrix = new char_ptr [m];
    matrix[0] = new char[m*n];
    for(int i = 1; i < m; ++i)
      matrix[i] = matrix[i-1] + n;
    for(int i = 0; i < m; i++)
      for(int j = 0; j < n; ++j)
    matrix[i][j] = 0;
  }
  void delete_face() {
    delete [] matrix[0];
    delete [] matrix;
  }
};

ostream& operator<< (ostream& stream, char** matrix) {
  // assuming the matrix has size m x n
  for(int i = 0; i < m; ++i) {
    for(int j = 0; j < n; ++j)
      stream << "[" << (int) matrix[i][j] << "]";
    stream << endl;
  }
  return stream;
}

typedef CGAL::Gmpz Integer;

typedef CGAL::Gmpq rat;

typedef CGAL::Cartesian<rat> Kernel;
typedef CGAL::Arr_segment_traits_2<Kernel> Traits_2;
typedef Traits_2::Point_2 Point_2;
typedef Traits_2::X_monotone_curve_2 Segment_2;

typedef CGAL::Arr_extended_dcel<Traits_2,
                VertexData,
                EdgeData,
                FaceData> Dcel;

typedef CGAL::Arrangement_2<Traits_2, Dcel> Arrangement_2;

typedef Arrangement_2::Face_handle Face_ptr;


void flip (int const& i, int const& j, char** matrix_src, char** matrix_dest) {
  // assuming i, j and k are in the correct interval, and matrices have
  // size m x n
  for(int r = 0; r < m; ++r)
    for(int s = 0; s < n; ++s)
      matrix_dest[r][s] = matrix_src[r][s];
  for (int r = int(1); r*i < int(m) && r*j < int(n); ++r)
    matrix_dest[r*i][r*j] = 1 - matrix_dest[r*i][r*j];
}

std::list <Face_ptr> dfs(Arrangement_2 const& graph, Face_ptr f) {
  std::list <Face_ptr> depth_first_result;
  std::stack <Face_ptr> vertices;

  vertices.push(f);
  f->data().discovered = true;
  while ( ! vertices.empty() ) {
    Face_ptr head = vertices.top();
    vertices.pop();
    
    depth_first_result.push_back(head);
    //    cout << head->data().matrix << endl;
    
    Arrangement_2::Ccb_halfedge_circulator first, curr, flip_op;
    first = curr = head->outer_ccb();
    do {
      Face_ptr son = curr->twin()->face();
      if( ! son->is_unbounded() && ! son->data().discovered) {
    vertices.push(son);
    son->data().discovered = true;
    flip(curr->data().i, curr->data().j, head->data().matrix, son->data().matrix);
      }
    } while(++curr != first);
  }
  
  return depth_first_result;
}

int main(int argc, char** argv) {
  if(argc != 3) cout << "usage : ./generator [m] [n]\nExit 1.\n", exit(1);

  string mm(*++argv); string nn(*++argv);
  (void) sscanf(nn.c_str(),"%d",&n);
  (void) sscanf(mm.c_str(),"%d",&m);
  
  cout << " -- " << m << "," << n << "-cubes generation --" << endl;

  Arrangement_2 arr;

  //general case
  for(int j(1); j < n; j++)
    for(int i(1); i < m; i++)
      for(int k(1); k < i+j; k++) {
    if(boost::math::gcd(j,boost::math::gcd(i,k)) == 1) {
      rat x0, y0, x1, y1;
      if(k > j) {
        y0 = rat(1);
        x0 = rat (k - j, i);
      }
      else if (k < j)
        x0 = rat(0), y0 = rat(k, j);
      else
        x0 = rat(0), y0 = rat(1);
      if(k < i) {
        y1 = rat(0);
        x1 = rat(k , i);
      }
      else if (k > i)
        x1 = rat(1), y1 = rat(k - i, j);
      else 
        x1 = rat(1), y1 = rat(0);
      
      insert(arr, Segment_2 (Point_2 (x0,y0), Point_2 (x1,y1)));
    }
      }


  //horizontal segments
  for(int j(1); j < n; j++)
    for(int k(1); k < j; k++)
      if(boost::math::gcd(j,k) == 1)
    insert(arr, Segment_2 (Point_2 (rat (0), rat (k, j)),
                   Point_2 (rat(1), rat(k, j))));

  insert(arr, Segment_2 (Point_2 (rat(0), rat(0)),
             Point_2 (rat(1), rat(0))));
  
  insert(arr, Segment_2 (Point_2 (rat(0), rat(1)),
             Point_2 (rat(1), rat(1)))); 


  //vectical segments
  for(int i(1); i < m; i++)
    for(int k(1); k < i; k++)
      if(boost::math::gcd(i,k) == 1)
    insert(arr, Segment_2 (Point_2 (rat(k, i), rat(0)),
                   Point_2 (rat(k, i), rat(1))));

  insert(arr, Segment_2 (Point_2 (rat(0), rat(0)),
             Point_2 (rat(0), rat(1))));

  insert(arr, Segment_2 (Point_2 (rat(1), rat(0)),
             Point_2 (rat(1), rat(1))));

  Arrangement_2::Face_handle start; //flat m,n-cube
  Arrangement_2::Edge_iterator eit;
  

  for(eit = arr.edges_begin(); eit != arr.edges_end(); ++eit) {
    Point_2 a = eit->source()->point(), b = eit->target()->point();

    if((a.x() == 0 && a.y() == 0) || (b.x() == 0 && b.y() == 0))
      if( ! eit->face()->is_unbounded())
    start = eit->face();
      else
    start = eit->twin()->face();

    if((a.x() - b.x()) != 0) { 
      rat alpha = (a.y() - b.y()) / (a.x() - b.x());
      rat beta = a.y() - alpha*a.x();

      Integer p = alpha.numerator(), q = alpha.denominator(),
    p_ = beta.numerator(), q_ = beta.denominator();
      Integer gcd_ = boost::math::gcd(q*q_, boost::math::gcd(-p*q_, p_*q));
      int i = (int) (-p*q_/gcd_).to_double(), j = (int) (q*q_/gcd_).to_double();
      eit->set_data(EdgeData(i,j));
      eit->twin()->set_data(EdgeData(i,j));
    }
    else {
      rat k = a.x();
      Integer p = k.numerator(), q = k.denominator();
      Integer gcd_ = boost::math::gcd(p,q);
      int i = (int) (q/gcd_).to_double();
      eit->set_data(EdgeData(i, 0));
      eit->twin()->set_data(EdgeData(i, 0));
    }
  }
 
  std::list <Face_ptr> rote_sequences = dfs(arr, start);
  for(list<Face_ptr>::iterator it = rote_sequences.begin();
      it != rote_sequences.end(); it++) {
    cout << (*it)->data().matrix << endl;
    (*it)->data().delete_face();
  }
  cout << "Number of m,n-cubes of plans through (0,0,0) : " << rote_sequences.size() << endl;

  return 0;
}