Floating Point Program gives Invalid Result

Question

RECENT EDIT

I am trying to run this floating point Quadratic Equation program on x86 MASM. This code is found in the Kip Irvine x86 textbook and I want to see how it works visually. The following code is below:

include irvine32.inc 
.DATA
a REAL4 3.0
b REAL4 7.0
cc REAL4 2.0
posx REAL4 0.0
negx REAL4 0.0

.CODE


main proc 
; Solve quadratic equation - no error checking
; The formula is: -b +/- squareroot(b2 - 4ac) / (2a)
fld1 ; Get constants 2 and 4
fadd st,st ; 2 at bottom
fld st ; Copy it
fmul a ; = 2a

fmul st(1),st ; = 4a
fxch ; Exchange
fmul cc ; = 4ac

fld b ; Load b
fmul st,st ; = b2
fsubr ; = b2 - 4ac
; Negative value here produces error
fsqrt ; = square root(b2 - 4ac)
fld b ; Load b
fchs ; Make it negative
fxch ; Exchange

fld st ; Copy square root
fadd st,st(2) ; Plus version = -b + root(b2 - 4ac)
fxch ; Exchange
fsubp st(2),st ; Minus version = -b - root(b2 - 4ac)

fdiv st,st(2) ; Divide plus version
fstp posx ; Store it
fdivr ; Divide minus version
fstp negx ; Store it

call writeint
exit 
main endp 
end main

So I was able to have my program compile, execute and work completely. However, Whenever I run the program, I get this as the result:

[phone]

Why is the result so large? Also I tried calling writefloat but it says this procedure is not in the Irvine32.inc or Macros.inc Library. Can someone show me why is it not working and what needs to be fixed? Thanks.

Dorran

You don't end the program in any way so the CPU will continue reading garbage data and executing it. Use a debugger to see what's happening always. – Aliphatic

I do not see the C part of this question. – Velarde

By the time you get down to fstp negx, st (top of FPU stack) has one of the roots of the quadratic equation. fstp negx takes the value in st and places it in memory at negx and pops the stack. You can simply remove fstp negx and print the top of the FPU stack by simply using the WriteFloat function with call WriteFloat.If you want the value both stored in negx and printed then you can change fstp negx to fst negx and follow it with call WriteFloat – Tattan

Of course remove call writeint altogether since it prints out the signed integer in EAX. It doesn't write floating point values. – Tattan

WriteFloat must be there except you've got a real old irvine32.inc and irvine32.lib. Get a newer link library from Irvine's homepage (Example programs and link library source code for the Seventh Edition) and install it. – Exert

@Exert But perhaps Gordon has a 5th edition textbook, or got his paths wrong? – Granulocyte

@IwillnotexistIdonotexist: The newer files are complete and compatible. I see no reason to stay behind. Wrong paths would cause the same error at call writeint. – Exert

Answer 1

Floating point numbers are processed in special registers in a special processor (FPU) and stored in a special format and cannot be treated as integers (WriteInt). Floating point numbers contain further information about the number like sign and exponent. The number itself is changed to a number between 1 and 2 with an appropriate exponent, where the leading 1 is normally hidden. Take a look here just for the double format: https://en.wikipedia.org/wiki/Double-precision_floating-point_format. These numbers are unlikely to be precise.

At least since 11 years the Irvine32 library provides the function WriteFloat to show the value of the FPU register ST0 in exponential form. It does not pop or free that register.

Change

call writeint

to

fld posx                ; Load floating point number into ST0
call WriteFloat         ; Write ST0
ffree st[0]             ; Free ST0 - don't forget it!
call Crlf               ; New line
fld negx                ; Load floating point number into ST0
call WriteFloat         ; Write ST0
ffree st[0]             ; Free ST0 - don't forget it!
call Crlf               ; New line

If your library don't have WriteFloat I suggest to download and install the newest files from Irvine's homepage: http://www.kipirvine.com/asm/examples/index.htm (Example programs and link library source code for the Seventh Edition). You can also use another library, e.g. the C-runtime library (msvcrt.inc and msvcrt.lib) or Raymond Filiatreault's FPU library.

If you can't use a library that provides floating point routines you have to convert the numbers by yourself:

INCLUDE irvine32.inc

.DATA
    a REAL4 3.0
    b REAL4 7.0
    cc REAL4 2.0
    posx REAL4 0.0
    negx REAL4 0.0

    buf BYTE 1024 DUP (?)

.CODE

