Commit 352dcf39926213b9d9b20760b7c1f803d28798e6 - eclib

+3

-1

configure.ac less more

2	2	# Process this file with autoconf to produce a configure script.
3	3
4	4	AC_PREREQ([2.65])
5		AC_INIT([eclib], [20141106], [john.cremona@gmail.com])
	5	AC_INIT([eclib], [20150228], [john.cremona@gmail.com])
6	6	AM_INIT_AUTOMAKE([-Wall])
7	7	AC_MSG_NOTICE([Configuring eclib...])
8	8	AC_CONFIG_SRCDIR([libsrc])
9	9	#AC_CONFIG_HEADERS([config.h])
10	10	AC_CONFIG_MACRO_DIR([m4])
	11	AM_PROG_AR
11	12
12	13	# Initialise libtools:
13	14	LT_INIT

20	21
21	22	# Checks for system features and libraries, via macros defined in m4/ directory
22	23	AX_CXX_COMPILE_STDCXX_0X
	24	# Uncomment this to check for C++-11 compliance (no C++-11 features are used)
23	25	#AX_CXX_COMPILE_STDCXX_11([],[optional])
24	26	AX_CXX_HEADER_TR1_UNORDERED_MAP
25	27	AX_CXX_HEADER_UNORDERED_MAP

+82

-40

libsrc/cubic.cc less more

78	78	}
79	79
80	80	// The quantity called C_1 in the paper, = Norm(h2-h0) and should be
81		// POSITIVE for a reduced form:
	81	// NON-NEGATIVE for a reduced form:
82	82
83	83	bigint cubic::j_c1() const
84	84	{

108	108	}
109	109
110	110	// MINUS the quantity called C_2 in the paper, = -Norm(h0-h1) and should be
111		// NEGATIVE for a reduced form:
	111	// NON-POSITIVE for a reduced form:
