Two Ulongs to Double

[deltarho[1859]]

I suspect that some of you will have seen this but I could not find it.

In my PRNG research I found this at Sebastiano Vigna's website here, towards the bottom.

Code: Select all

#include <stdint.h>
static inline double to_double(uint64_t x) {
  const union { uint64_t i; double d; } u = { .i = UINT64_C(0x3FF) << 52 | x >> 12 };
  return u.d - 1.0;
}

The above generates uniform doubles in the unit interval.

I am using it to convert 2 x 32-bit unsigned integers from a 32-bit output PRNG to a 64-bit Double.

To test it I total 10^7 random Ulong pairs and then determine the average, aiming at 0.5. On one run I got 0.5000907638144706 in 0.36 seconds, so it is fairly swift.

Code: Select all

Const DoubleBits = 1023ull Shl 52
Union UllToDbl
  i As Ulongint
  d As Double
End Union
 
' Adaptation of mkulongint by, I think, dodicat?
' In practice n1 And n2 will be from a PRNG Ulong Output
' Vigna's code within the following

Function TwoUlongsToDouble( n1 as Ulong, n2 as Ulong ) as Double
Dim Ans As UllToDbl
Dim As Ulongint res
  Cast(Ulong Ptr,@res)[0]=n1
  Cast(Ulong Ptr,@res)[1]=n2
  Ans.i = DoubleBits Or (res Shr 12)
  Return Ans.d - 1.0
End Function

Randomize 

Dim As Ulong i
Dim As Double t, tot

' CULng(Rnd*&hFFFFFFFFul) by MrSwiss
t = Timer
For i = 1 To 10^7
  tot += TwoUlongsToDouble( CULng(Rnd*&hFFFFFFFFul), CULng(Rnd*&hFFFFFFFFul) )
Next
t = Timer - t

Print tot/10^7, t

Sleep

[MrSwiss]

@deltarho[],

while you might be a talented mathematician, your code contains:
too many translation (and other) mistakes ... (36 Seconds -- NOT 0.36)

Corrected version (34 Seconds FBC 64 straight):

Code: Select all

'Const DoubleBits = 1023ull Shl 52   ' translation error = 2047 (dec)
Const As ULongInt DoubleBits = &h03FF Shl 52    ' make certain: NOT UInteger!!!

Union UllToDbl
    As Ulongint i
    As Double   d
End Union
 
' recoded: by MrSwiss

Function TwoUlongsToDouble( n1 as Ulong, n2 as Ulong ) as Double
    Dim Ans As UllToDbl
    Dim As Ulongint res = (n1 Shl 32) + n2    ' no ptr cast's needed (slow)
    Ans.i = DoubleBits Or (res Shr 12)
    Return Ans.d - 1.0
End Function

Randomize , 5

Dim As Ulong i
Dim As Double t, tot

' CULng(Rnd*&hFFFFFFFFul) by MrSwiss
t = Timer
For i = 1 To 10^7
    tot += TwoUlongsToDouble( CULng(Rnd * &hFFFFFFFFul), CULng(Rnd * &hFFFFFFFFul) )
Next
t = Timer - t

Print tot/10^7, t

Sleep

[deltarho[1859]]

@MrSwiss

Every correction you made is wrong.

Warning 33(0): Shift value greater than or equal to number of bits in data type. WinFBE will not continue on warnings.

&h03FF Shl 52

That is why I used 1023ull Shl 52

&h03FF => 1023 and not 2047

n1 As Ulong and you are doing a 'n1 Shl 32'. There is no room to Shl 32. Got another Warning 33(0)

Randomize , 5? Any Rnd will do for illustration purposes.

I am still getting 0.36 seconds.

The Cast(Ulong Ptr,@res)[0]=n1 and so on is Dodicat's code - I think.

I must admit that I do not feel comfortable questioning your code - I regard you as one of the FreeBASIC masters.

BTW, I am not a talented mathematician. I burnt out years ago. I am now just an old fart.

[MrSwiss]

@deltarho[],