double2dec PROC C USES edi              ; Args: ST(0): FPU-register to convert, EDI: pointer to string
LOCAL CONTROL_WORD:WORD, TEN:WORD, TEMP:WORD, DUMMY:QWORD

    ; modifying rounding mode
    fstcw CONTROL_WORD
    mov ax, CONTROL_WORD
    or ah, 00001100b            ; Set RC=11: truncating rounding mode
    mov TEMP, ax
    fldcw TEMP                  ; Load new rounding mode

    ; Check for negative
    ftst                        ; ST0 negative?
    fstsw ax
    test ah, 001b
    jz @F                       ; No: skip the next instructions
    mov byte ptr [edi], '-'     ; Negative sign
    add edi, 1
    @@:
    FABS                        ; Abs (upper case to differ from C-library)

    ; Separate integer and fractional part & convert integer part into ASCII
    fst st(1)                   ; Doubling ST(0) - ST(1)=ST(0)
    frndint                     ; ST(0) to integer
    fsub st(1), st(0)           ; Integral part in ST(0), fractional part in ST(1)

    ; Move 10 to st(1)
    mov TEN, 10
    fild TEN
    fxch

    xor ecx, ecx                ; Push/pop counter

    @@:                         ; First loop
    fst st(3)                   ; Preserve ST(0)
    fprem                       ; ST(0) = remainder ST(0)/ST(1)
    fistp word ptr TEMP         ; ST(3) -> ST(2) !
    push word ptr TEMP
    inc ecx
    fld st(2)                   ; Restore ST(0)
    fdiv st(0), st(1)
    frndint                     ; ST(0) to integer
    fxam                        ; ST0 == 0.0?
    fstsw ax
    sahf
    jnz @B                      ; No: loop

    fxch st(2)                  ; ST(0) <-> ST(2) (fractional part)
    ffree st(2)
    ffree st(3)

    @@:                         ; Second loop
    pop ax
    or al, '0'
    mov [edi], al
    inc edi
    loop @B                     ; Loop ECX times

    mov byte ptr [edi], '.'     ; Decimal point
    add edi, 1

    ; Isolate digits of fractional part and store ASCII
    get_fractional:
    fmul st(0), st(1)           ; Multiply by 10 (shift one decimal digit into integer part)
    fist word ptr TEMP          ; Store digit
    fisub word ptr TEMP         ; Clear integer part
    mov al, byte ptr TEMP       ; Load digit
    or al, 30h                  ; Convert digit to ASCII
    mov byte ptr [edi], al      ; Append it to string
    add edi, 1                  ; Increment pointer to string
    fxam                        ; ST0 == 0.0?
    fstsw ax
    sahf
    jnz get_fractional          ; No: once more
    mov byte ptr [edi], 0       ; Null-termination for ASCIIZ

    ; clean up FPU
    ffree st(0)                 ; Empty ST(0)
    ffree st(1)                 ; Empty ST(1)
    fldcw CONTROL_WORD          ; Restore old rounding mode

    ret                         ; Return: EDI points to the null-terminated string
double2dec ENDP


main proc
    ; Solve quadratic equation - no error checking
    ; The formula is: -b +/- squareroot(b2 - 4ac) / (2a)
    fld1 ; Get constants 2 and 4
    fadd st,st ; 2 at bottom
    fld st ; Copy it
    fmul a ; = 2a

    fmul st(1),st ; = 4a
    fxch ; Exchange
    fmul cc ; = 4ac

    fld b ; Load b
    fmul st,st ; = b2
    fsubr ; = b2 - 4ac
    ; Negative value here produces error
    fsqrt ; = square root(b2 - 4ac)
    fld b ; Load b
    fchs ; Make it negative
    fxch ; Exchange

    fld st ; Copy square root
    fadd st,st(2) ; Plus version = -b + root(b2 - 4ac)
    fxch ; Exchange
    fsubp st(2),st ; Minus version = -b - root(b2 - 4ac)

    fdiv st,st(2) ; Divide plus version
    fstp posx ; Store it
    fdivr ; Divide minus version
    fstp negx ; Store it

    ; Write the results

    fld posx            ; Load floating point number into ST0
    lea edi, buf        ; EDI: pointer to a buffer for a string
    call double2dec     ; Convert ST0 to buf and pop
    mov edx, edi        ; EDX: pointer to a null-terminated string
    call WriteString    ; Irvine32

    call Crlf           ; Irvine32: New line

    fld negx            ; Load floating point number into ST0
    lea edi, buf        ; EDI: pointer to a buffer for a string
    call double2dec     ; Convert ST0 to buf and pop
    mov edx, edi        ; EDX: pointer to a null-terminated string
    call WriteString    ; Irvine32

    call Crlf           ; Irvine32: New line

    exit                ; Irvine32: ExitProcess
main ENDP
end main

Exert

Answer 2

With thanks to Michael Petch

Your bug lies not in the computation itself, which in MASM is correct, but in your printing of the result. For printing floating-point numbers, writeint is incorrect; You should use WriteFloat, which has its own calling convention to be followed.

WriteFloat accepts a single float in st(0) and prints it to console [1]. It does not pop the value from the x87 stack.

Therefore, immediately after your FPU code, you should add

fld  posx
call WriteFloat
call Crlf
fld  negx
call WriteFloat
call Crlf
emms

You should also have the correct MASM includes at the beginning. Something akin to:

