Artifact efb81b11541234a1bc204976ccf2ba4a3d9ee926d51ec0e73a0c92700621654a:
- Executable file
r36/src/randpoly.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 15045) [annotate] [blame] [check-ins using] [more...]
module randpoly; % A random (generalized) polynomial generator % F.J.Wright@Maths.QMW.ac.uk, 14 July 1994 % Based on a port of the Maple randpoly function. % See RANDPOLY.TEX for documentation, and RANDPOLY.TST for examples. create!-package('(randpoly),'(contrib misc)); symbolic smacro procedure apply_c c; % Apply a coefficient generator function c that returns % a prefix form and convert it to standard quotient form. simp!* apply(c, nil); symbolic procedure apply_e e; % Apply an exponent generator function e % and check that it has returned an integer. if fixp(e := apply(e, nil)) then e else rederr "randpoly expons function must return an integer"; put('randpoly, 'simpfn, 'randpoly); flag('(randpoly), 'listargp); % allow single list argument symbolic procedure randpoly u; % Use: randpoly(vars, options) % vars: expression or list of expressions -- usually indeterminates % options: optional equations (of the form `option = value') % or keywords specifying properties: % % Option Use Default Value % ------ --- ------------- % coeffs = procedure to generate a coefficient rand(-99 .. 99) % expons = procedure to generate an exponent rand(6) % (must return an integer) % degree, deg or maxdeg = maximum total degree 5 % ord or mindeg = minimum total degree 0 (or 1) % (defaults to 1 if any "variable" is an equation) % terms = number of terms generated if sparse 6 % dense the polynomial is to be dense sparse % sparse the polynomial is to be sparse sparse begin scalar v, univar, c, e, trms, d, o, s, p, sublist; % Default values for options: %% c := rand {-99 .. 99}; % rand is a psopfn c := function(lambda(); random(199) - 99); d := 5; o := 0; trms := 6; s := 'sparse; % Evaluate arguments and process "variables": begin scalar wtl!*; % Ignore weights when evaluating args, including in revlis. v := car(u := revlis u); v := if eqcar(v, 'list) then cdr v else << univar := t; v . nil >>; v := for each vv in v collect begin scalar tmpvar; if eqexpr vv then << vv := !*eqn2a vv; o := 1 >> else if kernp simp!* vv then return vv; tmpvar := gensym(); sublist := {'equal, tmpvar, vv} . sublist; return tmpvar end; if univar then v := car v end; % Process options: for each x in cdr u do if x eq 'dense or x eq 'sparse then s := x else if not (eqexpr x and (if cadr x eq 'coeffs and functionp caddr x then c := caddr x else if cadr x eq 'expons and functionp caddr x then e := caddr x else if cadr x memq '(degree deg maxdeg) and natnump caddr x then d := caddr x else if cadr x memq '(ord mindeg) and natnump caddr x then o := caddr x else if cadr x eq 'terms and natnump caddr x then trms := caddr x)) then typerr(x, "optional randpoly argument"); % Generate the random polynomial: p := nil ./ 1; if o <= d then if s eq 'sparse then if null e then for each x in rand!-mons!-sparse(v,trms,d,o,univar) do p := addsq(p, multsq(apply_c c, x ./ 1)) else if univar then for i := 1 : trms do p := addsq(p, multsq(apply_c c, !*kp2q(v, apply_e e))) else for i := 1 : trms do (if numr cc then p := addsq(p, << for each vv in v do cc := multsq(cc, !*kp2q(vv, apply_e e)); cc >> )) where cc = apply_c c else % s eq 'dense if univar then << p := apply_c c; if o > 0 then p := multsq(p, mksq(v,o)); for i := o + 1 : d do p := addsq(p, multsq(apply_c c, mksq(v,i))) >> else for each x in rand!-mons!-dense(v,d,o) do p := addsq(p, multsq(apply_c c, x ./ 1)); return % Make any necessary substitutions for temporary variables: if sublist then simp!* subeval append(sublist, {mk!*sq p}) else p end; symbolic procedure functionp f; % Returns t if f can be applied as a function. getd f or eqcar(f,'lambda); symbolic procedure natnump n; % Returns t if n is a natural number. fixp n and n >= 0; symbolic smacro procedure kp2f(k, p); % k : unique kernel, p : natural number > 0 % Returns k^p as a standard form, taking account of % both asymptotic let rules and weightings. numr mksq(k, p); symbolic procedure !*kp2f(k, p); % k : unique kernel, p : natural number % Returns k^p as a standard form, taking account of % both asymptotic let rules and weightings. if p > 0 then kp2f(k, p) else 1; symbolic procedure !*kp2q(k, p); % k : unique kernel, p : any integer % Returns k^p as a standard quotient, taking account of % both asymptotic let rules and weightings. if p > 0 then mksq(k, p) else if zerop p then 1 ./ 1 else % Is this the right behaviour? % cf. part of procedure invsq in POLY.RED if null numr(k := mksq(k, -p)) then rederr "Zero divisor" else revpr k; symbolic procedure rand!-mons!-dense(v, d, o); % v : list of variables, % d : max total degree, o : min total degree. % Recursively returns a dense list of multivariate monomials % with total degree in [o, d] as STANDARD FORMS. begin scalar v_1; v_1 := car v; v := cdr v; return if null v then % single variable (if o > 0 then kp2f(v_1,o) else 1) . for i := o + 1 : d collect kp2f(v_1,i) else append( rand!-mons!-dense(v, d, o), for i := 1 : d join (for each x in rand!-mons!-dense(v, d - i, max(0, o - i)) collect multf(v_1!^i, x)) where v_1!^i = kp2f(v_1,i) ) end; symbolic procedure rand!-mons!-sparse(v, trms, d, o, univar); % v : (list of) variable(s), trms: number of terms, % d : max total degree, o : min total degree. % Returns a sparse list of at most trms monomials % with total degree in [o, d] as STANDARD FORMS. begin scalar n, v_1, maxtrms, otrms, s; if univar then maxtrms := d + 1 - o else << n := length v; v_1 := car v; otrms := if zerop o then 0 else binomial(n + o - 1, n); % max # terms to degree o-1: maxtrms := binomial(n + d, n) - otrms % max # terms in poly := max # terms to degree d - otrms >>; % Choose a random subset of the maxtrms terms by "index": s := rand!-comb(maxtrms, min(maxtrms,trms)); return if univar then for each ss in s collect !*kp2f(v, ss + o) else for each ss in s collect begin scalar p; p := 1; % Convert term "index" to exponent vector: ss := nil . inttovec(ss + otrms, n); for each vv in v do p := multf(!*kp2f(vv, car(ss := cdr ss)), p); return p end end; % Support procedures for randpoly % =============================== global '(!_binomialk !_binomialb !_binomialn); % binomial in the specfn package is implemented as an algebraic % operator, and I suspect is not very efficient. It will not clash % with the following implementation, which is about 50% faster on % my PC for binomial(200, 100) (with its caching disabled); symbolic procedure binomial(n, k); % Returns the binomial coefficient ASSUMING n, k integer >= 0. begin scalar n1, b; % Global !_BinomialK, !_BinomialB, !_BinomialN if k = 0 then return 1; if n < 2*k then return binomial(n,n-k); n1 := n+1; if !_binomialn = n then << % Partial result exits ... b := !_binomialb; if !_binomialk <= k then for i := !_binomialk+1 : k do b := quotient((n1-i)*b,i) else for i := !_binomialk step -1 until k+1 do b := quotient(i*b,n1-i) >> else << % First binomial computation b := 1; for i := 1 : k do b := quotient((n1-i)*b,i); !_binomialn := n >>; !_binomialk := k; return !_binomialb := b end; symbolic procedure rand!-comb(n, m); % Returns a list containing a random combination of m of the % first n NATURAL NUMBERS, ASSUMING integer n >= m >= 0. % (The values returned are 1 less than those % returned by the Maple randcomb function.) if m = n then for i := 0 : m - 1 collect i else begin scalar s; if n - m < m then begin scalar r; r := rand!-comb(n, n - m); for rr := 0 : n - 1 do if not(rr memq r) then s := rr . s end else for i := 0 : m - 1 do begin scalar rr; while (rr := random n) memq s do; % nothing s := rr . s end; return s end; symbolic procedure inttovec(m, n); % Returns the m'th (in lexicographic order) list of n % non-negative integers, ASSUMING integer m >= 0, n > 0. inttovec1(n, inttovec!-solve(n,m)); symbolic procedure inttovec1(n, dm); % n > 0 : integer; dm : dotted pair, d . m', where % d = norm of vector in N^n and m' = index of its tail in N^{n-1}. % Returns list representing vector in N^n, constructed recursively. % First vector element v_1 satisfies v_1 = d - norm of tail vector. if n = 1 then car dm . nil else ( (car dm - car dm1) . inttovec1(n - 1, dm1) ) where dm1 = inttovec!-solve(n - 1, cdr dm); % dotted pair symbolic procedure inttovec!-solve(n, m); % n > 0, m >= 0 : integer % Main subalgorithm to compute the vector in N^n with index m. % Returns as a dotted pair d . m' the norm (total degree) d in N^n % and the index m' of the tail sub-vector in N^{n-1}. % d is computed to satisfy ^{d-1+n}C_n <= m < ^{d+n}C_n, % where ^{d+n}C_n = number of n-vectors of norm d. if m = 0 or n = 1 then m . 0 else begin scalar d, c, cc; d := 0; cc := 1; % cc = ^{d+n}C_n repeat << c := cc; d := d + 1; % c = ^{d+n}C_n cc := quotient((n + d)*c, d); % cc = ^{d+1+n}C_n >> until cc > m; return d . (m - c) end; % Support for anonymous procedures (`proc's), % ========================================== % especially random number generators. % =================================== % Based partly on Maple's proc and rand function, and intended mainly % for use with the randpoly operator. % Interval code based on numeric.red and gnuplot.red by Herbert Melenk. % Create .. infix operator, avoiding warning if already defined. % (It is pre-defined via plot/plothook.sl at least in PSL-REDUCE.) newtok '( (!. !.) !*interval!*); precedence .., or; algebraic operator ..; put('!*interval!*, 'prtch, '! !.!.! ); put('rand, 'psopfn, 'rand); symbolic procedure rand u; % Returns a random number generator, and compiles it if COMP is on. % Optional second argument generates a named procedure. if null u or (cdr u and cddr u) then rederr "rand takes 1 or 2 arguments" else begin scalar fname, fn; if cdr u and not idp(fname := reval cadr u) then typerr(fname, "procedure name"); fn := if fixp(u := reval car u) and u > 0 then {'random, u} else if eqcar(u,'!*interval!*) then begin scalar a, b; if not(fixp(a := cadr u) and fixp(b := caddr u) and a<b) then rederr "rand range argument a .. b must have integer a,b with a < b"; return if zerop a then {'random, b + 1} else {'plus2, a, {'random, b - a + 1}} end else typerr(u, "integer or integer range"); fn := {'lambda, nil, fn}; %% putd compiles the function if !*comp is non-nil. return if fname then << flag({fname}, 'opfn); putd(fname, 'expr, fn) >> else if !*comp then putd(gensym(), 'expr, fn) else fn end; % Redefine the algebraic-mode random operator. remflag('(random), 'opfn); put('random, 'psopfn, 'evalrandom); symbolic procedure evalrandom u; % More flexible interface to the random function. if null u or cdr u then rederr "random takes a single argument" else if eqcar(u := reval car u,'!*interval!*) then begin scalar a, b; if not(fixp(a := cadr u) and fixp(b := caddr u) and a < b) then rederr "random range argument a .. b must have integer a,b with a < b"; return if zerop a then random(b + 1) else a + random(b - a + 1) end else if fixp u and u > 0 then random u % N.B. random also makes this check, but does not accept a range else typerr(u, "integer or integer range"); % Proc turns its argument expression into a lambda expression % and compiles it if COMP is ON. Provides a general version of rand. % Some such mechanism is necessary to prevent expressions containing % random number generators being evaluated too early. put('proc, 'psopfn, 'proc); symbolic procedure proc u; % Returns an anonymous procedure definition, % compiled if COMP is ON. if null u then rederr "proc requires at least a body argument" else << % aeval!* necessary instead of aeval here to avoid caching % and hence loss of possible randomness within loops: u := {'lambda, nil, {'aeval!*, mkquote car u}}; if !*comp then putd(gensym(), 'expr, u) else u >>; % User access procedures to evaluate and display procs etc. % ======================================================== % Not necessary -- provided only for test purposes and curious users! put('evalproc, 'psopfn, 'evalproc); symbolic procedure evalproc r; % r : proc, arg_1, ..., arg_n; args optional % Evaluates a proc applied to the subsequent arguments. apply(getproc car r, revlis cdr r); put('showproc, 'psopfn, 'showproc); symbolic procedure showproc r; % Displays a proc. (if codep rr then rederr "Argument is a compiled proc -- cannot display" else << terpri(); rprint subst('plus, 'plus2, rr); >>) where rr = getproc car r; symbolic procedure getproc r; % Support procedure: get a proc body. ( if idp r then getfnbody r or ( (r := get(r, 'avalue)) and car r eq 'scalar and eqcar(r := cadr r, 'lambda) and r ) or getfnbody r else if pairp r then if car r eq 'lambda then r % share variable else if eqcar(r := ( (if x then apply(x, {cdr r})) where x = get(car r, 'psopfn) ), 'lambda) then r ) or rederr "Argument is not a proc"; symbolic procedure getfnbody r; (r := getd r) and cdr r; endmodule; end;