File: polys.cc

package info (click to toggle)
eclib 2014-09-21-1
  • links: PTS, VCS
  • area: main
  • in suites: jessie, jessie-kfreebsd
  • size: 4,216 kB
  • ctags: 4,287
  • sloc: cpp: 45,827; makefile: 222; sh: 108
file content (115 lines) | stat: -rw-r--r-- 3,229 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
// polys.cc : implements interface to NTL polynomials
//////////////////////////////////////////////////////////////////////////
//
// Copyright 1990-2012 John Cremona
// 
// This file is part of the eclib package.
// 
// eclib is free software; you can redistribute it and/or modify it
// under the terms of the GNU General Public License as published by the
// Free Software Foundation; either version 2 of the License, or (at your
// option) any later version.
// 
// eclib 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 General Public License
// for more details.
// 
// You should have received a copy of the GNU General Public License
// along with eclib; if not, write to the Free Software Foundation,
// Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA
// 
//////////////////////////////////////////////////////////////////////////
 
#include <eclib/marith.h>
#include <eclib/polys.h>

FqPoly reduce(const ZPoly& f, const galois_field& Fq)
{
  NewFqPoly(Fq,fmodq);
  SetDegree(fmodq,Degree(f));
  for(int i=0; i<=Degree(f); i++)
    SetCoeff(fmodq,i,ZtoGF(Fq,PolyCoeff(f,i)));
  return fmodq;
}

vector<gf_element> roots(const FqPoly& f)
{
  // make f monic:
  FqPoly f1=f;
  MakeMonic(f1);
  // reduce to distinct roots case:
  ZZ_pX X; SetX(X); 
  ZZ_pX g = PowerXMod(ZZ_p::modulus(),f1)-X;
  vec_ZZ_p r; FindRoots(r,GCD(f1,g)); 
  vector<gf_element>ans;
  for(int i=0; i<r.length(); i++) ans.push_back(r[i]);
  return ans;
}


vector<bigint> rootsmod(const vector<bigint>& coeffs, bigint q)
{
  galois_field Fq(q);
  NewFqPoly(Fq,f);
  unsigned long i, deg = coeffs.size()-1;
  SetDegree(f,deg);
  for (i=0; i<=deg; i++) SetCoeff(f,i,ZtoGF(Fq,coeffs[i]));

  vector<gf_element> r = roots(f);
  vector<bigint>ans;
  for(i=0; i<r.size(); i++) ans.push_back(LiftGF(r[i]));

  sort(ans.begin(),ans.end());
  return ans;
}

//#define TRACE_ROOTS
vector<bigrational> roots(const vector<bigint>& coeffs)
{
#ifdef TRACE_ROOTS
  cout<<"Finding rational roots of polynomial f with coefficients "<<coeffs<<endl;
#endif
  ZZX f,g; ZZ c;
  vector<bigrational> ans;
  int i, d = coeffs.size()-1;  // degree
  if(d<1) return ans;
  for(i=0; i<=d; i++)
    SetCoeff(f,d-i,coeffs[i]);
#ifdef TRACE_ROOTS
  cout<<"f = "<<f<<endl;
#endif
  vec_pair_ZZX_long factors;
  factor(c,factors,f);
#ifdef TRACE_ROOTS
  cout<<"f has "<<factors.length()<<" factors"<<endl;
#endif
  for(i=0; i<factors.length(); i++)
    {
      g = factors[i].a;
#ifdef TRACE_ROOTS
      cout<<"factor "<<g<<" has degree "<<deg(g)<<endl;
#endif
      if(deg(g)==1)
        {
          bigrational root = bigrational(-coeff(g,0),coeff(g,1));
#ifdef TRACE_ROOTS
          cout<<"root "<<root<<endl;
#endif
          ans.push_back(root);
        }
    }
  sort(ans.begin(), ans.end());
  return ans;
}

// find the number of roots of X^3 + bX^2 + cX + d = 0 (mod p)
int nrootscubic(const bigint& b, const bigint& c, const bigint& d, const bigint& p)
{
  vector<bigint> coeffs;
  coeffs.push_back(d);
  coeffs.push_back(c);
  coeffs.push_back(b);
  coeffs.push_back(BIGINT(1));
  return rootsmod(coeffs,p).size();
}