1 // msubspace.cc: implementations of multiprecision subspace class
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 //#define CHECK_RESTRICT // define this to make restrict_mat and prestrict check
25 // invariance of subspaces.
26 #include <eclib/marith.h>
27 #include <eclib/msubspace.h>
28
29 // Definitions of nonmember, nonfriend operators and functions:
30
combine(const msubspace & s1,const msubspace & s2)31 msubspace combine(const msubspace& s1, const msubspace& s2)
32 {
33 bigint d = s1.denom * s2.denom;
34 mat_m b1=s1.basis, b2=s2.basis;
35 long nr = b1.nro, nc = b2.nco;
36 mat_m b = b1*b2;
37 bigint g; long n=nr*nc; bigint* bp=b.entries;
38 while ((n--)&&(!is_one(g))) g=gcd(g,*bp++);
39 if(!(is_zero(g)||is_one(g)))
40 {
41 d/=g; bp=b.entries; n=nr*nc; while(n--) (*bp++)/=g;
42 }
43 vec_i p = s1.pivots[s2.pivots];
44 return msubspace(b,p,d);
45 }
46
47 //This one is used a LOT
restrict_mat(const mat_m & m,const msubspace & s)48 mat_m restrict_mat(const mat_m& m, const msubspace& s)
49 { long i,j,k,d = dim(s), n=m.nro;
50 bigint dd = s.denom;
51 mat_m ans(d,d);
52 const mat_m& sb = s.basis;
53 bigint *ap, *a=m.entries, *b=sb.entries, *bp, *c=ans.entries, *cp;
54 int *pv=s.pivots.entries;
55 for(i=0; i<d; i++)
56 {
57 bp=b; k=n; ap=a+n*(pv[i]-1);
58 while(k--)
59 {
60 cp=c; j=d;
61 while(j--)
62 {
63 *cp++ += *ap * *bp++;
64 }
65 ap++;
66 }
67 c += d;
68 }
69 // N.B. The following check is strictly unnecessary and slows it down,
70 // but is advisable!
71 #ifdef CHECK_RESTRICT
72 int check = 1; n = sb.nrows();
73 for (i=1; (i<=n) && check; i++)
74 for (j=1; (j<=d) && check; j++)
75 check = (dd*m.row(i)*sb.col(j) == sb.row(i)*ans.col(j));
76 if (!check)
77 {
78 cerr<<"Error in restrict_mat: msubspace not invariant!"<<endl;
79 }
80 #endif
81 return ans;
82 }
83
kernel(const mat_m & mat,int method)84 msubspace kernel(const mat_m& mat, int method)
85 {
86 long rank, nullity, n, r, i, j;
87 bigint d, zero; zero=0;
88 vec_i pcols,npcols;
89 mat_m m = echelon(mat,pcols,npcols, rank, nullity, d, method);
90 long dim = m.ncols();
91 mat_m basis(dim,nullity);
92 for (n=1; n<=nullity; n++) basis.set(npcols[n],n,d);
93 for (r=1; r<=rank; r++)
94 { i = pcols[r];
95 for (j=1; j<=nullity; j++) basis.set(i,j, -m(r,npcols[j]));
96 }
97 msubspace ans(basis, npcols, d);
98 return ans;
99 }
100
image(const mat_m & mat,int method)101 msubspace image(const mat_m& mat, int method)
102 {
103 vec_i p,np;
104 bigint d;
105 long rank, nullity;
106 mat_m b = transpose(echelon(transpose(mat),p,np,rank,nullity,d,method));
107 msubspace ans(b,p,d);
108 return ans;
109 }
110
eigenspace(const mat_m & mat,const bigint & lambda,int method)111 msubspace eigenspace(const mat_m& mat, const bigint& lambda, int method)
112 {
113 mat_m m = addscalar(mat,-lambda);
114 msubspace ans = kernel(m,method);
115 return ans;
116 }
117
subeigenspace(const mat_m & mat,const bigint & l,const msubspace & s,int method)118 msubspace subeigenspace(const mat_m& mat, const bigint& l, const msubspace& s, int method)
119 {
120 mat_m m = restrict_mat(mat,s);
121 msubspace ss = eigenspace(m, l*(denom(s)),method);
122 msubspace ans = combine(s,ss );
123 return ans;
124 }
125
pcombine(const msubspace & s1,const msubspace & s2,const bigint & pr)126 msubspace pcombine(const msubspace& s1, const msubspace& s2, const bigint& pr)
127 {
128 bigint d = s1.denom * s2.denom; // redundant since both should be 1
129 mat_m b1=s1.basis, b2=s2.basis;
130 mat_m b = matmulmodp(b1,b2,pr);
131 vec_i p = s1.pivots[s2.pivots];
132 return msubspace(b,p,d);
133 }
134
prestrict(const mat_m & m,const msubspace & s,const bigint & pr)135 mat_m prestrict(const mat_m& m, const msubspace& s, const bigint& pr)
136 { long i,j,k,d = dim(s), n=m.nro;
137 bigint dd = s.denom; // will be 1 if s is a mod-p msubspace
138 mat_m ans(d,d);
139 const mat_m& sb = s.basis;
140 bigint *ap, *a=m.entries, *b=sb.entries, *bp, *c=ans.entries, *cp;
141 int *pv=s.pivots.entries;
142 for(i=0; i<d; i++)
143 {
144 bp=b; k=n; ap=a+n*(pv[i]-1);
145 while(k--)
146 {
147 cp=c; j=d;
148 while(j--)
149 {
150 *cp += mod(*ap * *bp++, pr);
151 *cp = mod(*cp, pr);
152 cp++;
153 }
154 ap++;
155 }
156 c += d;
157 }
158 #ifdef CHECK_RESTRICT
159 mat_m& left = matmulmodp(m,sb,pr), right = matmulmodp(sb,ans,pr);
160 if(dd!=1) left*=dd;
161 int check = (left==right);
162 if (!check)
163 {
164 cerr<<"Error in prestrict: msubspace not invariant!"<<endl;
165 }
166 #endif
167 return ans;
168 }
169
pkernel(const mat_m & mat,const bigint & pr)170 msubspace pkernel(const mat_m& mat, const bigint& pr)
171 {
172 long rank, nullity, n, r, i, j;
173 vec_i pcols,npcols;
174 mat_m m = echmodp(mat,pcols,npcols, rank, nullity, pr);
175 long dim = m.ncols();
176 mat_m basis(dim,nullity);
177 bigint one; one=1;
178 for (n=1; n<=nullity; n++) basis.set(npcols[n],n,one);
179 for (r=1; r<=rank; r++)
180 { i = pcols[r];
181 for (j=1; j<=nullity; j++) basis.set(i,j, -m(r,npcols[j]));
182 }
183 msubspace ans(basis, npcols, one);
184 return ans;
185 }
186
pimage(const mat_m & mat,const bigint & pr)187 msubspace pimage(const mat_m& mat, const bigint& pr)
188 {
189 vec_i p,np;
190 long rank, nullity;
191 mat_m b = transpose(echmodp(transpose(mat),p,np,rank,nullity,pr));
192 msubspace ans(b,p,BIGINT(1));
193 return ans;
194 }
195
peigenspace(const mat_m & mat,const bigint & lambda,const bigint & pr)196 msubspace peigenspace(const mat_m& mat, const bigint& lambda, const bigint& pr)
197 {
198 mat_m m = addscalar(mat,-lambda);
199 msubspace ans = pkernel(m,pr);
200 return ans;
201 }
202
psubeigenspace(const mat_m & mat,const bigint & l,const msubspace & s,const bigint & pr)203 msubspace psubeigenspace(const mat_m& mat, const bigint& l, const msubspace& s, const bigint& pr)
204 {
205 mat_m m = prestrict(mat,s,pr);
206 msubspace ss = peigenspace(m, l*(denom(s)),pr);
207 msubspace ans = pcombine(s,ss,pr);
208 return ans;
209 }
210
211
212 //Attempts to lift from a mod-p msubspace to a normal Q-msubspace by expressing
213 //basis as rational using modrat and clearing denominators
214 //
lift(const msubspace & s,const bigint & pr,int trace)215 msubspace lift(const msubspace& s, const bigint& pr, int trace)
216 {
217 bigint modulus=pr,dd,n,d; long nr,nc,nrc;
218 int succ,success=1;
219 bigint lim=sqrt(pr>>1);
220 mat_m m = s.basis; bigint *mp;
221 if(trace)
222 {
223 cout << "Lifting mod-p msubspace.\n basis mat_m mod "<<pr<<" is:\n";
224 cout << m;
225 cout << "Now lifting back to Q.\n";
226 cout << "lim = " << lim << "\n";
227 }
228 nr = m.nro; nc = m.nco; nrc = nr*nc; mp=m.entries;
229 dd=1;
230 while(nrc--)
231 {
232 succ = modrat(*mp++,modulus,lim,n,d);
233 success = success && succ;
234 dd=lcm(d,dd);
235 }
236 if(!success)
237 cout << "Problems encountered with modrat lifting of msubspace." << endl;
238 dd=abs(dd);
239 if(trace) cout << "Common denominator = " << dd << "\n";
240 nrc=nr*nc; mp=m.entries;
241 while(nrc--)
242 {
243 *mp=mod(dd*(*mp),pr);
244 mp++;
245 }
246 msubspace ans(m, pivots(s), dd);
247 return ans;
248 }
249