double trouble

[deltarho[1859]]

@BasicCoder2

There is no need for the test 'If w1(i)<>0 Then'. If w1(i) = 0 then 'w1(i)=Int(w1(i))/100' will still be zero.

[caseih]

Ain't floating point fun? Edited this post to make it shorter and maybe less presumptuous.

There is no getting away from the fact that 93/100 is 0.93. With FreeBASIC 'Print 93/100' gives 0.9300000000000001. With PowerBASIC we get 0.93.

Whatever the reason for the difference, FreeBASIC is getting it wrong.

The question is: Should FreeBASIC be corrected?

The only way to correct it is to use a different type that does not involve binary fractions. All numbers based on fractions can have this issue, regardless of the base. For example "1/3." In all cases no matter how many decimal places you work it out to, the answer will always be less than 1/3.

But let's get back to binary. Binary has similar problems, but with different numbers since it's variations of 1/2^n. For example, here is a number that cannot be represented accurately in IEEE floating point: 0.2. No matter how you work it out as a series of binary fractions, the binary representation will always be slightly less than 0.2 (0.0011001100110011.. etc). Let's consider different "roudings" of the binary number:
0.0011 = 0.1875
0.00110011 = 0.19921875
0.001100110011 = 0.199951171
etc.

Let's consider 0.93:
0.93 = 0.1110 1110 0001 0100 0111 1010...
There isn't a good way of truncating (rounding) the binary digits to get 0.93.

FreeBASIC isn't getting it wrong per se. It's IEEE floating point getting it wrong. For example even simple math doesn't "work" in binary: 0.1 + 0.2 can never quite be equal to 0.3. At certain precisions it's close enough to not matter of course.

This problem cannot be completely removed, but you can certainly hide it. There are many articles on inaccuracies in MS Excel and how to deal with them. And there are techniques to use when writing code to convert a binary number to decimal that can help, such as setting a threshhold on the number of zeros in a row, or number of 9s in a row.

[jevans4949]

You have the same problem with pocket calculators; the more sophisticated ones get round it by holding your intermediate result as an (improper) fraction. Try doing 1/3+1/3+1/3 on older calculators and you will get 0.999999. A modern one will have internally calculated 3/3 and will give you the correct answer. At some point the calculator will have to convert to decimal or binary to perform an operation, eg square root or something.

In banking systems, where accuracy is important, currency amounts will be defined (in binary) internally in minor units (cents), or else the cents will be held in a separate field, so that some calculations, (e.g. charges) can easily be done on whole major units. Often the minor units are held as positive, and the major units can go negative - so if your account is 5 cents overdrawn it is held as (-1+.95) dollars. Before British currency went decimal in 1971, there wore 240 pence in a pound, which could under this system be stored in a byte.

If you need calculation to a known number of decimal places, you can calculate in binary in hundredths (or whatever) and when you need to print, convert to character and edit in the decimal point where it applies. Of course you still have to decide about rounding following a division, by examining the remainder. Different countries have different rules about rounding when calculating the percetage due as tax as well.

[deltarho[1859]]

@caseih

OK, points taken.

If we look at BasicCoders' code the statements using arithmetic don't need double precision – single precision is more than adequate.

If we declared w1() as Single the problem ceases to exist. The underlying issues persist, but BasicCoders's concern doesn't.

[caseih]

Good point. There are lots of gotchas with IEEE floating point, and precision doesn't always give you accuracy. Perhaps with single precision, when printing out the number to the screen, PRINT stops far earlier, knowing that the number of decimal places possible are much less. So the extra .00000000000...1 is simply lost in the calculations, and the rules that PRINT uses (actually c library scanf) to generate the digits sees the zeros and truncates it. There are research papers published on how to best convert and print a floating point number.

Another thing I thought of was that maybe if you use double precision for the actual math, and then truncate the result to single precision, that might work better sometimes also. Calculators may use a certain number of bits precision, but truncate the decimal result to 10 decimal places.

[caseih]

the more sophisticated ones get round it by holding your intermediate result as an (improper) fraction.

IEEE actually defined a fractional type where numbers are stored as essentially two integers. I think this was sometimes called "floating slash." But it never caught on as far as hardware support goes. Theoretically it allows arbitrary decimal accuracy, obviously at the expense of needing a lot more bits to encode things.

