FPU Rounding

[HSE]

  JJ

I was making things complex but is a good idea to make a more simple code.

Code Select Expand

include \masm32\include64\masm64rt.inc
.data
    f1 REAL8 123.45678901234567890
    f2 real8 100.0
    f3 REAL8 0.01
    f4 REAL8 0.0
    buffer db 248 dup (0)
    pbuff qword offset buffer
   
.code
entry_point proc
    movsd xmm2,f1
    mulsd xmm2,f2
    cvtsd2si rcx,xmm2
    cvtsi2sd xmm2, rcx
    mulsd xmm2,f3
    movsd f4, xmm2

    invoke vc_sprintf, pbuff, chr$("%f",10), f4
    conout pbuff, lf
   
    waitkey
    invoke ExitProcess,0
entry_point endp
end

Code Select Expand

123.460000

Press any key to continue...

[JK]

Sorry for the delay (no time yesterday) and thanks for your answers.

My first thought was multiply, round, divide with code like this:

Code Select Expand

.data

n   qword 100                                         ;10^2 -> round to 2 decimal places
float1 real8 123.4567890
float2 real8 123.4567890E100


.code

  int 3

  finit                                               ;set max. precision
  FILD n                                              ;10^(# of places)
  FLD float1                                          ;float
  FMUL ST(0),ST(1)                                    ;x 10^(# of places)
  FRNDINT                                             ;round to integer
  FDIVRP                                              ;/ 10^(# of places)
  FSTP ST(0)


  FILD n                                              ;10^(# of places)
  FLD float2                                          ;float
  FMUL ST(0),ST(1)                                    ;x 10^(# of places)
  FRNDINT                                             ;round to integer
  FDIVRP                                              ;/ 10^(# of places)
  FSTP ST(0)



... works lika charm in qword range, but as you can see fails miserably with larger numbers. This is where i´m stuck.

Mathematically seen large numbers can be rounded just as small numbers, there is no qword limit. So maybe i must cut off the exponent part, round the mantissa and then add the exponent again?

... and what if don´t want to round only the decimal places after the decimal point, but round to a total number of decimal places. I would have to normalize the number before rounding (e.g. 100.234 -> 1.00234E2, round to 4 places -> 1.002E3 = 100.2) - how would i do that? 

[mineiro]

Hi mineiro,
I've tried to test your code, which looks like Masm64 SDK code, but I get different results, see attachment. Where are printf and CStr defined? I can't find them in the Masm64 SDK, so I used vc_printf and chr$().

Sorry sir jj2007, I was playing in Linux. Change in your attached file following lines and will work on win64:

Code Select Expand


.code
WinMain proc
LOCAL msg:MSG
movsd xmm0,f2           ;123.45678901234567890
cvtsi2sd xmm1,one       ;1
cvtsi2sd xmm2,hundred   ;100
divsd xmm1,xmm2         ;1/100
movsd _1_100,xmm1       
mulsd xmm0,xmm2         ;12345.6789...
cvttsd2si rdx,xmm0      ;12345              ;<<<<<<<
movsd temp,xmm0
invoke vc_printf, chr$("%d",13,10), rdx     ;<<<<<<<<
movsd xmm0,temp         ;12345.6789...
mulsd xmm0,_1_100       ;1/100=0.01
cvttsd2si rdx,xmm0      ;123                ;<<<<<<<
invoke vc_printf, chr$("%d",13,10), rdx     ;<<<<<<<
...


This is a win64 code:

Code Select Expand


;ml64 /c round.asm
;link /subsystem:console /entry:WinMainCrt round.obj

include \masm64\include64\masm64rt.inc ; *** Masm64 SDK window template ***
.data
align 16 ; [url=http://masm32.com/board/index.php?topic=10659.msg118268#msg118268]mineiro[/url]
f2 REAL8 123.45678901234567890
one dq 1
hundred dq 100
_1_100 dq 0 ;1/100
temp dq 0

.code
public WinMainCrt
WinMainCrt::
movsd xmm0,f2           ;123.45678901234567890
cvtsi2sd xmm1,one       ;1
cvtsi2sd xmm2,hundred   ;100
divsd xmm1,xmm2         ;1/100
movsd _1_100,xmm1       
mulsd xmm0,xmm2         ;12345.6789...
cvttsd2si rdx,xmm0      ;12345              ;<<<<<<<
movsd temp,xmm0
invoke vc_printf, chr$("%d",13,10), rdx     ;<<<<<<<<
movsd xmm0,temp         ;12345.6789...
mulsd xmm0,_1_100       ;1/100=0.01
cvttsd2si rdx,xmm0      ;123                ;<<<<<<<
mov rax,0
invoke vc_printf, chr$("%d",13,10), rdx     ;<<<<<<<
ret
end


[mineiro]

... works lika charm in qword range, but as you can see fails miserably with larger numbers. This is where i´m stuck.

Code Select Expand


Precision loss. Using Real4 we have 8 decimal digits of precision (lossless), with real8 we have 15 digits (maybe 16 or 17). An example is when we print pi, loss occurs.

pi real8    3.141592653589793 238462643383279502884197169399
after print 3.141592653589793 1159979634685441851615906

This occurs because the number of digits in a continuous division always tends to increase. If we remove the integer part at each step:
3.1415
3.1415-3=0.1415                       3
1/0.1415=7.067137809-7=               7
1/0.067137809=14.894736842-14=       14
1/0.894736842=1.117647059-1=          1
1/0.117647059=8.5-8=0.5               8
1/0.5=2                               2

The reverse process:
Using integers                             Using floats
1/(2) = 1/2                              | 1/2=0.5
2/(2*8+1) = 2/17                         | 1/(0.5+8)=0.117647059
17/(17*1+2) = 17/19                      | 1/(0.117647059+1)=0.894736842
19/(19*14+17) = 19/283                   | 1/(0.894736842+14)=0.067137809
283/(283*7+19)= 283/2000                 | 1/(0.06713780+7)=0.1415
2000/(2000*3+283) = 2000/6283            | 1/(0.1415+3)=0.318319274

6283/2000 = 3.1415
or 2000/6283 = 0.318319274 then 1/0.318319274 = 3.1415


Mathematically seen large numbers can be rounded just as small numbers, there is no qword limit.

Code Select Expand


Yes, but you may lose accuracy if you don't deal with some problems that will appear, it will depend on how you deal with the problem.
The answer is arithmetic coding.
We initially defined 3 variables: low, medium, high. (our range)
The low will start completely filled with zeros.
The high will start completely filled with the last digit of the base, in this case 9 for the decimal base.
The middle will be (low + ((high - low) / 2)), so we will echo digits in binary to remake the number again.
We write the 3 variables in order and see where our number fits:

low  = 000
mid  = 499
high = 999


step 1:
low  = 000
our number fits here *31415, we are before the middle, we echo 0
mid  = 499
high = 999
the number is before the middle, so the middle becomes high, we calculate new middle, see where our number fits:


step2:
low  = 000
mid  = 249
*31415, we are past the middle, we echo 1
high = 499
half turns low


step3:
low  = 249
*                                       0
mid  = 374
high = 499
half turns high


step 4:
low  = 249
mid  = 311
*                                       1
high = 374
half turns low
In this step a phenomenon happens, our number starts with 3, the middle starts with 3, the top starts with 3. This will never change in the future, hence we know that we have the initial number 3.

From here we have two options to continue:
-------------------------------------
In data compression they remove the number 3 from the high and low part (the bits are the same), insert a 0 in the low part, insert a 9 (infinity) in the high part, and continue the compression:

before
mid  = 311
high = 374

after
low  = 110
high = 749
The process is much more complex than that, just a notion. If you like this way, read about arithmetic coding + data compression.
It's necessary some adjusts (1/4) and (3/4) if numbers are getting so much closer but don't change (like 499 and 500).
One example is: Imagine that the next number is zero, and zero will not fit in that range!!! (will be before low)
The formula that I'm using also changes if dealing with data compression.
https://en.wikipedia.org/wiki/Arithmetic_coding
-------------------------------------
We are not compressing data, we can reset the process and start again from the starting point.
We just need to remove the 3 from our pi and advance to the next digit: 31415-"3" = 1415
step 1:
low = 000
our number fits here *1415, we are before the middle, we echo 0
mid = 499
high = 999
...

Having the bits echoed we can reconstruct the given number:
0101

We start the process by setting low, medium and high. We take a bit, if it's zero then it's between low and middle, if it's 1 then it's between middle and high.
step1: 0101
low = 000
mid = 499
high = 999
The first bit is 0, our number is between low and middle.

step2: 101
low = 000
mid = 249
hih = 499
The current bit is now 1, we are between the middle and high....

The process stops when we have the 3 upper digits equal.

So, you can play with floating number or unsigned integers number.
But same problem happens here:
The number 5 its the question. So, you need deal with (1/2) initially, but take into account (1/4) and (3/4) and do some verify.
Each radix base will have this problem. You need check numbers after our actual number.


--------------------------------
Other choice is, using 31415 as example, but subtracting 1 from high
set a range [00000,100000], divide this range by 10 (decimal base) (from 0 to 9 digits, equal probabilities).

00000 0 09999 1 19999 2 29999 3 39999 4 49999 5 59999 6 69999 7 79999 8 89999 9 999999

Where our number fits? Between [29999 3 39999]
and go on.
39999-29999=10000 / 10 = 1000

29999  30999  31999  32999  33999  34999  35999  36999  37999  38999   39999

Where our number 31415 fits? Between [30999  31999]
Ops, left most numbers are the same, we can throw them away and continue our process.
So:
before
[30999  31999]
after
[09990  19999]  ;insert a zero in low, a 9 (infinite) in high, continue.

...


[JK]

Thanks, i will need some time to digest that ...

[HSE]

Mathematically seen large numbers can be rounded just as small numbers, there is no qword limit.

interesting. FPU limit is 1.8×10¹⁹ (or 2⁶⁴ ). Apparently numbers equals or less than this is are named "integral values".

Then number to round with mul,round,div must be less than 1.8x10^{(19-decimal_positions)}.

What the purpose to round so big numbers?

[Gunther]

HSE,

What the purpose to round so big numbers?

sometimes this is necessary. It depends on the circumstances.

The SI unit of the amount of substance is the mole. It contains 6.02214076 ^. 10²³ elementary entities. That's a factor of 10 000 more than your 19 decimal places - for example.

[HSE]

The SI unit of the amount of substance is the mole. It contains 6.02214076 ^. 10²³ elementary entities.

But perhaps you must use big integers for that, no FPU (like is JK question).

[raymond]

interesting. FPU limit is 1.8×10¹⁹ (or 2⁶⁴ ). Apparently numbers equals or less than this is are named "integral values".

Then number to round with mul,round,div must be less than 1.8x10^{(19-decimal_positions)}.

What the purpose to round so big numbers?

Sorry to disappoint you but the largest number which can be represented with the FPU in REAL10 mode is 1.19x10⁴⁹³², but with only a maximum accuracy of 19 decimal digit or 64 binary bits mind you.

One should ALWAYS round the result of a computation to the accuracy of the least accurate component used for the computation. For example, if 6.02214076 is the least accurate component, ALL results should be rounded to AT MOST 9 total (integral and fractional) decimal digits.

[HSE]

Sorry to disappoint you but the largest number which can be represented with the FPU in REAL10 mode is 1.19x10⁴⁹³², but with only a maximum accuracy of 19 decimal digit or 64 binary bits mind you.

Thanks Raymond, I was not forgetting that. 

Just FRNDINT will fail if number in st(0) is bigger than 2⁶⁴, returning same number in st(0) and C1 will be set to 0.

One should ALWAYS round the result of a computation to the accuracy of the least accurate component used for the computation. For example, if 6.02214076 is the least accurate component, ALL results should be rounded to AT MOST 9 total (integral and fractional) decimal digits.

Good to remember. Perhaps JK idea is related to that.

[JK]

least accurate component

Yes, with scientific calculations you must have an eye on accuracy (aka digits of precision), otherwise you might calculate nonsense. Therefore rounding is a must in certain places. Of course this "integer" limit of FRNDINT comes from the fact that the FPU can hold integers up to 10^18. 64 bits don´t allow for larger integer values. On the other side much larger floating point values are possible (with a maximum accurracy of 18 digits). So how would i round such floating point number to a specific number of decimal digits?

My best bet currently is:

separate the fractional part (which should be possible with FPREM), subtract from original value, multiply fraction, round, divide again and add to the original value.

Rounding to a total number of digits:

save integers as float, load again, divide by ten (so it´s 0.xxx.....) and then do the same as above
added: is normalizing a number necessary at all?

I´m going to investigate this ...   

[HSE]

Sorry to disappoint you but ... but with only a maximum accuracy of 19 decimal digit or 64 binary bits mind you.

Raymond always make the point about things we often forget      .

If you number is bigger than 1.8x10¹⁹ last part of number is lost by FPU, and you have no fractional part, only integers

Congratulation JK, your problem is solved 

[raymond]

If you number is bigger than 1.8x10¹⁹ last part of number is lost by FPU, and you have no fractional part, only integers

... and possibly truncated (specially if you look only at the 64 bits of the significand field).

But, you would rarely want to retain an accuracy even close to 19 decimal digits. And, your rounding procedure should detect such situations where the integer part exceeds the required accuracy and treat it separately anyway.

[hutch--]

Here is the barbarian view from someone who counts with their fingers if it does the job. 

1. Perform the calculation in whatever level of precision you require.
2. Truncate the display result at your predetermined decimal point range.
3. Read the next integer after the decimal point truncation end.
4. Depending on your preferred rounding method, move the last digit in the truncated end either UP or DOWN.
5. Leave the original calculation alone so you can reuse it. 

[jj2007]

Purest Masm64 SDK code attached.

Code Select Expand

src1 REAL10 123.45678901234567890
src2 REAL10 123456789012345.678
...
Result=123.46
Result=123456789012345.69

Same but with Print Str$("Result=%Hf\n", dest), where dest is REAL10: Result=123456789012345.68

Beware of rounding errors ;-)

The FPU has its limits (and SIMD instructions are worse)