rint.h

//     Copyright (c) 2012 Vadym Kliuchnikov sqct(dot)software(at)gmail(dot)com, Dmitri Maslov, Michele Mosca
//
//     This file is part of SQCT.
// 
//     SQCT is free software: you can redistribute it and/or modify
//     it under the terms of the GNU Lesser General Public License as published by
//     the Free Software Foundation, either version 3 of the License, or
//     (at your option) any later version.
// 
//     SQCT is distributed in the hope that it will be useful,
//     but WITHOUT ANY WARRANTY; without even the implied warranty of
//     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//     GNU Lesser General Public License for more details.
// 
//     You should have received a copy of the GNU Lesser General Public License
//     along with SQCT.  If not, see <http://www.gnu.org/licenses/>.
// 

#ifndef RINT_OPT_PROOF_H
#define RINT_OPT_PROOF_H

#include <gmpxx.h>
#include <complex>
#include "hprhelpers.h"




/// \brief Finds sde given denominator exponent and gde, see formula (1) http://arxiv.org/abs/1206.5236
int sde( int denom_exponent, int gde );
/// \brief Returns canonical form of complex number with complex argument in \f$ [0,\pi/8] \f$
std::complex<double> canonical( const std::complex<double>& val, int& w_power, bool& conj );

/// \brief Specifies type of product of two elements of type T
template< class T>
struct product_type
{
    /// \brief Product type
    typedef T pr_type;
};

/// \brief Specifies that procut of two int number must be of type long
template<>
struct product_type<int>
{
    /// \brief Product type
    typedef long int pr_type;
};

template < class TInt = long int >
struct ring_int_real;

/// \brief Integer in the ring of the form  \f$ a + \omega b + \omega ^ 2 c + \omega ^3 d,\f$
/// for\f$ \omega = e^{ i \frac{\pi}{4} k }\f$ and
/// integers \f$ a,b,c,d \f$.
///
/// \tparam TInt Defines the underlying ring to use ( int, long int, arbitrary precision integers -- mpz_class )
/// \note Precompiled for Tint in {int, long int, mpz_class, resring <8>}
template < class TInt = long int >
struct ring_int
{
    /// \brief Type used for high precision floating point data
    typedef hprHelpers::hpr_real mpclass;
    /// \brief Type of underlying ring
    typedef TInt int_type;
    /// \brief Type of product two elements of int_type
    typedef typename product_type< int_type >::pr_type pr_type;

    /// \brief Convert between different underlying rings
    template<class T>
    explicit ring_int( const ring_int<T>& val );
    /// \brief Value of new object is \f$ 0 + \omega 0 + \omega ^ 2 0 + \omega ^3 0,\f$
    ring_int();
    /// \brief Value of new object is \f$ a + \omega b + \omega ^ 2 c + \omega ^3 d,\f$
    ring_int ( int_type a, int_type b, int_type c, int_type d);
    /// \brief Value of new object is \f$ a + \omega 0 + \omega ^ 2 0 + \omega ^3 0,\f$
    ring_int ( int_type a );
    /// \brief Sets value to \f$ a + \omega b + \omega ^ 2 c + \omega ^3 d,\f$
    void set( int_type a, int_type b, int_type c, int_type d);
    ring_int& operator =(const ring_int& b);
    ring_int (const ring_int& b);

    /// \brief Returns canonical form of the number \f$x\f$, such that \f$Arg(x)\f$ belongs to \f$[0, \pi/8 )\f$
    ring_int canonical() const;

