FPU compare real to integer problem

[jj2007]

SmplMath and MasmBasic use two different methods to accommodate rounding errors:

Code Select Expand

.if faEQ(100.01,100.02,1.0E-3) ; +- 1.0E-3
print "is equal"
.endif
Fcmp FP8(100.01), FP8(100.02), 6 ; ca six digits
.if Zero?
print "is equal"
.endif


In other words, SmplMath sets an absolute error range, while MasmBasic uses a relative precision. What is "better" depends on your app.
In both cases, the point here is that the FPU, when reporting "equal", means "equal relative to my internal 64-bit precision".

The problem is you won't notice the difference. Or, more precisely: most of the time, you won't notice... :bg

test4=12345678   ; try taking away the 8...

include \masm32\MasmBasic\MasmBasic.inc   ; download
include \macros\smplmath\math.inc

.data?
MyReal10   REAL10 ?

.code
main proc
LOCAL fSlvTLS()   ; otherwise smplmath says it's not thread-safe
   
   print "exact compare: ",13,10
   .if fEQ(1.01,1.02)
      print "    fEQ( 1.01, 1.02):           equal      (exact)",13,10
   .else
      print "    fEQ( 1.01, 1.02):           not equal  (exact)",13,10
   .endif

   print chr$(13,10)
   print "compare using a threshold factor: ",13,10

   .if faEQ(1.01,1.02)
      print "    faEQ( 1.01, 1.02, 1.0E-2):  equal      (+-1%) ",13,10
   .else
      print "    faEQ( 1.01, 1.02, 1.0E-2):  not equal  (+-1%)",13,10
   .endif

   .if faEQ(1.01,1.02,1.0E-3)
      print "    faEQ( 1.01, 1.02, 1.0E-3):  equal      (+-1 per mille) ",13,10
   .else
      print "    faEQ( 1.01, 1.02, 1.0E-3):  not equal  (+-1 per mille)",13,10
   .endif   
   
   .if faEQ(100.01,100.02,1.0E-3)
      print "    faEQ( 100.01, 100.02, 1.0E-3):  equal      (+-1 per mille) ",13,10
   .else
      print "    faEQ( 100.01, 100.02, 1.0E-3):  not equal  (+-1 per mille)",13,10
   .endif   

   print chr$(13,10), "Same with MasmBasic Fcmp:", 13, 10
   Fcmp FP8(100.01), FP8(100.02), high   ; ca REAL8
   .if Zero?
      print "    Fcmp 100.01, 100.02, high:  equal       ",13,10
   .else
      print "    Fcmp 100.01, 100.02, high:  not equal  ",13,10
   .endif
   Fcmp FP8(100.01), FP8(100.02), 5   ; ca. five digits
   .if Zero?
      print "    Fcmp 100.01, 100.02, 5:     equal       ",13,10
   .else
      print "    Fcmp 100.01, 100.02, 5:     not equal  ",13,10
   .endif   
   Fcmp FP8(100.01), FP8(100.02), 6   ; ca six digits
   .if Zero?
      print "    Fcmp 100.01, 100.02, 6:  equal       ",13,10
   .else
      print "    Fcmp 100.01, 100.02, 6:     not equal  ",13,10
   .endif   
   Fcmp FP8(100.01), FP8(100.02), low   ; ca REAL4
   .if Zero?
      print "    Fcmp 100.01, 100.02, low:   equal       ",13,10, 10
   .else
      print "    Fcmp 100.01, 100.02, low:  not equal  ",13,10, 10
   .endif
   test4=12345678   ; try taking away the 8...

   print "little test x=sqrt(x/7)^2*28/x (should be exactly 4):", 13, 10, "for x="
   print str$(test4), ":", 13, 10
   void Str$(0)   ; initialise some MB routines
   push test4
   push 7
   push 28
   push test4
   fild stack
   pop eax
   fild stack
   pop eax
   fild stack
   pop eax
   fild stack
   pop eax
   deb 4, "FPU", ST(0), ST(1), ST(2), ST(3)   ;
   fdivr   ; 4/7=0.57
   fsqrt   ; 0.7
   fld st
   fmul   ; 4/7
   fmul   ; 4/7*28=16
   fdivr   ; 16/4=4
   fstp MyReal10
   deb 4, "R10", MyReal10
   Fcmp MyReal10, 4, top
   .if Zero?
      print "4 is 4", 13, 10
   .else
      print "4 is not 4", 13, 10
   .endif
   inkey chr$(13, 10, 10)
   exit
   
main endp
end main

Output:

Code Select Expand

exact compare:
    fEQ( 1.01, 1.02):           not equal  (exact)

compare using a threshold factor:
    faEQ( 1.01, 1.02, 1.0E-2):  equal      (+-1%)
    faEQ( 1.01, 1.02, 1.0E-3):  not equal  (+-1 per mille)
    faEQ( 100.01, 100.02, 1.0E-3):  equal      (+-1 per mille)

