Two Ulongs to Double

[coderJeff]

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.

heh, may be. :)

You got a lot of interesting things going on in this topic...I will read again that post. Did PB have same behaviour 32bit/64bit?

[MrSwiss]

Did PB have same behaviour 32bit/64bit?

Afaik: PB didn't ever have a x64 version (16/32 bit only).

At the time I've used PB, it was 16 bits (on DOS), x86/x286 ... just after
TurboBASIC was sold back, to its developer Robert Zale (by Borland).
PB ver. 2.1x ...

[deltarho[1859]]

i.e. should expect to see 1 through 254 occurring twice as often as 0 or 255.

This mapping bias has been exercising folk for some years.

The problem with mapping using [0,1) is that [0,1) is a function of mapping itself from [0, 2^32-1] to [0,1) assuming that a generator is outputting Ulongs. A better way is to map [0, 2^32-1] to our desired range and by-pass [0,1).

However, to do that we need the generator's Ulong output and none of FB's generators gives us that. All of yours truly generators, from RNDMT, PCG32II, CryptoRndII and CMWC, do just that. More to the point they all have a range function. We supply the lower and upper bounds of the range and the range function requests a Ulong from the generator engine. After a bit of asm we then get a Long returned.

As a test, I looped through 102400000 CMWC, my latest offering, range requests incrementing each count(0 To 255) depending upon the returned value. The expected value for each 'cubby hole' is then 400,000. Well, we are talking random numbers here so we won't get that. I sorted the results and got count(0) = 398424 and count(255) = 401731; the first being 1576 below expectation and the last being 1731 above expectation. A sorted count(0)/count(255) will tend to 1 as the number of requests tends to infinity.

So, no bias.

Added: With Dim pcg As pcg32, pcg.range(0,255) clocks over 200MHz.

[coderJeff]

deltarho, nice summary of the problem, very well put. Yes, fb's generators are using ulong with the last step to divide the ulong by 2^32 to return a double. Can get back the exact ulong value with rnd()*cdbl(4294967296ull); not so appealing if you want the fastest.

Re SHL 32, etc. thank you for taking time to explain. wiki page Coercion and Conversion is supposed to help, but it is missing a few things. I see. SHIFT is a statement in PB. And SHL is binary operator in FB.

A few more details, to add to your info and MrSwiss info:

1) variable = expression, the expression is fully evaluated, then assigned (stored) to the variable.

2) fbc automatically coerces types, most commonly from a type that is smaller than the machine register size up to the machine register size (32/64bit depending on target), otherwise programmer would be required to explicitly convert everything to compatible types. Coercion is different depending on target 32/64bit.

3) expr SHL N
- expr gets loaded to a machine register 32bit or 64bit, depends on fbc target
- SHL 32, shifts the machine register, warning on 32bit. On 64bit, it's fine
- on 64bit, we would get a warning on SHL 64, or higher
- at least SHL Page warns us that the results are undefined if N is negative or more than bits in the register.
- FYI, the undefined behavior under the hood is that, fbc translates to expr SHL (N MOD BITS), where bits is 32 or 64, depends on fbc target.

4) Can use C### conversions where ever needed because fbc optimizes them away, if there is no change in data type size.

5) This code is the minimum that makes sense for both 32/64bit with result in 64bit variable:

Code: Select all

dim as ulong n1, n2
Dim As ulongint res = (culongint(n1) Shl 32) + n2

FYI

My speed timings, probably a little slow because I am running an un-optimized debug version of fb's rtlib:

Code: Select all

fbc -gen gas -target win32
TwoUlongsToDouble:  0.4999107723015807     1.476572121705772
TwoUlongsToDouble32:  0.5000478563958972   1.422788166046075
TwoUlongsToDouble64:  0.4999743510336273   1.452752335902858

fbc -gen gcc -target win64
TwoUlongsToDouble:  0.4999608822256317     2.615885008126497
TwoUlongsToDouble32:  0.4997475645844072   2.571060933754779
TwoUlongsToDouble64:  0.4999529535941112   2.547443617833778

[deltarho[1859]]

@coderJeff

Had T...32 and T...64 seriously outperformed T then I would adopt them but since they don't and would require a conditional compile then I'll stay with T.

It is worth noting that I have an asm/SSE version of T and your T...32 and T...64 are pretty much matching it.

It is also worth noting that passing a function millions of times in PowerBASIC would definitely see a performance boost when converting to asm. With FreeBASIC and the compiler optimized to the hilt then that may not be the case.

[sancho3]

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.

Ran it on Ubuntu 64:

Code: Select all

 0.4999994417868154          0.4158210754394531

[deltarho[1859]]

A bit more.

Here is a function which will map either two Ulongs, (0,2^32-1), as above, or a single Ulong into (0,1]. Using two Ulongs we get 53-bit granularity and with one Ulong we get 32-bit granularity.

We can also use the function to map a Ulongint, (0,2^64-1), into (0,1] if we break the Ulongint into two Ulongs first. The following uses GetUlongs by dodicat.

Code: Select all

Function UlongToDouble( n1 As Ulong, n2 As Ulong = 0 ) As Double
Asm
  movd xmm1, dword Ptr [n1]
  movd xmm0, dword Ptr [n2]
  punpckldq xmm0, xmm1
  psrlq xmm0, 12
  mov eax, 1
  cvtsi2sd xmm1, eax
  por xmm0, xmm1
  subsd xmm0, xmm1
  movq [Function], xmm0
End Asm
End Function

Sub GetUlongs( n As Ulongint, Byref L1 As Ulong, Byref L2 As Ulong )
  L1=Cast(Ulong Ptr,@n)[0]
  L2=Cast(Ulong Ptr,@n)[1]
End Sub

' Single Ulong To Double
Print UlongToDouble( 0 ) ' should give 0
Print UlongToDouble( 2^30 ) ' should give 2^30/2^32 = 0.25
Print UlongToDouble( 2^31 ) ' should give 2^31/2^32 = 0.5
Print UlongToDouble( 2^32-1 ) ' (2^31-1)/2^32 = 0.9999999997671693563461 calculated to 22 sig figs

' Ulongint And two Ulongs To Double
Dim As Ulong ul1, ul2
GetUlongs( &hffffffffffffffffull, ul1, ul2 )
Print UlongToDouble( ul1, ul2 )
Randomize
ul1 = Culng(Rnd*&hFFFFFFFFul) : ul2 = Culng(Rnd*&hFFFFFFFFul)
Print UlongToDouble( ul1, ul2 )

Sleep

[dafhi]

here's what i do

Code: Select all

type sng_x2
  as Single           a,b
END TYPE

type ulong_x2
  as ULong            a,b
END TYPE

union Union64
  as sng_x2           s
  as ulong_x2         u
  as double           d
  as ULongInt         uu
END UNION

but keep that asm coming

[deltarho[1859]]

@dafhi Hmmm.

Show me how 2^30 maps into 0.25 as I do with 'Print UlongToDouble( 2^30 ) ' should give 2^30/2^32 = 0.25'

[dafhi]

[edited]
so you're creating a new mapping that's human-understandable, where "1" is 1000000.. and 0.5 is 0100000.. ?

[old post]
here's what 0.25 looks like

Code: Select all

Sub ShowBits( p as any ptr, Newline As boolean = true)
  '' init gfx or u will see wrong color
  '' use ScreenRes x,y, 32
  Dim As ulong  I, AryFG(3), AryBG(3)
  var           oldcol = color()
  AryFG(0) = RGB(200,200,200)
  AryFG(1) = RGB(255,0,0)
  AryFG(2) = RGB(0,255,0)
  AryFG(3) = RGB(0,96,255)
  AryBG(0) = RGB(100,100,100)
  AryBG(1) = RGB(127,0,0)
  AryBG(2) = RGB(0,127,0)
  AryBG(3) = RGB(0,0,127)
  dim as ulong ptr  v = p
  For I = 0 To 3
    Color AryFG(I), AryBG(I)
    print bin(*v shr (24 - i*8), 8);
  Next
  If Newline Then Print
  color oldcol
End Sub
sub ShowBits64(p as any ptr, Newline As boolean = true)
  showbits p+4, false
  ShowBits p, newline
END SUB


type sng_x2
  as Single           a,b
END TYPE

type ulong_x2
  as ULong            a,b
END TYPE

union Union64
  as sng_x2           s
  as ulong_x2         u
  as double           d
  as ULongInt         uu
END UNION


screenres 800,600, 32

dim as Union64   n

n.d = .25
showbits64 @n

sleep

[deltarho[1859]]

@dafhi: I think that we are talking two about different subjects.

What I am talking about is:

( 0, 1, 2, ....., 2^32-3, 2^32-2, 2^32-1 ) being mapped into ( 0,1 ].

So, 2^30 maps into 0.25 from 2^30/2^32 = 0.25

[dafhi]

Nice. Thanks for the asm.

[deltarho[1859]]

The asm is an adaptation of code written by a guy at the PureBasic forums, Wilbert, who gave me permission to use it. I did not have to ask but that is what I do.