[deltarho[1859]]

Another thing I thought of was that maybe if you use double precision for the actual math, and then truncate the result to single precision, that might work better sometimes also.

That is what I was doing with my 'trick': Print Cast(Single,w1(i))

[hhr]

I would do the printout something like this:

Code: Select all

Sub sPrint(n As Double, decimals As Byte)
   Dim As Double a,b
   a = Fix(n)
   b = Abs(Frac(n))
   b = Int(b*(10^decimals)+0.5)
   Print a & "." & b 'Strangely enough, it works.
End Sub

Dim As Single a
Dim As Double b
b = 0.93
a = b
Print "single: ";a
Print "double: ";b
Print "shortened double: ";
sPrint b,2

Sleep

I have done a little research into Quasi Monte Carlo Integration.
There it was important for me to consider the numerical instabilities, where rounding errors lead to disturbances of the result.
The 'trick' is a good idea here.

I can't try PowerBasic, but I would find it interesting how PowerBasic reacts to a program like this:

Code: Select all

Dim As Double a,s,n
a = 1/3
Do
   s += a
   n += 1
   Print a,s/n
Loop

[srvaldez]

hhr
I changed your code a bit

Code: Select all

#COMPILE EXE
#DIM ALL

FUNCTION PBMAIN () AS LONG

    DIM a AS DOUBLE, s AS DOUBLE, n AS DOUBLE
    DIM i AS LONG
    a = 1/3
    FOR i=1 TO 100
       s += a
       n += 1
       PRINT i, a, s/n
    NEXT i

    WAITKEY$
END FUNCTION

only lines 42 thru 49 print .333333333333334 the rest print .333333333333333

Code: Select all

 41            .333333333333333            .333333333333333
 42            .333333333333333            .333333333333334
 43            .333333333333333            .333333333333334
 44            .333333333333333            .333333333333334
 45            .333333333333333            .333333333333334
 46            .333333333333333            .333333333333334
 47            .333333333333333            .333333333333334
 48            .333333333333333            .333333333333334
 49            .333333333333333            .333333333333334
 50            .333333333333333            .333333333333333
 

[deltarho[1859]]

@hhr

Your first code looked good, but we only need one failure to drag us back to the drawing board. I didn't examine the code.

Try b = 85.07. I get:

Code: Select all

single:  85.07
double:  85.06999999999999
shortened double: 85.7

[hhr]

@srvaldez
It seems PowerBasic prints out 15 digits, FreeBasic prints 16 digits.

@deltarho
Yes, you are right, the line with the comment is wrong. I have to use strings. I will try tomorrow.

[srvaldez]

I said in a previous post that since double only holds about 15.95 digits that the default print precision should be 15
single is better off since it holds about 7.2 digits, so for single printing 7 digits should be ok
[edit] single has a problem in the seventh digit, the output is virtually the same in FB and PB, so PB's only advantage is to limit the output of double to 15 digits

Code: Select all

 87            .3333333      .3333335
 88            .3333333      .3333335
 89            .3333333      .3333335
 90            .3333333      .3333335
 91            .3333333      .3333335
 92            .3333333      .3333336
 93            .3333333      .3333336
 94            .3333333      .3333336
 95            .3333333      .3333336
 96            .3333333      .3333336
 97            .3333333      .3333336
 98            .3333333      .3333335
 99            .3333333      .3333335
 100           .3333333      .3333335
 

but caseih already explained why, there's no point in this as it's part of using binary floating point

[caseih]

I suspect (know, actually) that PB implemented its own printing routine, and made a good choice to limit to 15 digits. Whereas FB's runtime uses the Standard C Library's sscanf() or similar to do the conversion to string, and I guess the C library decided on 16 digits for some reason. Would be interesting to see if the same choice is made in other C libraries. GLibC, uLibc, libc, newlib, the MSVCRT, etc. I'm not sure if the standard dictates this sort of thing or not.

[caseih]

@hhr

Your first code looked good, but we only need one failure to drag us back to the drawing board. I didn't examine the code.

I think the problem is that when looking at the fractional part of the number, leading zeros are important, but trailing zeros are not. However any INTEGER type is going to throwaway the leading zeros and keep the trailing zeros.

[deltarho[1859]]

This thread could rumble on for some time.

BasicCoder's opening issue is resolved by simply declaring w1() as Single.