Floating Point Arithmetics. 1

We recently solved a problem of Daniel with a rather simple computation error in thread184-1738661.

Such calculation errors pop up from time to time and are not a CPU error or VFP bug, it's simply the limitation of floating point arithmetic. I get to the general overview of that later, but you may read an introduction to the problems in Python documentation Floating Point Arithmetic: Issues and Limitation, which shows the general problem is language independent.

Daniel wanted a value truncated at 6 decimal places. The basic computation is x=a/b and to get the truncation of that simple division he used FLOOR(1000000*a/b)/1000000, because FLOOR() like INT() truncates all decimal places, you first shift the number 6 decimal places to the left, then truncate decimal places of that number (which would have been 7th and higher decimal place before the shift) and then shift back 6 places to the right. So you expect 6 decimal places, not rounded up by the 7th decimal place, simply truncated.

That failed with the rather harmless looking example of computing 265.96/20. If you SET DECIMALS TO 6 and let VFP execute ? 265.96/20 you get the correct and precise result of 13.298000. Using Floor(1000000*265.96/20)/1000000 you get 13.297999

Mike Lewis solved that by using ROUND() in a way, which truncates a number to a certain precision instead of rounding it. ROUND has the advantage of having built in the computation to a certain precision and we only need to "tame" it to truncate instead of round. The trick is to subtract the rounding offset of 5 at the next decimal place. In case of 6 decimal places that's ROUND(x - 0.0000005), 6) for positive x, and ROUND(x + 0.0000005), 6) for negative x, or ROUND(x - SIGN(x) * 0.0000005), 6) in general. Round(265.96/20-0.0000005,6) is 13.298000, problem solved.

We can generalize this as a TRNUCATE() function, if you like:

FUNCTION Truncate()
Lparameters tnNumber, tnPrecision

Return Round(tnNumber-0.5*Sign(tnNumber)*10^-tnPrecision,tnPrecision)
ENDFUNC

Later on I'd like to put this new function to the test and discuss other proposed solutions and why they work, too or fail. For now my time is up, though.

Bye, Olaf.

 

[Mike Lewis]

A vey useful summary, Olaf.

A key point is that the floating-point format represents an approximation of actual numbers. The simple reason is that the number of possible numbers in the relevant range is higher than the number of combinations of bits available. This is most often seen when you do a calculation involving numbers with several decimal places, and the result is slightly different than if you had used a fixed-point representation - for example, if you had done the calculation by hand or with a calculator. It's tempting to think of this as a bug, but of course it is merely a limitation of the system.

It will be interesting to examine your TRUNCATE() function further, and to see if there are any circumstances in which it fails.

Mike

__________________________________
Mike Lewis (Edinburgh, Scotland)

Visual FoxPro articles, tips and downloads

 

[Olaf Doschke]

Part2

Reviewing other proposed solutions:

my solutions:
I first proposed Round(265.96/20,6) and only later realized Daniel wanted to truncate, not round.
Then Mike, whoms solution I already explained came up with subtracting 0.0000005
I changed my solution to computing a string instead of the value:

vt=265.96
qq=20
lcNum =  Alltrim(Str(1000000*vt/qq,20))
lcResult = Left(lcNum,Len(lcNum)-6)+"."+Right(lcNum,6)

You can also get back the value by using VAL(lcResult), but indeed that solution fails, because STR rounds, simply try this case:

vt=4.000001
qq=2
lcNum =  Alltrim(Str(1000000*vt/qq,20))
lcResult = Left(lcNum,Len(lcNum)-6)+"."+Right(lcNum,6)

That results in 2.000001 while the exact value is 2.0000005 and truncated should instead be 2.000000

Also it fails for values of vt/qq lower than 1 because of not generating leading zeros.

Nasib Kalsi's solution:

x = 14.736256732 && try 265.96/20 instead
y = 1000000.000000
? int(x*y)/y

The Int() function differs from Floor() in the value range of negative numbers, but this is not the major change. The solution actually fails for x = 265.96/20, but still is worth looking at.
I think Nasib tried to enforce a certain precision by using 1000000.000000 instead of 1000000 and that has an effect, though not sufficient.