Same with MasmBasic Fcmp:
    Fcmp 100.01, 100.02, high:  not equal
    Fcmp 100.01, 100.02, 5:     equal
    Fcmp 100.01, 100.02, 6:     not equal
    Fcmp 100.01, 100.02, low:   equal

little test x=sqrt(x/7)^2*28/x (should be exactly 4):
for x=12345678:

FPU
ST(0)           12345678.0000000000
ST(1)           7.00000000000000000
ST(2)           28.0000000000000000
ST(3)           12345678.0000000000
R10     MyReal10        4.00000000000000000
4 is 4

[RuiLoureiro]

Thanks for your replies
I need to give answers to all of you. Wait please

Dave,
1) use FINIT once at the beginning of the program - it sets up the FPU for 80-bit precision
when you start a program under windows, it is set up for 64-bit precision

2) use REAL10 instead of DT or TBYTE PTR

3) this code leaves one of the FPU registers loaded - make sure you FFREE it

        1. The first instruction is exactly FINIT
        3. Every succesful calculation starts loading the FPU
           with the first operand, next the second operand,
           next the operation and end with store the result at a
           REAL10 variable.
           Example 1:

        X   dt 2
        Y   dt 5
        Z   dt ?
                        mov     eax, offset X
                        mov     ebx, offset Y
                        ;
                        fld     tbyte ptr [eax]     ; st(1)
                        fld     tbyte ptr [ebx]     ; st(0)
                        fadd                        ; new st(0) = st(1)+st(0)
                        mov     eax, offset Z                             
                        fstp    tbyte ptr [eax]     ; store st(0) and pop

            What to FFREE ?

           Example 2:

        X   dt 2
        Y   dt 5
        Z   dt ?
                        mov     eax, offset X
                        ;
                        fldpi                       ; st(2)
                        fld     st(0)               ; st(1) get a copy                       
                       
                        fld     tbyte ptr [eax]     ; st(0)
                        fadd                        ; new st(0) = st(1)+st(0)
                        mov     eax, offset Z                             
                        fstp    tbyte ptr [eax]     ; store st(0) and pop
                       
                        fstp    st                  ; remove source


            What to FFREE ?
            I think all this examples are correct and we dont need to free nothing,
            we dont need to do finit againg. Correct Dave ?

[RuiLoureiro]

qWord,

        This precision problem can be solved by using a threshold-factor:
        pseudo code:

            if abs(a-b) <= abs(a*factor) then equal(requires: a != 0)

    0. Thanks for your reply;
   
    1. I dont know exactly what to do; What factor to use

    2. I think you understood the problem;
       I repeat it (this is exactly one example)

        a) From keyboard we type a string "2^(2^2)";
        b) The first thing to do is to solve (2^2);
        c) We convert the first 2 to a real10 variable X
           and the second 2 to a real10 variable Y;
        d) We apply the first formula to get X^Y in a real10 variable Z;

            So far so good!

            Now i want to use the second formula X^Zi where Zi is an integer
            ( X is real10 and Zi dword )
            or the first formula X^Z ( X is real10 and Z is real10 )

            Well, if Z=X^Y is an integer than Z=Trunc(Z) or Z-Trunc(Z)=0
                  if  Z=X^Y is not an integer than Z != Trunc(Z) or Z-Trunc(Z) != 0

            Correct ?

            I write 2 or 4 if integer and 2.0 or 4.0 if real, ok ?
           
            Solving this by hand, we have 2.0^2.0 is a real 4.0.

            So Z=4.0 and Trunc(4.0) = 4 than «Z is equal to Trunc(Z)»

            FPU tell me «Z is lower than Trunc(Z)»

            I dont know what is the result of ABS(Z-Trunc(Z)) using the FPU till now

[RuiLoureiro]

Raymond,

        Thank you for your replies

        1. I know your fpulib: FpuAtoFL, FpuTrunc, FpuXexpY, FpuCos, etc. etc.
           I read all what is written there.

        2. I think that FpuXexpY solves X^Y if and only if X>0
           it doesn't matter at all whether the "Y" is a float or an integer.

           Am i correct ?           
       
        3. The problem is when X is negative. How do you solve it in that case ?
           This answer your question, i think

[RuiLoureiro]

jj,
        Thanks for your reply

        The problem is not with fcomip. If you read my answers maybe you
        understand the problem

vanjast,
          Thank you for your help

[RuiLoureiro]

raymond,
        Please, read the answer i gave to qWord.

        Yes, integer32 is a dword.
       
        I used the code you post and i got Z is not_a_dword_integer
        It is 4.0... but not a_dword_integer

        It doesnt solve my problem.
        In your last post it seems you have only one formula to solve
        the problem of (X)^(Y) where X and Y are integers or reals.
        I am waiting your answer.
        Thanks

[dedndave]

            What to FFREE ?

the FPU has 8 registers
they behave differently from CPU general registers in many ways
one important difference is that each register is always marked as empty or full
if all 8 registers are full, an error (exception) will occur if you attempt to load another value

