DiophantineSolver.h

Go to the documentation of this file.

/* This file is part of MODEL, the Mechanics Of Defect Evolution Library.
 *
 * Copyright (C) 2015 by Giacomo Po <gpo@ucla.edu>
 *
 * MODEL is distributed without any warranty under the
 * GNU General Public License (GPL) v2 <http://www.gnu.org/licenses/>.
 */
 
#ifndef gbLAB_DiophantineSolver_h_
#define gbLAB_DiophantineSolver_h_
 
#include <Eigen/Dense>
#include "IntegerMath.h"
#include <climits>
#include <cmath>
#include "../Utilities/sortIndices.h"
namespace oILAB {
template <typename IntScalarType> class DiophantineSolver {
public:
  // Finds such x and y, that a  x + b  y = gcd(a, b)
  // https://github.com/ADJA/algos/blob/master/NumberTheory/DiophantineEquation.cpp
  static IntScalarType extended_gcd(IntScalarType a, IntScalarType b,
                                    IntScalarType &x, IntScalarType &y) {
    if (b == 0) {
      x = 1;
      y = 0;
      return a;
    }
    IntScalarType x1, y1;
    IntScalarType g = extended_gcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return g;
  }
  /*
          // Find u such that a_1 u_1 + a_2 u_2 + .... + a_n u_n = 1
          // method 1 - does not work, the algorithm is incorrect
 
          // "A fast algorithm to find reduced hyperplane unit cells and solve
     N-dimensional Bézout's identities."
          //  Acta Crystallographica Section A: Foundations and Advances 77.5
     (2021).
 
          static Eigen::Vector<IntScalarType,Eigen::Dynamic> solveBezout(const
     Eigen::Vector<IntScalarType,Eigen::Dynamic>& p)
          {
              int n= p.size();
              Eigen::Vector<IntScalarType,Eigen::Dynamic> out(n);
              out= Eigen::Vector<IntScalarType,Eigen::Dynamic>::Zero(n);
 
              if (n==2)
              {
                  IntScalarType g= extended_gcd(p(0),p(1),out(0),out(1));
                  std::cout << p(0) << "   " << p(1) << std::endl;
                  if (g<0) out= -out;
                  return out;
              }
 
              int ind=0;
              for (auto pi: p)
              {
                  if (abs(pi)==1)
                  {
                      out(ind) = pi;
                      return out;
                  }
                  ind++;
              }
 
              Eigen::Vector<IntScalarType,Eigen::Dynamic> pS(n);
              ind= -1;
              for (auto i: sortIndices(p))
              {
                  ind++;
                  if (p(i) == 0) break;
                  pS(ind)= p(i);
              }
 
              IntScalarType pSi0= pS(ind);
              int numNonZeros= ind+1;
              Eigen::Vector<IntScalarType,Eigen::Dynamic>
     quotients(numNonZeros-1); Eigen::Vector<IntScalarType,Eigen::Dynamic>
     remainders(numNonZeros-1); for(int i=0; i<quotients.size(); i++)
              {
                  quotients(i)= pS(i)/pSi0;
                  remainders(i)= pS(i)%pSi0;
 
                  //if (pS(i)/pSi0 < 0)
                  //    quotients(i)= floor((double)pS(i)/pSi0);
                  //else
                  //    quotients(i)= pS(i)/pSi0;
                  //remainders(i)= pS(i)-pSi0*quotients(i);
              }
 
              Eigen::Vector<IntScalarType,Eigen::Dynamic> u(numNonZeros-1);
              u= solveBezout(remainders);
 
              out(Eigen::seq(0,numNonZeros-2))= u;
              out(numNonZeros-1)= -quotients.dot(u);
 
              Eigen::Vector<IntScalarType,Eigen::Dynamic> outRestored(n);
              ind= 0;
              for (auto i : sortIndices(p))
              {
                  outRestored(i) = out(ind);
                  ind++;
              }
              return outRestored;
          }
  */
  // Find u such that a_1 u_1 + a_2 u_2 + .... + a_n u_n = 1
  // Method 0:
  // "A fast algorithm to find reduced hyperplane unit cells and solve
  // N-dimensional Bézout's identities."
  //  Acta Crystallographica Section A: Foundations and Advances 77.5 (2021).
  template <typename ArrayType>
  static ArrayType solveBezout(const ArrayType &a) {
    int n = a.size();
    if (n < 2)
      throw std::runtime_error("the size of arrays should be at least two");
    if (abs(IntegerMath<IntScalarType>::gcd(a)) != 1)
      throw std::runtime_error("No solution since the gcd is not 1 or -1.");
 
    ArrayType u(n);
 
    IntScalarType p, q;
    IntScalarType g12 =
        extended_gcd(a(0), a(1), p, q); // g12 - gcd of a(0) and a(1)
 
    // form the n-1 sized vector na= {g,a2,...,a_{n-1}}
    Eigen::Vector<IntScalarType, Eigen::Dynamic> na(n - 1);
    na(0) = g12;
    na(Eigen::seq(1, n - 2)) = a(Eigen::seq(2, n - 1));
 
    Eigen::Vector<IntScalarType, Eigen::Dynamic> k(n - 1);
    if (n == 2) {
      IntScalarType g = extended_gcd(a(0), a(1), u(0), u(1));
      if (g < 0)
        u = -u;
      return u;
    } else {
      k = solveBezout(na);
      // {p*k0,q*k0,k1,..,k_{n-2}}
      u(0) = p * k(0);
      u(1) = q * k(0);
      u(Eigen::seq(2, n - 1)) = k(Eigen::seq(1, n - 2));
      return u;
    }
  }
 
  // Solves equation a  x + b  y = c, writes answer to x and y
  static void solveDiophantine2vars(IntScalarType a, IntScalarType b,
                                    IntScalarType c, IntScalarType &x,
                                    IntScalarType &y) {
 
    IntScalarType g = extended_gcd(a, b, x, y);
 
    if (c % g != 0) {
      std::cout << a << "  " << b << "  " << c << "   " << g << std::endl;
      puts("Impossible");
      exit(0);
    }
 
    c /= g;
 
    x = x * c;
    y = y * c;
  }
    };
    } // namespace oILAB
 
#endif