用单精度浮点模拟双精度的运算，实现扩展浮点精度的算法

实际上，任何精度的浮点运算也可以像BigInteger那样计算，可以用多个浮点数据表示更高精度的浮点数据。而这就要求实现带余数的浮点四则运算。

带余数的浮点四则运算的定义：

加法和减法，余数和结果的和能精确表示运算值。

乘法和除法，余数和结果必须同号，对于乘法，余数和结果的和能精确表示运算值，而对于除法，余数就是原来的意思。

 

不多解释，代码如下：（修复某些可能存在错误）

#include <iostream>
#include <xmmintrin.h>

using namespace std;


typedef struct {
    float hi;
    float lo;
} doublf;

/**
 * BigFloat的实际值为：
 * 使用这个数据类型的缺点是数据存储冗余，优点是便于计算。
**/
typedef struct {
    int32_t len;//data中有效数据的长度
    int32_t exp;//指数部分
    float data[0];//浮点数据，每个都是1.x的数据
} BigFloat;

inline int32_t ifabs(float a){
    int32_t p=*(int32_t*)&a;
    return p&0x7fffffff;
}

//o=5
doublf addf(float a,float b){//带余数加法
    doublf ret;
    ret.hi=a+b;
    if(ifabs(a)>ifabs(b))ret.lo=b-(ret.hi-a);
    else ret.lo=a-(ret.hi-b);
    return ret;
}

//o=5
doublf subf(float a,float b){//带余数减法
    doublf ret;
    ret.hi=a-b;
    if(ifabs(a)>ifabs(b))ret.lo=(a-ret.hi)-b;
    else ret.lo=a-(ret.hi+b);
    return ret;
}

//o=14
doublf mule(float a,float b){
    doublf ret;
    float m,n;//结合整数运算优化
    int32_t ap,bp,am,bm;
    uint32_t mp;
    ap=*(int32_t*)&a;
    bp=*(int32_t*)&b;
    ret.hi=n=a*b;
    am=ap&0xff800000;
    bm=bp&0xff800000;
    am=am+bm-0x4b000000;
    m=*(float*)&am;
    ap|=0x800000;
    bp|=0x800000;
    mp=(uint32_t)ap*(uint32_t)bp;
    ap=(mp&0x7fffff)|am;
    ret.lo=*(float*)&ap;
    if(mp&0x800000){
        if(n+m!=n)ret.lo-=m;
    }
    else ret.lo-=m;
    return ret;
}

//o=14
doublf mulf(float a,float b){//带余数乘法
    unsigned int csr=_mm_getcsr();
    _mm_setcsr((csr & ~_MM_ROUND_MASK) | _MM_ROUND_TOWARD_ZERO);
    doublf ret=mule(a,b);
    _mm_setcsr(csr);
    return ret;
}

//o=17
doublf divf(float a,float b){//带余数除法
    unsigned int csr=_mm_getcsr();
    _mm_setcsr((csr & ~_MM_ROUND_MASK) | _MM_ROUND_TOWARD_ZERO);
    doublf ret,tmp;
    ret.hi=a/b;
    tmp=mule(ret.hi,b);
    ret.lo=(a-tmp.hi)-tmp.lo;
    _mm_setcsr(csr);
    return ret;
}

//检验正确
//o=20,8
doublf add(doublf a,doublf b){//最後一位会产生误差
    doublf ret,tmp;
    /* 计算 */
    if(ifabs(a.hi)<ifabs(b.hi)){
        tmp=a;
        a=b;
        b=tmp;
    }
    ret=addf(b.lo,a.lo);
    tmp=addf(ret.hi,b.hi);
    ret.hi=tmp.hi;
    ret.lo+=tmp.lo;
    tmp=addf(tmp.hi,a.hi);
    ret.hi=tmp.hi;
    ret.lo+=tmp.lo;//即使高低位重叠也不会有影响
    return ret;
}

//检验正确
//o=22,9
doublf sub(doublf a,doublf b){//最後一位会产生误差
    doublf tmp;
    /* 计算 */
    tmp.hi=-b.hi;
    tmp.lo=-b.lo;
    return add(a,tmp);
}