You can experiment with adding decimal places to see how that changes VFPs results:

? 0.05 * 1 && prints 0.05
? 0.05 * 1.0 && prints 0.050
? 0.05 * 1.00 && prints 0.0500
? 0.05 * 1.000 && prints 0.05000

Why? VFP is not only converting a decimal number into the binary float representation and approximation, it also stores how many decimals the original decimal representation had. And then it does a side calculation about how many decimal places the result can have and displays that number of decimal places. If you multiply a number with N decimal places with another number of M decimal places, the result can have up to N+M decimal places. You simply have to think about an extreme case. Eg for N=1 and M=2 multiply 0.1 with 0.01 you get 0.001 having B+M=1+2=3 decimal places.

VFP isn't doing that very well, though, eg try ? 0.001^3, this will not be the same as 0.001*0.001*0.001, with the ^ exponent operator VP strays with the number of decimal places or precision of the base number. In division VFP sums the decimal places, too, though the rule doesn't fit for division.

And if the values vt and qq in the expression vt/qq come from a table this info about the decimal places is not stored in a float field.

To get back to Nasib Kalsi, it made sense to add the precision of the factor and divisor 1000000, but it doesn't work out in general, unfortunately.

Vilhelm-Ion Praisach's observation:
He did some experiments to find out the internal value of 265.96/20. It's a useful experiment. The final thought, that string operations word, eg ?VAL(TRANSFORM(value /qty*1000000)) is 13298000. In general TRANSFORM would fail like STR because of rounding and it also doesn't truncate results having additional places. At least the string resulting from TRANSFORM(value /qty*1000000) has to be processed by removing any decimal places and setting a new decimal point in it before using VAL on it.

So far for now. In the next part I'll show some experiments.

 

[vgulielmus]

The number of insignificant zeros affects the result of the ROUND() function
?ROUND(1.0050000000000,2) && Outputs 1

There are a lot of interesting behavior.
I've made a set of experiments that can be reproduced, here:

1) Regarding Inputmask

http://www.foxite.com/archives/interesting-behavior-of-inputmask-0000406050.htm

2) Regarding Transform()

http://www.foxite.com/archives/interesting-behavior-of-inputmask-0000406150.htm

3) Regarding CAST()

http://www.foxite.com/archives/interesting-behavior-of-inputmask-0000406287.htm

4) Regarding ROUND()

http://www.foxite.com/archives/interesting-behavior-of-inputmask-0000408963.htm

Respectfully,
Vilhelm-Ion Praisach
Resita, Romania

 

[vgulielmus]

