为什么这些相同的 QB 计算会产生略有不同的值？

Question

所以，我正在尝试将一些非常古老且古老的工程分析 QBasic 4.5 代码移植到 C 中。我正在尝试精确匹配结果，但我发现我不太明白 QB 是如何进行数学运算的。

比如这两行

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 变为 1.07920125E-4（浮点数 0x38e2532d）

e2 变为 1.0792013E-4（浮点数 0x38e2532e）

这两者略有不同。谁能帮我理解为什么？非常感谢。

Answer 1

对于d2 和e2，我得到了相同的输出，即使就值的原始字节表示而言也是如此。这是一些带注释的输出：

# 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

代码：

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

编辑

正如@BlackJack 在 cmets 中解释的那样，您注意到的效果似乎出现在编译文件时。由于这是给出的线索，我在 DOSBox 中使用了 CodeView 调试器，以下是删减的结果：

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

BASIC 编译器 (BC.EXE) 将 d2 的赋值减少为对浮点常量的简单赋值（即，它评估表达式本身并将您的代码优化为单个赋值，而不是执行所有操作你指定）。但是，对e2 的赋值并不是那么简单，因为它包含一个非常量表达式a!。

为了解决这个问题并尝试保持尽可能高的精度，它将1 / (2 * a * 144) 更改为数学上等效的(1 / 288) / a，并且1 / 288 的近似值存储在偏移量0x0050，这就是为什么@ 987654333@ 最终加载了该偏移量。在加载了SINGLE 值后，它将它除以a 的值（偏移量0x0036）并将结果存储在e2（偏移量0x003E）中。您可以使用CONST a = 32.174 将e2 赋值为与d2 相同，但此时不能更改其值。

现在可能有人想知道为什么这只发生在编译时而不是在 IDE 中，我真的不知道。我最好的猜测是 IDE 在 FP 堆栈上保留尽可能多的浮点数以保持精度，因此它不使用 a 的 32 位舍入值，而是使用存储在 FP 上的现有 80 位值堆栈已经，如果它仍然存储在那里。这样，精度损失就会减少，因为将 80 位值存储在 FP 堆栈之外需要四舍五入到最接近的 32 位或 64 位值，具体取决于您指定存储值的位置。当然，如果出于某种原因，FP 堆栈上需要超过 8 个值，则需要将其中一个换出为另一个腾出空间，最终会出现精度损失。

@BlackJack 还指出 IDE 是在解释代码而不是通过优化对其进行编译，这可能是代码在 IDE 中运行时字节表示相同但在编译版本中不同的原因。也就是说，d2 和 e2 的计算都以完全相同的方式执行，而不是像 BC.EXE 那样将d2 的计算优化为单个值。

无论如何，您可能不会在您的 C 代码中注意到它，因为您的现代编译器比 BC.EXE 更智能，并且在优化方面有更多内存可供使用，即使没有现代浮点的帮助点技术，例如 SSE2。

Answer 2

您使用的是哪个 QB 版本？您如何打印或输出变量d2 和e2？

当我在 DOSBox 0.74 中的 QuickBASIC 4.5 中尝试您的程序时，我没有得到不同的值，当我 PRINT 时，d2 和 e2 是相同的。

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

感叹号运算符会将其类型转换为 SINGLE（单精度，4 个字节），因此它与 AS SINGLE 相同。也许您的d2 = 1! /.. 行中的值32.174 正在以某种方式被类型转换为DOUBLE？