1 // htconst.h: declarations of functions for height bounds
2 //////////////////////////////////////////////////////////////////////////
3 //
4 // Copyright 1990-2012 John Cremona
5 //
6 // This file is part of the eclib package.
7 //
8 // eclib is free software; you can redistribute it and/or modify it
9 // under the terms of the GNU General Public License as published by the
10 // Free Software Foundation; either version 2 of the License, or (at your
11 // option) any later version.
12 //
13 // eclib is distributed in the hope that it will be useful, but WITHOUT
14 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 // FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 // for more details.
17 //
18 // You should have received a copy of the GNU General Public License
19 // along with eclib; if not, write to the Free Software Foundation,
20 // Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA
21 //
22 //////////////////////////////////////////////////////////////////////////
23
24 // allow for multiple includes
25 #ifndef _ECLIB_HTCONST_H_
26 #define _ECLIB_HTCONST_H_
27
28 #include <eclib/elog.h>
29 #include <eclib/egr.h>
30 #include <eclib/sieve_search.h>
31
32 double silverman_bound(const Curvedata& CD);
33 double cps_bound(const Curvedata& CD);
height_constant(const Curvedata & CD)34 inline double height_constant(const Curvedata& CD)
35 {
36 // return silverman_bound(CD);
37 double b1 = silverman_bound(CD), b2 = cps_bound(CD);
38 return min(b1,b2);
39 }
40 double egr_height_constant(const Curvedata& CD);
41
42 bigfloat lower_height_bound_alt(const Curvedata& CD);
43 bigfloat lower_height_bound_search(const Curvedata& CD, const bigfloat& reg);
44 bigfloat lattice_const(int r);
45
46 // returns lower bound for height of non-torsion points of good
47 // reduction, following Cremona & Siksek in ANTS7
48 bigfloat egr_lower_height_bound(const Curvedata& CD);
49
50 // Class to find point of minimal height by searching.
51
52 // If egr_flag is set it ignores points which do not have everywhere
53 // good reduction (at all finite primes)
54
55 class point_min_height_finder : public point_processor {
56 Curvedata *E;
57 ComponentGroups CG; // used if egr_flag to test for egr
58 bigint a1,a2,a3,a4,a6;
59 vector<bigint> c;
60 int iso, egr_flag, verbose;
61 bigfloat min_ht;
62 Point Pmin;
63 vector<Point> all_points; //store points found
64 public:
65 point_min_height_finder(Curvedata* EE, int egr=0, int verb=0);
~point_min_height_finder()66 ~point_min_height_finder() {};
67 int process(const bigint& x, const bigint& y, const bigint& z);
68 void search(bigfloat h_lim);
get_min_ht()69 bigfloat get_min_ht() const {return min_ht;}
get_min_ht_point()70 Point get_min_ht_point() const {return Pmin;}
points()71 vector<Point> points() const {return all_points;}
72 };
73
74 // class Interval represents a closed interval [lh,rh] where either
75 // empty=1; or empty=0 and lh <= rh; flags rhinf, lhinf denote
76 // rh=infty and lh=-infty resp.
77
78 class Interval {
79 bigfloat lh, rh;
80 bool empty, lhinf, rhinf;
81 public:
Interval()82 Interval() : empty(0), lhinf(1), rhinf(1) {show(1);}
Interval(const bigfloat & a,const bigfloat & b)83 Interval(const bigfloat& a, const bigfloat& b)
84 :lh(a), rh(b), empty(a>b), lhinf(0), rhinf(0) {show(2);}
Interval(const bigfloat & a)85 Interval(const bigfloat& a)
86 :lh(a), empty(0), lhinf(0), rhinf(1) {show(3);}
Interval(const Interval & I)87 Interval(const Interval& I)
88 :lh(I.lh), rh(I.rh), empty(I.empty), lhinf(I.lhinf), rhinf(I.rhinf) {show(4);}
89 void operator=(const Interval& I)
90 {lh=I.lh; rh=I.rh; empty=I.empty; lhinf=I.lhinf; rhinf=I.rhinf;show(5);}
91 friend ostream& operator<< (ostream&s, const Interval&);
92 // void show(int i) {cout<<i<<": "<<(*this)<<endl;}
show(int i)93 void show(int i) {;}
is_empty()94 bool is_empty() const {return empty;}
95 void intersect(const Interval& I);
96 friend Interval intersect(const Interval& I, const Interval& J);
97 };
98
intersect(const Interval & I,const Interval & J)99 inline Interval intersect(const Interval& I, const Interval& J)
100 {
101 Interval K(I); K.intersect(J); return K;
102 }
103
104 vector<Interval> intersect(const vector<Interval>& L1, const vector<Interval>& L2);
105 //void output(const vector<Interval>& L);
106
107 // class Interval01 represents a closed subinterval [lh,rh] of
108 // [0,1], where either empty=1; or empty=0 and lh <= rh.
109
110 class Interval01 {
111 bigfloat lh, rh;
112 bool empty;
113 public:
Interval01()114 Interval01() : lh(to_bigfloat(0)), rh(to_bigfloat(1)), empty(0) {;}
Interval01(const bigfloat & a,const bigfloat & b)115 Interval01(const bigfloat& a, const bigfloat& b)
116 :lh(a), rh(b), empty(a>b) {;}
117 friend ostream& operator<< (ostream&s, const Interval01&);
118 // void show(int i) {cout<<i<<": "<<(*this)<<endl;}
show(int i)119 void show(int i) {;}
is_empty()120 bool is_empty() const {return empty;}
121 void intersect(const Interval01& I);
122 friend Interval01 intersect(const Interval01& I, const Interval01& J);
123 friend Interval01 operator/(const Interval01& I, const long n);
124 friend Interval01 operator+(const Interval01& I, const bigfloat& shift);
125 };
126
intersect(const Interval01 & I,const Interval01 & J)127 inline Interval01 intersect(const Interval01& I, const Interval01& J)
128 {
129 Interval01 K(I); K.intersect(J); return K;
130 }
131
132 vector<Interval01> intersect(const vector<Interval01>& L1,
133 const vector<Interval01>& L2);
134
135 vector<long> annihilators(const Curvedata& CD, long n);
136
137 // Class to compute lower bound for height of non-torsion points of good
138 // reduction, following Cremona & Siksek in ANTS7
139
140 class CurveHeightConst : public Curvedata, Cperiods {
141 bigfloat c; // archimidean constribution
142 bigfloat e3; // largest (or only) real root
143 bigfloat lower, upper;
144 int n_max;
145 int n_ann;
146 vector<long> ann; // indices of E0/E1 for first few primes
147 bigfloat D(const long n) const; // "denomContrib"
Bnmu(const long n,const bigfloat & mu)148 bigfloat Bnmu(const long n, const bigfloat& mu) const // = B_n(mu)
149 { return exp(n*n*mu+c-D(n)); }
150 int test_target(const bigfloat& target, long k);
151 vector<Interval> canonicalHeightInterval(const bigfloat& target, long k);
152 vector<Interval01> canonicalHeightInterval01(const bigfloat& target, long k);
153 void compute_phase1();
154 void compute_phase2();
155 vector<Interval> solveLEQ(long n, const bigfloat& B);
156 vector<Interval> solveGEQ(long n, const bigfloat& B);
157 vector<Interval01> solveLEQ01(long n, const bigfloat& B);
158 vector<Interval01> solveGEQ01(long n, const bigfloat& B);
159 vector<bigfloat> solveEllNPower(long n, const bigfloat& x1);
160 vector<bigfloat> solveEllNPower01(long n, const bigfloat& x1);
161 bigfloat psi(const bigfloat& x);
162 vector<bigfloat> ordinates(const bigfloat& x);
ztopoint(const bigcomplex & z)163 vector<bigcomplex> ztopoint(const bigcomplex& z)
164 {return ellztopoint(z,I2bigfloat(a1),I2bigfloat(a2),I2bigfloat(a3));}
pointtoz(const bigfloat & x,const bigfloat & y)165 bigcomplex pointtoz(const bigfloat& x, const bigfloat& y)
166 {return ellpointtoz(*this,*this,x,y);}
167 public:
168 CurveHeightConst(const Curvedata& CD);
compute()169 void compute() {compute_phase1(); compute_phase2(); }
get_value()170 bigfloat get_value() const {return lower;}
171 };
172
173 #endif
174