    /// \brief Ring addition
    ring_int& operator +=(const ring_int& y);
    /// \brief Ring substruction
    ring_int& operator -=(const ring_int& y);
    /// \brief Ring addition
    ring_int operator +(const ring_int& y) const;
    /// \brief Ring substruction
    ring_int operator -(const ring_int& y) const;
    /// \brief Ring additive inverse
    ring_int operator -() const;
    /// \brief Ring multiplication
    ring_int operator *(const ring_int& y) const;
    /// \brief Integer division of all integer coefficients by e
    ring_int operator /( int_type e) const;
    /// \brief Integer division of all integer coefficients by e
    ring_int& operator /=( int_type e);
    /// \brief Multiplies by 2^e
    ring_int& operator <<=( long int e);
    /// \brief Divides by 2^e
    ring_int& operator >>=( long int e);
    /// \brief Coefficient near \f$ \omega ^ k, 0 \le k \le 3 \f$
    const int_type& operator[] ( int k ) const;
    /// \brief Coefficient near \f$ \omega ^ k, 0 \le k \le 3 \f$
    int_type& operator[] ( int k );

    /// \brief Divides by \f$ \sqrt{2}^k \f$
    void div_eq_sqrt2( int k );
    /// \brief Divides by \f$ \sqrt{2} \f$
    void div_eq_sqrt2();
    /// \brief Divides by \f$ 2 \f$
    void div_eq_2();
    /// \brief Multiplies by \f$ \omega^k \f$
    void mul_eq_w( int k );
    /// \brief Multiplies by \f$ \omega^k \f$
    void mul_w( int k, ring_int<TInt>& out ) const;
    /// \brief Multiplies by \f$ \omega \f$
    void mul_eq_w();
    /// \brief Multiplies by \f$ \sqrt{2} \f$
    void mul_eq_sqrt2();
    /// \brief Computes gde of a base \f$ \sqrt{2} \f$
    /// \see Section 4 of http://arxiv.org/abs/1206.5236 for defintion and discussion
    int gde() const;
    /// \brief Returns true if all entries divisible by 2, false otherwise
    bool is_div_2() const;
    /// \brief Computes \f$ P(x) \f$
    /// \see Section 5 of http://arxiv.org/abs/1206.5236 for defintion and discussion
    pr_type ipxx() const;
    /// \brief Computes \f$ Q(x) \f$
    /// \see Section 5 of http://arxiv.org/abs/1206.5236 for defintion and discussion
    pr_type ipQxx() const;
    /// \brief Computes absolute value squared
    ring_int_real< pr_type > abs2() const;
     /// \brief Computes Chebyshev distance between integer coordinates
    int_type max_dist(const ring_int& y) const;

    /// \brief Returns \f$ e^{ i \frac{\pi}{4} k } \f$
    static ring_int omega( int k);
    /// \brief Returns \f$ \sqrt{2} = \omega - \omega^3 \f$
    static ring_int sqrt2();
    /// \brief Returns Complex conjugate
    ring_int conjugate() const;

    ring_int g_conjugate() const;
    /// \brief Sets value to its complex conjugate
    void conjugate_eq();

    ring_int i_canonical() const;
    ring_int w_canonical() const;

    /// \brief Returns complex double precision number approximately equal to
    /// \f$ \frac{1}{\sqrt{2}^d }( a + \omega b + \omega ^ 2 c + \omega ^3 d)\f$
    /// \param d Power of \f$ \sqrt{2} \f$ in the denominator
    std::complex<double> toComplex( int d ) const;

    /// \brief Finds complex double precision number approximately equal to
    /// \f$ z = \frac{1}{\sqrt{2}^d }( a + \omega b + \omega ^ 2 c + \omega ^3 d)\f$
    /// \param d Power of \f$ \sqrt{2} \f$ in the denominator
    /// \param re Contains \f$ Re(z)\f$ after fuinction call
    /// \param im Contains \f$ Im(z)\f$ after fuinction call
    void toComplex( int d, double& re, double& im ) const;

    /// \brief Returns complex high precision number approximately equal to
    /// \f$ \frac{1}{\sqrt{2}^d }( a + \omega b + \omega ^ 2 c + \omega ^3 d)\f$
    /// \param d Power of \f$ \sqrt{2} \f$ in the denominator
    std::complex<mpclass> toHprComplex( int d ) const;

    /// \brief Lexicographical order. Returns 1 if \f$x < y\f$, 0 if \f$x = y\f$, -1 if \f$x > y.\f$
    int le( const ring_int& y ) const;
    /// \brief Lexicographical order
    bool operator <( const ring_int& y ) const;
    /// \brief Returns true if two numbers are equal
    bool operator ==( const ring_int& y ) const;
    /// \brief Returns true if two numbers are not equal
    bool operator !=( const ring_int& y ) const;