112	112
113	113	bigint cubic::j_c2() const
114	114	{

143	143	// POSITIVE for a reduced form:
144	144
145	145	bigint cubic::j_c3() const
146		{
	146	{
147	147	bigint a = coeffs[0], b=coeffs[1], c=coeffs[2], d=coeffs[3];
148	148	bigint b2=b*b;
149	149	bigint b3=b*b2;

169	169	+ 27a2b2c2 + 6b5c - 648b2ca2d + 162adb4 + 1458a3d2*b;
170	170	}
171	171
	172	// The quantity C_4 (not in the paper), = Norm(h1)/8 and should be
	173	// NON-NEGATIVE for a reduced form with C1=0 (i.e. when h0=h2 we want h1>=0).
	174
	175	bigint cubic::j_c4() const
	176	{
	177	bigint a = coeffs[0], b=coeffs[1], c=coeffs[2], d=coeffs[3];
	178	bigint b2=b*b;
	179	bigint b3=b*b2;
	180	bigint b4=b*b3;
	181	bigint a2=a*a;
	182	bigint c2=c*c;
	183	bigint c3=c*c2;
	184
	185	return 27dc3a2 + (27d^2b3 - 54dc2b2 + 9c^4b)a + 9dcb4 - 2c3b3;
	186	}
	187
172	188	//#define DEBUG
173	189
174	190	bigcomplex cubic::hess_root() const
175	191	{
176	192	bigfloat discr = I2bigfloat(disc());
177		if(!is_positive(disc()))
	193	if(!is_positive(disc()))
178	194	{
179	195	cout<<"Error: hess_root called with negative dicriminant!\n";
180	196	return to_bigfloat(0);

183	199	bigfloat Q = I2bigfloat(q_semi());
184	200	bigfloat delta = sqrt(3*discr);
185	201	bigcomplex gamma(-Q,delta); gamma/=(2*P);
186		return gamma;
	202	return gamma;
187	203	}
188	204
189	205	void cubic::hess_reduce(unimod& m)

194	210	while(s)
195	211	{
196	212	s=0;
	213	// NB roundover(a,b) returns c such that a/b=c+x and -1/2 < x <= 1/2,
	214	// so after the shift (when P>0) we have -P <= Q < P.
197	215	k = roundover(-q_semi(),2*p_semi());
198	216	if(!is_zero(k))
199	217	{

209	227	cout << "invert: " << (*this) << endl;
210	228	#endif
211	229	}
	230	}
	231	// Now we have -P <= Q < P <= R and test for boundary condition
	232	if((p_semi()==r_semi()) && (q_semi()<0))
	233	{
	234	invert(m);
	235	#ifdef DEBUG
	236	cout << "Final inversion: " << (*this) << endl;
	237	#endif
212	238	}
213	239	if(a()<0) {::negate(coeffs[0]); ::negate(coeffs[2]);}
214	240	}

220	246	while(s)
221	247	{
222	248	s=0;
223		if(mat_c1()<0)
	249	if(mat_c1()<0)
224	250	{
225	251	s=1; invert(m);
226	252	#ifdef DEBUG

238	264	#endif
239	265	}
240	266	bigint plus1, minus1; plus1=1; minus1=-1;
241		while(mat_c2()>0)
	267	while(mat_c2()>0)
242	268	{
243	269	s=1; x_shift(plus1,m);
244	270	#ifdef DEBUG
245	271	cout << "Shift by +1: "<<(*this)<<endl;
246	272	#endif
247	273	}
248		while(mat_c3()<0)
	274	while(mat_c3()<0)
249	275	{
250	276	s=1; x_shift(minus1,m);
251	277	#ifdef DEBUG

265	291	while(s)
266	292	{
267	293	s=0;
268		if(j_c1()<0)
	294	if(j_c1()<0)
269	295	{
270	296	s=1; invert(m);
271	297	#ifdef DEBUG
272	298	cout << "invert: " << (*this) << endl;
273	299	#endif
274	300	}
275		bigfloat alpha = real_root();
276		bigfloat ra = I2bigfloat(a());
277		bigfloat rb = I2bigfloat(b());
278		bigfloat rc = I2bigfloat(c());
279		bigfloat h0 = (9raraalpha + 6rarb)alpha + 6rarc-rb*rb;
280		bigfloat h1 = 6(rarbalpha + (rbrb-rarc))alpha + 2rbrc;
281		k = Iround(-h1/(2*h0));
282		if (k!=0)
	301	if ((j_c2()>0) \|\| (j_c3()<0))
283	302	{
284	303	s=1;
285		x_shift(k,m);
286		#ifdef DEBUG
287		cout << "Initial shift by "<<k<<": "<<(*this)<<endl;
288		#endif
	304	bigfloat alpha = real_root();
	305	bigfloat ra = I2bigfloat(a());
	306	bigfloat rb = I2bigfloat(b());
	307	bigfloat rc = I2bigfloat(c());
	308	bigfloat h0 = (9raraalpha + 6rarb)alpha + 6rarc-rb*rb;
	309	bigfloat h1 = 6(rarbalpha + (rbrb-rarc))alpha + 2rbrc;
	310	#ifdef DEBUG
	311	cout << "h0 = "<<h0<<endl;
	312	cout << "h1 = "<<h1<<endl;
	313	cout << "-h1/(2h0) = "<<(-h1/(2h0))<<endl;
	314	#endif
	315	k = Iround(-h1/(2*h0));
	316	if (k!=0)
	317	{
	318	x_shift(k,m);
	319	#ifdef DEBUG
	320	cout << "Initial shift by "<<k<<": "<<(*this)<<endl;
	321	#endif
	322	}
	323	// Two loops to guard against rounding error in computing k:
	324	while(j_c2()>0)
	325	{
	326	x_shift(minus1,m);
	327	#ifdef DEBUG
	328	cout << "Shift by -1: "<<(*this)<<endl;
	329	#endif
	330	}
	331	while(j_c3()<0)
	332	{
	333	x_shift(plus1,m);
	334	#ifdef DEBUG
	335	cout << "Shift by +1: "<<(*this)<<endl;
	336	#endif
	337	}
289	338	}
290		while(j_c2()>0)
291		{
292		s=1; x_shift(minus1,m);
293		#ifdef DEBUG
294		cout << "Shift by -1: "<<(*this)<<endl;
295		#endif
296		}
297		while(j_c3()<0)
298		{
299		s=1; x_shift(plus1,m);
300		#ifdef DEBUG
301		cout << "Shift by +1: "<<(*this)<<endl;
302		#endif
303		}
304		#ifdef DEBUG
305		cout<<"C1="<<j_c1()<<", C2="<<j_c2()<<", C3="<<j_c3()<<endl;
306		#endif
307		}
308		if(a()<0) {::negate(coeffs[0]); ::negate(coeffs[1]);}
	339	if ((j_c1()==0) && (j_c4()<0))
	340	{
	341	s=1; invert(m);
	342	#ifdef DEBUG
	343	cout << "invert: " << (*this) << endl;
	344	#endif
	345	}
	346	#ifdef DEBUG
	347	cout<<"C1="<<j_c1()<<", C2="<<j_c2()<<", C3="<<j_c3()<<", C4="<<j_c4()<<endl;
	348	#endif
	349	}
	350	if(a()<0) {::negate(coeffs[0]); ::negate(coeffs[2]);}
309	351	}
310	352
311	353	// Just shifts x:

+1

-0

libsrc/eclib/cubic.h less more

101	101	bigint j_c1() const;
102	102	bigint j_c2() const;
103	103	bigint j_c3() const;
	104	bigint j_c4() const;
104	105
105	106	bigcomplex hess_root() const;
106	107	bigfloat real_root() const; // requires disc<0

+7

-7

libsrc/eclib/mat.h less more

68	68	mat& operator*=(scalar);
69	69	mat& operator/=(scalar);
70	70	const scalar* get_entries()const{return entries;}
	71	long nrows() const {return nro;}
	72	long ncols() const {return nco;}
	73	long rank() const;
	74	long nullity() const;
	75	long trace() const;
	76	vector<long> charpoly() const;
	77	long determinant() const;
71	78	void output(ostream&s=cout) const;
72	79	void output_pari(ostream&s=cout) const;
73	80	void output_pretty(ostream&s=cout) const;

75	82	void read_from_file(string filename); // binary input
76	83
77	84	// non-member (friend) functions and operators
78		friend long nrows(const mat&);
79		friend long ncols(const mat&);
80	85	friend void add_row_to_vec(vec& v, const mat& m, long i);
81	86	friend void sub_row_to_vec(vec& v, const mat& m, long i);
82	87	friend mat operator*(const mat&, const mat&);

135	140	mat submat(const mat& m, const vec& iv, const vec& jv);
136	141	mat echelon(const mat& m, vec& pcols, vec& npcols,
137	142	long& rk, long& ny, scalar& d, int method=0); // default method 0: scalars
138		long rank(const mat&);
139		long nullity(const mat&);
140		long trace(const mat&);
141		vector<long> charpoly(const mat&);
142		long determinant(const mat&);
143	143	mat addscalar(const mat&, scalar);
144	144	vec apply(const mat&, const vec&);

+7

-7

libsrc/eclib/mmatrix.h less more

63	63	//the parameter here is a dummy just to distinguish these
64	64	mat_i shorten(int) const;
65	65	mat_l shorten(long) const;
	66	long nrows() const {return nro;}
	67	long ncols() const {return nco;}
	68	long rank() const;
	69	long nullity() const;
	70	bigint trace() const;
	71	vector<bigint> charpoly() const;
	72	bigint determinant() const;
66	73
67	74	// non-member (friend) functions and operators
68		friend long nrows(const mat_m&);
69		friend long ncols(const mat_m&);
70	75	friend mat_m operator*(const mat_m&, const mat_m&);
71	76	friend vec_m operator*(const mat_m&, const vec_m&);
72	77	friend int operator==(const mat_m&, const mat_m&);

114	119	long& rk, long& ny, bigint& d, int method=0);
115	120	mat_m echelon(const mat_m& m, vec_l& pcols, vec_l& npcols,
116	121	long& rk, long& ny, bigint& d, int method=0);
117		long rank(const mat_m&);
118		long nullity(const mat_m&);
119		bigint trace(const mat_m&);
120		vector<bigint> charpoly(const mat_m&);
121		bigint determinant(const mat_m&);
122	122	mat_m addscalar(const mat_m&, const bigint&);
123	123	vec_m apply(const mat_m&, const vec_m&);
124	124

+1

-1

libsrc/eclib/msubspace.h less more

80	80	msubspace peigenspace(const mat_m& mat, const bigint& lambda, const bigint& pr);
81	81	msubspace psubeigenspace(const mat_m& mat, const bigint& l, const msubspace& s, const bigint& pr);
82	82
83		inline int dim(const msubspace& s) {return ncols(s.basis);} // the dimension
	83	inline int dim(const msubspace& s) {return s.basis.ncols();} // the dimension
84	84	inline bigint denom(const msubspace& s) {return s.denom;} // the denominator
85	85	inline vec_i pivots(const msubspace& s) {return s.pivots;} // the pivot vector
86	86	inline mat_m basis(const msubspace& s) {return s.basis;} // the basis matrix

+5

-3

libsrc/eclib/smat.h less more

64	64	smat select_rows(const vec& rows) const;
65	65	void setrow ( int i, const svec& v); // i counts from 1
66	66	svec row(int) const; // extract row i as an svec
67
	67	int nrows() const {return nro;}
	68	int ncols() const {return nco;}
	69	long rank(scalar mod=DEFAULT_MODULUS);
	70	long nullity(const scalar& lambda, scalar mod=DEFAULT_MODULUS); // nullity of this-lambda*I
	71
68	72	// non-member (friend) functions and operators
69	73
70		friend inline int nrows(const smat& A) {return A.nro;}
71		friend inline int ncols(const smat& A) {return A.nco;}
72	74	friend inline vector<int> dim(const smat& A)
73	75	{vector<int>d; d.push_back(A.nro);d.push_back(A.nco);return d;}
74	76	friend vec operator* (smat& m, const vec& v);

+5

-11

libsrc/eclib/smat_elim.h less more

65	65	void remove( list& L ); // L must be ordered
66	66	ordlist( int m = 10) : list(m) {;}
67	67	};
68
	68
69	69	scalar modulus;
70	70	int rank;
	71
71	72	ordlist* column; // an array of lists, oner per col, of row nos
72	73	// which have nonzero entries in that col
73		int *position;
	74	int *position;
74	75	int *elim_col; // in which row was col eliminated;
75	76	int *elim_row; //order of elimination of rows;
76	77	void clear_col(int,int,list&, int fr = 0, int fc = 0,int M = 0,int *li =0);

82	83	long n_active_entries(); // number of active entries
83	84	double active_density(); // density of non-eliminated part
84	85	void report(); // report current state (density etc)
85
	86
86	87	public:
87	88	int get_rank( ) { return rank; }
88	89	void init( );

