ivy: add an operator, rand, for floating-point random numbers

APL (or maybe just Dyalog) uses ?0 for a float in the range [0,1),
which we could do but it's clunky.

We also can't just do ?1.0 because numbers are reduced before hitting
the operator: 0.5 becomes rational 1/2 and 1.0 becomes integer one.

Instead, we add a new operator, rand, that behaves exactly like the
roll operator ? but ranges over all floats in [0, n) instead of
just integers.

To get high-precision random floats, we create two integers with
the size of the mantissa and divide them, rand()/2^floatPrec, to
get a high-precision mantissa, then scale it appropriately. I doubt
it's mathematically perfect but it seems pretty good; I included a
test to show that it gives mean and standard deviations as expected.

Fixes #167

doc.go

@@ -69,6 +69,7 @@ Unary operators
 
 	Name              APL   Ivy     Meaning
 	Roll              ?B    ?       One integer selected randomly from the first B integers
+	Random            ?0    rand    Like ?, but floating point. (APL uses ?0 as rand in [0,1)).
 	Ceiling           ⌈B    ceil    Least integer greater than or equal to B
 	                                If B is complex, the complex ceiling, as defined by McDonnell
 	Floor             ⌊B    floor   Greatest integer less than or equal to B

testdata/unary_bigfloat.ivy

@@ -201,6 +201,24 @@ box sqrt 2
 first sqrt 2
 	1.41421356237
 
+# Random numbers
+)seed 0
+rand 3
+	1.96599645204
+
+# High precision, not float64s.
+)seed 0
+1e1000 * rand 1e-1000
+	0.655332150681
+
+# Test stats. For a uniform distribution, ideally: mean = .5,  𝛔 = 1/√12 = 0.2887
+)seed 0
+x = rand@ 0 != iota 10000
+mean = (+/ x) / 10000
+sigma = sqrt 1/10000 * +/(x - mean) ** 2
+mean sigma
+	0.500853881349 0.288396812059
+
 # There is a (possible?) but in big.Float.Copy that made this not work.
 # BigFloat.inverse would change the value of floatOne. This is a simple
 # check against that; when broken, the output is 0.5 0.5 (!).

value/unary.go

@@ -7,6 +7,7 @@ package value
 import (
 	"fmt"
 	"math/big"
+	"math/bits"
 	"unicode/utf8"
 )
 
@@ -127,6 +128,34 @@ func printValue(c Context, v Value) Value {
 	return v
 }
 
+// bigRand returns a uniformly distributed BigFloat in the range [0, f).
+// For [0, 1), the mean should be 0.5 and 𝛔 should be 1/√12, or 0.2887.
+// A test of a million values yielded 0.500161824174 0.288777488704.
+func bigRand(c Context, f *big.Float) Value {
+	if f.Sign() < 0 {
+		Errorf("rand of negative number")
+	}
+	// We want two integers FloatPrec bits long.
+	prec := c.Config().FloatPrec()
+	// Word is a uint. We build max by setting all FloatPrec bits.
+	// (If FloatPrec is not a multiple of UintSize, it's still OK, we'll
+	// just make integers longer than the precision will hold, but not much.)
+	maxBits := make([]big.Word, (prec+(bits.UintSize-1))/bits.UintSize)
+	for i := range len(maxBits) {
+		maxBits[i] = big.Word(^uint(0))
+	}
+	max := big.NewInt(0).SetBits(maxBits)
+	max.Add(max, bigIntOne.Int) // big.Int.Rand is [0,n)
+	// Now pick a random integer from [0 to max)
+	rand := big.NewInt(0).Rand(c.Config().Random(), max)
+	// Make it a float and normalize it.
+	x := big.NewFloat(0).SetInt(rand)
+	x.Quo(x, big.NewFloat(0).SetInt(max))
+	// Finally, scale it up to [0,v).
+	x.Mul(x, f)
+	return BigFloat{x}
+}
+
 func init() {
 	ops := []*unaryOp{
 		{
@@ -149,6 +178,31 @@ func init() {
 			},
 		},
 
+		{
+			name:        "rand",
+			elementwise: true,
+			fn: [numType]unaryFn{
+				intType: func(c Context, v Value) Value {
+					var x big.Float
+					x.SetInt64(int64(v.(Int)))
+					return bigRand(c, &x)
+				},
+				bigIntType: func(c Context, v Value) Value {
+					var x big.Float
+					x.SetInt(v.(BigInt).Int)
+					return bigRand(c, &x)
+				},
+				bigRatType: func(c Context, v Value) Value {
+					var x big.Float
+					x.SetRat(v.(BigRat).Rat)
+					return bigRand(c, &x)
+				},
+				bigFloatType: func(c Context, v Value) Value {
+					return bigRand(c, v.(BigFloat).Float)
+				},
+			},
+		},
+
 		{
 			name:        "j",
 			elementwise: true,