Advanced Search
All Posts :
Lately I have been experimenting with structured Gaussian elimination. This is a technique for reducing large sparse systems of linear equations to much smaller dense ones, which can then be solved using a modular method. Needless to say, I had to generate some large sparse linear systems. I wanted the equations to be written as polynomials, because that is the natural sparse representation in Maple and it makes programming structured Gaussian elimination easier (you can use has and indets, for example). So I tried my favorite randpoly command. This was me trying to generate one linear equation: 
> vars := [seq(x[i],i=1..10000)]:
> f := randpoly(vars,degree=1,terms=20);
Error, (in f) too many levels of recursion
When I complained about this problem before I was offered alternatives by many of the helpful people who read mailing lists, newsgroups, and bug reports. These alternatives work, but I like the randpoly command. It has a simple syntax and it gives me exactly what I want. So I wrote my own command. It's not the fastest method, but it is good for me. It supports all of the randpoly syntax, except the expons option, which doesn't make sense for linear equations. The "variables" can actually be anything: monomials, polynomials, whatever. Use the homogeneous option to avoid introducing an arbitrary constant term. Here's the code: 
randlinear := proc(vars1)
  local opts, cof, hom, ter, den, vars, d, r;

   opts := table(['coeffs'=rand(-100..100), 'homogeneous'=false, 'terms'=5,
     op(select(type, [args], 'equation'))]);
   hom := member('homogeneous', [args]):
   den := member('dense', [args]):
   cof := opts['coeffs']:
   ter := opts['terms']:
   if type(vars1,list) then
     vars := `if`(hom=true, vars1, [op(vars1),1])
   elif type(vars1, set) then
     vars := [op(vars1),`if`(hom=true, NULL, 1)]
   else 
     vars := [vars1,`if`(hom=true, NULL, 1)]
   end if;
   d := `if`(den, 1, max(min(1,evalf(ter/nops(vars))),0));
   r := LinearAlgebra:-RandomVector(nops(vars), density=d, generator=cof);
   inner(convert(r,list), vars)
end proc:
linear polynomial random

Please Wait...