1 // -*- mode:C++ ; compile-command: "g++ -I.. -g -c vecteur.cc" -*- 2 /* 3 * Copyright (C) 2000,2014 B. Parisse, Institut Fourier, 38402 St Martin d'Heres 4 * 5 * This program is free software; you can redistribute it and/or modify 6 * it under the terms of the GNU General Public License as published by 7 * the Free Software Foundation; either version 3 of the License, or 8 * (at your option) any later version. 9 * 10 * This program is distributed in the hope that it will be useful, 11 * but WITHOUT ANY WARRANTY; without even the implied warranty of 12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 * GNU General Public License for more details. 14 * 15 * You should have received a copy of the GNU General Public License 16 * along with this program. If not, see <http://www.gnu.org/licenses/>. 17 */ 18 #ifndef _GIAC_VECTEUR_H 19 #define _GIAC_VECTEUR_H 20 #include "first.h" 21 #include "index.h" 22 #include <complex> 23 #include <iostream> 24 #ifdef HAVE_LIBGSL 25 #include <gsl/gsl_vector.h> 26 #include <gsl/gsl_matrix.h> 27 #include <gsl/gsl_permutation.h> 28 #endif // HAVE_LIBGSL 29 30 #ifndef NO_NAMESPACE_GIAC 31 namespace giac { 32 #endif // ndef NO_NAMESPACE_GIAC 33 34 const double randrange = 100.0; 35 36 struct unary_function_ptr; 37 struct ref_vecteur; 38 typedef vecteur matrice; // same type but different name 39 typedef vecteur::const_iterator const_iterateur; 40 typedef vecteur::iterator iterateur; 41 typedef std::complex<double> complex_double; 42 #if defined(HAVE_LONG_DOUBLE) && !defined(PNACL) 43 typedef long double long_double; 44 #else 45 typedef double long_double; 46 #endif 47 typedef std::complex<long_double> complex_long_double; 48 double complex_abs(const complex_double & c); 49 double complex_long_abs(const complex_long_double & c); 50 51 // make a matrix with free rows 52 // (i.e. it is possible to modify the answer in place) 53 matrice makefreematrice(const matrice & m); 54 gen freecopy(const gen & g); // this one makes a free copy of a vector, not of a matrix 55 // vecteur related functions 56 vecteur makevecteur(const gen & a); 57 vecteur makevecteur(const gen & a,const gen & b); 58 vecteur makevecteur(const gen & a,const gen & b,const gen & c); 59 vecteur makevecteur(const gen & a,const gen & b,const gen & c,const gen & d); 60 vecteur makevecteur(const gen & a,const gen & b,const gen & c,const gen & d,const gen & e); 61 vecteur makevecteur(const gen & a,const gen & b,const gen & c,const gen & d,const gen & e,const gen & f); 62 vecteur makevecteur(const gen & a,const gen & b,const gen & c,const gen & d,const gen & e,const gen & f,const gen & g); 63 vecteur makevecteur(const gen & a,const gen & b,const gen & c,const gen & d,const gen & e,const gen & f,const gen & g,const gen & h); 64 vecteur makevecteur(const gen & a,const gen & b,const gen & c,const gen & d,const gen & e,const gen & f,const gen & g,const gen & h,const gen & i); 65 vecteur makevecteur(const gen & a,const gen & b,const gen & c,const gen & d,const gen & e,const gen & f,const gen & g,const gen & h,const gen & i,const gen &j); 66 vecteur makevecteur(const gen & a,const gen & b,const gen & c,const gen & d,const gen & e,const gen & f,const gen & g,const gen & h,const gen & i,const gen &j,const gen & k); 67 vecteur makevecteur(const gen & a,const gen & b,const gen & c,const gen & d,const gen & e,const gen & f,const gen & g,const gen & h,const gen & i,const gen &j,const gen & k,const gen & l); 68 vecteur makevecteur(const gen & a,const gen & b,const gen & c,const gen & d,const gen & e,const gen & f,const gen & g,const gen & h,const gen & i,const gen &j,const gen & k,const gen & l,const gen & m); 69 vecteur makevecteur(const gen & a,const gen & b,const gen & c,const gen & d,const gen & e,const gen & f,const gen & g,const gen & h,const gen & i,const gen &j,const gen & k,const gen & l,const gen & m,const gen &n); 70 71 gen makesequence(const gen & a); 72 gen makesequence(const gen & a,const gen & b); 73 gen makesequence(const gen & a,const gen & b,const gen & c); 74 gen makesequence(const gen & a,const gen & b,const gen & c,const gen & d); 75 gen makesequence(const gen & a,const gen & b,const gen & c,const gen & d,const gen & e); 76 gen makesequence(const gen & a,const gen & b,const gen & c,const gen & d,const gen & e,const gen & f); 77 gen makesequence(const gen & a,const gen & b,const gen & c,const gen & d,const gen & e,const gen & f,const gen & g); 78 gen makesequence(const gen & a,const gen & b,const gen & c,const gen & d,const gen & e,const gen & f,const gen & g,const gen & h); 79 gen makesequence(const gen & a,const gen & b,const gen & c,const gen & d,const gen & e,const gen & f,const gen & g,const gen & h,const gen & i); 80 81 ref_vecteur * makenewvecteur(const gen & a); 82 ref_vecteur * makenewvecteur(const gen & a,const gen & b); 83 ref_vecteur * makenewvecteur(const gen & a,const gen & b,const gen & c); 84 ref_vecteur * makenewvecteur(const gen & a,const gen & b,const gen & c,const gen & d); 85 ref_vecteur * makenewvecteur(const gen & a,const gen & b,const gen & c,const gen & d,const gen & e); 86 ref_vecteur * makenewvecteur(const gen & a,const gen & b,const gen & c,const gen & d,const gen & e,const gen & f); 87 ref_vecteur * makenewvecteur(const gen & a,const gen & b,const gen & c,const gen & d,const gen & e,const gen & f,const gen & g); 88 ref_vecteur * makenewvecteur(const gen & a,const gen & b,const gen & c,const gen & d,const gen & e,const gen & f,const gen & g,const gen & h); 89 ref_vecteur * makenewvecteur(const gen & a,const gen & b,const gen & c,const gen & d,const gen & e,const gen & f,const gen & g,const gen & h,const gen & i); 90 91 // numeric root utilities 92 bool francis_schur(std_matrix<gen> & H,int n1,int n2,std_matrix<gen> & P,int maxiter,double eps,bool is_hessenberg,bool complex_schur,bool compute_P,bool no_lapack,GIAC_CONTEXT); 93 94 class matrix_double:public std::vector< std::vector<giac_double> >{ 95 public: 96 // inherited constructors matrix_double()97 matrix_double() : std::vector< std::vector<giac_double> >() { }; matrix_double(int i)98 matrix_double(int i) : std::vector< std::vector<giac_double> >(i) { }; matrix_double(int i,const std::vector<giac_double> v)99 matrix_double(int i,const std::vector<giac_double> v) : std::vector< std::vector<giac_double> >(i,v) { }; matrix_double(const matrix_double::const_iterator b,const matrix_double::const_iterator e)100 matrix_double(const matrix_double::const_iterator b,const matrix_double::const_iterator e) : std::vector< std::vector<giac_double> >(b,e) { }; 101 void dbgprint() const ; 102 }; 103 104 class matrix_complex_double:public std::vector< std::vector<complex_double> >{ 105 public: 106 // inherited constructors matrix_complex_double()107 matrix_complex_double() : std::vector< std::vector<complex_double> >() { }; matrix_complex_double(int i)108 matrix_complex_double(int i) : std::vector< std::vector<complex_double> >(i) { }; matrix_complex_double(int i,const std::vector<complex_double> v)109 matrix_complex_double(int i,const std::vector<complex_double> v) : std::vector< std::vector<complex_double> >(i,v) { }; 110 void dbgprint() const ; 111 }; 112 113 bool balanced_eigenvalues(matrix_double & H,vecteur & res,int maxiter,double eps,bool is_hessenberg,GIAC_CONTEXT); 114 bool francis_schur(matrix_double & H,int n1,int n2,matrix_double & P,int maxiter,double eps,bool is_hessenberg,bool compute_P); 115 bool francis_schur(matrix_complex_double & H,int n1,int n2,matrix_complex_double & P,int maxiter,double eps,bool is_hessenberg,bool compute_P); 116 117 bool matrice2lapack(const matrice & m,double * A,GIAC_CONTEXT); 118 void lapack2matrice(double * A,unsigned rows,unsigned cols,matrice & R); 119 bool lapack_schur(std_matrix<gen> & H,std_matrix<gen> & P,bool compute_P,vecteur & eigenvalues,GIAC_CONTEXT); 120 bool lapack_schur(matrix_double & H,matrix_double & P,bool compute_P,vecteur & eigenvalues); 121 bool eigenval2(matrix_double & H,int n2,giac_double & l1, giac_double & l2); 122 bool std_matrix_gen2std_matrix_giac_double(const std_matrix<gen> & H,matrix_double & H1,bool crunch); 123 bool std_matrix_giac_double2std_matrix_gen(const matrix_double & H,std_matrix<gen> & H1); 124 125 matrice companion(const vecteur & w); 126 gen a_root(const vecteur & v,const std::complex<double> & c0,double eps); 127 vecteur proot(const vecteur & v); 128 vecteur proot(const vecteur & v,double eps); 129 vecteur proot(const vecteur & v,double & eps,int & rprec); 130 vecteur real_proot(const vecteur & v,double eps,GIAC_CONTEXT); 131 gen symb_proot(const gen & e) ; 132 gen _proot(const gen & e,GIAC_CONTEXT); 133 extern const unary_function_ptr * const at_proot ; 134 gen symb_pcoeff(const gen & e) ; 135 gen _pcoeff(const gen & e,GIAC_CONTEXT); 136 vecteur pcoeff(const vecteur & v); 137 extern const unary_function_ptr * const at_pcoeff ; 138 gen symb_peval(const gen & arg1,const gen & arg2) ; 139 extern const unary_function_ptr * const at_peval; 140 gen _peval(const gen & e,GIAC_CONTEXT); 141 142 gen spread_convert(const gen & g,int g_row,int g_col,GIAC_CONTEXT); 143 vecteur lcell(const gen & g); 144 // these int are used to translate relative cells to spreadsheet names 145 bool iscell(const gen & g,int & r,int & c,GIAC_CONTEXT); 146 std::string printcell(const vecteur & v,GIAC_CONTEXT); 147 // given g=cell() or its argument at row i, column j 148 // return 0 if not a cell, 1 if a cell, then compute r and c s.t. g refers to (r,c), 149 // return 2 if g is e.g. A1:B4 compute ref of A1 and B4 150 int cell2pos(const gen & g,int i,int j,int & r,int & c,int & r2,int & c2); 151 // return cell(r,c) argument at (i,j) with same absolute/relative addressing 152 // as g 153 gen pos2cell(const gen & g,int i,int j,int r,int c,int r2,int c2); 154 // insert nrows/ncols of fill in m, e.g. fill= [0,0,2] for a spreadsheet 155 // or ["","",2] or 0 for a matrix 156 matrice matrice_insert(const matrice & m,int insert_row,int insert_col,int nrows,int ncols,const gen & fill,GIAC_CONTEXT); 157 // erase nrows/ncols 158 matrice matrice_erase(const matrice & m,int insert_row,int insert_col,int nrows,int ncols,GIAC_CONTEXT); 159 // extract submatrix 160 matrice matrice_extract(const matrice & m,int insert_row,int insert_col,int nrows,int ncols); 161 void makespreadsheetmatrice(matrice & m,GIAC_CONTEXT); 162 matrice extractmatricefromsheet(const matrice & m,bool value=true); 163 // eval spreadsheet, compute list of dependances in lc 164 void spread_eval(matrice & m,GIAC_CONTEXT); 165 166 gen makesuite(const gen & a); 167 gen makesuite(const gen & a,const gen & b); 168 gen makesuite_inplace(const gen & a,const gen & b); 169 vecteur mergevecteur(const vecteur & v,const vecteur & w); 170 vecteur mergeset(const vecteur & v,const vecteur & w); 171 int vrows(const vecteur & a); 172 void addvecteur(const vecteur & a,const vecteur & b,vecteur & res); 173 vecteur addvecteur(const vecteur & a,const vecteur & b); 174 void subvecteur(const vecteur & a,const vecteur & b,vecteur & res); 175 vecteur subvecteur(const vecteur & a,const vecteur & b); 176 vecteur negvecteur(const vecteur & a); // calls negmodpoly 177 // Multiplication by a scalar 178 void multvecteur(const gen & a,const vecteur & b,vecteur & res); 179 vecteur multvecteur(const gen & a,const vecteur & b); 180 void divvecteur(const vecteur & b,const gen & a,vecteur & res); 181 vecteur divvecteur(const vecteur & b,const gen & a); 182 // Matrix times vecteur 183 void multmatvecteur(const matrice & a,const vecteur & b,vecteur & res); 184 vecteur multmatvecteur(const matrice & a,const vecteur & b); 185 void multvecteurmat(const vecteur & a,const matrice & b,vecteur & res); 186 vecteur multvecteurmat(const vecteur & a,const matrice & b); 187 // Scalar product (does not conjugate, see scalarproduct in misc.h) 188 gen dotvecteur(const vecteur & a,const vecteur & b); 189 giac_double dotvecteur_double(const std::vector<giac_double> & a,const std::vector<giac_double> & c); 190 giac_double dotvecteur(const std::vector<giac_double> & a,const std::vector<giac_double> & c); 191 gen dotvecteur(const gen & a,const gen & b); 192 gen generalized_dotvecteur(const vecteur & a,const vecteur & b,int pos); 193 vecteur generalized_multmatvecteur(const matrice & a,const vecteur & b); 194 // Vect product (3-d) 195 gen complex2vecteur(const gen & g,GIAC_CONTEXT); 196 vecteur cross(const vecteur & v,const vecteur & w,GIAC_CONTEXT); 197 gen cross(const gen & g1,const gen & g2,GIAC_CONTEXT); 198 // ckmvmult check a and b 199 // if a and b are matrices returns matrix product 200 // if a is a matrix and b a vecteur return matr*vect 201 // otherwise returns the dot product of a and b 202 gen ckmultmatvecteur(const vecteur & a,const vecteur & b,GIAC_CONTEXT); 203 void multmatvecteur(const matrix_double & H,const std::vector<giac_double> & w,std::vector<giac_double> & v); 204 205 void vecteur2vector_int(const vecteur & v,int modulo,std::vector<int> & res); 206 bool vecteur2vectvector_int(const vecteur & v,int modulo,std::vector< std::vector<int> > & res); 207 void vector_int2vecteur(const std::vector<int> & v,vecteur & res); 208 void vectvector_int2vecteur(const std::vector< std::vector<int> > & v,vecteur & res); 209 bool iszero(const std::vector<int> & p); 210 211 // matrice related functions 212 bool ckmatrix(const matrice & a,bool allow_embedded_vect); 213 bool ckmatrix(const matrice & a); 214 bool ckmatrix(const gen & a); 215 bool ckmatrix(const gen & a,bool); 216 bool is_squarematrix(const matrice & a); 217 bool is_squarematrix(const gen & a); 218 bool is_fully_numeric(const vecteur & v, int withfracint = 0); 219 bool is_fully_numeric(const gen & a, int withfracint = 0); 220 bool is_exact(const vecteur & v); 221 bool is_exact(const gen & g); 222 // the following functions do not check that a is indeed a matrix 223 int mrows(const matrice & a); 224 int mcols(const matrice & a); 225 void mdims(const matrice &m,int & r,int & c); 226 void mtran(const matrice & a,matrice & res,int ncolres=0); 227 matrice mtran(const matrice & a); 228 gen _tran(const gen & a,GIAC_CONTEXT); 229 extern const unary_function_ptr * const at_tran ; 230 void transpose_double(const matrix_double & a,int r0,int r1,int c0,int c1,matrix_double & at); 231 232 // mmult assumes dimensions are correct 233 void mmult(const matrice & a,const matrice & b,matrice &res); 234 void mmult_atranb(const matrice & a,const matrice & btran,matrice & res); 235 matrice mmult(const matrice & a,const matrice & b); 236 bool mmultck(const matrice & a,const matrice & b,matrice & res); 237 matrice mmultck(const matrice & a,const matrice & b); 238 239 extern int strassen_limit; 240 void strassen_mod(bool skip_reduce,bool add,const std::vector< std::vector<int> > & A,const std::vector< std::vector<int> > & Btran,std::vector< std::vector<int> > & C,int p,int arbeg=0,int arend=0,int acbeg=0,int acend=0,int brbeg=0,int brend=0,int bcbeg=0); 241 void mmult_mod(const std::vector< std::vector<int> > & A,const std::vector< std::vector<int> > & Btran,std::vector< std::vector<int> > & C,int p,int Ar0=0,int Ar1=0,int Ac0=0,int Ac1=0,int Brbeg=0,int Brend=0,int Bcbeg=0,int Crbeg=0,int Ccbeg=0,bool add=false); 242 243 gen mtrace(const matrice & a); 244 gen ckmtrace(const gen & a,GIAC_CONTEXT); 245 extern const unary_function_ptr * const at_trace ; 246 247 gen common_deno(const vecteur & v); 248 249 bool convert(const vecteur & v,std::vector<giac_double> & v1,bool crunch); 250 bool convert(const vecteur & v,std::vector< complex_double > & v1,bool crunch); 251 bool convert(const std::vector<giac_double> & v,vecteur & v1); 252 void matrice2std_matrix_gen(const matrice & m,std_matrix<gen> & M); 253 void std_matrix_gen2matrice(const std_matrix<gen> & M,matrice & m); 254 bool vecteur2index(const vecteur & v,index_t & i); 255 void vect_vector_int_2_vect_vecteur(const std::vector< std::vector<int> > & N,std_matrix<gen> & M); 256 void vect_vecteur_2_vect_vector_int(const std_matrix<gen> & M,int modulo,std::vector< std::vector<int> > & N); 257 bool vecteur2vectvector_int(const vecteur & v,int modulo,std::vector< std::vector<int> > & res); 258 void vectvector_int2vecteur(const std::vector< std::vector<int> > & v,vecteur & res); 259 int dotvector_int(const std::vector<int> & v,const std::vector<int> & w,int modulo); 260 bool multvectvector_int_vector_int(const std::vector< std::vector<int> > & M,const std::vector<int> & v,int modulo,std::vector<int> & Mv); 261 void tran_vect_vector_int(const std::vector< std::vector<int> > & N,std::vector< std::vector<int> > & tN); 262 void apply_permutation(const std::vector<int> & permutation,const std::vector<int> &x,std::vector<int> & y); 263 void vecteur2vector_int(const vecteur & v,int modulo,std::vector<int> & res); 264 265 enum matrix_algorithms { 266 RREF_GAUSS_JORDAN=0, 267 RREF_GUESS=1, 268 RREF_BAREISS=2, 269 RREF_MODULAR=3, 270 RREF_PADIC=4, 271 RREF_LAGRANGE=5 272 }; 273 // For approx linear combination, anything ||<eps will be replaced by 0 274 void linear_combination(const gen & c1,const vecteur & v1,const gen & c2,const vecteur & v2,const gen & c,const gen & cinv,vecteur & v,double eps,int cstart=0); 275 gen exact_div(const gen & a,const gen & b); 276 277 // row reduction from line l and column c to line lmax and column cmax 278 // lmax and cmax are not included 279 // line are numbered starting from 0 280 // if fullreduction is false, reduction occurs under the diagonal only 281 // if dont_swap_below !=0, for line numers < dont_swap_below 282 // the pivot is searched in the line instead of the column 283 // hence no line swap occur 284 // convert_internal = false if we do not want conversion to rational fractions 285 // algorithm=0 Gauss-Jordan, 1 guess, 2 Bareiss, 3 modular, 4 p-adic 286 // rref_or_det_or_lu = 0 for rref, 1 for det, 2 for lu, 3 lu w/o permutation 287 // if unsure fullreduction=true, dont_swap_below=0, convert_internal=true 288 // algorithm=1, rref_or_det_or_lu=0 289 // Returns 0 on failure, 1 on success, 2 if success inverting and no need to remove identity 290 int mrref(const matrice & a, matrice & res, vecteur & pivots, gen & det,int l, int lmax, int c,int cmax, 291 int fullreduction,int dont_swap_below,bool convert_internal,int algorithm,int rref_or_det_or_lu,GIAC_CONTEXT); 292 // holds temporary work storage for block operation 293 struct smallmodrref_temp_t { 294 std::vector< std::vector<int> > Ainvtran,Ainv,CAinv; 295 std::vector<int> permblock,maxrankblock; 296 vecteur pivblock; 297 std::vector<int> y,y1,y2,y3; 298 std::vector<longlong> z,z1,z2,z3; 299 }; 300 bool in_modrref(const matrice & a, std::vector< std::vector<int> > & N,matrice & res, vecteur & pivots, gen & det,int l, int lmax, int c,int cmax,int fullreduction,int dont_swap_below,int modulo,int carac,int rref_or_det_or_lu,const gen & mult_by_det_mod_p,bool inverting,bool no_initial_mod,smallmodrref_temp_t * workptr); 301 bool modrref(const matrice & a, matrice & res, vecteur & pivots, gen & det,int l, int lmax, int c,int cmax, 302 int fullreduction,int dont_swap_below,const gen & modulo,bool ckprime,int rref_or_det_or_lu); 303 bool mrref(const matrice & a, matrice & res, vecteur & pivots, gen & det,GIAC_CONTEXT); 304 bool modrref(const matrice & a, matrice & res, vecteur & pivots, gen & det,const gen& modulo); 305 306 307 // finish full row reduction to echelon form if N is upper triangular 308 // this is done from lmax-1 to l 309 void smallmodrref_upper(std::vector< std::vector<int> > & N,int l,int lmax,int c,int cmax,int modulo); 310 // finish row reduction for matrices with much more columns than rows 311 // version adapted for threads parallelization 312 // assumes that all columns are reduced in parallel, pivots are searched 313 // starting at column 0 314 void in_thread_smallmodrref_upper(std::vector< std::vector<int> > & N,int l,int lpivot,int lmax,int c,int cmax,int modulo,int parallel); 315 void thread_smallmodrref_upper(std::vector< std::vector<int> > & N,int l,int lmax,int c,int cmax,int modulo,int parallel); 316 void free_null_lines(std::vector< std::vector<int> > & N,int l,int lmax,int c,int cmax); 317 int smallmodrref_lastpivotcol(const std::vector< std::vector<int> > & K,int lmax); 318 319 void smallmodrref(int nthreads,std::vector< std::vector<int> > & N,vecteur & pivots,std::vector<int> & permutation,std::vector<int> & maxrankcols,longlong & idet,int l, int lmax, int c,int cmax,int fullreduction,int dont_swap_below,int modulo,int rref_or_det_or_lu,bool reset,smallmodrref_temp_t * workptr,bool allow_block,int carac); 320 void doublerref(matrix_double & N,vecteur & pivots,std::vector<int> & permutation,std::vector<int> & maxrankcols,double & idet,int l, int lmax, int c,int cmax,int fullreduction,int dont_swap_below,int rref_or_det_or_lu,double eps); 321 void modlinear_combination(vecteur & v1,const gen & c2,const vecteur & v2,const gen & modulo,int cstart,int cend=0); 322 void modlinear_combination(std::vector<int> & v1,int c2,const std::vector<int> & v2,int modulo,int cstart,int cend,bool pseudo); 323 vecteur fracmod(const vecteur & v,const gen & modulo,gen * den=0,int prealloc=128); 324 gen modproduct(const vecteur & v, const gen & modulo); 325 matrice mrref(const matrice & a,GIAC_CONTEXT); 326 gen _rref(const gen & a,GIAC_CONTEXT); // first non 0 elem in row is 1 327 extern const unary_function_ptr * const at_rref ; 328 void add_identity(matrice & arref); 329 void add_identity(std::vector< std::vector<int> > & arref); 330 bool remove_identity(matrice & res); 331 bool remove_identity(std::vector< std::vector<int> > & res,int modulo); 332 333 void mdividebypivot(matrice & a,int lastcol,GIAC_CONTEXT); // in-place div by pivots 334 void mdividebypivot(matrice & a,int lastcol=-1); // in-place div by pivots 335 // if lastcol==-1, divide last col, if lastcol==-2 do not divide last col 336 // if lastcol>=0 stop dividing at lastcol 337 void midn(int n,matrice & res); 338 matrice midn(int n); 339 gen _idn(const gen & e,GIAC_CONTEXT); 340 extern const unary_function_ptr * const at_idn ; 341 342 gen fieldcoeff(const gen &F); 343 vecteur vranm(int n,const gen & f,GIAC_CONTEXT); 344 matrice mranm(int n,int m,const gen & f,GIAC_CONTEXT); // random matrix using f 345 gen _ranm(const gen & e,GIAC_CONTEXT); 346 extern const unary_function_ptr * const at_ranm ; 347 gen _randvector(const gen & e,GIAC_CONTEXT); 348 349 bool read_reduction_options(const gen & a_orig,matrice & a,bool & convert_internal,int & algorithm,bool & minor_det,bool & keep_pivot,int & last_col); 350 bool minv(const matrice & a,matrice & res,bool convert_internal,int algorithm,GIAC_CONTEXT); /* if unsure, use true for convert_internal and 1 for algorithm */ 351 bool smallmodinv(const std::vector< std::vector<int> > & a,std::vector< std::vector<int> > & res,int modulo,longlong & det_mod_p); 352 bool modinv(const matrice & a,matrice & res,const gen & modulo,gen & det_mod_p); 353 // solve a*x=b where a and b have integer coeffs 354 // using a p-adic algorithm, n is the precision required 355 // c must be the inverse of a mod p 356 vecteur padic_linsolve_c(const matrice & a,const vecteur & b,const matrice & c,unsigned n,const gen & p,unsigned reconstruct=0); 357 // solve a*x=b where a and b have integer coeffs using a p-adic algorithm 358 // lcmdeno of the answer may be used to give an estimate of the 359 // least divisor element of a if b is random 360 // returns 0 if no invertible found, -1 if det==0, 1 otherwise 361 int padic_linsolve(const matrice & a,const vecteur & b,vecteur & res,gen & p,gen & det_mod_p,gen & h2,unsigned reconstruct=0,int maxtry=4); 362 // a is a matrix with integer coeffs 363 // find p such that a mod p has the same rank 364 // rankline and rankcols are the lines/cols used for the submatrix 365 // asub of max rank, ainv is the inverse of asub mod p 366 // return -1 or the rank 367 int padic_linsolve_prepare(const matrice & a,gen & p,std::vector<int> & ranklines, std::vector<int> & rankcols,matrice & asub,matrice & ainv,vecteur & compat,vecteur & kernel); 368 // solve a prepared non Cramer linear system 369 bool padic_linsolve_solve(const matrice & a,const gen & p,const std::vector<int> & ranklines,const std::vector<int> & rankcols,const matrice & asub,const matrice & ainv,const vecteur & compat,const vecteur & b,vecteur & sol); 370 gen _padic_linsolve(const gen & g,GIAC_CONTEXT); 371 372 matrice minv(const matrice & a,GIAC_CONTEXT); 373 gen mdet(const matrice & a,GIAC_CONTEXT); 374 gen _det(const gen & a,GIAC_CONTEXT); 375 extern const unary_function_ptr * const at_det ; 376 gen _det_minor(const gen & g,bool convert_internal,GIAC_CONTEXT); 377 extern const unary_function_ptr * const at_det_minor ; 378 gen det_minor(const matrice & a,vecteur lv,bool convert_internal,GIAC_CONTEXT); 379 void sylvester(const vecteur & v1,const vecteur & v2,matrice & res); 380 matrice sylvester(const vecteur & v1,const vecteur & v2); 381 gen _sylvester(const gen & a,GIAC_CONTEXT); 382 extern const unary_function_ptr * const at_sylvester ; 383 384 gen _hessenberg(const gen & g,GIAC_CONTEXT); 385 extern const unary_function_ptr * const at_hessenberg ; 386 387 bool balance_krylov(const matrix_double & H,std::vector<giac_double> & d,int niter=5,double cutoff=1e-8); 388 bool probabilistic_pmin(const matrice & m,vecteur & w,bool check,GIAC_CONTEXT); 389 vecteur mpcar_int(const matrice & A,bool krylov,GIAC_CONTEXT,bool compute_pmin); 390 391 void mod_pcar(std_matrix<gen> & N,vecteur & res,bool compute_pmin); 392 bool mod_pcar(const matrice & A,std::vector< std::vector<int> > & N,int modulo,bool & krylov,std::vector<int> & res,GIAC_CONTEXT,bool compute_pmin); 393 bool mod_pcar(std::vector< std::vector<int> > & N,int modulo,bool & krylov,std::vector<int> & res,GIAC_CONTEXT,bool compute_pmin); 394 vecteur mpcar_hessenberg(const matrice & A,int modulo,GIAC_CONTEXT); 395 gen _pcar_hessenberg(const gen & g,GIAC_CONTEXT); 396 extern const unary_function_ptr * const at_pcar_hessenberg ; 397 void mhessenberg(std::vector< std::vector<int> > & H,std::vector< std::vector<int> > & P,int modulo,bool compute_P); 398 void hessenberg(std_matrix<gen> & H,std_matrix<gen> & P,GIAC_CONTEXT); 399 void hessenberg_ortho(std_matrix<gen> & H,std_matrix<gen> & P,GIAC_CONTEXT); 400 void hessenberg_ortho(std_matrix<gen> & H,std_matrix<gen> & P,int firstrow,int n,bool compute_P,int already_zero,double eps,GIAC_CONTEXT); 401 void qr_ortho(std_matrix<gen> & H,std_matrix<gen> & P,bool computeP,GIAC_CONTEXT); 402 void hessenberg_schur(std_matrix<gen> & H,std_matrix<gen> & P,int maxiter,double eps,GIAC_CONTEXT); 403 gen _hessenberg(const gen & g0,GIAC_CONTEXT); 404 405 vecteur mpcar(const matrice & a,vecteur & B,bool compute_B,GIAC_CONTEXT); 406 vecteur mpcar(const matrice & a,vecteur & Bv,bool compute_Bv,bool convert_internal,GIAC_CONTEXT); 407 gen _pcar(const gen & a,GIAC_CONTEXT); 408 extern const unary_function_ptr * const at_pcar ; 409 410 // if jordan is false, errors for non diagonalizable matrices 411 // if jordan is true, d is a matrix, not a vecteur 412 bool egv(const matrice & m,matrice & p,vecteur & d, GIAC_CONTEXT, bool jordan,bool rational_jordan_form,bool eigenvalues_only); 413 gen symb_egv(const gen & a); 414 matrice megv(const matrice & a,GIAC_CONTEXT); 415 gen _egv(const gen & a,GIAC_CONTEXT); 416 extern const unary_function_ptr * const at_egv ; 417 gen _svd(const gen & a,GIAC_CONTEXT); 418 419 gen symb_egvl(const gen & a); 420 vecteur megvl(const matrice & a,GIAC_CONTEXT); 421 gen _egvl(const gen & a,GIAC_CONTEXT); 422 extern const unary_function_ptr * const at_egvl ; 423 424 gen symb_jordan(const gen & a); 425 vecteur mjordan(const matrice & a,bool rational_jordan,GIAC_CONTEXT); 426 gen _jordan(const gen & a,GIAC_CONTEXT); 427 gen jordan(const gen & a,bool rational_jordan,GIAC_CONTEXT); 428 gen _rat_jordan(const gen & a,GIAC_CONTEXT); 429 extern const unary_function_ptr * const at_jordan ; 430 matrice rat_jordan_block(const vecteur & v,int n,bool pseudo); 431 gen _rat_jordan_block(const gen &args,GIAC_CONTEXT); 432 matrice pseudo_rat_to_rat(const vecteur & v,int n); 433 434 // if g is an expression, replace x by m in g (m assumed to be diagonal) 435 matrice diagonal_apply(const gen & g,const gen & x,const matrice & m,GIAC_CONTEXT); 436 // apply an analytic function to a square matrix 437 matrice analytic_apply(const unary_function_ptr * u,const matrice & m,GIAC_CONTEXT); 438 matrice analytic_apply(const gen &ux,const gen & x,const matrice & m,GIAC_CONTEXT); 439 matrice matpow(const matrice & m,const gen & n,GIAC_CONTEXT); 440 gen _matpow(const gen & a,GIAC_CONTEXT); 441 442 bool mker(const matrice & a,vecteur & v,int algorithm,GIAC_CONTEXT); 443 bool mker(const matrice & a,vecteur & v,GIAC_CONTEXT); // algorithm=0 444 vecteur mker(const matrice & a,GIAC_CONTEXT); 445 gen _ker(const gen & a,GIAC_CONTEXT); 446 extern const unary_function_ptr * const at_ker ; 447 448 bool mimage(const matrice & a, vecteur & v,GIAC_CONTEXT); 449 vecteur mimage(const matrice & a,GIAC_CONTEXT); 450 gen _image(const gen & a,GIAC_CONTEXT); 451 extern const unary_function_ptr * const at_image ; 452 453 gen symb_cross(const gen & arg1,const gen & arg2); 454 gen symb_cross(const gen & a); 455 gen _cross(const gen & a,GIAC_CONTEXT); 456 extern const unary_function_ptr * const at_cross ; 457 458 gen symb_size(const gen & a); 459 gen _size(const gen & a,GIAC_CONTEXT); 460 extern const unary_function_ptr * const at_size ; 461 bool vecteur2index(const vecteur & v,index_t & i); 462 #ifdef HAVE_LIBGSL 463 gsl_vector * vecteur2gsl_vector(const vecteur & v,GIAC_CONTEXT); // allocate 464 int vecteur2gsl_vector(const vecteur & v,gsl_vector * w,GIAC_CONTEXT); // no alloc 465 int vecteur2gsl_vector(const_iterateur it,const_iterateur itend,gsl_vector * w,GIAC_CONTEXT); 466 vecteur gsl_vector2vecteur(const gsl_vector * v); 467 int matrice2gsl_matrix(const matrice & m,gsl_matrix * w,GIAC_CONTEXT); 468 gsl_matrix * matrice2gsl_matrix(const matrice & m,GIAC_CONTEXT); 469 matrice gsl_matrix2matrice(const gsl_matrix * v); 470 vecteur gsl_permutation2vecteur(const gsl_permutation * p,GIAC_CONTEXT); 471 #endif // HAVE_LIBGSL 472 473 bool mlu(const matrice & a0,vecteur & P,matrice & L,matrice & U,GIAC_CONTEXT); 474 gen lu(const gen & a,GIAC_CONTEXT); 475 extern const unary_function_ptr * const at_lu ; 476 gen qr(const gen & a,GIAC_CONTEXT); 477 extern const unary_function_ptr * const at_qr ; 478 matrice thrownulllines(const matrice & res); 479 480 extern const unary_function_ptr * const at_lll_reduce ; 481 extern const unary_function_ptr * const at_lll ; 482 483 gen _cholesky(const gen & a,GIAC_CONTEXT); 484 extern const unary_function_ptr * const at_cholesky ; 485 gen _svd(const gen & a,GIAC_CONTEXT); 486 extern const unary_function_ptr * const at_svd ; 487 gen _basis(const gen & a,GIAC_CONTEXT); 488 extern const unary_function_ptr * const at_basis ; 489 gen _ibasis(const gen & a,GIAC_CONTEXT); 490 extern const unary_function_ptr * const at_ibasis ; 491 492 // inert function used for internal representation in spreadsheets 493 gen _cell(const gen & a,GIAC_CONTEXT); 494 extern const unary_function_ptr * const at_cell ; 495 496 int alphaposcell(const std::string & s,int & r); 497 gen l2norm(const vecteur & v,GIAC_CONTEXT); 498 matrice gramschmidt(const matrice & m,bool normalize,GIAC_CONTEXT); 499 matrice gramschmidt(const matrice & m,matrice & r,bool normalize,GIAC_CONTEXT); 500 // lll decomposition of m 501 matrice lll(const matrice & M,matrice & L,matrice & O,matrice &A,GIAC_CONTEXT); 502 matrice lll(const matrice & m,GIAC_CONTEXT); 503 gen _lll(const gen & g,GIAC_CONTEXT); 504 505 506 void matrice2std_matrix_gen(const matrice & m,std_matrix<gen> & M); 507 void std_matrix_gen2matrice(const std_matrix<gen> & M,matrice & m); 508 void std_matrix_gen2matrice_destroy(std_matrix<gen> & M,matrice & m); 509 510 bool is_integer_vecteur(const vecteur & m,bool intonly=false); 511 bool is_integer_matrice(const matrice & m,bool intonly=false); 512 bool is_mod_vecteur(const vecteur & m,std::vector<int> & v,int & p); 513 bool is_mod_matrice(const matrice & m,std::vector< std::vector<int> > & M,int & p); 514 515 typedef void (*bezout_fonction)(const gen & a,const gen & b,gen & u,gen &v,gen &d); 516 517 struct environment; 518 bool ismith(const matrice & Aorig, matrice & U,matrice & A,matrice & V,GIAC_CONTEXT); 519 bool smith(const std_matrix<gen> & Aorig,std_matrix<gen> & U,std_matrix<gen> & A,std_matrix<gen> & V,environment * env,GIAC_CONTEXT); 520 bool ihermite(const matrice & Aorig, matrice & U,matrice & A,GIAC_CONTEXT); 521 bool hermite(const std_matrix<gen> & Aorig,std_matrix<gen> & U,std_matrix<gen> & A,environment * env,GIAC_CONTEXT); 522 gen _ihermite(const gen & g,GIAC_CONTEXT); 523 gen _ismith(const gen & g,GIAC_CONTEXT); 524 #if !defined NSPIRE && !defined FXCG 525 gen _csv2gen(const gen & g,GIAC_CONTEXT); 526 matrice csv2gen(std::istream & i,char sep,char nl,char decsep,char eof,GIAC_CONTEXT); 527 #endif 528 529 // find index i of x in v that is i such that v[i] <= x < v[i+1] 530 // where v[-1]=-inf, and v[v.size()]=+inf 531 int dichotomy(const std::vector<giac_double> & v,double x); 532 533 #ifndef NO_NAMESPACE_GIAC 534 } // namespace giac 535 #endif // ndef NO_NAMESPACE_GIAC 536 537 #endif // _GIAC_VECTEUR_H 538