2

So, I'm trying to port some very old and venerable engineering analysis QBasic 4.5 code into C. I'm trying to match results exactly, and I find that I can't quite understand how QB does its math.

For example, these two lines

DIM a AS SINGLE
DIM d2 AS SINGLE
DIM e2 AS SINGLE

a = 32.174
d2 = 1! / (2! * 32.174 * 144!)
e2 = 1! / (2! * a! * 144!)

d2 becomes 1.07920125E-4 (floating point 0x38e2532d)

e2 becomes 1.0792013E-4 (floating point 0x38e2532e)

which are ever so slightly different. Can anyone help me understand why? Thanks very much.

Chris
  • 629
  • 1
  • 10
  • 28

2 Answers2

5

I'm getting the same output for both d2 and e2, even in terms of the raw byte representation of the values. Here's some annotated output:

# Calculation results
d2: 38 E2 53 2E
e2: 38 E2 53 2E
 1.079201E-04 =  1.079201E-04

# Result of changing the last byte (some mantissa bits) to alter the value,
# proving they're not equal
d2: 38 E2 53 2F
e2: 38 E2 53 2E
 1.079201E-04 <> 1.079201E-04

# Result above may just be luck. This result alters the first byte
# (some exponent bits) to prove that the intended bits were altered.
d2: 39 E2 53 2E
e2: 38 E2 53 2E
 4.316805E-04 <> 1.079201E-04

Code:

DIM a AS SINGLE
DIM SHARED d2 AS SINGLE
DIM SHARED e2 AS SINGLE

a = 32.174
d2 = 1! / (2! * 32.174 * 144!)
e2 = 1! / (2! * a! * 144!)

' Print the hex representation of the bytes
' and show they're initially equal.
CALL printHex
PRINT

' Change the last byte of the mantissa by 1 bit.
' Show that doing this makes the two values unequal.
DEF SEG = VARSEG(d2)
    POKE VARPTR(d2), PEEK(VARPTR(d2)) + 1
DEF SEG
CALL printHex
PRINT

' Show that the correct byte was poked by reverting mantissa change and
' altering exponent.
DEF SEG = VARSEG(d2)
    POKE VARPTR(d2), PEEK(VARPTR(d2)) - 1
    POKE VARPTR(d2) + 3, PEEK(VARPTR(d2) + 3) + 1
DEF SEG
CALL printHex

SUB printHex
    'SHARED variables used:
    '  d2, e2

    DIM d2h AS STRING * 8, e2h AS STRING * 8

    ' Get bytes of d2 and e2, storing them as hexadecimal values
    ' in d2h and e2h.
    DEF SEG = VARSEG(d2)
        MID$(d2h, 1) = hexByte$(PEEK(VARPTR(d2) + 3))
        MID$(d2h, 3) = hexByte$(PEEK(VARPTR(d2) + 2))
        MID$(d2h, 5) = hexByte$(PEEK(VARPTR(d2) + 1))
        MID$(d2h, 7) = hexByte$(PEEK(VARPTR(d2)))
    DEF SEG = VARSEG(e2)
        MID$(e2h, 1) = hexByte$(PEEK(VARPTR(e2) + 3))
        MID$(e2h, 3) = hexByte$(PEEK(VARPTR(e2) + 2))
        MID$(e2h, 5) = hexByte$(PEEK(VARPTR(e2) + 1))
        MID$(e2h, 7) = hexByte$(PEEK(VARPTR(e2)))
    DEF SEG

    ' Print the bytes, separating them using spaces.
    PRINT "d2: "; MID$(d2h, 1, 2); " "; MID$(d2h, 3, 2); " ";
    PRINT MID$(d2h, 5, 2); " "; MID$(d2h, 7, 2)
    PRINT "e2: "; MID$(e2h, 1, 2); " "; MID$(e2h, 3, 2); " ";
    PRINT MID$(e2h, 5, 2); " "; MID$(e2h, 7, 2)

    ' Print whether d2 is equal to e2.
    IF d2 = e2 THEN
        PRINT d2; "= "; e2
    ELSE
        PRINT d2; "<>"; e2
    END IF
END SUB

FUNCTION hexByte$ (b%)
    ' Error 5 is "Illegal function call".
    ' This can only happen if b% is outside the range 0..255.
    IF b% < 0 OR b% > 255 THEN ERROR 5

    ' MID$("0" + HEX$(15), 2 + (-1)) => MID$("0F",  1) => "0F"
    ' MID$("0" + HEX$(16), 2 + ( 0)) => MID$("010", 2) => "10"
    hexByte$ = MID$("0" + HEX$(b%), 2 + (b% < 16))
END FUNCTION

EDIT

As @BlackJack explained in the comments, the effects you're noticing appear to occur when the file is compiled. Since that was the clue given, I used the CodeView debugger in DOSBox, and here's the abridged result:

5:      a = 32.174
057D:0030 C70636002DB2   MOV       Word Ptr [0036],B22D
057D:0036 C70638000042   MOV       Word Ptr [0038],4200
6:      d2 = 1! / (2! * 32.174 * 144!)
057D:003C C7063A002D53   MOV       Word Ptr [003A],532D
057D:0042 C7063C00E238   MOV       Word Ptr [003C],38E2
7:      e2 = 1! / (2! * a! * 144!)
057D:0048 CD35065000     FLD       DWord Ptr [0050]; 00 CB 21 CD
057D:004D CD34363600     FDIV      DWord Ptr [0036]; 42 00 B2 2D
057D:0052 CD351E3E00     FSTP      DWord Ptr [003E]; e2 = result
057D:0057 CD3D           FWAIT

