Hi,
after realizing that Steffen and I share an office at RHUL I discussed this
thing with him today for some time.
The improved implementation (#980) still doesn't produce random polynomials
for two reasons: The code does not seem to produce monomials of degree up to
$d$ uniformly at random.
deg = degree
for i in range(n):
exp = ZZ.random_element(0,deg+1)
if exp <= deg:
exps.append(exp)
deg -= exp
else:
exps.append(0)
E.g. in this implementation the probability that a (0,0,0....,0,0,0,0)
exponent tuple is produced is (1/(d+1))^n because we have to hit 0 out of d+1
possibilities n times.
However, if we write down a list M of all possible monomials of degree <= d
this L has length
#M = \sum_{i=0}^d ( binomial(n + i -1, i ) )
and the probability to hit the constant monomial is (1/#M). Plugging in some
numbers: n = 10, d = 2
implemented: (1 / 3)^10
uniformly at random: 1 / ( 1 + 10 + 55)
To fix this we need to consider choosing monomials of a degree D at random
which can be done by finding a bijection between the integers up to some N =
binomial( n + d -1 , d ) and the set of compositions. Mike (mhansen)
provided this code to do it:
def random_monomial(n, degree):
total = binomial(n+degree-1,degree)
random_index = randint(0, total)
comb = choose_nk.from_rank(random_index, n+degree-1, n-1)
monomial = [ comb[0] ]
monomial += map(lambda i: comb[i+1]-comb[i]-1, range(n-2))
monomial.append( n+degree-1-comb[-1]-1 )
and this explanation:
http://www.cse.ucsd.edu/classes/sp07/cse21/Ron_Notes/starbin.pdf
Side note: This code does not correctly for me as it can run into an infinite
loop which probably can be tracked down to _comb_largest with a parameter 'w'
< 0.
This construction is random if the random number generator ("randint") is
random. Btw. how random is randint?
Now we need to choose a degree d at random. We cannot do this uniformly at
random because the classes (associated to the degree) are not of the same
size (e.g., 1, 10, 55 above) We don't want to choose classes uniformly but
elements and again mhansen beat me to an implementation:
def random_monomial_less_than_degree(n, degree):
#Get the counts of the total number of monomials
#of degree d
counts = [1] #degree 0
for d in range(1, degree):
counts.append(binomial(n+d-1, d))
total = sum(counts)
#Select a random one
random_index = randint(0, total-1)
#Figure out which degree it corresponds to
d = 0
while random_index >= counts[d]:
random_index -= counts[d]
d += 1
#Generate the corresponding monomial
comb = choose_nk.from_rank(random_index, n+d-1, n-1)
monomial = [ comb[0] ]
monomial += map(lambda i: comb[i+1]-comb[i]-1, range(n-2))
monomial.append( n+d-1-comb[-1]-1 )
return monomial
However, an implementation without the need to explicitly construct the counts
array would be good.
With this in place one can start thinking about constructing random
polynomials from random monomials. For this we don't have an elegant solution
yet. We (=Steffen and me) distinguish two cases:
Let #T be the number of terms requested.
If $MM is the number of all possible monomials of all allowed degree <= d and
#T < $MM/2 then we generate #T _distinct_ monomials as above and assign
random coefficients (possibly zero). The overhead of requesting distinct
monomials should be at most 2 because if we are looking for 1/2 of the search
space, every second choice should be a double.
If #T > $MM/2 we generate all possible monomials, pick #T and assign
coefficients. This should be more effective for this case because we don't
generate any doubles (but memory is more precious than cpu cycles, so we
might have to move that barrier up a bit).
Does this sound like a sensible approach?
Martin
--
name: Martin Albrecht
_pgp: http://pgp.mit.edu:11371/pks/lookup?op=get&search=0x8EF0DC99
_www: http://www.informatik.uni-bremen.de/~malb
_jab: [EMAIL PROTECTED]
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/
-~----------~----~----~----~------~----~------~--~---
