qideal.h

// FILE qideal.h
#ifndef __QIDEAL_H__
#define __QIDEAL_H__

#include "quads.h"

class Factorization;

class Qideal  {
protected:
  INT a,b,c;   // HNF Z-basis: I=c[a,b+w]
  INT ac;      // ac = a*c = smallest integer
  INT nm;      // nm = a*c*c = norm
  long iclass;  // 0 if principal, 1 if not, -1 if undetermined
  Quad g0,g1;   // Reduced Z-basis (if iclass!=-1):
                // I=<g0,g1>, g0 non-zero of minimal norm,
                // and iclass=0 iff I=<g0> is principal
  long index;   // Index (from 1) of this ideal in the standard
                // sorting of ideals of the same norm, or -1 if not known
  Factorization *F; // pointer to the ideal's factorization (only set when needed)
  vector<Quad> the_residues;
  void make_residues();

public:
//constructors
  Qideal();                          // base initialization for class `Qideal'
  Qideal(const Qideal& );            // the copy constructor
  Qideal(const INT&, const INT&, const INT&); // construct given Z-basis abc
  explicit Qideal(const INT&);               // princ ideal gen by integer
  explicit Qideal(const Quad&);               // princ ideal gen by Quad
  explicit Qideal(const vector<Quad>&);       // ideal generated by list of Quads */
  explicit Qideal(const pair< vector<INT>, vector<INT>> &);  // ideal from Z-gens */
  explicit Qideal(const string& s);           // ideal from label N.i
  Qideal& operator= (const Qideal& I)
  { a=I.a;
    b=I.b;
    c=I.c;
    ac=I.ac;
    nm=I.nm;
    iclass=I.iclass;
    g0=I.g0;
    g1=I.g1;
    index=I.index;
    F = 0;  // the reason why we cannot just use the default copy
    return *this;
  }

  // destructor // needs to free up F
  ~Qideal();

  // member access
  INT get_a() const {return a;}
  INT get_b() const {return b;}
  INT content() const {return c;}
  INT smallest_integer() const {return ac;}
  INT norm() const {return nm;}
  Quad zgen(int i) const {return (i? Quad(b*c,c): Quad(ac));} // HNF Z-module gens
  Quad gen();          // smallest element, so a generator iff principal
  vector<Quad> gens(); // reduced Z-module ggens
  vector<INT> get_rv() const {return {ac, b*c};} // real parts of Z-module gens
  vector<INT> get_iv() const {return { ZERO, c};}   // imag parts of Z-module gens

  void set_index(int ind=0); // if 0 (default) computes the correct index
  long get_index()
  {
    if (index<1) set_index();
    return index;
  }

  Factorization factorization(); // sets F if necessary then returns F

  //
  //operators
  //
  int operator== (const Qideal& f) const {return (a==f.a)&&(b==f.b)&&(c==f.c);}
  int operator!= (const Qideal& f) const {return (a!=f.a)||(b!=f.b)||(c!=f.c);}
  // just for sorting:
  int operator<(const Qideal& J) const {return (nm<J.nm) || ((nm==J.nm) && (index<J.index)); }
//
  Qideal operator+(const INT&) const;
  Qideal operator+(const Quad&) const;
  Qideal operator+(const Qideal&) const;
  friend Qideal operator+(const INT&, const Qideal&);
  friend Qideal operator+(const Quad&, const Qideal&);
  void operator+=(const INT&);
  void operator+=(const Quad&);
  void operator+=(const Qideal&);
//
  Qideal operator*(const INT&) const;
  Qideal operator*(const Quad&) const;
  Qideal operator*(const Qideal&) const;
  friend Qideal operator*(const INT&, const Qideal&);
  friend Qideal operator*(const Quad&, const Qideal&);
  void operator*=(const INT&);
  void operator*=(const Quad&);
  void operator*=(Qideal);

  Qideal intersection(const Qideal& I);
  Quad second_generator(const Quad& a); // with nonzero a in this, return b such that this=(a,b)

  // Assuming this*J is principal, sets g to a generator and returns a
  // 2x2 matrix of determinant g whose columns generate this and J,
  // the first column being (g0,g1)
  mat22 AB_matrix(Qideal&J, Quad&g);
  // as above with J = conj(this), g=norm(this)
  mat22 AB_matrix();
  // As above with (2,1)-entry in N
  mat22 AB_matrix_of_level(const Qideal&N, Quad&g);
  mat22 AB_matrix_of_level(const Qideal&J, const Qideal&N, Quad&g);