//检验正确
doublf mul(doublf a,doublf b){
#if 0//最後一位会产生误差
//o=137,25
    unsigned int csr=_mm_getcsr();
    _mm_setcsr((csr & ~_MM_ROUND_MASK) | _MM_ROUND_TOWARD_ZERO);
    doublf ret,tmp,tmp2,tmp3,tmp4;
    /* 计算 */
    tmp=mule(a.lo,b.lo);
    tmp2=mule(a.lo,b.hi);
    tmp3=mule(a.hi,b.lo);
    tmp4=mule(a.hi,b.hi);
    _mm_setcsr(csr);
    ret=add(tmp,tmp2);//应允许误差存在
    ret=add(ret,tmp3);
    ret=add(ret,tmp4);
    return ret;
#else//最後两位会产生误差
//o=25,8
    doublf ret=mulf(a.hi,b.hi);
    ret.lo+=(a.lo*b.lo+a.lo*b.hi)+a.hi*b.lo;
    ret=addf(ret.hi,ret.lo);//克服高低位重叠造成的影响
    return ret;
#endif
}

//检验正确
doublf div(doublf a,doublf b){//最後两位会产生误差
#if 1//o=27,10
    doublf ret,tmp;
    float d=b.hi+b.lo;
    ret.hi=a.hi/d;
    unsigned int csr=_mm_getcsr();
    _mm_setcsr((csr & ~_MM_ROUND_MASK) | _MM_ROUND_TOWARD_ZERO);
    /* 计算 */
    tmp=mule(b.hi,ret.hi);
    //tmp2=mule(b.lo,ret.hi);
    _mm_setcsr(csr);
    ret.lo=((((a.hi-tmp.hi)-tmp.lo)-ret.hi*b.lo)+a.lo)/d;//精度低但更快
    //a=sub(a,tmp);
    //a=sub(a,tmp2);
    //ret.lo=a.hi/d;
    ret=addf(ret.hi,ret.lo);//还是应该修正，以克服高低位重叠的影响
    return ret;
#else//o=24,8
    doublf ret;
    float d=b.hi+b.lo;
    ret=divf(a.hi,d);
    ret.lo=(ret.lo-ret.hi*((d-b.hi)-b.lo)+a.lo)/d;
    return ret;
#endif
}

inline double to_double(doublf a){
    return (double)a.hi+a.lo;
}

inline doublf to_doublf(double a){
    doublf ret;
    ret.hi=a;
    ret.lo=a-ret.hi;
    return ret;
}

/**
 * 这样就能有办法用单精度模拟双精度（可以有2^-46的精度），甚至模拟任何精度的浮点数据。
 * 关键是要支持带余数的浮点四则运算。
**/

int main()
{
    int i;
    uint64_t tmp;
    doublf df,af,tf;
    double dd=45671.154545,ad=1.001245648454,td,tt;
    do{
        tf=df=to_doublf(dd);
        af=to_doublf(ad);
        cin>>i;
        while(i--){
            tt=to_double(df=add(df,af));
            tmp=*(uint64_t*)&tt;
            cout<<"add="<<hex<<tmp<<endl;
            td=dd+ad;
            tmp=*(uint64_t*)&td;
            cout<<"chk="<<hex<<tmp<<endl;
            dd=tt;
        };
        df=tf;
        dd=to_double(tf);
        cin>>i;
        while(i--){
            tt=to_double(df=sub(df,af));
            tmp=*(uint64_t*)&tt;
            cout<<"sub="<<hex<<tmp<<endl;
            td=dd-ad;
            tmp=*(uint64_t*)&td;
            cout<<"chk="<<hex<<tmp<<endl;
            dd=tt;
        };
        df=tf;
        dd=to_double(tf);
        cin>>i;
        while(i--){
            tt=to_double(df=mul(df,af));
            tmp=*(uint64_t*)&tt;
            cout<<"mul="<<hex<<tmp<<endl;
            td=dd*ad;
            tmp=*(uint64_t*)&td;
            cout<<"chk="<<hex<<tmp<<endl;
            dd=tt;
        };
        df=tf;
        dd=to_double(tf);
        cin>>i;
        while(i--){
            tt=to_double(df=div(df,af));
            tmp=*(uint64_t*)&tt;
            cout<<"div="<<hex<<tmp<<endl;
            td=dd/ad;
            tmp=*(uint64_t*)&td;
            cout<<"chk="<<hex<<tmp<<endl;
            dd=tt;
        };
        cin>>dd;
        cin>>ad;
    }while(ad!=0.0);
    return 0;
}