INCLUDE     irvine32.inc
INCLUDE     floatio.inc
INCLUDE     macros.inc
INCLUDELIB  kernel32.lib
INCLUDELIB  user32.lib
INCLUDELIB  Irvine32.lib

Lacking MASM on my machine, I rewrote your program as a GNU C program with inline assembly, and besides each instruction put in a comment the state of the floating-point stack at that moment.

#include <stdio.h>


int main(void){
    asm(
    ".intel_syntax\n"
    ".data\n"
    "a:     .single 3.0\n"
    "b:     .single 7.0\n"
    "cc:    .single 2.0\n"
    "posx:  .single 0.0\n"
    "negx:  .single 0.0\n"

    ".text\n"
    "fld1\n"                     // [1]
    "fadd    %st, %st\n"         // [2]
    "fld     %st\n"              // [2,                     2]
    "fmul    dword ptr a\n"      // [2a,                    2]

    "fmul    %st(1), %st\n"      // [2a,                    4a]
    "fxch\n"                     // [4a,                    2a]
    "fmul    dword ptr cc\n"     // [4ac,                   2a]

    "fld     dword ptr b\n"      // [b,                     4ac,              2a]
    "fmul    %st, %st\n"         // [b^2,                   4ac,              2a]
    "fsubrp\n"                   // [b^2-4ac,               2a]
    "fsqrt\n"                    // [sqrt(b^2-4ac),         2a]
    "fld     dword ptr b\n"      // [b,                     sqrt(b^2-4ac),    2a]
    "fchs\n"                     // [-b,                    sqrt(b^2-4ac),    2a]
    "fxch\n"                     // [sqrt(b^2-4ac),         -b,               2a]

    "fld     %st\n"              // [sqrt(b^2-4ac),            sqrt(b^2-4ac), -b, 2a]
    "fadd    %st, %st(2)\n"      // [-b+sqrt(b^2-4ac),         sqrt(b^2-4ac), -b, 2a]
    "fxch\n"                     // [   sqrt(b^2-4ac),      -b+sqrt(b^2-4ac), -b, 2a]
    "fsubp   %st(2), %st\n"      // [-b+sqrt(b^2-4ac),      -b-sqrt(b^2-4ac), 2a]

    "fdiv    %st, %st(2)\n"      // [(-b+sqrt(b^2-4ac))/2a, -b-sqrt(b^2-4ac), 2a]
    "fstp    dword ptr posx\n"   // [ -b-sqrt(b^2-4ac),     2a]
    "fdivrp\n"                   // [(-b-sqrt(b^2-4ac))/2a]
    "fstp    dword ptr negx\n"   // []

    ".att_syntax\n"
    );

    extern float posx, negx;

    printf("posx: %+0.17f\nnegx: %+0.17f\n", posx, negx);

    return 0;
}

Printing

posx: -0.33333334326744080
negx: -2.00000000000000000

Which is correct:

• 3*(-1/3)^2 + 7*(-1/3) + 2 = 3/9 - 7/3 + 2 = 1/3-7/3+2 = -6/3+2 = -2+2 = 0
• 3*(-2)^2 + 7*(-2) + 2 = 3*4 -14 + 2 = 12-14+2 = -2+2 = 0

[1] § 12.2.7 Reading and Writing Floating-Point Values

Granulocyte

You are using a different assembler than the OP. You are using GCC inline assembler using Intel syntax. The OPs code does in fact work (except for the writeint and is proper syntax. It is using 32-bit MASM (Microsoft Assembler) with the Irvine32 library (it is a library from a well known Assembly language book).MASM unlike other assemblers has limited type checking. The size of each variable is known from the declarations which all use REAL4 which is 32-bits (size of a DWORD). MASM will automatically determine (where it can) the size of memory operands from the declarations – Tattan

Upon looking at the way the OPs code is structured it does appear the FPU stack is returned to a clear state after the final fstp negx. That is because MASM will encode fsubr to fsubrp st(1),st and fdivr to fdivrp st(1),st – Tattan

@MichaelPetch OK, I see. But in that case my own program is a faithful GAS reproduction of the MASM, and my commenting of activity on the x87 stack is exactly what OP thought he should see. There is no bug in the FPU code, only the display code. – Granulocyte

As I said in the code actually does the computation right. But the OP is also asking how to write the float to the display. He originally tagged this question MASM not GNU Assembler/GCC. He is using the well known Irvine32 library(someone kindly came in and tagged it later). His failure is pretty much about the incorrect use of WriteInt instead of WriteFloat (He just has to put WriteFloat twice in his code to print both roots). His undefined error is likely because he hasn't linked the Irvine32.lib library in. – Tattan

So although you provide an answer it really doesn't answer the question as posed by the OP. The OPs code was already minimal complete example that was syntactically correct and usable in MASM. – Tattan

@MichaelPetch I made the changes, with credit to you. – Granulocyte