you are obviously using: FBC 32 (I've clearly stated: FBC 64).
There are absolutely *no warnings/errors* ...
(FBC 32 is known, to have such issues, because of: Integer = Long and,
all those "behind the scenes" use of Integer)

Yes, the HEX to DEC conversion of mine was wrong, but in order to avoid
such issues, I'm simply using HEX notation ...

[deltarho[1859]]

@MrSwiss

I've clearly stated: FBC 64

Yes, you did - I rarely use FBC 64.

The difference in speed between

Code: Select all

Dim As Ulongint res = (n1 Shl 32) + n2    ' no ptr cast's needed (slow)

and

Code: Select all

Dim As Ulongint res
  Cast(Ulong Ptr,@res)[0]=n1
  Cast(Ulong Ptr,@res)[1]=n2

is negligible.

Perhaps the compiler is doing clever tricks.

Your code comes in at about 109 seconds, with either, giving 0.4999147802485142.

However, there is something wrong somewhere because with my code, FBC 32, I am still getting 0.36 seconds.

Yes, the HEX to DEC conversion of mine was wrong, but in order to avoid
such issues, I'm simply using HEX notation ...

I am from the north east of England - we still use decimal up there. <smile>

[MrSwiss]

Perhaps the compiler is doing clever tricks.

No, don't think so.

The difference is: Integer = LongInt (FBC 64), which avoids the: *shift value ...* issues.
In FBC 32, another TypeCast such as: (CULngInt(n1) shl 32) + n2, would do the trick, too.
(but, every one of those, is costing speed)

However, there is something wrong somewhere because with my code, FBC 32, I am still getting 0.36 seconds.

Code adapted for FBC 32, clock's at 48 Seconds ... seems to be a problem, on your side.
(unless you have the latest: i9++ from Intel, currently "under test" <grin>)

[deltarho[1859]]

Code adapted for FBC 32, clock's at 48 Seconds ... seems to be a problem, on your side.

48 Seconds!

Running from my IDE there is a slight hesitation before the console pops up, anti-virus checking, and then before I can blink, I get the result. Running the exe from my file manager the console pops up almost instantaneously and a split second later I get the result. I don't need Timer to tell me the result is coming in easily within a second. I am getting 0.42 seconds with gas.

unless you have the latest: i9++ from Intel, currently "under test" <grin>)

Intel i7-3770K @3.5GHz, 3.9GHz with turbo. Turbo always kicks in and stays in - high performance liquid cooler on board.

Will someone else do a speed check?

[deltarho[1859]]

Revised code. One of MrSwiss' suggestions taken on board ie Dim As Ulongint res - the original was convoluted, to say the least. Still 0.36 seconds.

FBC 32, -gen gcc -Wc -O3

Code: Select all

Const DoubleBits = 1023ull Shl 52

Union UllToDbl
  i As Ulongint
  d As Double
End Union
 
' In practice, n1 And n2 will be from a PRNG Ulong Output
' Vigna's code within the following

Function TwoUlongsToDouble( n1 As Ulong, n2 As Ulong ) As Double
Dim Ans As UllToDbl
Dim As Ulongint res = (CULngint(n1) Shl 32) + n2 ' MrSwiss
  Ans.i = DoubleBits Or (res Shr 12)
  Return Ans.d - 1.0
End Function

Randomize 

Dim As UInteger i
Dim As Double t, tot

' Culng(Rnd*&hFFFFFFFFul) by MrSwiss
t = Timer
For i = 1 To 10^7
  tot += TwoUlongsToDouble( Culng(Rnd*&hFFFFFFFFul), Culng(Rnd*&hFFFFFFFFul) )
Next
t = Timer - t

Print tot/10^7, t

Sleep

[MrSwiss]

Will someone else do a speed check?

I've figured it out:
Randomize , 5 <-- is the culprit (is using up, all the time ...)

Another improvement: use UInteger for "i" (as loop-iterator), aka:
ULong (FBC 32) or ULongInt (FBC 64)

[deltarho[1859]]

I've figured it out:

I should have spotted that being a random number nutter.

Randomize , 5 is great for a few numbers but not going past it 2x10^7 times it isn't.

Good catch!

Added: " use UInteger for "i" ". Done, thanks

[deltarho[1859]]

In the following, I am not trying to tell anyone how to suck eggs but some may benefit. In my second post above I wrote

n1 As Ulong and you are doing a 'n1 Shl 32'. There is no room to Shl 32. Got another Warning 33(0)

Somebody, including MrSwiss, should have jumped on that. Out of context, it is correct but in context it is wrong.

In assembly, we have shl working on registers. In PowerBASIC we have shl working on integral variables as SHIFT [SIGNED] {LEFT|RIGHT} ivar, countexpr. In both cases, the object being shifted is updated, for a non-zero shift. There is no 'in context' here.

From MrSwiss' first post we have this:

Code: Select all

Dim As Ulongint res = (n1 Shl 32) + n2

Assume we are using FBC 64. res is being dimensioned as 64 bit. n1 Shl 32 will see n1 being placed into res and then shifted to the left by 32 bits. res starts as empty but now has n1 at the 'big end'. The 'little end' is empty and we now place n2 there. FBC 64 does not say a word. FBC 32, on the other hand, issues a warning. It should not because n1 is not being altered. It's logic is the same as mine above.

With FBC 32 we could use (CULngInt(n1) shl 32) + n2 to avoid a compilation warning.

Suppose we have, again assuming we are using FBC 64

Code: Select all

Dim As Ulong res = (n1 Shl 32) + n2

Working through this n1 gets placed into res. However, because res is now 32 bit the shl sees n1 disappear into the ether and res is empty again. n2 now gets placed into res.

