Cfloat — 实现¶

文件结构¶

cfloat/
├── tiny_cfloat.h        (24 行 — 类型定义 + 声明)
└── tiny_cfloat.c        (70 行 — 实现)

依赖：tiny_math_config.h（通过头文件）、<math.h>（sqrtf, cosf, sinf, atan2f, logf, fabsf）。

实现概览¶

所有函数都是纯 C，按值传递、按值返回。没有动态分配或平台特定的分发——模块完全可移植。

函数	行数	关键实现细节
`tiny_cf_add`	3	逐元素加法
`tiny_cf_sub`	3	逐元素减法
`tiny_cf_mul`	3	标准 \((a+bi)(c+di) = (ac-bd) + i(ad+bc)\)
`tiny_cf_div`	20	Smith 鲁棒算法（见下方）
`tiny_cf_conj`	2	虚部取反
`tiny_cf_abs`	2	`sqrtf(re² + im²)`，零值防护
`tiny_cf_arg`	2	`atan2f(im, re)`
`tiny_cf_from_polar`	3	\(r\cdot\cos\theta + i\cdot r\cdot\sin\theta\)
`tiny_cf_sqrt`	8	主分支
`tiny_cf_log`	7	\(\ln\\|z\\| + i \cdot \arg(z)\)

关键算法¶

除法：Smith 鲁棒方法¶

朴素除法计算分母 \(|b|^2 = b_r^2 + b_i^2\)，对系数较大时可能溢出。Smith 方法通过缩放计算避免中间值溢出：

tiny_cfloat_t tiny_cf_div(tiny_cfloat_t a, tiny_cfloat_t b)
{
    float abs_re = fabsf(b.re), abs_im = fabsf(b.im);
    float mag = (abs_re > abs_im) ? abs_re : abs_im;

    // 近零分母防护
    if (mag < TINY_MATH_MIN_POSITIVE_INPUT_F32) {
        return sentinel_large_value(a);
    }

    if (abs_re >= abs_im) {
        float k = b.im / b.re;
        float den = b.re + k * b.im;
        r.re = (a.re + k * a.im) / den;
        r.im = (a.im - k * a.re) / den;
    } else {
        float k = b.re / b.im;
        float den = b.im + k * b.re;
        r.re = (a.re * k + a.im) / den;
        r.im = (a.im * k - a.re) / den;
    }
    return r;
}

平方根：主分支¶

tiny_cfloat_t tiny_cf_sqrt(tiny_cfloat_t a)
{
    if (a.re == 0.0f && a.im == 0.0f) return zero;
    float mag = tiny_cf_abs(a);
    r.re = sqrtf((mag + a.re) * 0.5f);
    r.im = (a.im >= 0.0f)
        ?  sqrtf((mag - a.re) * 0.5f)
        : -sqrtf((mag - a.re) * 0.5f);
    return r;
}

虚部的符号遵循 C99 约定：\(\operatorname{sgn}(\Im(\sqrt{z})) = \operatorname{sgn}(\Im(z))\)，分支切割在负实轴上。

在高层模块中的使用¶

tiny_cfloat_t 是被以下模块使用的基本类型：

decomp — Householder 反射需要复数运算来提取特征值
eigen — 特征值存储为 tiny_cfloat_t[]（实特征值的 \(\Im(\lambda)=0\)）
modal — 模态分析对复数值的特征值和振型进行操作
iterative — Arnoldi/Lanczos 迭代产生复数 Ritz 值

// 示例：Francis QR 返回的特征值为复数
tiny_cfloat_t eigs[10];
tiny_eig_general_f32(A, 10, eigs, workspace, ws_size);

for (int i = 0; i < 10; i++) {
    printf("λ%d = %f + i·%f\n", i, eigs[i].re, eigs[i].im);
}