  // test if this ideal is coprime to another ideal or a Quad:
  int is_coprime_to(const Qideal& I) const
  {
    INT g = gcd(ac, I.ac);
    INT one(1);
    return (g==one? 1: gcd(g, c*I.c*(b-I.b))==one);
  }
  int is_coprime_to(const Quad& alpha) const
  {
    return is_coprime_to(Qideal(alpha));
  }

  // versions returning more data

  // return 1 iff this is coprime to J; if so, set r in this and s in J with r+s=1
  int is_coprime_to(Qideal& J, Quad& r, Quad& s);

  // return 1 iff this is coprime to alpha; if so, set inverse so an inverse of alpha modulo this
  int is_coprime_to(const Quad& alpha, Quad& inverse);

  // Return 1 iff this is coprime to (c,d); if so, set x,y so c*x+d*y =1
  // modulo this ideal.  If fix=1, ensure that y is coprime to this.
  int is_coprime_to(const Quad& c, const Quad& d, Quad& x, Quad& y, int fix=0);

  // return J = (c/d)*this coprime to N, or (if anti=1) J such that J*this=(c) and d=1
  // (implemented in primes.cc)
  // These are not const since they call is_principal() which may fill data
  Qideal equivalent_coprime_to(const Qideal& N, Quad& c, Quad& d, int anti=0);
  // Same again if you don't need c,d
  Qideal equivalent_coprime_to(const Qideal& N, int anti=0)
  {
    Quad cc, dd;
    return equivalent_coprime_to(N, cc, dd, anti);
  }
  // return J coprime to N such that (J/this)^2, or (J*this)^2 if anti, is principal
  Qideal equivalent_mod_2_coprime_to(const Qideal& N, int anti=0);

  int genus_character(const INT& D); // one unram char value in {+1,-1}
  vector<int> genus_character(); // vector of values in {0,1} of length r = #prime factors of field disc
                                 // adding to 0 mod 2
  long genus_class(int contract=0); // integer in [0,2^r-1] whose bits are genus_character() when contract=0
                                    // or same reduced mod 2^{r-1} when contract=1
  int has_square_class() {return (genus_class()==0);}

  vector<INT> possible_unramified_twists();  // sublist of Quad::all_disc_factors() consisting of those D not 1
                                               // for which chi_D(Q)=+1 for all prime powers Q||N

  // return J such that J^2 is equivalent to this (or J^2*this is
  // principal if anti==1), or if no such J exists (i.e., if the ideal
  // class is not a square), return the zero ideal.  (implemented in
  // primes.cc)
  Qideal sqrt_class(int anti=0);

  // return J coprime to N such that J^2*this is principal; if no such
  // J exists (i.e., if the ideal class is not a square, return the
  // zero ideal.  (implemented in primes.cc)
  Qideal sqrt_coprime_to(const Qideal& N);
  //
  Qideal operator/(const INT&) const;
  Qideal operator/(const Quad&) const;
  Qideal operator/(const Qideal&) const;
  friend Qideal operator/(const INT&, const Qideal&);
  friend Qideal operator/(const Quad&, const Qideal&);
  void operator/=(const INT&);
  void operator/=(const Quad&);
  void operator/=(const Qideal&);

//functions defined in qideal.cc unless inline
  int is_zero() const {return c.is_zero();}
  int is_nonzero() const {return c.is_nonzero();}
  int is_principal();          // fills iclass if necessary
  int is_principal(Quad& g);   // same but puts generator into g
  int ideal_class();           // i s.t. this is equivalent to the i'th ideal in the class group
  int is_primitive() const {return c.is_one();}
  int is_square();
  int is_Galois_stable() const {return ::divides(a, (2*b+Quad::t));}
  int is_prime();
  int is_prime_power();
  int is_equivalent(const Qideal& I) const;
  int is_anti_equivalent(const Qideal& I) const;
  int contains(const INT& n) const  {return ::divides(ac, n);}
  int contains(const Quad& alpha) const;
  int contains(const Qideal& I) const  {return ::divides(c,I.c) && ::divides(ac,I.ac) && contains(I.zgen(1));}

