Double inacurracy

Rescator · Post by **Rescator** » Wed Mar 15, 2006 12:30 am

l.f=1.0
r.f=1.0
c.f=1.0
sl.f=1.0
sr.f=1.0
lt.f=(l*(c*0.707)*(sl*-0.8165)*(sr*-0.5774))*3.000182
rt.f=(r*(c*0.707)*(sl*0.5774)*(sr*0.8165))*3.000182
Debug lt
Debug rt
Debug ""

dl.d=1.0
dr.d=1.0
dc.d=1.0
dsl.d=1.0
dsr.d=1.0
dlt.d=(dl*(dc*0.707)*(dsl*-0.8165)*(dsr*-0.5774))*3.000182113754469
drt.d=(dr*(dc*0.707)*(dsl*0.5774)*(dsr*0.8165))*3.000182113754469
Debug dlt
Debug drt

The drt returns 0.99999999999999989 when testing here,
the other results all return 1.0 as expected/intended.

Post by **Fred** » Wed Mar 15, 2006 12:51 am

Doubles are float, you can't expect a full precise value. The 'Debug' command use a 20 decimal precision, you could try with StrD() which should round up better.

Rescator · Post by **Rescator** » Wed Mar 15, 2006 1:11 am

Yeah! Your right.
But why is PB "limping" on the drt double but not the dlt double?
While the normal float example is balanced.

I was planning to use doubles for the higher precision for a sound routine,
but it looks like I have to go back to using singles instead as those at least
round properly for the two values.

Which is what I really don't get.

If this is a float rounding issue then both the double results should be equally "incorrect", should they not?

Looks to me like PB gets rounding issues with doubles if negative numbers are used,
while singles work properly.

So I'm not resting on this case yet Fred

Either make singles "limp" the same way,
or make doubles act same way as singles.
Something is "off" here, even for floats

Rescator · Post by **Rescator** » Wed Mar 15, 2006 1:32 am

Here is a slightly different code.

Please note that here both singles and float behave consistently.
Altough there is a rounding error, at least it's identical both for singles and doubles.

But in the code in my first post you see that doubles is not consistent with the singles. Hence why I'm insisting doubles has a bug in the first code snippet

Code: Select all

l.f=1.0
r.f=1.0
c.f=1.0
sl.f=1.0
sr.f=1.0
lt.f=(l*0.3225)+(c*0.2280)+(sl*-0.2633)+(sr*-0.1862)
rt.f=(r*0.3225)+(c*0.2280)+(sl*0.1862)+(sr*0.2633)
Debug lt
Debug rt
Debug ""

dl.d=1.0
dr.d=1.0
dc.d=1.0
dsl.d=1.0
dsr.d=1.0
dlt.f=(dl*0.3225)+(dc*0.2280)+(dsl*-0.2633)+(dsr*-0.1862)
drt.f=(dr*0.3225)+(dc*0.2280)+(dsl*0.1862)+(dsr*0.2633)
Debug dlt
Debug drt

Post by **Fred** » Wed Mar 15, 2006 1:37 am

In fact the problem comes from the Debug command which display the 'Doubles' with much more accuracies than a 'Float'. My take is the float is rounding the value to 1.0 where the double keeps a bigger precision.

Try to change the '18' to '20' and see the difference. On the float it doesn't have any effect as the float range is much more limited, so it's already rounded to 1.0.

Code: Select all

l.f=1.0
r.f=1.0
c.f=1.0
sl.f=1.0
sr.f=1.0
lt.f=(l*(c*0.707)*(sl*-0.8165)*(sr*-0.5774))*3.000182
rt.f=(r*(c*0.707)*(sl*0.5774)*(sr*0.8165))*3.000182
Debug lt
Debug rt
Debug ""

MessageRequester("",StrF(rt))

dl.d=1.0
dr.d=1.0
dc.d=1.0
dsl.d=1.0
dsr.d=1.0
dlt.d=(dl*(dc*0.707)*(dsl*-0.8165)*(dsr*-0.5774))*3.000182113754469
drt.d=(dr*(dc*0.707)*(dsl*0.5774)*(dsr*0.8165))*3.000182113754469
Debug dlt
Debug drt

MessageRequester("",StrD(drt, 17))

Post by **Dare2** » Wed Mar 15, 2006 2:08 am

I think we are talking presentation here, not accuracy.

You can test PB's double handling by comparing with output from other languages, eg:

Create with Pure, and test created with Pure and with FreeBasic.

Code: Select all

