random number generation

49 views
Skip to first unread message

Craig Citro

unread,
Feb 25, 2007, 8:16:58 PM2/25/07
to sage-...@googlegroups.com
Hey all,

So I tried to generate a random polynomial today, and ran into some trouble. Here's what I did:

sage: R.<x> = ZZ['x']
sage: R.random_element(3)
<sage crashes>

I traced back the problem, and it's not clear what the right fix is. So R.random_element makes a list of the appropriate length and calls ZZ.random_element(0) to fill it up. In the comments, it clearly explains why 0: it assumes that this is the default argument to pass to get a choice that is as spread across the ring as possible. Now, the problem is that ZZ.random_element(0) calls GMP's mpz_urandomm() function, which throws a division by zero exception when you call it with 0; it only knows how to generate a random element between two bounds (at least, according to the comments in the source). So I'm not sure which of the following things is the "right" fix for this:

1) Fix ZZ.random_element so that it calls off to mpz_urandomm with a large bound; if so, I'm not sure what bound to use.

2) Change R.random_element so that it doesn't try to pass off a 0 to ZZ.random_element , which seems like the wrong solution -- this code is in PolynomialRing.random_element, which really shouldn't know in general how to ask for a random element in its base_ring. However, this means that we should document somewhere that in any ring, random_element(0) is the request for "as general a random element as possible."

3) Decide to think about it later, and throw an exception. This was apparently GMP's solution. ;)

I think (1) is the right solution, but I'm not sure what bound to use. But maybe one of the other choices is better? What does everyone else think?

-cc

didier deshommes

unread,
Mar 1, 2007, 10:30:11 PM3/1/07
to sage-...@googlegroups.com
On 2/25/07, Craig Citro <craig...@gmail.com> wrote:
> Hey all,
>
> So I tried to generate a random polynomial today, and ran into some trouble.
> Here's what I did:
>
> sage: R.<x> = ZZ['x']
> sage: R.random_element(3)
> <sage crashes>

That is a nice edge case. I would say that or you just return 0
everytime. Other rings seem to do that:
{{{
sage: RR.random_element(0)
0.000000000000000
sage: QQ.random_element(0)
0
sage: RDF.random_element(0)
0.143951483848
}}}

didier

William Stein

unread,
Mar 2, 2007, 12:57:29 PM3/2/07
to sage-...@googlegroups.com
On 3/1/07, didier deshommes <dfde...@gmail.com> wrote:
> On 2/25/07, Craig Citro <craig...@gmail.com> wrote:
> > Hey all,
> >
> > So I tried to generate a random polynomial today, and ran into some trouble.
> > Here's what I did:
> >
> > sage: R.<x> = ZZ['x']
> > sage: R.random_element(3)
> > <sage crashes>
>
> That is a nice edge case. I would say that or you just return 0
> everytime. Other rings seem to do that:
> {{{
> sage: RR.random_element(0)
> 0.000000000000000
> sage: QQ.random_element(0)
> 0
> sage: RDF.random_element(0)
> 0.143951483848
> }}}

> > I traced back the problem, and it's not clear what the right fix is. So


> > R.random_element makes a list of the appropriate length and calls
> > ZZ.random_element(0) to fill it up. In the comments, it clearly explains why

I've fixed this for sage > 2.2. The patch is attached in case you're
interested.

3251.patch

Robert Bradshaw

unread,
Mar 2, 2007, 5:15:55 PM3/2/07
to sage-...@googlegroups.com
It's always bugged me that the default distribution for integers (and
rationals) is just a uniform distribution over some small range. What
if instead we chose the distribution ZZ.random_element() = floor(1/r)
where r is uniformly distributed in (-1,1). Then P(n) = 1 / (2 |n| (|
n| + 1)) for all n in Z-{0}. This gives mostly small numbers with the
occasional large ones thrown in at ever decreasing probabilities.

A random rational could then be the ratio of two such integers.

- Robert

> <3251.patch>

William Stein

unread,
Mar 2, 2007, 5:17:38 PM3/2/07
to sage-...@googlegroups.com
Sure, go for it!


--
William Stein
Associate Professor of Mathematics
University of Washington

Reply all
Reply to author
Forward
0 new messages