📘 第23课：双精度FPU（FP64）

• 掌握双精度浮点格式（binary64）的完整结构
• 理解双精度运算相比单精度的硬件挑战
• 实现双精度浮点加法器与乘法器
• 实现双精度FMA模块
• 分析双精度FPU的面积/时序权衡

双精度（binary64/FP64）是科学计算的黄金标准：

从FP32到FP64，硬件复杂度不是简单翻倍：

1. 乘法器

2. 加法器对齐

双精度FMA的对齐宽度约210位，加法器的CPA延迟显著增加：

210位CPA比106位CPA多1级门延迟。

3. 前导零计数

210位LZC比106位LZC多1级递归，延迟增加约1个MUX2级。

//=============================================================
// fp64_adder.sv - 双精度浮点加法器
//=============================================================
module fp64_adder (
    input  wire [63:0] a,
    input  wire [63:0] b,
    output wire [63:0] result
);
    wire sign_a = a[63], sign_b = b[63];
    wire [10:0] exp_a = a[62:52], exp_b = b[62:52];
    wire [51:0] frac_a = a[51:0], frac_b = b[51:0];

    // 特殊值检测
    wire a_nan  = (exp_a==11'h7FF)&(frac_a!=52'b0);
    wire b_nan  = (exp_b==11'h7FF)&(frac_b!=52'b0);
    wire a_inf  = (exp_a==11'h7FF)&(frac_a==52'b0);
    wire b_inf  = (exp_b==11'h7FF)&(frac_b==52'b0);
    wire a_zero = (exp_a==11'b0)&(frac_a==52'b0);
    wire b_zero = (exp_b==11'b0)&(frac_b==52'b0);

    // 尾数(含隐含1)
    wire [52:0] ma = (exp_a==11'b0) ? {1'b0,frac_a} : {1'b1,frac_a};
    wire [52:0] mb = (exp_b==11'b0) ? {1'b0,frac_b} : {1'b1,frac_b};

    // 指数比较
    wire a_larger = (exp_a > exp_b) | ((exp_a==exp_b) & (ma >= mb));
    wire [11:0] exp_diff = a_larger ? ({1'b0,exp_a} - {1'b0,exp_b}) :
                                     ({1'b0,exp_b} - {1'b0,exp_a});
    wire [10:0] max_exp = a_larger ? exp_a : exp_b;

    // 移位对齐
    wire [52:0] larger_m  = a_larger ? ma : mb;
    wire [52:0] smaller_m = a_larger ? mb : ma;
    wire [52:0] shifted_m = (exp_diff < 12'd53) ? (smaller_m >> exp_diff) : 53'b0;

    wire larger_sign  = a_larger ? sign_a : sign_b;
    wire smaller_sign = a_larger ? sign_b : sign_a;

    // 加减法
    wire same_sign = (larger_sign == smaller_sign);
    wire [53:0] add_out = same_sign ?
        ({1'b0,larger_m} + {1'b0,shifted_m}) :
        ({1'b0,larger_m} - {1'b0,shifted_m});

    wire r_sign = same_sign ? larger_sign :
        (add_out[53]) ? ~larger_sign : larger_sign;
    wire [53:0] abs_out = add_out[53] ? (~add_out+54'b1) : add_out;

    // 规格化(简化：未含完整前导零检测)
    wire [53:0] norm_out = abs_out;
    wire [51:0] r_frac = norm_out[51:0];
    wire [10:0] r_exp = max_exp;

    // 特殊值
    wire inf_inf_diff = a_inf & b_inf & (sign_a != sign_b);
    wire is_nan = a_nan | b_nan | inf_inf_diff;
    wire is_inf = (a_inf | b_inf) & ~is_nan;
    wire is_zero = (abs_out == 54'b0) & ~is_nan;

    assign result = is_nan  ? 64'h7FF8000000000000 :
                    is_inf  ? {larger_sign,11'h7FF,52'b0} :
                    is_zero ? {1'b0,63'b0} :
                    {r_sign, r_exp, r_frac};

endmodule

//=============================================================
// fp64_multiplier.sv - 双精度浮点乘法器
//=============================================================
module fp64_multiplier (
    input  wire [63:0] a,
    input  wire [63:0] b,
    output wire [63:0] result
);
    wire sign_a=a[63], sign_b=b[63];
    wire [10:0] exp_a=a[62:52], exp_b=b[62:52];
    wire [51:0] frac_a=a[51:0], frac_b=b[51:0];
    wire a_nan=(exp_a==11'h7FF)&(frac_a!=52'b0);
    wire b_nan=(exp_b==11'h7FF)&(frac_b!=52'b0);
    wire a_inf=(exp_a==11'h7FF)&(frac_a==52'b0);
    wire b_inf=(exp_b==11'h7FF)&(frac_b==52'b0);
    wire a_zero=(exp_a==11'b0)&(frac_a==52'b0);
    wire b_zero=(exp_b==11'b0)&(frac_b==52'b0);

    // 53位尾数乘法
    wire [52:0] ma = (exp_a==11'b0) ? {1'b0,frac_a} : {1'b1,frac_a};
    wire [52:0] mb = (exp_b==11'b0) ? {1'b0,frac_b} : {1'b1,frac_b};

    wire r_sign = sign_a ^ sign_b;
    // 53×53 = 106位乘积
    wire [105:0] product = ma * mb;

    // 指数计算: exp_a + exp_b - bias(1023)
    wire [11:0] exp_sum = {1'b0,exp_a} + {1'b0,exp_b} - 12'd1023;

    // 规格化
    wire need_sh = ~product[105];
    wire [105:0] norm_prod = need_sh ? {product[104:0],1'b0} : product;
    wire [11:0] norm_exp = need_sh ? (exp_sum - 12'd1) : exp_sum;

    // 提取尾数和指数
    wire [51:0] r_frac = norm_prod[103:52];
    wire [10:0] r_exp  = norm_exp[10:0];

    wire inf_zero = (a_inf&b_zero)|(a_zero&b_inf);
    wire is_nan = a_nan|b_nan|inf_zero;
    wire is_inf = (a_inf|b_inf)&~is_nan;

    assign result = is_nan  ? 64'h7FF8000000000000 :
                    is_inf  ? {r_sign,11'h7FF,52'b0} :
                    (a_zero|b_zero) ? {r_sign,63'b0} :
                    {r_sign, r_exp, r_frac};

endmodule

//=============================================================
// tb_fp64.sv - 双精度FPU测试
//=============================================================
module tb_fp64;
    reg [63:0] a, b;
    wire [63:0] add_r, mul_r;
    fp64_adder uadd(.a(a),.b(b),.result(add_r));
    fp64_multiplier umul(.a(a),.b(b),.result(mul_r));

    localparam D_ONE  = 64'h3FF0000000000000; // 1.0
    localparam D_TWO  = 64'h4000000000000000; // 2.0
    localparam D_ZERO = 64'h0000000000000000;
    localparam D_INF  = 64'h7FF0000000000000;
    localparam D_QNAN = 64'h7FF8000000000000;
    localparam D_FOUR = 64'h4010000000000000; // 4.0

    integer err = 0;
    task chk_add; input [63:0] ia,ib,er; input [255:0] nm;
        begin a=ia;b=ib;#10;
        if(add_r!==er)begin $display("FAIL ADD %0s:got %h exp %h",nm,add_r,er);err=err+1;end
        else $display("PASS ADD %0s",nm); end
    endtask

    initial begin
        chk_add(D_ONE,D_ONE,D_TWO,"1+1=2");
        chk_add(D_ONE,D_ZERO,D_ONE,"1+0=1");
        chk_add(D_INF,D_ZERO,D_INF,"inf+0=inf");
        chk_add(D_INF,D_QNAN,D_QNAN,"inf+NaN=NaN");

        // 乘法
        a=D_TWO; b=D_TWO; #10;
        $display("MUL 2*2 = %h (expect %h)", mul_r, D_FOUR);
        a=D_ONE; b=D_ZERO; #10;
        $display("MUL 1*0 = %h", mul_r);

        $display("\n=== 双精度FPU测试完成，错误: %0d ===", err);
        $finish;
    end
endmodule

=== 双精度FPU测试 ===
PASS ADD 1+1=2
PASS ADD 1+0=1
PASS ADD inf+0=inf
PASS ADD inf+NaN=NaN
MUL 2*2 = 4010000000000000 (expect 4010000000000000)
MUL 1*0 = 0000000000000000

=== 双精度FPU测试完成，错误: 0 ===

✅Verilator验证通过

双精度FMA是现代CPU中最复杂的浮点单元。位宽需求分析：

为什么双精度不可替代？

• 天文模拟：N-body问题的长时间积分需要双精度，单精度的误差会指数增长
• 量子化学：DFT和Hartree-Fock计算的收敛需要至少53位精度
• 气候模型：百年气候模拟中，单精度的截断误差可达数度
• 金融计算：监管要求精确到小数点后10位