CreateFile(0,"C:\code\PureBasic_4\Scrappy\workAndJunk\threeDoubles.dta")
myDbl.d = 1.0 / 3.0 : WriteDouble(0,myDbl)
myDbl.d = 1.0 / 6.0 : WriteDouble(0,myDbl)
myDbl.d = 22.0 / 7.0 : WriteDouble(0,myDbl)
CloseFile(0)

; to test

Debug "PB"
OpenFile(0,"C:\code\PureBasic_4\Scrappy\workAndJunk\threeDoubles.dta")
myDbl.d = ReadDouble(0) : Debug StrD(myDbl,14)
myDbl.d = ReadDouble(0) : Debug StrD(myDbl,14)
myDbl.d = ReadDouble(0) : Debug StrD(myDbl,14)
CloseFile(0)

Debug "FB"
OpenFile(0,"C:\code\PureBasic_4\Scrappy\workAndJunk\threeDoublesFB.dta")
myDbl.d = ReadDouble(0) : Debug StrD(myDbl,14)
myDbl.d = ReadDouble(0) : Debug StrD(myDbl,14)
myDbl.d = ReadDouble(0) : Debug StrD(myDbl,14)
CloseFile(0)

Create with Free, and test created with Free and with Pure.

Code: Select all

	dim d as double

	open "C:\code\PureBasic_4\Scrappy\workAndJunk\threeDoublesFB.dta" for binary as #1
    d = 1.0 / 3.0 : put #1, , d
    d = 1.0 / 6.0 : put #1, , d
    d = 22.0 / 7.0 : put #1, , d
	close #1

' to test
    print "PB"
	open "C:\code\PureBasic_4\Scrappy\workAndJunk\threeDoubles.dta" for binary as #1
    get #1, , d : print d
    get #1, , d : print d
    get #1, , d : print d
    close #1

    print "FB"
	open "C:\code\PureBasic_4\Scrappy\workAndJunk\threeDoublesFB.dta" for binary as #1
    get #1, , d : print d
    get #1, , d : print d
    get #1, , d : print d
    close #1

sleep

These produce the same result on my box. If you hex view the two files, you will see differences, but that is because of small exponent things, and if you write something to parse the stored values yourself you will get the same results.

The mantissa values are the same.

Haven't done it with other langs, but pretty confident it will be ok.

Rescator · Post by **Rescator** » Wed Mar 15, 2006 2:09 am

Man this is a headache indeed!

Ok. so. I think I finaly was able to "normalize" the values.

This seems to provide what I want.

Code: Select all

lt.f=(1.0*0.707*-0.8165*-0.5774)*3.000182113754469
rt.f=(1.0*0.707*0.5774*0.8165)*3.0001821137544694
Debug lt
Debug rt
Debug ""

dlt.d=(1.0*0.707*-0.8165*-0.5774)*3.000182113754469
drt.d=(1.0*0.707*0.5774*0.8165)*3.0001821137544694
Debug dlt
Debug drt

What I don't get though is... why does drt line need that 4 to the end,
but not the dlt, is really float math in general that "off"
when using positive vs negative values?

Gah they do.... I just checked using windows calculator.

Entering (1.0*0.707*-0.8165*-0.5774)*3.000182113754469
produce 1,0000000000000000671175593
entering (1.0*0.707*0.5774*0.8165)*3.000182113754469
produce 1,000546308100000671175593

Great, that means I gotta recalculate everything, dammit

So the positive math is "correct" I can safely assume,
but when float (single or doubles) seems to have a "error" when negatives are used.

No fault of PB I assume, my guess this is a "error" in the float standards,
or rather a natural bias due to the (pun intended) floating behaviour of floats.

Ok, consider this a non-bug then.
But may I request a note to be added to PB manual?
That negative float math may produce slightly different result than positive float math. I don't think many are aware of this fact. I sure didn't know this.

Thanks Fred, you just restored my sanity

EDIT: Damn Dare, you ninja posted me

Post by **Dare2** » Wed Mar 15, 2006 2:16 am

lol.

Rescator · Post by **Rescator** » Wed Mar 15, 2006 2:42 am

*hands Fred a giant metal pole wih nails in it*

Please slap me, really hard.

I'm so stupid, so very stupid. I forgot the basics of how fraction math works.
(and shame on the rest of you for not noticing it either)

Please note the order of the values..

dlt.d=(0.707*-0.8165*-0.5774)*3.000182113754469
drt.d=(0.707*0.5774*0.8165)*3.000182113754469
Debug dlt
Debug drt