[dodicat]

You can do a straight map using 0.. 2^31 -> 0...5
I am not sure how to incorporate the second parameter n2.
The asm looks like an or then a subtraction.

Code: Select all

Function UlongToDouble( n1 As Ulong, n2 As Ulong = 0 ) As Double
Asm
  movd xmm1, dword Ptr [n1]
  movd xmm0, dword Ptr [n2]
  punpckldq xmm0, xmm1
  psrlq xmm0, 12
  mov eax, 1
  cvtsi2sd xmm1, eax
  por xmm0, xmm1
  subsd xmm0, xmm1
  movq [Function], xmm0
End Asm
End Function

Sub GetUlongs( n As Ulongint, Byref L1 As Ulong, Byref L2 As Ulong )
  L1=Cast(Ulong Ptr,@n)[0]
  L2=Cast(Ulong Ptr,@n)[1]
End Sub

' Single Ulong To Double
Print UlongToDouble( 0 ) ' should give 0
Print UlongToDouble( 2^30 ) ' should give 2^30/2^32 = 0.25
Print UlongToDouble( 2^31 ) ' should give 2^31/2^32 = 0.5
Print UlongToDouble( 2^32-1 ) ' (2^31-1)/2^32 = 0.9999999997671693563461 calculated to 22 sig figs
Print UlongToDouble( 2018 )
Print UlongToDouble( 100000000 )
'left out using n2
print
print

'======================


function map(a as double,b as double,x as double,c as double,d as double) as double
    return ((d)-(c))*((x)-(a))/((b)-(a))+(c)
end function

Function UlongToDouble2(n1 as ulong,n2 as ulong=0) as double
    'n2??  not ure how to incorporate
    return map(0,2^31,n1,0,.5) 
end function

' Single Ulong To Double
Print UlongToDouble2( 0 ) ' should give 0
Print UlongToDouble2( 2^30 ) ' should give 2^30/2^32 = 0.25
Print UlongToDouble2( 2^31 ) ' should give 2^31/2^32 = 0.5
Print UlongToDouble2( 2^32-1 ) ' (2^31-1)/2^32 = 0.9999999997671693563461 calculated to 22 sig figs
Print UlongToDouble2( 2018 )
Print UlongToDouble2( 100000000 )

'left out using n2

Sleep  

[deltarho[1859]]

@dodicat

At an initial glance, I did not believe that your function map() would give identical results to mine and was astonished when it did.

The asm is faster and, of course, the n2 allows mapping 64 bits of information to give us a 53-bit granularity.

Wow!

@dafhi

With some applications we need to map into ( -1, 1 ] as opposed to ( 0, 1 ]. We could use UlongToDouble(), multiply by 2 and then subtract 1. However, that may be costly so, since you are an asm fan, I have stayed in the asm domain with UlongToDoubleS() where S is short for signed.

Code: Select all

Function UlongToDoubleS( n1 As Ulong, n2 As Ulong = 0 ) As Double
Asm
  movd xmm1, dword Ptr [n1]
  movd xmm0, dword Ptr [n2]
  punpckldq xmm0, xmm1
  psrlq xmm0, 12
  mov eax, 2
  cvtsi2sd xmm1, eax
  por xmm0, xmm1
  subsd xmm0, xmm1
  mov eax, 1
  cvtsi2sd xmm1, eax
  subsd xmm0, xmm1
  movq [Function], xmm0
End Asm
End Function
 
Sub GetUlongs( n As Ulongint, Byref L1 As Ulong, Byref L2 As Ulong )
  L1=Cast(Ulong Ptr,@n)[0]
  L2=Cast(Ulong Ptr,@n)[1]
End Sub
 
' Single Ulong to Double
Print UlongToDoubleS( 0 ) ' should give -1
Print UlongToDoubleS( 2^30 ) ' should give (2^30/2^32)*2 - 1 = -0.5
Print UlongToDoubleS( 2^31 ) ' should give (2^31/2^32)*2 - 1  = 0
Print UlongToDoubleS( 2^32-1 ) ' ((2^32-1)/2^32)*2 - 1 = 0.9999999995343387126923 calculated to 22 sig figs
 
' Ulongint And two Ulongs To Double
Dim As Ulong ul1, ul2
GetUlongs( &hffffffffffffffffull, ul1, ul2 )
Print UlongToDoubleS( ul1, ul2 )
Randomize
ul1 = Culng(Rnd*&hFFFFFFFFul) : ul2 = Culng(Rnd*&hFFFFFFFFul)
Print UlongToDoubleS( ul1, ul2 )
 
Sleep