Some findings from the experiments with ROUND() (extracted from

http://www.foxite.com/archives/interesting-behavior-of-inputmask-0000408963.htm)

The errors begins with certain powers of 2 and takes up to a certain value, much less than the next power of 2.
The difference between the two extremes reveals a very interesting pattern.
The edges are practically the same (for example, if the errors begins for 128 then very likely are spreading until to 163). Clearly at least one of the causes is the internal representation

For i.065
The errors appears to the numbers 16384.065 16385.065 16386.065 ... 20971.065

20971-16384=1000111101011 (base 2)
5242-4096=10001111010 (please compare it with the previous and the next number)

For i.165
10-8=10
83885-65536=100011110101101

For i.265
4
163-128=100011

For i.365
20-16=100
81-64=10001
2621-2048=1000111101

For i.465
327-256=1000111
1310-1024=100011110
41942-32768=10001111010110

For i.565
0
4
5242-4096=10001111010

For i.665
9-8=1
83885-65536=100011110101101

For i.765
4
163-128=100011

For i.865
19-16=11
81-64=10001
2620-2048=1000111100

For i.965
326-256=1000110
1309-1024=100011101
41942-32768=10001111010110


Respectfully,
Vilhelm-Ion Praisach
Resita, Romania

 

[DSummZZZ]

I spent the better part of a week toying with single vs double precision, conversion, etc., for Fox apps that may require it. If you have some spare time and are bored, peruse this FAQ I eventually came up with:

Conversion Functions: Dec,Hex,Binary,Oct Back and Forth
faq184-4461


-Dave Summers-

Even more Fox stuff at:

http://www.dave-summers.net/foxprolinks.htm

 

[dbMark]

This is what I came up with to address the issue a year and a half ago. I'd be happy to hear whether this is a viable solution.

FLOOR, as well as related CEILING and INT, work with the integer portions of numbers. What my example does is shift the decimal point right (positive parm2) or left (negative parm2), applies FLOOR and then shifts the decimal point back. It uses exponents of 10 (10, 100, 1000, etc) when multiplying and then dividing the number in parm1. By using exponents to shift the decimal point the issue of rounding is sidestepped.

FUNCTION FloorDec
   PARAMETERS in_num, in_dec  && both parameters required
   RETURN FLOOR(in_num*10^in_dec)/10^in_dec  && shift decimal point back and forth so FLOOR() preserves desired segment
ENDFUNC

Or that line could be split into separate steps to visualize it better.

RetVal = in_num*10^in_dec
RetVal = FLOOR(retVal)
RetVal = retVal/10^in_dec

Note: We can’t use INT here since it ignores the pos/neg sign status. INT returns the next lower integer number, FLOOR returns the next lower integer based on the value and CEILING returns the next higher integer based on the value. A behavior to keep in mind if the function will be used on negative numbers.

 

[Olaf Doschke]

This is what Daniel used in thread184-1738661, besides - as you point out - it does not use 10^6 but 1000000.

Just testing a special case is no proof, unless it also doesn't work, so I tested:

Set decimals to 6

value = 265.96
qty = 20
? FloorDec(value/qty,6)

FUNCTION FloorDec
   PARAMETERS in_num, in_dec  && both parameters required
   RETURN FLOOR(in_num*10^in_dec)/10^in_dec  && shift decimal point back and forth so FLOOR() preserves desired segment
ENDFUNC

This results in 13.297999 here, so it still has the same problem, using exponents didn't help.

Even testing more directly ? FLOOR(265.96/20*10^6)/10^6 results in the same wrong decimal places.
So definitely this has a counter example proving it doesn't work out. Using exponents is not the solution.

It's always a possibility some different computation of the same value solves a problem with certain numbers, but has problems with other numbers. That's why I'd like to test drive this more and report back. I just don't have the time to dive into this single problem. I keep answering other questions as a higher priority, besides working of course

But this isn't forgotten. To have a reference I need a computation based on decimal math, computing with strings of numbers without putting them into numeric variables, so you can judge if certain egde cases of floating point binary numbers work out with a certain function or not. Decimal math libraries exist in other languages, eg a quick google found a library for TCL here: [URL unfurl="true"]http://docs.activestate.com/activetcl/8.4/tcllib/math/decimal.html[/url]

Bye, Olaf.

 

[dbMark]

Back when dBase III was introduced in the mid-1980s the app had issues with handling numbers. If you added, subtracted, etc enough numbers the could not be compared to equal what was expected. As I recall it was explained that the floating point arithmetic used was accurate only to about 15 digits and they were comparing too far out, maybe to 18 digits. So they released dBase III Plus which reduced the sensitivity of comparisons to 15 digits and perhaps other tweaks under the hood. I believe Visual FoxPro handles floating point numbers similarly, discarding minute values at the far right than cannot be relied upon to be accurate.

So the code I posted above has similar limitations. Use a number with enough digits where the 9's extend too far to the right and it could probably still hit the limitations of floating point numbers as implemented by VFP.

 

[Olaf Doschke]

dbMark,

this example I gave has the correct result of 13.298 and that's only 3 of the 6 decimal places accuracy wanted.

It's wrong to assume 15 decimal places accuracy mean you can rely on 15 decimal places of any computation. The edge cases I talked about can mean a cascade of rounding you see at earlier places, already, as this example shows.

The problem is simply, that you can't even store 1/20 exact, needless to speak of 265.96/20.

The accuracy of 15 decimal places is just a comparison of the binary accuracy of the mantisse part of a floating point number. A double float has 53 bits mantissa, which means the lowest bit has a value of 2^-53, which compares to about 10^-16. computational 53*Log(2)/Log(10), but it doesn't mean the accuracy of all computations is that high, you get errors simply letting the complirer convert a decimal number in your source code to the floating point value in memory and computation errors have error propagation, so in the end your value can even get wrong in higher decimal places than 10^-16 or 10^-15. Sometimes you get the exact result in decimal simply because the reconversion of the float to a decimal display cancels out the inexactness. So eg ? 0.1 prints 0.1, while in memory 0.1 can't be stored exactly. ? 0.1 printing out 0.1 is no proof, 0.1 can be stored exact, we know by definition of the float format it can't be stored exactly. And 0.05 or 1/20 is just shifted one bit.

You get correct results most of the time because of errors cancelling out each other instead of accumulating and you can easily show a simple effect: Create a million random numbers with high variance, that means with very high and very low exponents, eg ? RAND()*10^(RAND()*80-40), store them in a double field. Sum them in the order created and then sum them sorted from lowest to highest value. You get a more accurate sum by summing the sorted values. The reasoning is simple. If you sum ~10^-40 and ~*10^32 you get 20^32, the lower value is below the accuracy of the hight value and so doesn't contribute to the sum. If you sum many 10^-40 values they sum up to higher values having a chance to contribute to the sum.

Bye, Olaf.

 

[Olaf Doschke]

corrected: ...If you sum ~10^-40 and ~10^32 you get 10^32... (as if you added 0 - that also already works with less difference in the exponents, 16 is enough)

Bye, Olaf.

 

[vgulielmus]

I guess VFP has some specific issues.

?ROUND(1.0050000000000,2)==1             && .T.
?ROUND(1.005000000000,2)==1.01           && .T.

* The previous is and exception
?ROUND(2.0050000000000,2)==2.01          && .T.

* A way to get this results
CREATE CURSOR xx (q N(20,7),p N(20,6))
INSERT INTO xx VALUES (1,1.005)
?ROUND(xx.q*xx.p,2)                      && 1
CREATE CURSOR yy (q N(20,6),p N(20,6))
INSERT INTO yy VALUES (1,1.005)
?ROUND(yy.q*yy.p,2)                      && 1.01

* The previous is and exception
INSERT INTO xx VALUES (2,1.005)
?ROUND(xx.q*xx.p,2)                      && 2.01

Here

http://www.foxite.com/archives/interesting-behavior-of-inputmask-0000408963.htm

is a little more about this

Respectfully,
Vilhelm-Ion Praisach
Resita, Romania

 

[dbMark]

Hmm, I thought I had good code, we've used it without issues until now. Here is try #2, at your convenience? I don't know if I should use 18 decimals, but that's the maximum allowed so I did.

My working theory is, if the value is off but I can see it as I want it, then use what I see - and don't place it back into a variable until I finish calculating! Likely my code had worked for me since my data was picked up from a field and not a variable and that's how I avoided the floating point issue?

FUNCTION FloorDec
   PARAMETERS in_num, in_dec  && both parameters required
   * Clean up in_num in case it is floating point data by placing it into a string
   * Note: Do not store intermediate value into a variable or floating point issue returns
   RETURN FLOOR( VAL( LEFT( STR(in_num, 50, 18) , 50) ) *10^in_dec)/10^in_dec
ENDFUNC

 

[Olaf Doschke]

Picking up a value from a field only differs from a variable, if it's a field type different from double float, eg numeric fields really store strings of numbers, but while calculating everything is double floats.

Well, you can try out your function yourself, why don't you do that?

I already said

It's always a possibility some different computation of the same value solves a problem with certain numbers, but has problems with other numbers.

Which means: If some change of code solves a problem with eg 265.96/20, it can still have precision problems with other numbers and vice versa. Your code can be fine for you and you may stay with it. If you want precise math even in simple calculations, there is no perfect way to work with floating point numbers, you'd need to work out a library working with strings instead of doubles and even that won't allow you to store numbers like 1/3 exact. The problem of endless numbers is limiting you in any system. Calulating with decimals would only less often lead to errors in the domain of decimal places, unless you have such things as 1/3.

Think about it this way: How many instances of other errors did you have, typos, syntax error in your code, hardware defects. Compared to the harm floating point imperfection can have, it's neglectable besiodes other errors in terms of costss and time invested to fix them. We just often need to justify computational math, as end users have no idea how computers compute, even though it's the core application they were invented for.

Bye, Olaf.

 

[Olaf Doschke]

One thing about your new FloorDec() function: It uses STR and that rounds. I proposed that and already discussed, why it wouldn't fit to the problem of truncation with other values than 265.96/20, simply reread here, try ? FloorDec(10000000000.000001/2,6) it should be 5000000000.000000, but it's 5000000000.000001 rounded up.

Bye, Olaf.

 

[dbMark]

And I got inconsistent results.
SET DECIMALS TO 6
? 10000000000.000001
10000000000.000000
? FloorDec(10000000000.000001/2,6)
5000000000.000001

* One digit less:
? 1000000000.000001
1000000000.000001
? FloorDec(1000000000.000001/2,6)
500000000.000000

* One digit more:
? 100000000000.000001
100000000000.000001
? FloorDec(1000000000.000001/2,6)
50000000000.000000

 

[Olaf Doschke]

Yes, that's the nature of such things, you can't predict it, you can't even rely on a certain range of "normal" values working exact, sometimes the cancelling out of errors makes longer numbers work again, as in your last sample. You might find random numbers with only 6 digits also not working out.

Bye, Olaf.

 

[NasibKalsi]

Why? VFP is not only converting a decimal number into the binary float representation and approximation, it also stores how many decimals the original decimal representation had. And then it does a side calculation about how many decimal places the result can have and displays that number of decimal places. If you multiply a number with N decimal places with another number of M decimal places, the result can have up to N+M decimal places. You simply have to think about an extreme case. Eg for N=1 and M=2 multiply 0.1 with 0.01 you get 0.001 having B+M=1+2=3 decimal places.

VFP isn't doing that very well, though, eg try ? 0.001^3, this will not be the same as 0.001*0.001*0.001, with the ^ exponent operator VP strays with the number of decimal places or precision of the base number. In division VFP sums the decimal places, too, though the rule doesn't fit for division.

I agreed very much so. Taking it further, I created c# (2012) dll and then tried the solution. It always seem to work so far.

such as

objMath = CreateObject("FoxCOM.FoxLib")
x = 265.96
y = 20
z = 1000000  && Accuracy to 6 decimal places
? int(objMath.Divide(x,y)*z)/z  && Not sure divide by z would fail for some numbers

If anyone require Fox Library (FoxLib) code, I can post it here. The example, how to create dll was taken from posts in this forum and google.

 

[Olaf Doschke]

The difficulty is deciding what results are correct. Correcting the one case, which doesn't work with VFP floating point arithmetics is no proof, but it's obvious any other language has other cases of failing with floats. It's simply in the nature of these numbers.

I'm still building up a string math library and am now at the division, essential for this case of 265.96/20. What I would like to test then is whether Mike Lewis' idea of using a modified ROUND is enough, then we can use that as ideal way for truncation, working better than INT or FLOOR and won't need an external library.

My string math library currently compute sums and products with exact precision only limited by string length. Of course it computes much slower. Division is something not letting you deduct the number of needed decimal places by the number of decimal places of the operands. Think alone of the simple case 1/3. You get infinite decimal places, but it's a repeating decimal, which can be notated as finite number. If you divide by such numbers it can get messy, though.

In the end you have to live with a finite precision in division, too. Even if you think of multiplying with the inverse number, that means dividing 1 by it, so a division is unavoidable. Fortunately the case of 265.96/20 is a finite result, so that will not be an edge case, but of course I want to test very many floating point numbers to find other non working cases.

I expect to have a result like ROUND works correct up to a certain number of decimal places or the percentage of miscalculation cases grows in a certain way when increasing the decimal places you want.

And as a side solution the string math library can always compute a certain precision by the definition of the division algorithm used. It'll be something we all learned in school as long division or handwritten division, see

http://en.wikipedia.org/wiki/Long_division

Bye, Olaf.