World: r3wp

suggestion: R3

The reason, I prefer your first version is:
- I always know the endpoints exactly.
- The mean is well defined.
- The function is simpler and faster.

REBOL [
    Author: "Ladislav Mecir"
    Purpose: {Converts decimal to ordinal and back}
]

make object! [
    int-min: to integer! #{8000000000000000}
    set 'dec-as-int func [x [decimal!]] [
        x: to integer! to binary! x
        either x < 0 [int-min - x] [x]
    ]
    set 'int-as-dec func [x [integer!]] [
        if x < 0 [x: int-min - x]
        to decimal! to binary! x
    ]
]

your "difference in ulps" is just the difference between the ordinal 
numbers

So an R3 ulps could be:

ulps: func [
	lo	[decimal!]
	hi	[decimal!]
][
	(to-integer to-binary hi) - to-integer to-binary lo
]

>> ulps 0.0 0.1
== 4591870180066957722
>> ulps 0.1 0.2
== 4503599627370496

(as-int hi) - (as-int lo)

It seems, numbers lie much closer around zero than around 0.1.

sorry, (dec-as-int hi) - (dec-as-int lo)

>> (dec-as-int 0.1) - (dec-as-int 0.0)
== 4591870180066957722
>> (dec-as-int 0.2) - (dec-as-int 0.1)
== 4503599627370496

Same result. Now I'm confused! :-)

If this is correct, then my suggestion to do
random 1.0
many times will actually show the problem.

yes, you get different result only when trying (dec-as-int 0.0) - 
(dec-as-int -0.0)

etc.

The highest decimal is
>> to-decimal #{7fef ffff ffff ffff}
== 1.79769313486232e+308


There is the same number of ulps between 0.0 and 2.0 as between 1.0 
and the highest decimal:

>> ulps 0.0 2.0
== 4611686018427387904
>> ulps 1.0 to-decimal #{7fef ffff ffff ffff}
== 4611686018427387903


I think, the space between number are different all the way from 
zero to the highest decimal, but I'm not sure.

space between *numbers* are ...

Funnily enough, 1.5 is the middle value when talking ulps:

>> ulps 0.0 1.5
== 4609434218613702656
>> ulps 1.5 hi
== 4609434218613702655

where hi is:
hi: to-decimal #{7fef ffff ffff ffff}

This mean, there is as many different possible decimals between 0.0 
and 1.5 as between 1.5 and the highest possible 64-bit decimal, when 
talking floating points in a computer using the 64-bit IEEE 754 standard.

how about this, do you think, that there is a problem with it too?
two**62: to integer! 2 ** 62
two**9: to integer! 2 ** 9
two**53: to integer! 2 ** 53
r-uni: func [x [decimal!] /local d] [
	d: random two**62
	; transform to 0..2**53-1
	d: d // two**53
	; transform to 0.0 .. x (0.0 inclusive, x exclusive)
	d / two**53 * x
]

a modification:

two**62: to integer! 2 ** 62
two**53: to integer! 2 ** 53
two**53+1: two**53 + 1
r-uni: func [x [decimal!] /local d] [
	d: random two**62
	; transform to 0..2**53
	d: d // two**53+1
	; transform to 0.0 .. x (0.0 inclusive, x inclusive)
	d / two**53 * x
]

When you two are done, could you look at CureCode ticket #1027 and 
R3's alpha 66 behavior and see if it is correct? I can't judge.

LOL, that is what we are just discussing

anyway, this is always a problem, when we try to generate numbers 
in an "exotic" interval, not in the usual [0.0 .. 1.0]

I'll just grab some food, then I look at your r-uni.

hmm, seems, that the former may still produce X when rounding?

In any case, this chat should be the basis for the docs that go with 
it. People aren't aware of how much thought goes into these things.

>> two**62 / two**53+1
== 512.0

So two**53+1 * 512 should equal two**62, but it doesn't:

>> two**62
== 4611686018427387904
>> two**53+1 * 512
== 4611686018427388416

Why is that?

(I wanted to check , if d // two**53+1 gave an even distribution, 
but there seem to be something wrong with calculations of large numbers.)

Why is that?
 - rounding

This should give an math overflow, I think:

>> to integer! 2 ** 62
== 4611686018427387904

Like this does in R2:

>> to integer! 2 ** 31
** Math Error: Math or number overflow

no, R3 uses 64-bit integers

no wait, it's 2 ** 63, that should give overflow.

yeah

But!? :-)

>> to integer! 2 ** 62 - 1
== 4611686018427387904
>> to integer! 2 ** 62
== 4611686018427387904

Isn't that a bug?

no, it is the precision limit of IEEE754

I see.

(I mean, it exceeds the precision limit)

Yes, I understand.

After you do:

d: d // two**53+1


Are we sure, the new d is even distributed? That every value, d can 
now take, is equal possible?

hmm, to make sure, I should use rejection

as follows:

rnd: func [
	{random range with rejection, not using last bits}
	value [integer!]
	/local r s
] [
	s: two-to-62 - (two-to-62 // value)
	until [
		r: random two-to-62
		r <= s
	]
	s: s / value
	r - (r // s) / s + 1
]

If all possible values of d is equal possible, and d now can have 
values from 0 to two**53, then

d / two**53

will return values from 0.0 to 1.0 both inclusive, and 

d / two**53 * x

will return values from 0.0 to x both inclusive. Right?

yes

...multiplication may still trigger rounding, though

It give same result as your first RANDOM:

d: to-decimal to-binary 1023
blk: []
insert/dup blk 0 1024
random/seed now

loop 1000000 [i: to-integer to-binary r-uni d blk/(i + 1): blk/(i 
+ 1) + 1]

>> print [blk/1 blk/512 blk/1024]
462 1036 481

that is caused by the multiplication rounding

yes

I like it better than the version, where max is changed to 0.0

another variant would be to generate the uniform deviates in the 
definite interval not allowing multiplication to "mess things"

Is r-uni better than the version, you gave earlier?

tt62: to integer! 2 ** 62
r: func [x [decimal!]] [(random tt62) - 1 / tt62 * x]