  // for alpha in this ideal, return its Z-coeffs {x,y} w.r.t.Z-gens.
  vector<INT> zcoeffs(const Quad& alpha) const {return {(alpha.r-b*alpha.i)/ac, alpha.i/c}; }
  Quad zcombo(const INT& x, const INT& y) const { return Quad(x*ac+y*b*c, y*c); }
  int divides(const INT& n) const {return contains(n);}
  int divides(const Quad& alpha) const  {return contains(alpha);}
  int divides(const Qideal& I) const  {return contains(I);}

  Qideal divide_out(const Qideal& I); // largest factor of this coprime to I
  void make_primitive();              // divide out content
  Qideal primitive_part() const;      // largest primitive factor of this (=this/content)

  Qideal conj() const;            // returns the conjugate ideal

  Quad reduce(const Quad& alpha); // reduction of alpha modulo this ideal (not const; may fill())
  Quad resnum(long i); // the i'the residue mod this, in standard order (0'th is 0)
  long numres(const Quad& alpha) const; // the index of a residue mod this, in standard order (0'th is 0)
  // return a list of (reduced) residues modulo this ideal:
  vector<Quad> residues();
  // return a list of (reduced) invertible residues modulo this ideal
  vector<Quad> invertible_residues();
  // return a list of (reduced) invertible residues modulo this ideal, and a list of their inverses
  pair<vector<Quad>, vector<Quad>> invertible_residues_and_inverses();

// i/o
  friend string gens_string(Qideal& I); // not const as it calls I.fill()
  friend ostream& operator<<(ostream& s, const Qideal& x);
  friend istream& operator>>(istream& s, Qideal& x);

  friend class Quadprime;

private:
  int ok() const;                 // checks that [a,b+w] *is* an ideal
  void fill();                    // determines iclass, g0, g1
  void abc_from_HNF(const vector<INT>&);
};

// return i if I is equivalent to the i'th ideal in Jlist, else -1
int find_ideal_class(const Qideal& I, const vector<Qideal>& Jlist);
// return i if I is equivalent mod squares to the i'th ideal in Jlist, else -1
int find_ideal_class_mod_squares(const Qideal& I, const vector<Qideal>& Jlist);
// return i where I is equivalent to the i'th ideal in Quad::class_group
inline int find_ideal_class(const Qideal& I) {return find_ideal_class(I, Quad::class_group);}
// return the equivalent ideal in Quad::class_group
Qideal class_representative(Qideal I);
// return the ideal J in Quad::class_group for which I*J = <g> is principal and set g
Qideal class_anti_representative(Qideal I, Quad& g);

// An AB-matrix with given first column
mat22 AB_matrix(const Quad& a, const Quad& c);

Qideal Qideal_from_norm_index(INT N, int i); // i'th ideal of norm N

long val(const Qideal& factor, const Qideal& dividend);

// If a is in (c,d) return 1 and x,y such that a=c*x+d*y, else return 0
int express2gens(const Quad& a, const Quad& c, const Quad& d, Quad& x, Quad& y);

// These three functions return lists which are not sorted in the standard way
vector<Qideal> primitive_ideals_with_norm(INT N, int both_conj=1);
vector<Qideal> ideals_with_norm(INT N, int both_conj=1);
vector<Qideal> ideals_with_bounded_norm(INT maxnorm, int both_conj=1);

string ideal_label(Qideal& I); // returns label of ideal I
string gens_string(Qideal& I);  // returns string of gens, of the form (x) if principal or (x,y) ideal I

// Class to hold sorted lists of ideals of given norm
class Qideal_lists {
  static map<INT, vector<Qideal>> N_to_Ilist;
public:
  static vector<Qideal> ideals_with_norm(INT N, int both_conj=1);
  static vector<Qideal> ideals_with_bounded_norm(INT maxnorm, int both_conj=1);

  static void clear() //  clear the static list
  {
    N_to_Ilist.clear();
  }
};

// compute a list of ideals coprime to N whose classes generate the 2-torsion
vector<Qideal> make_nulist(Qideal& N);

// return 1 iff a is the square mod M of some r in reslist
int squaremod(const Quad& a, const Qideal& M, const vector<Quad>& reslist);
vector<int> makechitable(const Qideal& L, const vector<Quad>& reslist);

// test functions used in qidltest.cc

void residuetest(Qideal& I);

// return a pair {rv,iv} of vectors of the real and imagainary parts
// of a, a*w for a in gens; so that the [rv[i],iv[i]] are Z-module
// generators of the ideal generated by the gens
pair<vector<INT>, vector<INT>> restrict_scalars(const vector<Quad>& gens);


#endif

// END OF FILE qideal.h