The BASIC Compiler (BC.EXE) reduced the assignment to d2 to a simple assignment of a floating-point constant (i.e. it evaluated the expression itself and optimized your code to that single assignment rather than performing all of the operations you specified). However, the assignment to e2 isn't that simple since it contains a non-constant expression a!.

To combat the problem and attempt to keep as much precision as possible, it changed 1 / (2 * a * 144) to the mathematically equivalent (1 / 288) / a, and the approximate value of 1 / 288 was stored at offset 0x0050, which is why FLD ended up loading that offset. After loading that SINGLE value, it divided it by the value of a (offset 0x0036) and stored the result in e2 (offset 0x003E). You can make the assignment to e2 the same as d2 by using CONST a = 32.174, but you cannot change its value then.

Someone is probably wondering by now why this only happens when compiled and not in the IDE, and I honestly don't know. My best guess is the IDE keeps as many floats on the FP stack as much as it can to retain precision, so instead of using the 32-bit rounded value of a, it uses the existing 80-bit value stored on the FP stack already, if it is still stored there. This way, there's less precision loss since storing the 80-bit value outside of the FP stack requires rounding to the nearest 32-bit or 64-bit value, depending on where you specify to store the value. Of course, if more than 8 values are needed on the FP stack for some reason, one will need to be swapped out to make room for another, and the precision loss will manifest eventually.

@BlackJack also pointed out that the IDE is interpreting the code rather than compiling it with optimizations, which could be the reason that the byte representations are the same when the code is run in the IDE but different in a compiled version. That is, both the computations of d2 and e2 are executed in the exact same way rather than optimizing the calculation of d2 to a single value as BC.EXE does.

In any case, you likely don't notice it in your C code because your modern compiler is a lot smarter and has more memory to work with when it comes to optimizing than BC.EXE, even without the help of modern floating-point technologies like SSE2.

  • 2
    Nice, that is some thorough investigation and explanation. My only comment is that the repeating `MID$+PRINT` part could be put in a `SUB`. – BdR Jun 29 '16 at 07:43
  • Wow, that is nice work. You are obviously one of the world's remaining experts in QB. :) I guess the takeaway is that I am experiencing the problem and you are not... – Chris Jun 29 '16 at 15:15
  • @Chris Not really an expert. Just converting some knowledge of C and DOS to QB. –  Jun 29 '16 at 19:26
  • @BdR Agreed. I've updated the code, moving those common bits of code into a `SUB` and changing the `DEF FNhex$` definition to `FUNCTION hexByte$` since `DEF FN` and `SUB` don't get along. –  Jun 29 '16 at 22:19
  • 1
    I get different output depending on wether I run the program from within the QB4.5 IDE or the compiled EXE. The standalone executable outputs what the OP showed in the question (2D and 2E as last byte). – BlackJack Jan 31 '18 at 21:29
  • @BlackJack I never thought about that honestly; I assumed the IDE and the compiler behaved the same, but apparently that isn't true! Anyway, I just edited my answer with some info about what I found. Thanks for the help. :-) –  Feb 01 '18 at 04:04
  • I'm just guessing too: Most probably running in the IDE doesn't compile but interpret what's there in the code without optimizing/constant folding. – BlackJack Feb 01 '18 at 08:20
  • @BlackJack Makes perfect sense. I always miss the simple stuff when I'm working with things that seem more complex to me. :-P –  Feb 02 '18 at 00:05
2

Which QB version are you using? And how are you printing or outputing your variables d2 and e2?

When I try your program in QuickBASIC 4.5 in DOSBox 0.74 I don't get different values, d2 and e2 are the same when I PRINT them.

a =  32.174
d2 =  1.079201E-04 
e2 =  1.079201E-04 

The exlamation mark operator will typecast it to SINGLE (single precision, 4 bytes), so it's same as AS SINGLE. Maybe the value 32.174 in your line d2 = 1! /.. is being typecast to DOUBLE somehow?

BdR
  • 2,770
  • 2
  • 17
  • 36
  • Looks like when you PRINT them there is rounding going on... the difference is past that "201". I am comparing them by actually using "PUT#" to write them to file (I was doing that anyway), and then looking at the hex. Also, when I put the suffix '!' on that 32.174, The IDE deletes it. From which I assume that the SINGLE is implicit and therefore the '!' is superfluous. – Chris Jun 28 '16 at 21:12
  • Oh sorry... using QB4.5 on a 32-bit machine. I have DOSBOX on 64-bit but haven't tried that. – Chris Jun 28 '16 at 21:13
  • I'm not familiar with `PUT#` to output and examine SINGLE variables. But anyway my QuickBasic IDE also removes the `!` automatically after editing the line, so the variable is already implicitly SINGLE like you say. Btw when you add a `#` behind it to indicate it's a DOUBLE that is not removed. – BdR Jun 28 '16 at 21:34