117	118	return s;
118	119	}
119	120
120		long rank(smat& sm, scalar mod=DEFAULT_MODULUS);
121
122		inline long nullity(const smat& sm, const scalar& lambda, scalar mod=DEFAULT_MODULUS) // nullity of sm-lambda*I
123		{
124		smat sma(sm); sma-=lambda; return ncols(sm)-rank(sma,mod);
125		}
126
127	121	class ssubspace {
128	122
129	123	public:

143	137	inline scalar mod() const {return modulus;} // the (prime) modulus
144	138
145	139	// non-member (friend) functions and operators
146		friend int dim(const ssubspace& s) {return ncols(s.basis);}
	140	friend int dim(const ssubspace& s) {return s.basis.ncols();}
147	141	friend vec pivots(const ssubspace& s) {return s.pivots;}
148	142	friend smat basis(const ssubspace& s) {return s.basis;}
149	143	friend ssubspace combine(const ssubspace& s1, const ssubspace& s2);

+1

-1

libsrc/eclib/sub.h less more

79	79	subspace psubeigenspace(const mat& m, scalar l, const subspace& s, scalar pr);
80	80
81	81
82		inline int dim(const subspace& s) {return ncols(s.basis);} // the dimension
	82	inline int dim(const subspace& s) {return s.basis.ncols();} // the dimension
83	83	inline scalar denom(const subspace& s) {return s.denom;} // the denominator
84	84	inline vec pivots(const subspace& s) {return s.pivots;} // the pivot vector
85	85	inline mat basis(const subspace& s) {return s.basis;} // the basis matrix

+1

-1

libsrc/egr.cc less more

794	794	cout<<"In comp_map_image, m="<<m;
795	795	cout<<"moduli = "<<moduli<<endl;
796	796	#endif
797		int npts=nrows(m), np=ncols(m);
	797	int npts=m.nrows(), np=m.ncols();
798	798	int i, j, jj;
799	799	if(np==0) return ans;
800	800	for(j=1; j<=np; j++)

+178

-0

libsrc/gpslave.cc less more

	0	// parislave.cc: class for starting up a "slave" background gp process
	1	//////////////////////////////////////////////////////////////////////////
	2	//
	3	// Copyright 1990-2012 John Cremona
	4	//
	5	// This file is part of the eclib package.
	6	//
	7	// eclib is free software; you can redistribute it and/or modify it
	8	// under the terms of the GNU General Public License as published by the
	9	// Free Software Foundation; either version 2 of the License, or (at your
	10	// option) any later version.
	11	//
	12	// eclib is distributed in the hope that it will be useful, but WITHOUT
	13	// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
	14	// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
	15	// for more details.
	16	//
	17	// You should have received a copy of the GNU General Public License
	18	// along with eclib; if not, write to the Free Software Foundation,
	19	// Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA
	20	//
	21	//////////////////////////////////////////////////////////////////////////
	22
	23	#ifdef USE_GP_FACTORING
	24
	25	#include <unistd.h> // for unlink() (not needed on linux)
	26	#include <sstream>
	27	#include <eclib/gpslave.h>
	28	#include <eclib/marith.h>
	29
	30	//#define DEBUG_PARI_SLAVE
	31
	32	parislave the_pari_slave; // The one and only instance
	33
	34	// test if we have an executable gp to use...
	35	int do_we_have_pari(string& gp_path)
	36	{
	37	string path_to_gp(getenv("PATH_TO_GP"));
	38	if (path_to_gp.str().size()==0)
	39	gp_path = string("gp");
	40	else
	41	gp_path = string(path_to_gp)+string("gp");
	42	string comm = string("[ -x ") + gp_path + string(" ]");
	43	#ifdef DEBUG_PARI_SLAVE
	44	cout<<"In do_we_have_pari, command="<<comm<<endl;
	45	#endif
	46	return !system(comm.c_str());
	47	}
	48
	49	parislave::parislave()
	50	{
	51	// test if we have an executable gp to use...
	52	string gpcommand;
	53	dummy = !do_we_have_pari(gpcommand);
	54	if(dummy) return;
	55
	56	strcpy(gpinfilename,"/tmp/gp_in_XXXXXX");
	57	mkstemp(gpinfilename);
	58	string comm=string("rm ")+gpinfilename; system(comm.c_str());
	59	#ifdef DEBUG_PARI_SLAVE
	60	cout<<"Using file "<<gpinfilename<<" for gp input"<<endl;
	61	#endif
	62
	63	comm = string("mkfifo -m 600 ")+gpinfilename;
	64	#ifdef DEBUG_PARI_SLAVE
	65	cout<<"About to run command "<<comm<< " via system() call"<<endl;
	66	#endif
	67	system(comm.c_str());
	68
	69	comm = gpcommand+" -q -f < "+gpinfilename;
	70	#ifdef DEBUG_PARI_SLAVE
	71	cout<<"About to call popen( "<<comm<< " )"<<endl;
	72	#endif
	73	gp = popen(comm.c_str(),"r");
	74	gpin.open(gpinfilename);
	75	}
	76
	77	parislave::~parislave()
	78	{
	79	if(dummy) return;
	80	#ifdef DEBUG_PARI_SLAVE
	81	cout<< "Shutting down gp process..."<<endl;
	82	#endif
	83	ofstream gpin(gpinfilename);
	84	gpin <<"\\q\n"<<flush; // to shut down the gp process
	85	#ifdef DEBUG_PARI_SLAVE
	86	cout<< "Deleting temp file "<<gpinfilename<<endl;
	87	#endif
	88	/*
	89	string comm("rm ");
	90	comm+=gpinfilename;
	91	system(comm.c_str());
	92	*/
	93	unlink(gpinfilename);
	94	}
	95
	96	int parislave::is_prime(const bigint& n)
	97	{
	98	if(dummy)
	99	{
	100	#ifdef DEBUG_PARI_SLAVE
	101	cout<<"dummy call to pari_slave::is_prime(): no gp running!"<<endl;
	102	#endif
	103	return 1;
	104	}
	105	#ifdef DEBUG_PARI_SLAVE
	106	cout << "Writing to gp's input stream: "
	107	<< "print(isprime(" << n << "))\n"<<flush;
	108	#endif
	109	gpin << "print(isprime(" << n << "))\n"<<flush;
	110	#ifdef DEBUG_PARI_SLAVE
	111	cout << "Reading from gp's output stream"<<endl;
	112	#endif
	113	int ans;
	114	fscanf(gp,"%s",gpoutput);
	115	#ifdef DEBUG_PARI_SLAVE
	116	cout<<"just read "<<gpoutput<<" from gp output"<<endl;
	117	#endif
	118	istringstream gpout(gpoutput);
	119	gpout>>ans;
	120	#ifdef DEBUG_PARI_SLAVE
	121	cout<<"Read ans="<<ans<<" from gp output"<<endl;
	122	cout<<"Finished reading from gp output"<<endl;
	123	#endif
	124	return ans;
	125	}
	126
	127	vector<bigint> parislave::factor(const bigint& n, int proof)
	128	{
	129	vector<bigint> plist; bigint p;
	130	if(dummy)
	131	{
	132	#ifdef DEBUG_PARI_SLAVE
	133	cout<<"Calling pdivs_trial()"<<endl;
	134	#endif
	135	return pdivs_trial(n);
	136	}
	137	#ifdef DEBUG_PARI_SLAVE
	138	cout << "Writing to gp's input stream: "
	139	<< "print(Vec(factor(" << n << ")[,1]))\n"<<flush;
	140	#endif
	141	gpin << "print(Vec(factor(abs(" << n << "))[,1]))\n"<<flush;
	142	#ifdef DEBUG_PARI_SLAVE
	143	cout << "Reading from gp's output stream"<<endl;
	144	#endif
	145	int first=1; char c = ',';
	146	while(c!=']')
	147	{
	148	fscanf(gp,"%s",gpoutput);
	149	#ifdef DEBUG_PARI_SLAVE
	150	cout<<"just read "<<gpoutput<<" from gp output"<<endl;
	151	#endif
	152	istringstream gpout(gpoutput);
	153	if(first) gpout>>skipws>>c; // swallow leading "["
	154	first=0;
	155	gpout>>p;
	156	#ifdef DEBUG_PARI_SLAVE
	157	cout<<"Reading p="<<p<<" from gp output"<<endl;
	158	#endif
	159	plist.push_back(p);
	160	gpout>>skipws>>c>>skipws; // swallow ",", but it might turn out to be "]"
	161	}
	162	#ifdef DEBUG_PARI_SLAVE
	163	cout<<"Finished reading from gp output"<<endl;
	164	if(proof) cout<<"Proving primality of factors found..."<<flush;
	165	#endif
	166	if(proof)
	167	for(vector<bigint>::const_iterator pi=plist.begin(); pi!=plist.end(); pi++)
	168	{
	169	p =*pi;
	170	if(!is_prime(p))
	171	{
	172	cout<<"WARNING: pari's factor() returned p="<<p<<" for which isprim(p) FAILS!! Please report.";
	173	}
	174	}
	175	return plist;
	176	}
	177	#endif // USE_GP_FACTORING

+5

-5

libsrc/homspace.cc less more

299	299	cout << "pivots = "<<pivs << endl;
300	300	cout << "denom = "<<d1 << endl;
301	301	}
302		rk = ncols(sp);
	302	rk = sp.ncols();
