With probability p the routine takes the Alexander dual of I; the default value of p is .05, and it can be set with the option AlexanderProbility.
Otherwise uses the Metropolis algorithm to produce a random walk on the space of square-free ideals. Note that there are a LOT of square-free ideals; these are the Dedekind numbers, and the sequence (with 1,2,3,4,5,6,7,8 variables) begins 3,6,20,168,7581, 7828354, 2414682040998, 56130437228687557907788. (see the Online Encyclopedia of Integer Sequences for more information). Given I in a polynomial ring S, we make a list ISocgens of the square-free minimal monomial generators of the socle of S/(squares+I) and a list of minimal generators Igens of I. A candidate "next" ideal J is formed as follows: We choose randomly from the union of these lists; if a socle element is chosen, it's added to I; if a minimal generator is chosen, it's replaced by the square-free part of the maximal ideal times it. the chance of making the given move is then 1/(#ISocgens+#Igens), and the chance of making the move back would be the similar quantity for J, so we make the move or not depending on whether random RR < (nJ+nSocJ)/(nI+nSocI) or not; here random RR is a random number in [0,1].
i1 : setRandomSeed(currentTime()) o1 = 1638138489 |
i2 : S=ZZ/2[vars(0..3)] o2 = S o2 : PolynomialRing |
i3 : J = monomialIdeal"ab,ad, bcd" o3 = monomialIdeal (a*b, a*d, b*c*d) o3 : MonomialIdeal of S |
i4 : randomSquareFreeStep J
o4 = {monomialIdeal (a*b*c, a*d, b*c*d), {a*b*c, a*d, b*c*d}, {c*d, b*d, b*c,
------------------------------------------------------------------------
a*c, a*b}}
o4 : List
|
With 4 variables and 168 possible monomial ideals, a run of 5000 takes less than 6 seconds on a reasonably fast machine. With 10 variables a run of 1000 takes about 2 seconds.
i5 : setRandomSeed(1) o5 = 1 |
i6 : rsfs = randomSquareFreeStep o6 = randomSquareFreeStep o6 : MethodFunctionWithOptions |
i7 : J = monomialIdeal 0_S o7 = monomialIdeal () o7 : MonomialIdeal of S |
i8 : time T=tally for t from 1 to 5000 list first (J=rsfs(J,AlexanderProbability => .01));
-- used 6.75788 seconds
|
i9 : #T o9 = 168 |
i10 : T
o10 = Tally{monomialIdeal () => 45 }
monomialIdeal (a*b*c, a*b*d) => 33
monomialIdeal (a*b*c, a*b*d, a*c*d) => 22
monomialIdeal (a*b*c, a*b*d, a*c*d, b*c*d) => 26
monomialIdeal (a*b*c, a*b*d, b*c*d) => 37
monomialIdeal (a*b*c, a*b*d, c*d) => 31
monomialIdeal (a*b*c, a*c*d) => 25
monomialIdeal (a*b*c, a*c*d, b*c*d) => 25
monomialIdeal (a*b*c, a*d) => 31
monomialIdeal (a*b*c, a*d, b*c*d) => 33
monomialIdeal (a*b*c, a*d, b*d) => 36
monomialIdeal (a*b*c, a*d, b*d, c*d) => 35
monomialIdeal (a*b*c, a*d, c*d) => 47
monomialIdeal (a*b*c, b*c*d) => 30
monomialIdeal (a*b*c, b*d) => 39
monomialIdeal (a*b*c, b*d, a*c*d) => 36
monomialIdeal (a*b*c, b*d, c*d) => 39
monomialIdeal (a*b*c, c*d) => 24
monomialIdeal (a*b*c, d) => 20
monomialIdeal (a*b*d, a*c*d) => 35
monomialIdeal (a*b*d, a*c*d, b*c*d) => 32
monomialIdeal (a*b*d, b*c*d) => 33
monomialIdeal (a*b*d, c*d) => 45
monomialIdeal (a*b, a*c) => 37
monomialIdeal (a*b, a*c*d) => 42
monomialIdeal (a*b, a*c*d, b*c*d) => 31
monomialIdeal (a*b, a*c, a*d) => 29
monomialIdeal (a*b, a*c, a*d, b*c*d) => 38
monomialIdeal (a*b, a*c, a*d, b*d) => 28
monomialIdeal (a*b, a*c, a*d, b*d, c*d) => 22
monomialIdeal (a*b, a*c, a*d, c*d) => 30
monomialIdeal (a*b, a*c, b*c) => 37
monomialIdeal (a*b, a*c, b*c*d) => 35
monomialIdeal (a*b, a*c, b*c, a*d) => 26
monomialIdeal (a*b, a*c, b*c, a*d, b*d) => 33
monomialIdeal (a*b, a*c, b*c, a*d, b*d, c*d) => 31
monomialIdeal (a*b, a*c, b*c, a*d, c*d) => 28
monomialIdeal (a*b, a*c, b*c, b*d) => 29
monomialIdeal (a*b, a*c, b*c, b*d, c*d) => 23
monomialIdeal (a*b, a*c, b*c, c*d) => 26
monomialIdeal (a*b, a*c, b*c, d) => 29
monomialIdeal (a*b, a*c, b*d) => 28
monomialIdeal (a*b, a*c, b*d, c*d) => 25
monomialIdeal (a*b, a*c, c*d) => 20
monomialIdeal (a*b, a*c, d) => 27
monomialIdeal (a*b, a*d) => 30
monomialIdeal (a*b, a*d, b*c*d) => 31
monomialIdeal (a*b, a*d, b*d) => 30
monomialIdeal (a*b, a*d, b*d, c*d) => 25
monomialIdeal (a*b, a*d, c*d) => 34
monomialIdeal (a*b, b*c) => 25
monomialIdeal (a*b, b*c*d) => 30
monomialIdeal (a*b, b*c, a*c*d) => 30
monomialIdeal (a*b, b*c, a*d) => 27
monomialIdeal (a*b, b*c, a*d, b*d) => 32
monomialIdeal (a*b, b*c, a*d, b*d, c*d) => 32
monomialIdeal (a*b, b*c, a*d, c*d) => 32
monomialIdeal (a*b, b*c, b*d) => 42
monomialIdeal (a*b, b*c, b*d, a*c*d) => 35
monomialIdeal (a*b, b*c, b*d, c*d) => 31
monomialIdeal (a*b, b*c, c*d) => 24
monomialIdeal (a*b, b*c, d) => 34
monomialIdeal (a*b, b*d) => 41
monomialIdeal (a*b, b*d, a*c*d) => 43
monomialIdeal (a*b, b*d, c*d) => 31
monomialIdeal (a*b, c) => 32
monomialIdeal (a*b, c*d) => 26
monomialIdeal (a*b, c, a*d) => 36
monomialIdeal (a*b, c, a*d, b*d) => 36
monomialIdeal (a*b, c, b*d) => 38
monomialIdeal (a*b, c, d) => 38
monomialIdeal (a*b, d) => 16
monomialIdeal (a*c*d, b*c*d) => 27
monomialIdeal (a*c, a*b*d) => 30
monomialIdeal (a*c, a*b*d, b*c*d) => 34
monomialIdeal (a*c, a*b*d, c*d) => 26
monomialIdeal (a*c, a*d) => 32
monomialIdeal (a*c, a*d, b*c*d) => 39
monomialIdeal (a*c, a*d, b*d) => 16
monomialIdeal (a*c, a*d, b*d, c*d) => 23
monomialIdeal (a*c, a*d, c*d) => 27
monomialIdeal (a*c, b*c) => 21
monomialIdeal (a*c, b*c*d) => 24
monomialIdeal (a*c, b*c, a*b*d) => 24
monomialIdeal (a*c, b*c, a*b*d, c*d) => 25
monomialIdeal (a*c, b*c, a*d) => 37
monomialIdeal (a*c, b*c, a*d, b*d) => 28
monomialIdeal (a*c, b*c, a*d, b*d, c*d) => 35
monomialIdeal (a*c, b*c, a*d, c*d) => 28
monomialIdeal (a*c, b*c, b*d) => 28
monomialIdeal (a*c, b*c, b*d, c*d) => 18
monomialIdeal (a*c, b*c, c*d) => 17
monomialIdeal (a*c, b*c, d) => 26
monomialIdeal (a*c, b*d) => 30
monomialIdeal (a*c, b*d, c*d) => 23
monomialIdeal (a*c, c*d) => 20
monomialIdeal (a*c, d) => 27
monomialIdeal (a*d, b*c*d) => 30
monomialIdeal (a*d, b*d) => 32
monomialIdeal (a*d, b*d, c*d) => 49
monomialIdeal (a*d, c*d) => 52
monomialIdeal (a, b) => 27
monomialIdeal (a, b*c) => 13
monomialIdeal (a, b*c*d) => 17
monomialIdeal (a, b*c, b*d) => 11
monomialIdeal (a, b*c, b*d, c*d) => 21
monomialIdeal (a, b*c, c*d) => 12
monomialIdeal (a, b*c, d) => 26
monomialIdeal (a, b*d) => 10
monomialIdeal (a, b*d, c*d) => 15
monomialIdeal (a, b, c) => 38
monomialIdeal (a, b, c*d) => 20
monomialIdeal (a, b, c, d) => 31
monomialIdeal (a, b, d) => 34
monomialIdeal (a, c) => 30
monomialIdeal (a, c*d) => 21
monomialIdeal (a, c, b*d) => 26
monomialIdeal (a, c, d) => 37
monomialIdeal (a, d) => 15
monomialIdeal (b*c, a*b*d) => 26
monomialIdeal (b*c, a*b*d, a*c*d) => 26
monomialIdeal (b*c, a*b*d, c*d) => 29
monomialIdeal (b*c, a*c*d) => 21
monomialIdeal (b*c, a*d) => 29
monomialIdeal (b*c, a*d, b*d) => 31
monomialIdeal (b*c, a*d, b*d, c*d) => 29
monomialIdeal (b*c, a*d, c*d) => 47
monomialIdeal (b*c, b*d) => 50
monomialIdeal (b*c, b*d, a*c*d) => 35
monomialIdeal (b*c, b*d, c*d) => 30
monomialIdeal (b*c, c*d) => 14
monomialIdeal (b*c, d) => 23
monomialIdeal (b*d, a*c*d) => 33
monomialIdeal (b*d, c*d) => 42
monomialIdeal (b, a*c) => 35
monomialIdeal (b, a*c*d) => 43
monomialIdeal (b, a*c, a*d) => 39
monomialIdeal (b, a*c, a*d, c*d) => 27
monomialIdeal (b, a*c, c*d) => 22
monomialIdeal (b, a*c, d) => 30
monomialIdeal (b, a*d) => 37
monomialIdeal (b, a*d, c*d) => 32
monomialIdeal (b, c) => 26
monomialIdeal (b, c*d) => 26
monomialIdeal (b, c, a*d) => 35
monomialIdeal (b, c, d) => 34
monomialIdeal (b, d) => 33
monomialIdeal (c, a*b*d) => 28
monomialIdeal (c, a*d) => 35
monomialIdeal (c, a*d, b*d) => 25
monomialIdeal (c, b*d) => 23
monomialIdeal (c, d) => 30
monomialIdeal 1 => 21
monomialIdeal a => 10
monomialIdeal b => 39
monomialIdeal c => 20
monomialIdeal d => 40
monomialIdeal(a*b) => 35
monomialIdeal(a*b*c) => 26
monomialIdeal(a*b*c*d) => 32
monomialIdeal(a*b*d) => 32
monomialIdeal(a*c) => 26
monomialIdeal(a*c*d) => 38
monomialIdeal(a*d) => 42
monomialIdeal(b*c) => 14
monomialIdeal(b*c*d) => 29
monomialIdeal(b*d) => 22
monomialIdeal(c*d) => 43
o10 : Tally
|
i11 : J
o11 = {monomialIdeal (a*b, a*c, b*c, c*d), {a*b, a*c, b*c, c*d}, {b*d, a*d,
-----------------------------------------------------------------------
c}}
o11 : List
|
The object randomSquareFreeStep is a method function with options.