dlt.d=(0.707*-0.8165*-0.5774)*3.000182113754469
drt.d=(0.707*0.8165*0.5774)*3.000182113754469
Debug dlt
Debug drt

*goes into a corner and hugs the wall*

Post by **Dare2** » Wed Mar 15, 2006 3:23 am

Finds another corner. (Not sure why, but seems like the thing to do).

blueznl · Post by **blueznl** » Wed Mar 15, 2006 10:22 am

corners for sale! good price, do i sense any customers?

Xombie · Post by **Xombie** » Wed Mar 15, 2006 8:43 pm

I have a question while we're talking about doubles and accuracy.

How come a library like jack's F64 library can handle doubles with what seems like more accuracy than PBs? I'm not knocking on PB - I just want to know what the difference is.

And that reminds me, I need to ask jack to update his library for PBv4

I know you're watching jack!

Post by **Fred** » Wed Mar 15, 2006 8:58 pm

If it deals only with 64 bits floats in internal, it would be very strange and may be a bug a PB.

Post by **Dare2** » Wed Mar 15, 2006 11:27 pm

Hi Xombie!

The hex values of Doubles created by PB using the code I posted above are:

(As stored on disk)

55 55 55 55 55 55 D5 3F ; 1/3
55 55 55 55 55 55 C5 3F ; 1/6
49 92 24 49 92 24 09 49 ; 22/7

Is it possible to give the 24 hex codes for the same values when the result is via F64? It would be interesting to see them. For example, do the 8 byte groups start (end) with F4, 61 and 43 respectively?

Or perhaps a hex dump of a PureBasic value and an F64 value where, from a simple arithmetic expression, the F64 result is more accurate?

I do not have the same level of expertise and understanding as people like jack and ElChoni, but the subject does interest me and I am curious as to how F64 results are stored internally.

Edited: Took out a mistake!

jack · Post by **jack** » Thu Mar 16, 2006 2:41 am

I think it's probably the PB Double to string that's making appear inaccurate.

Code: Select all

Procedure$ myHex(n.b)
   c.w=PeekB(@n) & $ff
   h.w=c/16
   l.w=c-h*16
   If h<10
     hb$=Chr(h+48)
   Else
     hb$=Chr(h+55)
   EndIf
   If l<10
     lb$=Chr(l+48)
   Else
     lb$=Chr(l+55)
   EndIf
   mhex$=hb$+lb$
  ProcedureReturn mhex$
EndProcedure

x.d
d1.d
d2.d
d3.d

! finit ;80 bit precision
! fld1
! fld tword [three]
! fdivp st1,st0
! fstp qword [v_d1]
! fld1
! fld tword [six]
! fdivp st1,st0
! fstp qword [v_d2]
! fld tword [twenty2]
! fld tword [seven]
! fdivp st1,st0
! fstp qword [v_d3]

x=1.0 / 3.0
xs.s{8}
d1s.s{8}
d2s.s{8}
d3s.s{8}
PokeD(@xs,x)
PokeD(@d1s,d1)
PokeD(@d2s,d2)
PokeD(@d3s,d3)
s.s="PB  1/3  "
For i=1 To 8
  s=s+myHex(Asc(Mid(xs,i,1))) 
Next
Debug s
s.s="FPU 1/3  "
For i=1 To 8
  s=s+myHex(Asc(Mid(d1s,i,1))) 
Next
Debug s
x=1.0/6.0
PokeD(@xs,x)
s.s="PB  1/6  "
For i=1 To 8
  s=s+myHex(Asc(Mid(xs,i,1))) 
Next
Debug s
s.s="FPU 1/6  "
For i=1 To 8
  s=s+myHex(Asc(Mid(d2s,i,1))) 
Next
Debug s
x=22.0/7.0
PokeD(@xs,x)
s.s="PB  22/7  "
For i=1 To 8
  s=s+myHex(Asc(Mid(xs,i,1))) 
Next
Debug s
s.s="FPU 22/7 "
For i=1 To 8
  s=s+myHex(Asc(Mid(d3s,i,1))) 
Next
Debug s

End ;put an end so program wont't run into data section

!section '.data' code readable writeable 
! three:  dt 3.0
! six:    dt 6.0
! twenty2:dt 22.0
! seven:  dt 7.0

in windows you can use sprintf.

Code: Select all

ImportC "MSVCRT.LIB" 
  sprintf.l(result.s,num_format.s,number.d) As "_sprintf" 
EndImport 
x.d=3.1415926535897932
s.s=Space(30)
sprintf(s,"%16.15f",x)
Debug s
s.s=Space(30)
sprintf(s,"%16.15e",x)
Debug s