303	303	coord_vecs.resize(ngens+1); // 0'th is unused
304		for(i=1; i<=ngens; i++)
	304	for(i=1; i<=ngens; i++)
305	305	coord_vecs[i]=sp.row(i);
306	306	//sp=smat(0,0); // clear space
307	307	#else

404	404	}
405	405	if (verbose)
406	406	{
407		cout << "delta matrix done: size "<<nrows(deltamat)<<"x"<<ncols(deltamat)<<". "<<endl;
	407	cout << "delta matrix done: size "<<deltamat.nrows()<<"x"<<deltamat.ncols()<<". "<<endl;
408	408	if(verbose>1)
409	409	cout<<"deltamat = "<<deltamat<<endl;
410	410	cout << "About to compute kernel of delta"<<endl;

514	514	long i= coordindex[ind];
515	515	if (i>0) return m.row(i);
516	516	if (i<0) return -m.row(-i);
517		return vec(ncols(m));
	517	return vec(m.ncols());
518	518	}
519	519
520	520	long homspace::nfproj_coords_from_index(int ind, const vec& bas) const

609	609
610	610	vec homspace::proj_coords(long nn, long dd, const mat& m) const
611	611	{
612		vec ans = vec(ncols(m));
	612	vec ans = vec(m.ncols());
613	613	add_proj_coords(ans, nn, dd, m);
614	614	return ans;
615	615	}

+35

-117

libsrc/mat.cc less more