so - if you load (or push) 2 values, you should pop 2 to keep balance
this allows you to use the function repeatedly without filling all the registers

in the code i posted

Code Select Expand

        finit

        fld real10 ptr Y
        fild dword ptr Yi
        fcomip   st,st(1)
        jz       equal

FINIT empties all 8 registers
we then load 2 values with FLD and FILD
the FCOMIP instruction pops one of those values
but, we have left the other value in the FPU
if we want to use the code over and over without FINIT, we should FFREE st(0)
FINIT is a nice instruction   :P
but, it takes a lot of clock cycles

[jj2007]

        The problem is not with fcomip.

if  Z=X^Y is not an integer than Z != Trunc(Z)

The FPU will tell you Z != Trunc(Z) because of its precision limits. Your examples are confused, it would really help if you could isolate the problem and post the code. Can't be so difficult for a console app...

[qWord]

In other words, SmplMath sets an absolute error range, while MasmBasic uses a relative precision.

no, SmplMath error range is also relative: e.g. if abs(a-b) <= abs(a*1.E-2) ...
This method fails, if a=0.
Interesting to see, that MB is analysing the fraction bits - vieleicht klau ich die Idee :bg

[jj2007]

Interesting to see, that MB is analysing the fraction bits - vieleicht klau ich die Idee :bg

Finger weg! I just downloaded SmplMath, looks good :U I am in an advanced stage of rolling my own, so I will deliberately avoid looking at your sources. My tokeniser is working perfectly, but i got a little bit stuck with the preference order, and would need a) a full night of sleep, b) a full day of quiet coding, both of which won't happen anytime soon :(
Anyways, compliments for SmplMath :thumbu

[qWord]

I am in an advanced stage of rolling my own, so I will deliberately avoid looking at your sources. My tokeniser is working perfectly, but i got a little bit stuck with the preference order, and would need a) a full night of sleep, b) a full day of quiet coding, both of which won't happen anytime soon :(
Anyways, compliments for SmplMath :thumbu

It will work out in the end - I needed 3 attempts and 5 years to get it.

:bg

[raymond]

3. The problem is when X is negative. How do you solve it in that case ?

With X being negative, the result is a real number when Y is an integer. However, if Y is not an integer, the result is a complex number.
Can you handle complex numbers?
If yes, you could use the Zlib (similar to the Fpulib but for complex numbers) from http://www.ray.masmcode.com/complex.html
It would return the result with the correct sign for all cases, whether Y is integer or not.

Note: the downloadable file contains a dialog box which you can use to test any of the functions with any input parameters you want.

[raymond]

Here are a few examples of the previous post.

(-3.0)² = 9 + 0(i)

(-3.0)^2.2 =  9.070357505108625 + 6.590000471674018(i)
(3.0)^2.2 =  11.21157845653966 + 0(i) = sqrt(9.070357505108625² + 6.590000471674018(i)²)

(-3.0)^2.5 = 0 + 15.58845726811990(i)
(3.0)^2.5 = 15.58845726811990 + 0(i)

(-3.0)^2.8 = -17.53465226998805 + 12.73967058793206(i)
(3.0)^2.8 = 21.67402216752623 + 0(i) = sqrt((-17.53465226998805)² + 12.73967058793206(i)²)

(-3.0)^-2.8 = -0.037326574095097 + -0.027119343504831(i)
(3.0)^-2.8 = 0.046138182948723 + 0(i) = sqrt((-0.037326574095097)² + (-0.027119343504831(i))²)

[RuiLoureiro]

raymond,

Can you handle complex numbers?

        Yes, but in my present work i dont want to use it.
        In any case, thank you for your help.

With X being negative, the result is a real number when Y is an integer

        This is why we need 2 formulas:

                    1)  (X)^(Y)     X, Y real10
                    2)  (X)^(Yi)    X real10,  Yi integer32           

        In general, X and Y are the results of expressions.
        Because i am working with real10 numbers, i have a calculator
        (a proc that does this) that gives me the result in real10.
               
        For that reason, i want to get Yi from Y, making Yi=Trunc(Y)
        and not from a particular procedure that calculate Yi directly.
        The problem is to decide if Yi is the correct integer and if so
        we use formula 2).
        If Y is not an integer we use formula 1).
       
        Have you any opinion about this ?

About FpuSin

        Could you tell me why it gives Sin( pi ) = -0.00000... ?
        (But FpuCos doesnt give it with pi/2)

        And why not 0.00000... ?

        Thanks

[raymond]

Quote
Can you handle complex numbers?


        Yes, but in my present work i dont want to use it.

If Y is not an integer we use formula 1).

The only conclusion I can reach from the above quotes is that "formula 1" uses the absolute value of X and returns the length of the vector (without the angle) in the complex plane when X is negative. Otherwise, you don't have any choice but to use complex numbers to compute such powers and return both cartesian (or polar) coordinates in the complex plane.