    /// \brief Returns true if 2 divides x
    static bool is_div_2( pr_type e );
    /// \brief Computes gde of base 2 of integer
    static int gde2 ( pr_type e );
    /// \brief Returns true if gde of a base \f$ \sqrt{2} \f$ of \f$ a + \sqrt{2} b \f$ is 1
    static int isGde1 ( pr_type a, pr_type b );
    /// \brief Returns true if gde of a base \f$ \sqrt{2} \f$ of \f$ a + \sqrt{2} b \f$ is 0
    static int isGde0 ( pr_type a, pr_type b );
    /// \brief Returns true if gde of a base \f$ \sqrt{2} \f$ of \f$ a + \sqrt{2} b \f$ is 2
    static int isGde2 ( pr_type a, pr_type b );
    /// \brief Returns true if two integers can form a column
    bool is_compl( const ring_int& v) const;
    /// \brief Returns true if number is real
    bool is_im_eq0() const;

    bool divides(  const TInt& val ) const;

    /// \brief Stores integer coefficients of the ring integer
    int_type v[4];
};



/// \brief Represents real integers in the ring.
template < class TInt >
struct ring_int_real : public ring_int<TInt>
{
    /// \brief Basic type
    typedef ring_int<TInt> base;
    /// \brief Type of underlying ring
    typedef TInt int_type;
    /// \brief Sets value to \f$0\f$
    ring_int_real() = default;
    /// \brief Sets value to \f$a+\sqrt{2}b\f$
    ring_int_real ( int_type a, int_type b );
    ring_int_real ( const ring_int<TInt>& r );
    /// \brief Fast algorithm for computation of gde. See discussion in Section 5 of http://arxiv.org/abs/1206.5236
    int gde() const;
    ring_int_real g_conjugate() const;



    bool non_negative() const;
    ring_int_real& operator /= ( const ring_int_real& val );
    ring_int_real& operator /= ( const TInt& val );
    bool divides(  const ring_int_real& val ) const;
    bool divides(  const TInt& val ) const;

    void make_positive();

    int_type norm() const;
};

template < class TInt >
std::pair<int,long> unit_log( const ring_int_real<TInt>& val)
{
  int sign = 1; long power = 0; int power_sign = 1;
  ring_int_real<TInt> u = val, gu;
  if( val[0] < 0 && val[1] <=0 )
  {
    sign = -1; power_sign = -1;
  } else if ( val[0] > 0 && val[1] >= 0 )
  {
    sign = 1; power_sign = -1;
  }
  else if ( val[0] * val[1] < 0 )
  {
    if( (val[1] % 2 == 1 || val[1] % 2 == -1 ) != (val[1] < 0) )
    {
      sign = 1; power_sign = 1;
    }
    else
    {
      sign = -1; power_sign = 1;
    }
  }

  auto rr = ring_int_real<TInt>(sign,0);

  if( power_sign == 1 )
    gu = ring_int_real<TInt>(1,1);
  else
    gu = ring_int_real<TInt>(-1,1);

  while( u != rr )
  {
    u = u * gu;
    power += power_sign;
  }
  return std::make_pair(sign,power);
}

template < class TInt >
ring_int_real<TInt> unit_power( const std::pair<int,long>& v )
{
   ring_int_real<TInt> u(v.first,0), gu;
   long inc = 1;
   if( v.second > 0 )
   {
     inc = 1;
     gu = ring_int_real<TInt>(-1,1);
   }
   else
   {
     inc = -1;
     gu = ring_int_real<TInt>(1,1);
   }

   for( long i = 0; i != v.second; i+= inc )
     u = u * gu;

   return u;
}



typedef ring_int_real<mpz_class> zs2type;
typedef mpz_class ztype;
typedef ring_int<mpz_class> zwt;

#endif // RINT_OPT_PROOF_H