diff --git a/libsrc/eclib/mwprocs.h b/libsrc/eclib/mwprocs.h index b318452..7868f73 100644 --- a/libsrc/eclib/mwprocs.h +++ b/libsrc/eclib/mwprocs.h @@ -40,7 +40,7 @@ class mw : public point_processor { Curvedata *E; vector basis; int rank, maxrank; - bigfloat *height_pairs; + vector> height_pairs; bigfloat reg; int verbose, process_points; bigfloat& mat_entry(int i, int j); @@ -49,7 +49,6 @@ class mw : public point_processor { saturator satsieve; public: mw(Curvedata*, int verb=0, int pp=1, int maxr=999); - ~mw(); // processing of new points, with saturation at primes up to sat // (default MAXSATPRIME, none if sat==0) @@ -70,7 +69,7 @@ class mw : public point_processor { inline bigfloat& mw::mat_entry(int i, int j) { - return *(height_pairs + (i*MAXRANK) + j); + return height_pairs[i][j]; } class sieve { diff --git a/libsrc/mwprocs.cc b/libsrc/mwprocs.cc index ad0f65e..e611f6f 100644 --- a/libsrc/mwprocs.cc +++ b/libsrc/mwprocs.cc @@ -79,7 +79,7 @@ bigfloat min_real(vector array) //cout<<"minr finally " << minr << "\n"; return minr; } - + int order_real_roots(vector& bnd, vector roots); //checks (and returns) how many roots are actually real, and puts those in //bnd, in increasing order, by calling set_the_bound @@ -89,76 +89,70 @@ int set_the_bounds(vector& bnd, bigfloat x0, bigfloat x1, bigfloat x2); //is on [x0,infty]. The function returns 3 in the first case, 1 in the second. //If x0 overflows, it returns 0. A warning is printed out. -#define matentry(m,i,j) *((m)+((i)*MAXRANK)+(j)) - -bigfloat det(bigfloat *m, long m_size); +bigfloat det(vector>& m, long m_size); // fwd declaration: det and detminor jointly recursive -bigfloat* get_minor(bigfloat *m, long m_size, long i0, long j0) +vector> get_minor(vector>& m, long m_size, long i0, long j0) { long i, j, ii, jj; - bigfloat *minor = new bigfloat[MAXRANK*MAXRANK]; + vector> minor(MAXRANK, vector(MAXRANK)); for (i=0; i=i0)ii++; for (j=0; j=j0) jj++; - matentry(minor,i,j) = matentry(m,ii,jj); + minor[i][j] = m[ii][jj]; } } return minor; } -bigfloat det_minor(bigfloat *m, long m_size, long i0, long j0) +bigfloat det_minor(vector>& m, long m_size, long i0, long j0) { - bigfloat *minor = get_minor(m,m_size,i0,j0); - bigfloat det_return = det(minor, m_size-1); - delete [] minor; - return det_return; + vector> minor = get_minor(m,m_size,i0,j0); + return det(minor, m_size-1); } -bigfloat det(bigfloat *m, long m_size) +bigfloat det(vector>& m, long m_size) { + bigfloat one = to_bigfloat(1); switch (m_size) { - case 0: - return to_bigfloat(1); break; + case 0: + return one; break; case 1: - return matentry(m,0,0); break; + return m[0][0]; break; case 2: - return matentry(m,0,0)*matentry(m,1,1) - matentry(m,1,0)*matentry(m,0,1); + return m[0][0]*m[1][1] - m[0][1]*m[1][0]; break; default: - long i,j,i0; - bigfloat ans=to_bigfloat(1), pivot=matentry(m,0,0), piv, temp, - eps=to_bigfloat(1.0e-6); - for(i0=0; i00) // swap rows 0, i0: { - ans=to_bigfloat(-1); + negate = 1; for(j=0; j(MAXRANK)); } // NB We cannot use the default parameter mechanism as this must fit @@ -332,7 +321,7 @@ int mw::process(const Point& PP, int sat) if (rank==0) // first non-torsion point { reg = height(P); - mat_entry(0,0) = reg; // ie height_pairs[0][0] = reg + height_pairs[0][0] = reg; basis.push_back(P); rank=1; if (verbose) cout<<" is generator number 1\n"; #ifdef DEBUG