328	328
329	329	// Definitions of non-member, friend operators and functions
330	330
331		long nrows(const mat& m) {return m.nro;}
332		long ncols(const mat& m) {return m.nco;}
333
334	331	// add/sub row i of mat to v (implemented in mat.cc)
335	332	void add_row_to_vec(vec& v, const mat& m, long i)
336	333	{

672	669	mat transpose(const mat& m)
673	670	{
674	671	long i,j,nr,nc;
675		nr=ncols(m); nc=nrows(m);
	672	nr=m.ncols(); nc=m.nrows();
676	673	mat ans(nr, nc);
677	674	for (i=1; i<=nr; i++)
678	675	for (j=1; j<=nc; j++)

953	950	return ans;
954	951	}
955	952
956		// Original version
957
958		// mat echelon0(const mat& entries, vec& pcols, vec& npcols,
959		// long& rk, long& ny, long& d)
960		// {
961		// long nc, nr;
962		// long r,c,r2,r3,rmin;
963		// long min, mr2c,lastpivot;
964		// rk=0; ny=0; r=1; lastpivot=1;
965		//
966		// mat m(entries);
967		// nc=ncols(m); nr=nrows(m);
968		// pcols.init(nc);
969		// npcols.init(nc);
970		// for (c=1; (c<=nc)&&(r<=nr); c++)
971		// {min = abs(m(r,c));
972		// rmin = r;
973		// for (r2=r+1; (r2<=nr)&&(min!=1); r2++)
974		// { mr2c = abs(m(r2,c));
975		// if ((0<mr2c) && ((mr2c<min) \|\| (min==0))) { min=mr2c; rmin=r2 ;}
976		// }
977		// if (min==0) npcols[++ny] = c;
978		// else
979		// {pcols[++rk] = c;
980		// if (rmin>r) m.swaprows(r,rmin);
981		// for (r3 = r+1 ; r3<=nr; r3++)
982		// elimrows2(m,r,r3,c,lastpivot);
983		// lastpivot=min;
984		// r++;
985		// }
986		// }
987		// for (c = rk+ny+1; c<=nc; c++) npcols[++ny] = c;
988		// pcols = pcols.slice(rk);
989		// npcols = npcols.slice(ny); // truncate index vectors
990		// // cout << "In echelon:\n";
991		// // cout << "pcols = " << pcols << "\n";
992		// // cout << "npcols = " << npcols << "\n";
993		// d=1;
994		// lastpivot=1;
995		// if (ny>0)
996		// {for (r=1; r<=rk; r++) m.clearrow(r);
997		// for (r=1;r<=rk; r++)
998		// for (r2=r+1; r2<=rk; r2++)
999		// elimrows1(m,r2,r,pcols[r2]);
1000		// for (r=1; r<=rk; r++)
1001		// d = (d*m(r,pcols[r]))/gcd(d,m(r,pcols[r]));
1002		// d = abs(d);
1003		// for (r=1; r<=rk; r++)
1004		// {long fac = d/m(r,pcols[r]);
1005		// m.multrow(r,fac);
1006		// }
1007		// }
1008		// else
1009		// for (r=1; r<=rk; r++)
1010		// for (c=1; c<=nc; c++)
1011		// m.set(r,c, (c==pcols[r])); // 0 or 1 !
1012		// return m.slice(rk,nc);
1013		// }
1014
1015		long rank(const mat& entries)
	953	long mat::rank() const
1016	954	{
1017	955	long rk,nr,nc,r,c,r2,r3,rmin;
1018	956	long min, mr2c,lastpivot;
1019	957	rk=0; r=1; lastpivot=1;
1020		mat m(entries);
1021		nc=ncols(m); nr=nrows(m);
	958	mat m(*this);
	959	nc=m.ncols(); nr=m.nrows();
1022	960	for (c=1; (c<=nc)&&(r<=nr); c++)
1023	961	{ min = abs(m(r,c));
1024	962	rmin = r;

1038	976	return rk;
1039	977	}
1040	978
1041		long nullity(const mat& m)
1042		{
1043		return ncols(m)-::rank(m);
1044		}
1045
1046		long trace(const mat& a)
1047		{ long i; long ans=0;
1048		for (i=1; i<=nrows(a); i++) ans += a(i,i);
	979	long mat::nullity() const
	980	{
	981	return nco-rank();
	982	}
	983
	984	long mat::trace() const
	985	{ long i=0; scalar* aii=entries; long ans=0;
	986	for (; i<nro; i++, aii+=(nco+1))
	987	ans += *aii;
1049	988	return ans;
1050	989	}
1051
1052		/* OLD VERSION
1053		vector<long> charpoly(const mat& a)
1054		{ long i,k,r,n = nrows(a);
1055		vec tlist = vec(n);
1056		mat apower(a);
1057		vector<long> clist(n+1);
1058		tlist[1] = trace(a);
1059		for (i=2; i<=n; i++)
1060		{ apower*=a;
1061		tlist[i] = trace(apower);
1062		}
1063		clist[n]=1;
1064		for (k=1; k<=n; k++)
1065		{ long temp = 0;
1066		for (r=1; r<=k; r++) temp+= tlist[r]*clist[n+r-k];
1067		clist[n-k]= -temp/k;
1068		}
1069		return clist;
1070		}
1071		*/
1072
1073		// NEW VERSION -- FADEEV'S METHOD
1074
1075		vector<long> charpoly(const mat& a)
1076		{ long n = nrows(a);
1077		mat b(a);
	990
	991	// FADEEV'S METHOD
	992
	993	vector<long> mat::charpoly() const
	994	{ long n = nrows();
	995	mat b(*this);
1078	996	mat id(idmat((scalar)n));
1079	997	vector<long> clist(n+1);
1080		long t = trace(a);
	998	long t = trace();
1081	999	clist[n] = 1;
1082	1000	clist[n-1] = -t;
1083	1001	for (long i=2; i<=n; i++)
1084		{ b=a(b-tid); // cout << b; // (for testing only)
1085		t=trace(b)/i;
	1002	{ b=(this)(b-t*id); // cout << b; // (for testing only)
	1003	t=b.trace()/i;
1086	1004	clist[n-i] = -t;
1087	1005	}
1088	1006	if (!(b==t*id))

1092	1010	return clist;
1093	1011	}
1094	1012
1095
1096		long determinant(const mat& m)
1097		{
1098		vector<long> cp = charpoly(m);
1099		long det = cp[0];
1100		if (nrows(m)%2==1) det=-det;
1101		return det;
	1013	long mat::determinant() const
	1014	{
	1015	long det = charpoly()[0];
	1016	if (nrows()%2==1)
	1017	return -det;
	1018	else
	1019	return det;
1102	1020	}
1103	1021
1104	1022	mat addscalar(const mat& m, scalar c)
1105	1023	{
1106		return m+(c*idmat((scalar)nrows(m)));
	1024	return m+(c*idmat((scalar)m.nrows()));
1107	1025	}
1108	1026
1109	1027	vec apply(const mat& m, const vec& v) // same as *(mat, vec)
1110	1028	{
1111		long nr=nrows(m), nc=ncols(m);
	1029	long nr=m.nrows(), nc=m.ncols();
1112	1030	vec ans(nr);
1113	1031	if (nc==dim(v))
1114	1032	for (long i=1; i<=nr; i++)

1373	1291	#endif /* TRACE */
1374	1292	long nc,nr,r,c,r2,r3,rmin;
1375	1293	scalar min, mr2c;
1376		nr=nrows(entries), nc=ncols(entries);
	1294	nr=entries.nrows(), nc=entries.ncols();
1377	1295	mat m(nr,nc);
1378	1296	scalar mij=m.entries, entriesij=entries.entries;
1379	1297	long n=nr*nc;

1599	1517
1600	1518	mat ref_via_flint(const mat& M, scalar pr)
1601	1519	{
1602		long nr=nrows(M), nc=ncols(M);
	1520	long nr=M.nrows(), nc=M.ncols();
1603	1521	long i, j;
1604	1522
1605	1523	// copy of the modulus for FLINT

1640	1558	mat ref_via_flint(const mat& M, vec& pcols, vec& npcols,
1641	1559	long& rk, long& ny, scalar pr)
1642	1560	{
1643		long nr=nrows(M), nc=ncols(M);
	1561	long nr=M.nrows(), nc=M.ncols();
1644	1562	long i, j, k;
1645	1563
1646	1564	#ifdef TRACE_FLINT_RREF

1836	1754
1837	1755	double sparsity(const mat& m)
1838	1756	{
1839		long nr=nrows(m), nc=ncols(m);
	1757	long nr=m.nrows(), nc=m.ncols();
1840	1758	if(nr==0) return 1;
1841	1759	if(nc==0) return 1;
1842	1760	scalar* matij=m.entries;

+32

-34

libsrc/mmatrix.cc less more

386	386
387	387	// Definitions of non-member, friend operators and functions
388	388
389		long nrows(const mat_m& m) {return m.nro;}
390		long ncols(const mat_m& m) {return m.nco;}
391
392	389	mat_m operator*(const mat_m& m1, const mat_m& m2)
393	390	{
394	391	long j,k, m=m1.nro, n=m1.nco, p=m2.nco;

597	594	mat_m transpose(const mat_m& m)
598	595	{
599	596	long i,j,nr,nc;
600		nr=ncols(m); nc=nrows(m);
	597	nr=m.ncols(); nc=m.nrows();
601	598	mat_m ans(nr, nc);
602	599	for (i=1; i<=nr; i++)
603	600	for (j=1; j<=nc; j++)

746	743	return ans;
747	744	}
748	745
749		long rank(const mat_m& m1)
	746	long mat_m::rank() const
750	747	{
751	748	long rk,nr,nc,r,c,r2,r3,rmin;
752	749	bigint min, mr2c,lastpivot;
753	750	rk=0; r=1; lastpivot=1;
754		mat_m m(m1);
755		nc=ncols(m); nr=nrows(m);
	751	mat_m m(*this);
	752	nc=m.ncols(); nr=m.nrows();
756	753	for (c=1; (c<=nc)&&(r<=nr); c++)
757	754	{ min = abs(m(r,c));
758	755	rmin = r;

776	773	return rk;
777	774	}
778	775
779		long nullity(const mat_m& m)
780		{
781		return ncols(m)-rank(m);
782		}
783
784		bigint trace(const mat_m& a)
785		{ long i; bigint ans;
786		for (i=1; i<=nrows(a); i++) ans += a(i,i);
	776	long mat_m::nullity() const
	777	{
	778	return ncols()-rank();
	779	}
	780
	781	bigint mat_m::trace() const
	782	{ long i=0; bigint* aii=entries; bigint ans=BIGINT(0);
	783	for (; i<nro; i++, aii+=(nco+1))
	784	ans += *aii;
787	785	return ans;
788	786	}
789
790		// CHARPOLY -- FADEEV'S METHOD
791
792		vector<bigint> charpoly(const mat_m& a)
793		{ long n = nrows(a);
794		mat_m b(a);
	787
	788	// FADEEV'S METHOD
	789
	790	vector<bigint> mat_m::charpoly() const
	791	{ long n = nrows();
	792	mat_m b(*this);
795	793	mat_m id(midmat(n)), tid;
796	794	vector<bigint> clist(n+1);
797		bigint t = trace(a), ii;
	795	bigint t = trace(), ii;
798	796	clist[n] = 1;
799	797	clist[n-1] = -t;
800	798	for (long i=2; i<=n; i++)
801	799	{ tid=t*id;
802	800	b-=tid;
803		b=b*a; // cout << b; // (for testing only)
	801	b=b(this); // cout << b; // (for testing only)
804	802	ii=i;
805		t=trace(b)/ii;
	803	t=b.trace()/ii;
806	804	clist[n-i] = -t;
807	805	}
808	806	tid=t*id;

814	812	return clist;
815	813	}
816	814
817
818		bigint determinant(const mat_m& m)
819		{
820		vector<bigint> cp = charpoly(m);
821		bigint det = cp[0];
822		if (nrows(m)%2==1) det=-det;
823		return det;
	815	bigint mat_m::determinant() const
	816	{
	817	bigint det = charpoly()[0];
	818	if (nrows()%2==1)
	819	return -det;
	820	else
	821	return det;
824	822	}
825	823
826	824	mat_m addscalar(const mat_m& m, const bigint& c)
827	825	{
828		mat_m ans(midmat(nrows(m)));
	826	mat_m ans(midmat(m.nrows()));
829	827	ans*=c;
830	828	ans+=m;
831	829	return ans;

833	831
834	832	vec_m apply(const mat_m& m, const vec_m& v) // same as *(mat_m, vec_m)
835	833	{
836		long nr=nrows(m), nc=ncols(m);
	834	long nr=m.nrows(), nc=m.ncols();
837	835	vec_m ans(nr);
838	836	if (nc==dim(v))
839	837	for (long i=1; i<=nr; i++)

868	866	#endif /* TRACE */
869	867	long nc,nr,r,c,r2,r3,rmin;
870	868	bigint min, mr2c,lastpivot, temp;
871		nr=nrows(m1), nc=ncols(m1);
	869	nr=m1.nrows(), nc=m1.ncols();
872	870	mat_m m(nr,nc);
873	871	for (c=1; c<=nc; c++)
874	872	for (r=1; r<=nr; r++) m(r,c)=mod(m1(r,c),pr);

+3

-3

libsrc/msubspace.cc less more

68	68	// N.B. The following check is strictly unnecessary and slows it down,
69	69	// but is advisable!
70	70	#ifdef CHECK_RESTRICT
71		int check = 1; n = nrows(sb);
	71	int check = 1; n = sb.nrows();
72	72	for (i=1; (i<=n) && check; i++)
73	73	for (j=1; (j<=d) && check; j++)
74	74	check = (ddm.row(i)sb.col(j) == sb.row(i)*ans.col(j));

87	87	bigint d, zero; zero=0;
88	88	vec_i pcols,npcols;
89	89	mat_m m = echelon(mat,pcols,npcols, rank, nullity, d, method);
90		long dim = ncols(m);
	90	long dim = m.ncols();
91	91	mat_m basis(dim,nullity);
92	92	for (n=1; n<=nullity; n++) basis.set(npcols[n],n,d);
93	93	for (r=1; r<=rank; r++)

173	173	long rank, nullity, n, r, i, j;
174	174	vec_i pcols,npcols;
175	175	mat_m m = echmodp(mat,pcols,npcols, rank, nullity, pr);
176		long dim = ncols(m);
	176	long dim = m.ncols();
177	177	mat_m basis(dim,nullity);
178	178	bigint one; one=1;
179	179	for (n=1; n<=nullity; n++) basis.set(npcols[n],n,one);

+4

-4

libsrc/nfd.cc less more

329	329	vec_m nfd::ap(long p)
330	330	{
331	331	mat K = basis(h1->kern).as_mat();
332		long rk = nrows(K);
	332	long rk = K.nrows();
333	333	matop *matlist;
334	334	long k,l,n = h1->modulus, dims=dim(S);
335	335	vec_m apvec(dims);

372	372	mat_m nfd::heckeop(long p)
373	373	{
374	374	mat K = basis(h1->kern).as_mat();
375		long rk = nrows(K);
	375	long rk = K.nrows();
376	376	matop *matlist;
377	377	long j,k,l,n = h1->modulus, dimh=h1->h1dim(), dims=dim(S);
378	378	int bad = ::divides(p,n);

423	423
424	424	bigint inverse(const mat_m& a, mat_m& ainv)
425	425	{
426		long d = nrows(a);
	426	long d = a.nrows();
427	427	mat_m aug=colcat(a,midmat(d));
428	428	long rk, ny; vec pc,npc; bigint denom;
429	429	mat_m ref = echelon(aug, pc, npc, rk, ny, denom);

435	435	void showmatrix(const mat_m& m)
436	436	{
437	437	#ifdef OUTPUT_PARI_STYLE
438		long i,j, nc=ncols(m),nr=nrows(m);
	438	long i,j, nc=m.ncols(),nr=m.nrows();
439	439	cout << "[";
440	440	for(i=0; i<nr; i++)
441	441	{

+12

-8

libsrc/smat.cc less more

474	474
475	475	mat smat::operator*( const mat& m )
476	476	{
477		if( nco != nrows(m) )
	477	if( nco != m.nrows() )
478	478	{
479	479	cerr << "incompatible smat & mat in operator*\n";
480	480	abort();
481	481	}
482		mat product( nro, ncols(m) );
	482	mat product( nro, m.ncols() );
483	483	int i, j, d, t;
484	484	scalar ans;
485	485	for( i = 1; i <= nro; i++ )
486	486	{
487	487	d = col[i-1][0];
488		for( j = 1; j <= ncols(m); j++ )
	488	for( j = 1; j <= m.ncols(); j++ )
489	489	{
490	490	ans = 0;
491	491	for( t = 0; t < d; t++ ) ans += val[i-1][t]*m(col[i-1][t+1],j);

495	495	return product;
496	496	}
497	497
	498	long smat::nullity(const scalar& lambda, scalar mod) // nullity of this-lambda*I
	499	{
	500	smat sma(*this); sma-=lambda; return sma.ncols()-sma.rank(mod);
	501	}
498	502
499	503	// Definitions of non-member, friend operators and functions
500	504

518	522
519	523	vec operator* (smat& m, const vec& v)
520	524	{
521		int r = nrows(m), c=ncols(m);
	525	int r = m.nrows(), c=m.ncols();
522	526	if(c!=dim(v))
523	527	{
524	528	cout<<"Error in smatvec: wrong dimensions ("<<r<<"x"<<c<<")"<<dim(v)<<endl;

533	537
534	538	svec operator* ( const svec& v, const smat& A )
535	539	{
536		if( v.d != nrows(A) )
	540	if( v.d != A.nrows() )
537	541	{
538	542	cout << "incompatible sizes in v*A\n";
539	543	cout << "Dimensions "<<v.d<<" and "<<dim(A)<<endl;
540	544	abort();
541	545	}
542		svec prod(ncols(A));
	546	svec prod(A.ncols());
543	547	map<int,scalar>::const_iterator vi;
544	548	for(vi=v.entries.begin(); vi!=v.entries.end(); vi++)
545	549	prod += (vi->second)*(A.row(vi->first));

548	552
549	553	svec mult_mod_p( const svec& v, const smat& A, const scalar& p )
550	554	{
551		if( v.d != nrows(A) )
	555	if( v.d != A.nrows() )
552	556	{
553	557	cout << "incompatible sizes in v*A\n";
554	558	cout << "Dimensions "<<v.d<<" and "<<dim(A)<<endl;
555	559	abort();
556	560	}
557		svec prod(ncols(A));
	561	svec prod(A.ncols());
558	562	map<int,scalar>::const_iterator vi;
559	563	for(vi=v.entries.begin(); vi!=v.entries.end(); vi++)
560	564	prod.add_scalar_times_mod_p(A.row(vi->first), vi->second,p);

+4

-4

libsrc/smat_elim.cc less more

1040	1040	// the (i,j) entry of dmat goes in the remaining_rows[i-1]'th row,
1041	1041	// remaining_cols[j-1] column. For simplicity of coding, we create
1042	1042	// the new rows as svecs and the use setrow().
1043		int nrd = nrows(dmat); // may be less than nrr since 0 rows are trimmed
	1043	int nrd = dmat.nrows(); // may be less than nrr since 0 rows are trimmed
1044	1044	svec rowi(nco);
1045	1045	for(i=1; i<=nrd; i++)
1046	1046	{

1074	1074	}
1075	1075
1076	1076
1077		long rank(smat& sm, scalar mod)
1078		{
1079		smat_elim sme(sm,mod);
	1077	long smat::rank(scalar mod)
	1078	{
	1079	smat_elim sme(*this,mod);
1080	1080	vec pivs, npivs;
1081	1081	(void) sme.kernel(npivs,pivs);
1082	1082	return sme.get_rank();

+5

-5

libsrc/sub.cc less more

70	70	mat expressvectors(const mat& m, const subspace& s)
71	71	{ vec p = pivots(s);
72	72	long n = dim(s);
73		mat ans(n,ncols(m));
	73	mat ans(n,m.ncols());
74	74	for (int i=1; i<=n; i++) ans.setrow(i, m.row(p[i]));
75	75	return ans;
76	76	}

100	100	// N.B. The following check is strictly unnecessary and slows it down,
101	101	// but is advisable!
102	102	if(cr) {
103		// int check = 1, n = nrows(b);
	103	// int check = 1, n = b.nrows();
104	104	// for (i=1; (i<=n) && check; i++)
105	105	// for (j=1; (j<=d) && check; j++)
106	106	// check = (ddm.row(i)b.col(j) == b.row(i)*ans.col(j));

120	120	scalar d;
121	121	vec pcols,npcols;
122	122	mat m = echelon(m1,pcols,npcols, rank, nullity, d, method);
123		int dim = ncols(m);
	123	int dim = m.ncols();
124	124	mat basis(dim,nullity);
125	125	for (n=1; n<=nullity; n++) basis.set(npcols[n],n,d);
126	126	for (r=1; r<=rank; r++)

211	211	long rank, nullity, n, r, i, j;
212	212	vec pcols,npcols;
213	213	mat m = echmodp(m1,pcols,npcols, rank, nullity, pr);
214		int dim = ncols(m);
	214	int dim = m.ncols();
215	215	mat basis(dim,nullity);
216	216	for (n=1; n<=nullity; n++) basis.set(npcols[n],n,1);
217	217	for (r=1; r<=rank; r++)

228	228	long rank, nullity, i, j, jj, t, tt;
229	229	vec pcols,npcols;
230	230	mat m = echmodp_uptri(m1,pcols,npcols, rank, nullity, pr);
231		int dim = ncols(m);
	231	int dim = m.ncols();
232	232	mat basis(dim,nullity);
233	233	for(j=nullity; j>0; j--)
234	234	{

+4

-4

libsrc/xsplit.cc less more

91	91	data.the_opmat_ = smat(0,0); // releases its space
92	92	}
93	93	else {
94		if( nrows(data.submat_) == 0 ) {
	94	if( data.submat_.nrows() == 0 ) {
95	95	if( depth == 0 ) data.submat_ = h -> s_opmat(depth,1,verbose);
96	96	else data.submat_ = h -> s_opmat_restricted(depth,*(data.nest_),1,verbose);
97	97	}

142	142	// Save space (will recompute when needed)
143	143	//if( ( depth == 0 )
144	144	// && ( dim(s) > 0 )
145		// && ( nrows(data.submat_) > 1000 )
	145	// && ( data.submat_.nrows() > 1000 )
146	146	// && ( data.submatUsage_ == data.numChildren_ ) ) {
147	147	// data.submat_ = smat(0,0);
148	148	//}

541	541
542	542	mat sparse_restrict(const mat& m, const subspace& s)
543	543	{
544		if(dim(s)==nrows(m)) return m; // trivial special case, s is whole space
	544	if(dim(s)==m.nrows()) return m; // trivial special case, s is whole space
545	545	scalar dd = denom(s); // will be 1 if s is a mod-p subspace
546	546	mat b(basis(s));
547	547	smat sm(m), sb(b);

575	575
576	576	smat restrict_mat(const smat& m, const subspace& s)
577	577	{
578		if(dim(s)==nrows(m)) return m; // trivial special case, s is whole space
	578	if(dim(s)==m.nrows()) return m; // trivial special case, s is whole space
579	579	return mult_mod_p(m.select_rows(pivots(s)),smat(basis(s)),MODULUS);
580	580	}
581	581

+10

-4

progs/list_cubics.cc less more

0		// list_cubicx.cc: Program for listing integer cubics with given discriminant bound
	0	// list_cubics.cc: Program for listing integer cubics with given discriminant bound
1	1	//////////////////////////////////////////////////////////////////////////
2	2	//
3	3	// Copyright 1990-2012 John Cremona

87	87	{
88	88	cout<<disc<<"\t";
89	89	cubic g(a,b,c,d);
90		cout<<g<<"\t----------->\t";
	90	cout<<g;
	91	if(disc!=g.disc())
	92	cout<<" [WRONG DISC]";
	93	cout<<"\t---(reduces to)--->\t";
91	94	if(neg)
92	95	g.jc_reduce(m);
93	96	else
94		g.hess_reduce(m);
95		cout<<g<<endl;
	97	g.hess_reduce(m);
	98	cout<<g;
	99	if(disc!=g.disc())
	100	cout<<" [WRONG DISC]";
	101	cout<<endl;
96	102	}
97	103	}
98	104	}

+31

-31

progs/out/list_cubics.out less more

0	0	Enter discriminant bound (positive or negative): Positive discriminants up to 200
1		49 [1,-1,-2,1] -----------> [1,-1,-2,1]
2		81 [1,0,-3,1] -----------> [1,0,-3,1]
3		108 [1,0,-3,0] -----------> [1,0,-3,0]
4		117 [1,1,-3,0] -----------> [0,-3,-1,1]
5		125 [1,-1,-3,2] -----------> [1,2,-2,-1]
6		148 [1,0,-4,2] -----------> [1,-1,-3,1]
7		148 [1,1,-3,-1] -----------> [1,1,-3,-1]
8		169 [1,1,-4,1] -----------> [1,1,-4,1]
	1	49 [1,-1,-2,1] ---(reduces to)---> [1,-2,-1,1]
	2	81 [1,0,-3,1] ---(reduces to)---> [1,-3,0,1]
	3	108 [1,0,-3,0] ---(reduces to)---> [1,0,-3,0]
	4	117 [1,1,-3,0] ---(reduces to)---> [0,-3,-1,1]
	5	125 [1,-1,-3,2] ---(reduces to)---> [1,2,-2,-1]
	6	148 [1,0,-4,2] ---(reduces to)---> [1,-1,-3,1]
	7	148 [1,1,-3,-1] ---(reduces to)---> [1,1,-3,-1]
	8	169 [1,1,-4,1] ---(reduces to)---> [1,-4,1,1]
9	9	Enter discriminant bound (positive or negative): Negative discriminants down to 200
10		-23 [1,-1,0,1] -----------> [1,-1,0,1]
11		-23 [1,0,-1,1] -----------> [1,1,0,1]
12		-27 [1,0,0,1] -----------> [1,0,0,1]
13		-31 [1,1,-2,1] -----------> [1,-1,0,1]
14		-31 [1,1,0,1] -----------> [1,1,0,1]
15		-44 [1,1,-1,1] -----------> [1,1,-1,1]
16		-59 [1,-1,-1,2] -----------> [1,2,0,1]
17		-76 [1,0,-2,2] -----------> [1,-1,-3,1]
18		-83 [1,-1,-3,4] -----------> [1,-1,-1,2]
19		-87 [1,-1,-2,3] -----------> [1,1,-2,1]
20		-100 [1,-1,0,2] -----------> [1,-1,0,2]
21		-104 [1,0,-1,2] -----------> [1,0,-1,2]
22		-107 [1,1,-3,2] -----------> [1,-2,-4,1]
23		-108 [1,0,0,2] -----------> [1,0,0,2]
24		-116 [1,1,0,2] -----------> [1,1,0,2]
25		-135 [1,0,-3,3] -----------> [1,0,-3,1]
26		-147 [1,1,-1,2] -----------> [1,1,-1,2]
27		-152 [1,1,-2,2] -----------> [1,1,-2,2]
28		-172 [1,-1,-1,3] -----------> [1,2,0,2]
29		-175 [1,0,-5,5] -----------> [1,2,-3,1]
30		-176 [1,0,-4,4] -----------> [1,1,-3,1]
31		-199 [1,1,-4,3] -----------> [1,-1,-4,1]
32		-200 [2,2,-1,1] -----------> [1,1,-2,2]
	10	-23 [1,-1,0,1] ---(reduces to)---> [1,-1,0,1]
	11	-23 [1,0,-1,1] ---(reduces to)---> [1,-1,0,1]
	12	-27 [1,0,0,1] ---(reduces to)---> [1,0,0,1]
	13	-31 [1,1,-2,1] ---(reduces to)---> [1,1,0,1]
	14	-31 [1,1,0,1] ---(reduces to)---> [1,1,0,1]
	15	-44 [1,1,-1,1] ---(reduces to)---> [1,-1,1,1]
	16	-59 [1,-1,-1,2] ---(reduces to)---> [1,2,0,1]
	17	-76 [1,0,-2,2] ---(reduces to)---> [1,1,3,1]
	18	-83 [1,-1,-3,4] ---(reduces to)---> [1,1,1,2]
	19	-87 [1,-1,-2,3] ---(reduces to)---> [1,-1,2,1]
	20	-100 [1,-1,0,2] ---(reduces to)---> [1,-1,0,2]
	21	-104 [1,0,-1,2] ---(reduces to)---> [1,0,-1,2]
	22	-107 [1,1,-3,2] ---(reduces to)---> [1,2,4,1]
	23	-108 [1,0,0,2] ---(reduces to)---> [1,0,0,2]
	24	-116 [1,1,0,2] ---(reduces to)---> [1,1,0,2]
	25	-135 [1,0,-3,3] ---(reduces to)---> [1,0,3,1]
	26	-147 [1,1,-1,2] ---(reduces to)---> [1,1,-1,2]
	27	-152 [1,1,-2,2] ---(reduces to)---> [1,1,-2,2]
	28	-172 [1,-1,-1,3] ---(reduces to)---> [1,2,0,2]
	29	-175 [1,0,-5,5] ---(reduces to)---> [1,-2,3,1]
	30	-176 [1,0,-4,4] ---(reduces to)---> [1,-1,3,1]
	31	-199 [1,1,-4,3] ---(reduces to)---> [1,1,4,1]
	32	-200 [2,2,-1,1] ---(reduces to)---> [1,-1,2,2]
33	33	Enter discriminant bound (positive or negative):

+1

-1

progs/out/reduce_cubics.out less more

15	15	Mathews reduced cubic = [83,-149,294,-878]
16	16	after transform by [14,23;3,5]
17	17	Using JC/Julia reduction ...
18		JC/Julia reduced cubic = [83,-100,-245,-650]
	18	JC/Julia reduced cubic = [83,100,245,-650]
19	19	after transform by [14,9;3,2]
20	20	Enter cubic coeffs a, b, c, d:

+1

-1

tests/hecketest.cc less more

380	380	double sparsity(const mat_m& m)
381	381	{
382	382	double count=0;
383		long i,j,nr=nrows(m), nc=ncols(m);
	383	long i,j,nr=m.nrows(), nc=m.ncols();
384	384	for(i=0; i<nr; i++)
385	385	for(j=0; j<nc; j++)
386	386	if(!is_zero(m(i+1,j+1))) count=count+1;

+2

-2

tests/mattest.cc less more

122	122	cout << "Enter any number "; cin >> i;
123	123	}
124	124	{
125		vector<long> cp = charpoly(a);
	125	vector<long> cp = a.charpoly();
126	126	cout << "char. poly. of A has coefficients " << cp << endl;
127		cout << "det(A) = " << determinant(a) << endl;
	127	cout << "det(A) = " << a.determinant() << endl;
128	128	}
129	129	{
130	130	aug = colcat(a,idmat(r));

+2

-2

tests/mmattest.cc less more

109	109	cout << "Enter any number "; cin >> i;
110	110	}
111	111	{
112		vector<bigint> cp = charpoly(a);
	112	vector<bigint> cp = a.charpoly();
113	113	cout << "char. poly. of A has coefficients " << cp << endl;
114		cout << "det(A) = " << determinant(a) << endl;
	114	cout << "det(A) = " << a.determinant() << endl;
115	115	}
116	116	{
117	117	aug = colcat(a,midmat(r));

+6

-6

tests/mspace.cc less more

39	39	cout << "Enter entries of M: ";
40	40	cin >> m;
41	41	cout << " M = " << m;
42		cout << "Trace(M) = " << trace(m) << endl;
	42	cout << "Trace(M) = " << m.trace() << endl;
43	43	//
44	44	mat_m mpower=m;
45	45	for (i=2; i<=r; i++)
46	46	{mpower=mpower*m;
47	47	cout << "m^" << i << " = " << mpower;
48		cout << "Trace(m^" << i << ") = " << trace(mpower) << endl;
	48	cout << "Trace(m^" << i << ") = " << mpower.trace() << endl;
49	49	}
50	50	//
51	51	{
52		vector<bigint> cp = charpoly(m);
	52	vector<bigint> cp = m.charpoly();
53	53	cout << "char. poly. of m has coefficients " << cp << endl;
54	54	}
55		cout << "det(M) = " << determinant(m) << endl;
56		cout << "rank(M) = " << rank(m) << endl;
57		cout << "nullity(M) = " << nullity(m) << endl;
	55	cout << "det(M) = " << m.determinant() << endl;
	56	cout << "rank(M) = " << m.rank() << endl;
	57	cout << "nullity(M) = " << m.nullity() << endl;
58	58	{
59	59	msubspace ker = kernel(m);
60	60	mat_m kerbasis = basis(ker);

+2

-2

tests/smattest.cc less more

444	444	mat ker_mat = echmodp( m, pc, npc, rk, ny, pr);
445	445	cout << " rank using echmodp : " << rk;
446	446	int pop = 0;
447		int nro = nrows(ker_mat);
448		int nco = ncols(ker_mat);
	447	int nro = ker_mat.nrows();
	448	int nco = ker_mat.ncols();
449	449	for( int r = 1; r <= nro; r++ ) {
450	450	for( int c = 1; c <= nco; c++ ) {
451	451	pop += ( ker_mat( r, c ) != 0 );

+6

-6

tests/space.cc less more

35	35	cout << "Enter entries of M: ";
36	36	cin >> m;
37	37	cout << " M = " << m;
38		cout << "Trace(M) = " << trace(m) << endl;
	38	cout << "Trace(M) = " << m.trace() << endl;
39	39	/*
40	40	int i;
41	41	mat mpower=m;
42	42	for (i=2; i<=r; i++)
43	43	{mpower=mpower*m;
44	44	cout << "m^" << i << " = " << mpower;
45		cout << "Trace(m^" << i << ") = " << trace(mpower) << endl;
	45	cout << "Trace(m^" << i << ") = " << mpower.trace() << endl;
46	46	}
47	47	{
48		vector<long> cp = charpoly(m);
	48	vector<long> cp = m.charpoly();
49	49	cout << "char. poly. of m has coefficients " << cp << endl;
50	50	}
51		cout << "det(M) = " << determinant(m) << endl;
52		cout << "rank(M) = " << rank(m) << endl;
53		cout << "nullity(M) = " << nullity(m) << endl;
	51	cout << "det(M) = " << m.determinant() << endl;
	52	cout << "rank(M) = " << m.rank() << endl;
	53	cout << "nullity(M) = " << m.nullity() << endl;
54	54	*/
55	55	cout << endl << "Enter number of times to repeat kernel tests: ";
56	56	cin >> ntimes;

+1

-1

tests/tnfd.cc less more

83	83	if(bad) cout<<"W("; else cout<<"T(";
84	84	cout <<p<<") = ";
85	85	showmatrix(tp);
86		vector<bigint> cptp = charpoly(tp);
	86	vector<bigint> cptp = tp.charpoly();
87	87	for(i=0; i<dims; i++)
88	88	{
89	89	bigint temp = cptp[i];