From dc4c1c0c6b501a656827441791abb5cdc09b1f19 Mon Sep 17 00:00:00 2001 From: Doug Torrance Date: Fri, 13 Sep 2024 22:28:18 -0400 Subject: [PATCH 01/14] Simplify rawSetRandomQQ Have it call rawSetRandomQQ to de-duplicate code --- M2/Macaulay2/e/interface/random.cpp | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/M2/Macaulay2/e/interface/random.cpp b/M2/Macaulay2/e/interface/random.cpp index 2836c47c48..2c37b332cb 100644 --- a/M2/Macaulay2/e/interface/random.cpp +++ b/M2/Macaulay2/e/interface/random.cpp @@ -67,31 +67,26 @@ gmp_ZZ rawRandomInteger(gmp_ZZ maxN) return result; } -gmp_QQ rawRandomQQ(gmp_ZZ height) +void rawSetRandomQQ(mpq_ptr result, gmp_ZZ height) /* returns random a/b, where 1 <= b <= height, 1 <= a <= height */ /* if height is the null pointer, use the default height */ { - mpq_ptr result = getmemstructtype(mpq_ptr); - mpq_init(result); if (height == nullptr) height = maxHeight; mpz_urandomm(mpq_numref(result), state, height); mpz_urandomm(mpq_denref(result), state, height); mpz_add_ui(mpq_numref(result), mpq_numref(result), 1); mpz_add_ui(mpq_denref(result), mpq_denref(result), 1); mpq_canonicalize(result); - return moveTo_gmpQQ(result); } -void rawSetRandomQQ(mpq_ptr result, gmp_ZZ height) +gmp_QQ rawRandomQQ(gmp_ZZ height) /* returns random a/b, where 1 <= b <= height, 1 <= a <= height */ /* if height is the null pointer, use the default height */ { - if (height == nullptr) height = maxHeight; - mpz_urandomm(mpq_numref(result), state, height); - mpz_urandomm(mpq_denref(result), state, height); - mpz_add_ui(mpq_numref(result), mpq_numref(result), 1); - mpz_add_ui(mpq_denref(result), mpq_denref(result), 1); - mpq_canonicalize(result); + mpq_ptr result = getmemstructtype(mpq_ptr); + mpq_init(result); + rawSetRandomQQ(result, height); + return moveTo_gmpQQ(result); } gmp_RR rawRandomRR(unsigned long precision) From 59cf05486d5c8643d0fedb2a27422592fe2f254d Mon Sep 17 00:00:00 2001 From: Doug Torrance Date: Sat, 14 Sep 2024 13:29:31 -0400 Subject: [PATCH 02/14] Add routine for computing Farey approximations to engine This gives us the closest rational number to a given real number with bounded denominator. It's available at top level using the unexported rawFareyApproximation. --- M2/Macaulay2/d/interface.dd | 13 +++++ M2/Macaulay2/e/interface/random.cpp | 78 +++++++++++++++++++++++++++++ M2/Macaulay2/e/interface/random.h | 6 +++ 3 files changed, 97 insertions(+) diff --git a/M2/Macaulay2/d/interface.dd b/M2/Macaulay2/d/interface.dd index 921c6c7ddc..0785ae007a 100644 --- a/M2/Macaulay2/d/interface.dd +++ b/M2/Macaulay2/d/interface.dd @@ -51,6 +51,19 @@ export rawRandomZZ(e:Expr):Expr := ( is maxN:ZZcell do toExpr(Ccode(ZZ, "rawRandomInteger(", maxN.v, ")")) else WrongArgZZ()); setupfun("rawRandomZZ",rawRandomZZ); +export rawFareyApproximation(e:Expr):Expr := ( + when e + is a:Sequence do ( + if length(a) != 2 then return WrongNumArgs(2); + when a.0 + is x:RRcell do ( + when a.1 + is hgt:ZZcell do toExpr(Ccode(QQ, + "rawFareyApproximation(", x.v, ", ", hgt.v, ")")) + else WrongArgZZ(2)) + else WrongArgRR(1)) + else WrongNumArgs(2)); +setupfun("rawFareyApproximation", rawFareyApproximation); export rawRandomQQ(e:Expr):Expr := ( when e is Nothing do toExpr(Ccode(QQ, "rawRandomQQ(", "0)")) diff --git a/M2/Macaulay2/e/interface/random.cpp b/M2/Macaulay2/e/interface/random.cpp index 2c37b332cb..e1a02c2a41 100644 --- a/M2/Macaulay2/e/interface/random.cpp +++ b/M2/Macaulay2/e/interface/random.cpp @@ -67,6 +67,84 @@ gmp_ZZ rawRandomInteger(gmp_ZZ maxN) return result; } +void rawSetFareyApproximation(mpq_ptr result, gmp_RR x, gmp_ZZ height) +/* sets result = the nearest rational to x w/ denominator <= height */ +{ + int sgn; + mpfr_t fracpart, tmp1, tmp2; + mpz_t intpart, a, b, c, d, q, p; + mpq_t tmp3; + + sgn = mpfr_sgn(x); + mpfr_init2(fracpart, mpfr_get_prec(x)); + mpfr_abs(fracpart, x, MPFR_RNDN); + mpz_init(intpart); + mpfr_get_z(intpart, fracpart, MPFR_RNDD); + mpfr_frac(fracpart, fracpart, MPFR_RNDN); + + mpfr_init2(tmp1, mpfr_get_prec(x)); + mpfr_init2(tmp2, mpfr_get_prec(x)); + mpz_init_set_ui(a, 0); + mpz_init_set_ui(b, 1); + mpz_init_set_ui(c, 1); + mpz_init_set_ui(d, 1); + mpz_init_set_ui(p, 1); + mpz_init_set_ui(q, 2); + mpq_init(tmp3); + + while (mpz_cmp(q, height) <= 0) { + mpfr_mul_z(tmp1, fracpart, q, MPFR_RNDN); + if (mpfr_cmp_z(tmp1, p) <= 0) { + mpz_set(c, p); + mpz_set(d, q); + } else { + mpz_set(a, p); + mpz_set(b, q); + } + mpz_add(p, a, c); + mpz_add(q, b, d); + } + + // tmp1 = fracpart - a/b + mpfr_set_z(tmp1, a, MPFR_RNDN); + mpfr_neg(tmp1, tmp1, MPFR_RNDN); + mpfr_div_z(tmp1, tmp1, b, MPFR_RNDN); + mpfr_add(tmp1, tmp1, fracpart, MPFR_RNDN); + + // tmp2 = c/d - fracpart + mpfr_set_z(tmp2, c, MPFR_RNDN); + mpfr_div_z(tmp2, tmp2, d, MPFR_RNDN); + mpfr_sub(tmp2, tmp2, fracpart, MPFR_RNDN); + + if (mpfr_cmp(tmp1, tmp2) <= 0) { + mpq_set_z(result, a); + mpq_set_z(tmp3, b); + } else { + mpq_set_z(result, c); + mpq_set_z(tmp3, d); + } + + mpq_div(result, result, tmp3); + mpq_set_z(tmp3, intpart); + mpq_add(result, result, tmp3); + mpq_set_si(tmp3, sgn, 1); + mpq_mul(result, result, tmp3); + mpq_canonicalize(result); + + mpz_clears(intpart, a, b, c, d, p, q, nullptr); + mpfr_clears(fracpart, tmp1, tmp2, nullptr); + mpq_clear(tmp3); +} + +gmp_QQ rawFareyApproximation(gmp_RR x, gmp_ZZ height) +/* returns the nearest rational to x w/ denominator <= height */ +{ + mpq_ptr result = getmemstructtype(mpq_ptr); + mpq_init(result); + rawSetFareyApproximation(result, x, height); + return moveTo_gmpQQ(result); +} + void rawSetRandomQQ(mpq_ptr result, gmp_ZZ height) /* returns random a/b, where 1 <= b <= height, 1 <= a <= height */ /* if height is the null pointer, use the default height */ diff --git a/M2/Macaulay2/e/interface/random.h b/M2/Macaulay2/e/interface/random.h index dda424f6b6..b7ba77409c 100644 --- a/M2/Macaulay2/e/interface/random.h +++ b/M2/Macaulay2/e/interface/random.h @@ -28,6 +28,12 @@ int32_t rawRandomInt(int32_t max); gmp_ZZ rawRandomInteger(gmp_ZZ maxN); /* if height is the null pointer, use the default height */ +void rawSetFareyApproximation(mpq_ptr result, gmp_RR x, gmp_ZZ height); +/* sets result = the nearest rational to x w/ denominator <= height */ + +gmp_QQ rawFareyApproximation(gmp_RR x, gmp_ZZ height); +/* returns the nearest rational to x w/ denominator <= height */ + gmp_QQ rawRandomQQ(gmp_ZZ height); /* returns random a/b, where 1 <= b <= height, 1 <= a <= height */ /* if height is the null pointer, use the default height */ From 6fcb951181d4bec9f3d0b14618e0e24fc1d767a4 Mon Sep 17 00:00:00 2001 From: Doug Torrance Date: Sat, 14 Sep 2024 13:41:45 -0400 Subject: [PATCH 03/14] Update random(QQ) to use uniform distribution on [0, height] We then round the result to the nearest rational number with denominator bounded by the Height option using Farey approximation --- M2/Macaulay2/e/interface/random.cpp | 20 ++++++++++--------- M2/Macaulay2/e/interface/random.h | 8 ++++---- .../Macaulay2Doc/functions/random-doc.m2 | 6 ++++-- 3 files changed, 19 insertions(+), 15 deletions(-) diff --git a/M2/Macaulay2/e/interface/random.cpp b/M2/Macaulay2/e/interface/random.cpp index e1a02c2a41..05816170be 100644 --- a/M2/Macaulay2/e/interface/random.cpp +++ b/M2/Macaulay2/e/interface/random.cpp @@ -146,20 +146,22 @@ gmp_QQ rawFareyApproximation(gmp_RR x, gmp_ZZ height) } void rawSetRandomQQ(mpq_ptr result, gmp_ZZ height) -/* returns random a/b, where 1 <= b <= height, 1 <= a <= height */ -/* if height is the null pointer, use the default height */ +/* sets result = a sample from the uniform distribution on [0, height], */ +/* rounded to the nearest rational number with denominator bounded by height */ { + mpfr_t x; + if (height == nullptr) height = maxHeight; - mpz_urandomm(mpq_numref(result), state, height); - mpz_urandomm(mpq_denref(result), state, height); - mpz_add_ui(mpq_numref(result), mpq_numref(result), 1); - mpz_add_ui(mpq_denref(result), mpq_denref(result), 1); - mpq_canonicalize(result); + mpfr_init2(x, gmp_defaultPrecision); + mpfr_urandomb(x, state); + mpfr_mul_z(x, x, height, MPFR_RNDN); + rawSetFareyApproximation(result, x, height); + mpfr_clear(x); } gmp_QQ rawRandomQQ(gmp_ZZ height) -/* returns random a/b, where 1 <= b <= height, 1 <= a <= height */ -/* if height is the null pointer, use the default height */ +/* returns a sample from the uniform distribution on [0, height], */ +/* rounded to the nearest rational number with denominator bounded by height */ { mpq_ptr result = getmemstructtype(mpq_ptr); mpq_init(result); diff --git a/M2/Macaulay2/e/interface/random.h b/M2/Macaulay2/e/interface/random.h index b7ba77409c..18a923aa00 100644 --- a/M2/Macaulay2/e/interface/random.h +++ b/M2/Macaulay2/e/interface/random.h @@ -35,12 +35,12 @@ gmp_QQ rawFareyApproximation(gmp_RR x, gmp_ZZ height); /* returns the nearest rational to x w/ denominator <= height */ gmp_QQ rawRandomQQ(gmp_ZZ height); -/* returns random a/b, where 1 <= b <= height, 1 <= a <= height */ -/* if height is the null pointer, use the default height */ +/* returns a sample from the uniform distribution on [0, height], */ +/* rounded to the nearest rational number with denominator bounded by height */ void rawSetRandomQQ(mpq_ptr result, gmp_ZZ height); -/* sets result = random a/b, where 1 <= b <= height, 1 <= a <= height */ -/* if height is the null pointer, use the default height */ +/* sets result = a sample from the uniform distribution on [0, height], */ +/* rounded to the nearest rational number with denominator bounded by height */ gmp_RR rawRandomRR(unsigned long prec); /* returns a uniformly distributed random real with the given precision, in diff --git a/M2/Macaulay2/packages/Macaulay2Doc/functions/random-doc.m2 b/M2/Macaulay2/packages/Macaulay2Doc/functions/random-doc.m2 index c26c379701..6cf8261531 100644 --- a/M2/Macaulay2/packages/Macaulay2Doc/functions/random-doc.m2 +++ b/M2/Macaulay2/packages/Macaulay2Doc/functions/random-doc.m2 @@ -77,8 +77,10 @@ Node Description Text If the @TT "Height"@ option specifies a number @TT "h"@ and @TT "T"@ - is @TO "ZZ"@ then the integers returned are in the range @TT "[0, h)"@; - for @TO "QQ"@ the numerator and denominator are in the range @TT "[1, h]"@. + is @TO "ZZ"@ then the integers returned are in the range @TT "[0, h)"@. + If @TT "T"@ is @TO QQ@, then the results are drawn from the uniform + distribution on @TT "[0, h]"@ and rounded to the nearest rational number + with denominator bounded by @TT "h"@. Example random RR random CC_100 From 7a1f5b62ac9fdc2a13c0a8552d3a3e4314285ae5 Mon Sep 17 00:00:00 2001 From: Doug Torrance Date: Sat, 14 Sep 2024 14:07:52 -0400 Subject: [PATCH 04/14] Add unit tests for Farey approximation --- M2/Macaulay2/tests/normal/numbers.m2 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/M2/Macaulay2/tests/normal/numbers.m2 b/M2/Macaulay2/tests/normal/numbers.m2 index e5cd2738cd..f5397472be 100644 --- a/M2/Macaulay2/tests/normal/numbers.m2 +++ b/M2/Macaulay2/tests/normal/numbers.m2 @@ -1064,6 +1064,10 @@ assert( not isANumber ((1/0.-1/0.) + 1*ii) ) assert( not isANumber (1 + (ii/0.-ii/0.) ) ) assert( not isANumber ((1 + ii/0.) - ii/0. ) ) +importFrom(Core, "rawFareyApproximation") +assert( rawFareyApproximation(numeric pi, 10) == 22/7 ) +assert( rawFareyApproximation(numeric pi, 200) == 355/113 ) +assert( rawFareyApproximation(-pi, 10) == -22/7 ) -- Local Variables: -- compile-command: "make -C $M2BUILDDIR/Macaulay2/packages/Macaulay2Doc/test numbers.out" From 08088328149e8697feb29a9c4417e853ab912fa1 Mon Sep 17 00:00:00 2001 From: Doug Torrance Date: Sat, 14 Sep 2024 23:46:16 -0400 Subject: [PATCH 05/14] Define gmp_defaultPrecision for engine unit tests It's defined in the interpreter, which we don't link against here, and we need it for rawSetRandomQQ. --- M2/Macaulay2/e/unit-tests/testMain.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/M2/Macaulay2/e/unit-tests/testMain.cpp b/M2/Macaulay2/e/unit-tests/testMain.cpp index f102914ff9..b375f88ac1 100644 --- a/M2/Macaulay2/e/unit-tests/testMain.cpp +++ b/M2/Macaulay2/e/unit-tests/testMain.cpp @@ -2,6 +2,8 @@ #include #include +unsigned long gmp_defaultPrecision = 53; + int main(int argc, char **argv) { IM2_initialize(); From de2eb7537c749821f1325683e27e3d17a93f4d65 Mon Sep 17 00:00:00 2001 From: Doug Torrance Date: Wed, 18 Sep 2024 07:28:17 -0400 Subject: [PATCH 06/14] Reword division by zero error message in engine Previously misspelled the word "attempt". We use the same message as the same error in the interpreter. --- M2/Macaulay2/e/exceptions.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/M2/Macaulay2/e/exceptions.hpp b/M2/Macaulay2/e/exceptions.hpp index ec8ffa63a2..4471fb44cf 100644 --- a/M2/Macaulay2/e/exceptions.hpp +++ b/M2/Macaulay2/e/exceptions.hpp @@ -16,7 +16,7 @@ struct overflow_exception : public engine_error struct division_by_zero_error : public engine_error { explicit division_by_zero_error(const std::string &msg) : engine_error(msg) {} - explicit division_by_zero_error() : engine_error(std::string{"atttempt to divide by zero"}) {} + explicit division_by_zero_error() : engine_error(std::string{"division by zero"}) {} }; struct internal_error : public engine_error { From 51f614bb12ac000a5eb20ee5c4385690c9c875b5 Mon Sep 17 00:00:00 2001 From: Doug Torrance Date: Wed, 18 Sep 2024 07:32:56 -0400 Subject: [PATCH 07/14] Rename rawRandomRR -> rawRandomRRUniform For consistency w/ rawRandomRRNormal --- M2/Macaulay2/d/interface.dd | 6 +++--- M2/Macaulay2/e/interface/random.cpp | 6 +++--- M2/Macaulay2/e/interface/random.h | 2 +- M2/Macaulay2/m2/reals.m2 | 4 ++-- 4 files changed, 9 insertions(+), 9 deletions(-) diff --git a/M2/Macaulay2/d/interface.dd b/M2/Macaulay2/d/interface.dd index 0785ae007a..ba77bc9b17 100644 --- a/M2/Macaulay2/d/interface.dd +++ b/M2/Macaulay2/d/interface.dd @@ -70,12 +70,12 @@ export rawRandomQQ(e:Expr):Expr := ( is ht:ZZcell do toExpr(Ccode(QQ, "rawRandomQQ(", ht.v, ")")) else WrongArgZZ()); setupfun("rawRandomQQ",rawRandomQQ); -export rawRandomRR(e:Expr):Expr := ( +export rawRandomRRUniform(e:Expr):Expr := ( when e is prec:ZZcell do if !isULong(prec.v) then WrongArgSmallInteger() - else toExpr(Ccode(RR, "rawRandomRR(", toULong(prec.v), ")")) + else toExpr(Ccode(RR, "rawRandomRRUniform(", toULong(prec.v), ")")) else WrongArgZZ()); -setupfun("rawRandomRR",rawRandomRR); +setupfun("rawRandomRRUniform",rawRandomRRUniform); export rawRandomRRNormal(e:Expr):Expr := ( when e is prec:ZZcell do if !isULong(prec.v) then WrongArgSmallInteger() diff --git a/M2/Macaulay2/e/interface/random.cpp b/M2/Macaulay2/e/interface/random.cpp index 05816170be..6f19b71cfa 100644 --- a/M2/Macaulay2/e/interface/random.cpp +++ b/M2/Macaulay2/e/interface/random.cpp @@ -169,7 +169,7 @@ gmp_QQ rawRandomQQ(gmp_ZZ height) return moveTo_gmpQQ(result); } -gmp_RR rawRandomRR(unsigned long precision) +gmp_RR rawRandomRRUniform(unsigned long precision) /* returns a uniformly distributed random real with the given precision, in * range [0.0,1.0] */ { @@ -193,8 +193,8 @@ gmp_CC rawRandomCC(unsigned long precision) * [1.0,1.0] */ { gmp_CCmutable result = getmemstructtype(gmp_CCmutable); - result->re = const_cast(rawRandomRR(precision)); - result->im = const_cast(rawRandomRR(precision)); + result->re = const_cast(rawRandomRRUniform(precision)); + result->im = const_cast(rawRandomRRUniform(precision)); return reinterpret_cast(result); } diff --git a/M2/Macaulay2/e/interface/random.h b/M2/Macaulay2/e/interface/random.h index 18a923aa00..541b244661 100644 --- a/M2/Macaulay2/e/interface/random.h +++ b/M2/Macaulay2/e/interface/random.h @@ -42,7 +42,7 @@ void rawSetRandomQQ(mpq_ptr result, gmp_ZZ height); /* sets result = a sample from the uniform distribution on [0, height], */ /* rounded to the nearest rational number with denominator bounded by height */ -gmp_RR rawRandomRR(unsigned long prec); +gmp_RR rawRandomRRUniform(unsigned long prec); /* returns a uniformly distributed random real with the given precision, in * range [0.0,1.0] */ diff --git a/M2/Macaulay2/m2/reals.m2 b/M2/Macaulay2/m2/reals.m2 index 67d905d9c8..20fd49f6ed 100644 --- a/M2/Macaulay2/m2/reals.m2 +++ b/M2/Macaulay2/m2/reals.m2 @@ -243,9 +243,9 @@ truncate Number := {} >> o -> x -> ( else if x < 0 then ceiling x else 0) -- e.g., RRi's containing 0 as interior pt -random RR := RR => opts -> x -> x * rawRandomRR precision x +random RR := RR => opts -> x -> x * rawRandomRRUniform precision x random(RR,RR) := opts -> (x,y) -> x + random(y-x) -RR'.random = opts -> R -> rawRandomRR R.precision +RR'.random = opts -> R -> rawRandomRRUniform R.precision CC'.random = opts -> C -> rawRandomCC C.precision random RingFamily := opts -> R -> random(default R,opts) From 0471672b94c3d96be80d04de33c4c4814637e269 Mon Sep 17 00:00:00 2001 From: Doug Torrance Date: Thu, 19 Sep 2024 12:55:35 -0400 Subject: [PATCH 08/14] Memoize GF(ZZ,ZZ) This way, calling GF(p^n) 100 times won't give us 100 different GaloisField objects. This fixes a strange example in the "random(Ring)" docs where we called "tally for i to 100 list random GF 11" and got a Tally object with 100 different key-value pairs. --- M2/Macaulay2/m2/galois.m2 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/M2/Macaulay2/m2/galois.m2 b/M2/Macaulay2/m2/galois.m2 index b92bf30487..87861ea535 100644 --- a/M2/Macaulay2/m2/galois.m2 +++ b/M2/Macaulay2/m2/galois.m2 @@ -117,7 +117,9 @@ findGalois(ZZ,ZZ) := RingElement => opts -> (p,n) -> ( else error (opts.Strategy | " is not a valid argument for' Strategy' option") ) -GF(ZZ,ZZ) := GaloisField => opts -> (p,n) -> ( +GF(ZZ,ZZ) := GaloisField => ( + memo := new CacheTable; + opts -> (p,n) -> memo#(p, n, opts) ??= ( if not isPrime p then error "expected a prime number as base"; if n <= 0 then error "expected positive exponent"; if n == 1 and isMember(opts.Strategy, {"Old", "Aring"}) @@ -126,7 +128,7 @@ GF(ZZ,ZZ) := GaloisField => opts -> (p,n) -> ( primelem := findGalois(p,n,opts); GF(ring primelem, PrimitiveElement=>primelem, Strategy=>opts.Strategy, SizeLimit=>opts.SizeLimit, Variable=>opts.Variable) - ) + )) -- x = if x === null then getSymbol "a" else baseName x; -- R := (ZZ/p) (monoid [x]); From 09f97485f4e0b675e682cc518d375458c829f826 Mon Sep 17 00:00:00 2001 From: Doug Torrance Date: Fri, 20 Sep 2024 17:31:41 -0400 Subject: [PATCH 09/14] Add random(QQ) --- M2/Macaulay2/m2/reals.m2 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/M2/Macaulay2/m2/reals.m2 b/M2/Macaulay2/m2/reals.m2 index 20fd49f6ed..e2cae3468e 100644 --- a/M2/Macaulay2/m2/reals.m2 +++ b/M2/Macaulay2/m2/reals.m2 @@ -249,6 +249,9 @@ RR'.random = opts -> R -> rawRandomRRUniform R.precision CC'.random = opts -> C -> rawRandomCC C.precision random RingFamily := opts -> R -> random(default R,opts) +random QQ := QQ => opts -> x -> rawFareyApproximation( + random numeric x, opts.Height) + -- algebraic operations and functions RR.isBasic = CC.isBasic = RRi.isBasic = true From 97fd04eba859c674f9393101d58664c2b9cc3e20 Mon Sep 17 00:00:00 2001 From: Doug Torrance Date: Sat, 21 Sep 2024 17:58:26 -0400 Subject: [PATCH 10/14] Document random(QQ) --- .../Macaulay2Doc/functions/random-doc.m2 | 23 +++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/M2/Macaulay2/packages/Macaulay2Doc/functions/random-doc.m2 b/M2/Macaulay2/packages/Macaulay2Doc/functions/random-doc.m2 index 6cf8261531..5af2881cdc 100644 --- a/M2/Macaulay2/packages/Macaulay2Doc/functions/random-doc.m2 +++ b/M2/Macaulay2/packages/Macaulay2Doc/functions/random-doc.m2 @@ -61,6 +61,29 @@ Node SeeAlso setRandomSeed +Node + Key + (random, QQ) + Headline + get a random rational number + Usage + random x + Inputs + x:QQ + Height => ZZ + Outputs + :QQ -- randomly chosen from the interval $[0, x]$ + Description + Text + A random number is chosen from the uniform distribution on the interval + $[0, x]$ and then rounded (using the @wikipedia "Farey sequence"@) to the + nearest rational number with denominator bounded by the @CODE "Height"@ + option. + Example + apply(10, i -> random(7_QQ, Height => 5)) + SeeAlso + setRandomSeed + Node Key (random, Type) From a09c8fab6c3b59ab4a5c287fa8f12e8fe84ce69c Mon Sep 17 00:00:00 2001 From: Doug Torrance Date: Sun, 22 Sep 2024 00:44:01 -0400 Subject: [PATCH 11/14] Update random matrix Core test for new random(QQ) behavior --- M2/Macaulay2/tests/normal/randommat.m2 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/M2/Macaulay2/tests/normal/randommat.m2 b/M2/Macaulay2/tests/normal/randommat.m2 index 40ac9db11c..103bf9d1e4 100644 --- a/M2/Macaulay2/tests/normal/randommat.m2 +++ b/M2/Macaulay2/tests/normal/randommat.m2 @@ -17,8 +17,7 @@ assert( f == matrix {{47*x-16*y, 21*x-42*y}, {-17*x-23*y, 16*x-34*y}}) S = QQ[x,y]/(x^2,y^2) g = random(S^{2:0},S^{2:-1}) toString g -assert( g == matrix {{10/3*x+7/4*y, 6/7*x+9/8*y}, {5/6*x+8*y, 3/4*x+8/7*y}}) --- old version: assert( g == matrix {{x+5/6*y, 3/7*x+4/5*y}, {8*x+2/5*y, 3/4*x+1/6*y}} ) +assert( g == matrix {{59/10*x+74/9*y, 39/10*x+34/9*y}, {11/3*x+42/5*y, 35/6*x+1/2*y}}) -- check random isomorphisms are isomorphisms From 25efd13cafe044bb01c66ebb410918528494bb95 Mon Sep 17 00:00:00 2001 From: Doug Torrance Date: Sun, 22 Sep 2024 00:47:47 -0400 Subject: [PATCH 12/14] Update ComputationsBook test for new random(QQ) behavior --- .../tests/ComputationsBook/monomialIdeals/test.out.expected | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/M2/Macaulay2/tests/ComputationsBook/monomialIdeals/test.out.expected b/M2/Macaulay2/tests/ComputationsBook/monomialIdeals/test.out.expected index 4c6015d898..d508628dc3 100644 --- a/M2/Macaulay2/tests/ComputationsBook/monomialIdeals/test.out.expected +++ b/M2/Macaulay2/tests/ComputationsBook/monomialIdeals/test.out.expected @@ -570,9 +570,9 @@ o84 : MonomialIdeal of QQ[z , z , z , z , z , z , i85 : distraction I - 2 2 2 2 2 2 2 3 2 2 3 2 2 2 2 2 3 2 2 3 4 3 2 2 3 4 3 3 2 2 2 2 2 3 2 2 3 4 3 2 2 3 4 -o85 = ideal (x - 281x x - 658x x + 18870x + 102007x x + 81997x , - 162x x x - 355x x x + 116154x x x + 366639x x x x + 245660x x x - 20706840x x - 86836112x x x - 106699274x x x - 34723260x x + 27540x x + 87404x x x + 59285x x - 19746180x x - 81726348x x x - 102990913x x x - 41025220x x + 3520162800x + 18220181320x x + 32640507284x x + 23721732958x x + 5798784420x , - 170x x - 154x x + 172550x x + 330050x x x + 157388x x - 58052620x x - 171376684x x x - 163057048x x x - 50230488x x + 6475361200x + 26001828600x x + 37553372912x x + 22982200824x x + 4970805840x ) - 0 0 3 0 4 3 3 4 4 0 1 3 0 1 4 0 1 3 0 1 3 4 0 1 4 0 3 0 3 4 0 3 4 0 4 1 3 1 3 4 1 4 1 3 1 3 4 1 3 4 1 4 3 3 4 3 4 3 4 4 1 3 1 4 1 3 1 3 4 1 4 1 3 1 3 4 1 3 4 1 4 3 3 4 3 4 3 4 4 + 2 2 2 2 2 2 2 3 2 2 3 2 2 2 2 2 3 2 2 3 4 3 2 2 3 4 3 3 2 2 2 2 2 3 2 2 3 4 3 2 2 3 4 +o85 = ideal (x - 498x x - 338x x + 52785x + 83202x x + 28536x , - 159x x x - 196x x x + 118296x x x + 198930x x x x + 65464x x x - 20148480x x - 46550160x x x - 30683520x x x - 4829440x x + 24327x x + 56064x x x + 32144x x - 18099288x x - 49836834x x x - 42640512x x x - 10736096x x + 3082717440x + 10426525200x x + 12328804800x x + 5771001600x x + 792028160x , - 132x x - 417x x + 113520x x + 458016x x x + 314001x x - 28119168x x - 153120288x x x - 224821032x x x - 68632362x x + 1940336640x + 15229336320x x + 36676727040x x + 26414848080x x + 4305174720x ) + 0 0 3 0 4 3 3 4 4 0 1 3 0 1 4 0 1 3 0 1 3 4 0 1 4 0 3 0 3 4 0 3 4 0 4 1 3 1 3 4 1 4 1 3 1 3 4 1 3 4 1 4 3 3 4 3 4 3 4 4 1 3 1 4 1 3 1 3 4 1 4 1 3 1 3 4 1 3 4 1 4 3 3 4 3 4 3 4 4 o85 : Ideal of S From 82aabfe76f4e30a9843317b149d29798d2291233 Mon Sep 17 00:00:00 2001 From: Doug Torrance Date: Thu, 19 Sep 2024 13:10:08 -0400 Subject: [PATCH 13/14] Use a different random seed for randomBinaryForm example Otherwise, we end up calling QEPCAD, which will fail if it's not available (e.g., on the macOS GitHub builds). --- M2/Macaulay2/packages/CoincidentRootLoci/documentation.m2 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/M2/Macaulay2/packages/CoincidentRootLoci/documentation.m2 b/M2/Macaulay2/packages/CoincidentRootLoci/documentation.m2 index f7e7f45d92..08c6dbfbef 100644 --- a/M2/Macaulay2/packages/CoincidentRootLoci/documentation.m2 +++ b/M2/Macaulay2/packages/CoincidentRootLoci/documentation.m2 @@ -330,7 +330,7 @@ randomBinaryForm(d,r,) randomBinaryForm(d,,c)", Inputs => {"d" => ZZ,"r" => ZZ,"c" => ZZ}, Outputs => {RingElement => {"a random binary form of degree ",TT"d",", real rank ",TT "r"," and complex rank ",TT "c"}}, -EXAMPLE {"F = randomBinaryForm 5","F = randomBinaryForm(5,4,3)","(realrank F,complexrank F)","F = randomBinaryForm(6,4,4)","(realrank F,complexrank F)",}, +EXAMPLE {"setRandomSeed 2", "F = randomBinaryForm 5","F = randomBinaryForm(5,4,3)","(realrank F,complexrank F)","F = randomBinaryForm(6,4,4)","(realrank F,complexrank F)",}, SeeAlso => {realrank,complexrank}} document { Key => {realroots,(realroots,RingElement)}, From 7fabcd33146d69c8f1e80ab6cd433ff62c9eb91e Mon Sep 17 00:00:00 2001 From: Doug Torrance Date: Sun, 22 Sep 2024 09:04:57 -0400 Subject: [PATCH 14/14] Raise an error in random(QQ) if height is nonpositive --- M2/Macaulay2/d/interface.dd | 3 +-- M2/Macaulay2/e/interface/random.cpp | 15 ++++++++++++++- M2/Macaulay2/e/matrix.cpp | 14 ++++++++++++-- 3 files changed, 27 insertions(+), 5 deletions(-) diff --git a/M2/Macaulay2/d/interface.dd b/M2/Macaulay2/d/interface.dd index ba77bc9b17..7834562ffa 100644 --- a/M2/Macaulay2/d/interface.dd +++ b/M2/Macaulay2/d/interface.dd @@ -66,8 +66,7 @@ export rawFareyApproximation(e:Expr):Expr := ( setupfun("rawFareyApproximation", rawFareyApproximation); export rawRandomQQ(e:Expr):Expr := ( when e - is Nothing do toExpr(Ccode(QQ, "rawRandomQQ(", "0)")) - is ht:ZZcell do toExpr(Ccode(QQ, "rawRandomQQ(", ht.v, ")")) + is ht:ZZcell do toExpr(Ccode(QQorNull, "rawRandomQQ(", ht.v, ")")) else WrongArgZZ()); setupfun("rawRandomQQ",rawRandomQQ); export rawRandomRRUniform(e:Expr):Expr := ( diff --git a/M2/Macaulay2/e/interface/random.cpp b/M2/Macaulay2/e/interface/random.cpp index 6f19b71cfa..f361152948 100644 --- a/M2/Macaulay2/e/interface/random.cpp +++ b/M2/Macaulay2/e/interface/random.cpp @@ -3,6 +3,9 @@ #include "interface/random.h" #include "interface/gmp-util.h" +#include "exceptions.hpp" +#include "error.h" + #define INITIALMAXINT 10 #define IA 16807 @@ -152,6 +155,9 @@ void rawSetRandomQQ(mpq_ptr result, gmp_ZZ height) mpfr_t x; if (height == nullptr) height = maxHeight; + if (mpz_cmp_si(height, 0) <= 0) + throw exc::engine_error("expected a positive height"); + mpfr_init2(x, gmp_defaultPrecision); mpfr_urandomb(x, state); mpfr_mul_z(x, x, height, MPFR_RNDN); @@ -165,7 +171,14 @@ gmp_QQ rawRandomQQ(gmp_ZZ height) { mpq_ptr result = getmemstructtype(mpq_ptr); mpq_init(result); - rawSetRandomQQ(result, height); + + try { + rawSetRandomQQ(result, height); + } catch (const exc::engine_error& e) { + ERROR(e.what()); + return nullptr; + } + return moveTo_gmpQQ(result); } diff --git a/M2/Macaulay2/e/matrix.cpp b/M2/Macaulay2/e/matrix.cpp index 3e41f222c6..2832d6ec33 100644 --- a/M2/Macaulay2/e/matrix.cpp +++ b/M2/Macaulay2/e/matrix.cpp @@ -745,7 +745,12 @@ Matrix *Matrix::random( int32_t u = rawRandomInt((int32_t)10000); if (u > threshold) continue; } - mat.set_entry(j, i, R->random()); + try { + mat.set_entry(j, i, R->random()); + } catch (const exc::engine_error& e) { + ERROR(e.what()); + return nullptr; + } } } else if (special_type == 1) @@ -760,7 +765,12 @@ Matrix *Matrix::random( int32_t u = rawRandomInt((int32_t)10000); if (u > threshold) continue; } - mat.set_entry(j, i, R->random()); + try { + mat.set_entry(j, i, R->random()); + } catch (const exc::engine_error& e) { + ERROR(e.what()); + return nullptr; + } } } }