Update random(QQ) to use standard uniform distribution

M2/Macaulay2/d/interface.dd

@@ -51,18 +51,30 @@ 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)"))
-     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 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()

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
 {

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
@@ -67,34 +70,119 @@ gmp_ZZ rawRandomInteger(gmp_ZZ maxN)
   return result;
 }
 
-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 */
+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);
-  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);
+  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 */
+/* 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);
+  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);
+  rawSetFareyApproximation(result, x, height);
+  mpfr_clear(x);
+}
+
+gmp_QQ rawRandomQQ(gmp_ZZ 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);
+
+  try {
+    rawSetRandomQQ(result, height);
+  } catch (const exc::engine_error& e) {
+    ERROR(e.what());
+    return nullptr;
+  }
+
+  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] */
 {
@@ -118,8 +206,8 @@ gmp_CC rawRandomCC(unsigned long precision)
  * [1.0,1.0] */
 {
   gmp_CCmutable result = getmemstructtype(gmp_CCmutable);
-  result->re = const_cast<gmp_RRmutable>(rawRandomRR(precision));
-  result->im = const_cast<gmp_RRmutable>(rawRandomRR(precision));
+  result->re = const_cast<gmp_RRmutable>(rawRandomRRUniform(precision));
+  result->im = const_cast<gmp_RRmutable>(rawRandomRRUniform(precision));
   return reinterpret_cast<gmp_CC>(result);
 }
 

M2/Macaulay2/e/interface/random.h

@@ -28,15 +28,21 @@ 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 */
+/* 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);
+gmp_RR rawRandomRRUniform(unsigned long prec);
 /* returns a uniformly distributed random real with the given precision, in
  * range [0.0,1.0] */
 

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;
+              }
             }
         }
     }

M2/Macaulay2/e/unit-tests/testMain.cpp

@@ -2,6 +2,8 @@
 #include <M2/gc-include.h>
 #include <engine.h>
 
+unsigned long gmp_defaultPrecision = 53;
+
 int main(int argc, char **argv)
 {
   IM2_initialize();

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]);

M2/Macaulay2/m2/reals.m2

@@ -243,12 +243,15 @@ 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)
 
+random QQ := QQ => opts -> x -> rawFareyApproximation(
+    random numeric x, opts.Height)
+
 -- algebraic operations and functions
 
 RR.isBasic = CC.isBasic = RRi.isBasic = true

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)}, 

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)
@@ -77,8 +100,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

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
 

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"

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