So

Code: Select all

Dim As ulong n1, n2
n1 = 99
n2 = 45
Dim As Ulong res = (n1 Shl 32) + n2
Print res, n1
Sleep

will see 45 99 being printed as if n1 never existed. Notice that n1 has not changed.

With res as Ulongint we will see 425201762349 99 being printed

If your background is asm plus PowerBASIC then you may fall into the same trap as me. If your background is FreeBASIC with some C, perhaps, then you may be wondering what all the fuss is about.

[MrSwiss]

Assume we are using FBC 64. res is being dimensioned as 64 bit. n1 Shl 32 will see n1 being placed into res and then shifted to the left by 32 bits. res starts as empty but now has n1 at the 'big end'. The 'little end' is empty and we now place n2 there. FBC 64 does not say a word. FBC 32, on the other hand, issues a warning. It should not because n1 is not being altered. It's logic is the same as mine above.

No, don't agree with this reasoning, because:

• 1. operation: solve out the bracket = n1 shl 32 (load n1)
  2. operation: shift left 32 (in a 64 bit register)
  3. operation: load n2
  4. operation: add them up
  5. operation: assign (result) to res

In Assembly (simplified):

Code: Select all

ASM ' 64 bit assembly
    xor rax, rax
    xor rdx, rdx
    mov eax, dword ptr [n1]
    shl rax, 32
    mov edx, dword ptr [n2] 
    add rax, rdx
    mov qword ptr [res], rax
END ASM

The same operation in 32 bit ASM, is quite different ...
(just correct assignment needs to be done, to 64 bit variable, but limited to: [shl/shr 32])

[deltarho[1859]]

The difference between what I wrote and what MrSwiss wrote is similar to the difference between a thought experiment and a mathematical paper. Albert Einstein and Richard Feynman did not get Nobel prizes for their communication skills and thought experiments but for their mathematical papers. Stephen Hawking was told that the expected sales for his 'A Brief History of Time' would be halved for each formula he uses in his book. He heeded that advice and his book sold more than 10 million copies in 20 years and was translated into 35 languages by 2001.

[coderJeff]

' CULng(Rnd*&hFFFFFFFFul) by MrSwiss

I started answering this in the useful defines topic, then saw it here too.

I've come across this before in OpenGL, mapping floating point values to integer colours and pixel locations (and the reverse mapping).

Taking N = cubyte(rnd()*255), as example:
- 255 is the largest value the ubyte can hold; 1 less than the number of values it could contain.
- A ubyte has 256 possible values.

Assuming rnd() gives even distribution over interval [0, 1):

With rounding:
- might get 0 or 255, but only half as likely as any other number.
- 0 can only round down from 0.5.
- 255 can only round up from 254.5
- 1 can round up from 0.5 or down from 1.5.
- ditto for 2 through 254
- i.e. should expect to see 1 though 254 occurring twice as often as 0 or 255.

With truncating:
- will never get 255.

So, want to map the range [0, 1) to [0, 256) evenly:
R = rnd() * 256

Then, truncating the floating point gives an even distribution of integers
I = cubyte( fix( rnd() * 256 ) )

It's similar issue for rnd()*&hfffffffful, where the multiplier is one less the number of values possible.

If that's important or not, entirely up to the programmer. A simple cubyte(rnd()*255), might be good enough for the purpose.

[coderJeff]

In my PRNG research ...

Neat.

You have 2 x 32 bits to stuff in to 52 bits. Does it matter which bits are thrown away? This one throws away the high end bits, so no shifts in the 32 bit version.

Code: Select all

union UllToDbl
	type
		lo32 as ulong
		hi32 as ulong		
	end type
	d as double
end union

function TwoUlongsToDouble32( n1 as ulong, n2 as ulong ) as double
	dim ans as UllToDbl = any
	ans.lo32 = n1 
	ans.hi32 = n2 and &hffffful or &h3ff00000ul
	return ans.d - 1.0
end function

function TwoUlongsToDouble64( n1 as ulong, n2 as ulong ) as double
	dim ans as ulongint = (culngint(n2) shl 32 or n1) and &hfffffffffffffull or &h3ff0000000000000ull 
	return *cast(double ptr, @ans) - 1.0
end function

'' 32-bit
print TwoUlongsToDouble32(0, 0) 
print TwoUlongsToDouble32(0, 1)
print TwoUlongsToDouble32(1, 0)
print TwoUlongsToDouble32(&hfffffffful, &hfffffffful)

'' 64-bit
print TwoUlongsToDouble64(0, 0) 
print TwoUlongsToDouble64(0, 1)
print TwoUlongsToDouble64(1, 0)
print TwoUlongsToDouble64(&hfffffffful, &hfffffffful)

A couple variations, should work in 32/64 bit modes.