找回密码
 欢迎注册
楼主: Ickiverar

[分享] 高精度除法中,试商法的 3/2 算法

[复制链接]
发表于 2022-1-6 11:35:01 | 显示全部楼层
当m+n很大时候,是用b倒数乘以a做的,所以,严格来说,m+n比较小的时候才是传统除法,这时候,不同的试商策略,对运算时间的影响是很大的。
所以这个问题其实还和m,n取值有关系
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2024-3-26 10:10:20 | 显示全部楼层
本帖最后由 knate 于 2024-3-26 10:12 编辑

这里有两个估商实现
//要求 iRadix <= 2^31,
//估商时 除数 + 1 结果作为估算值.
uint32_t div_3_2_int_new(uint32_t a, uint32_t b, uint32_t c, uint32_t iDivsor_up, uint32_t iDivsor_down) {
        //这个不是真正的除法,是估算.
        //而且有个前提,结果一定是一位的.
        if (a == 0) {
                return ((__int64)b * cStatus.iRADIX + c) / ((__int64)iDivsor_up * cStatus.iRADIX + iDivsor_down + 1);
        }

        if (iDivsor_down == cStatus.iRADIX - 1) {
                iDivsor_up += 1;
                return (unsigned int)(((__int64)a * cStatus.iRADIX + b) / iDivsor_up);
        }
        {
                uint64_t divsor = iDivsor_up; divsor *= cStatus.iRADIX;
                divsor += iDivsor_down;
                //aR^2 + bR+c / dR+f
                // a*up * 2^32 + a * down + bR + c;
                uint64_t first, second, third;        //2^32
                uint64_t up, down;
                up = cStatus.iRADIX; up *= up; // R^2
                down = up & 0XFFFFFFFF;
                up >>= 32;

                first = a; first *= up;                //a *up

                second = b; second *= cStatus.iRADIX;
                down *= a;
                second += down;
                second += c;
                //a * down + bR + c;

                third = second & 0XFFFFFFFF;
                second >>= 32;

                first += second;

                uint32_t val = first >> 32;
                if (val) {
                        uint32_t num = asm_bsr(val);//查找first 最高位位置.
                        first <<= 31 - num;
                        first += third >> (num + 1);
                        divsor >>= num + 1;
                        divsor += 1;
                        return (uint32_t)(first / divsor);
                }
                first <<= 32;
                first += third;
                return (uint32_t)(first / (divsor + 1));
        }
        return 0;
}

uint32_t div_3_2_int_bin(uint32_t a, uint32_t b, uint32_t c, uint32_t iDivsor_up, uint32_t iDivsor_down) {
        //这个不是真正的除法,是估算.
        //而且有个前提,结果一定是一位的.
        if (a == 0) {
                return ((__int64)b * cStatus.iRADIX + c) / ((__int64)iDivsor_up * cStatus.iRADIX + iDivsor_down + 1);
        }

        if (iDivsor_down == cStatus.iRADIX - 1) {
                iDivsor_up += 1;
                return (unsigned int)(((__int64)a * cStatus.iRADIX + b) / iDivsor_up);
        }
        {
                uint64_t divsor = iDivsor_up; divsor *= cStatus.iRADIX;
                divsor += iDivsor_down + 1;
                iDivsor_up = divsor >> 32; iDivsor_down = divsor & 0XFFFFFFFF;
                //除数转换为2进制显示.

                unsigned int iQuotient = 0;
                uint64_t iRemainder_up, iRemainder_down;
                uint32_t iMultiple;
                // 2^64 / divsor  商: multiple  ,余数 高32位 iRemainder_up  ,余数 低32位 iRemainder_down

                uint64_t val;
                val = 0XFFFFFFFF; val <<= 32; val += 0XFFFFFFFF;
                iMultiple = val / divsor; val = val % divsor;
                if (val + 1 == divsor) {
                        //2^64 mod divsor 商和余数
                        iMultiple += 1;
                        iRemainder_up = 0;
                        iRemainder_down = 0;
                }
                else {
                        val += 1;
                        iRemainder_up = val >> 32;
                        iRemainder_down = val & 0XFFFFFFFF;
                }

                __int64 up, down;  // down 低32位,  up 高64位, 合计96位.  2进制表示.
                val = cStatus.iRADIX; val *= val;
                up = val >> 32;  down = val & 0XFFFFFFFF;
                up *= a;  down *= a;        //aR^2

                val = cStatus.iRADIX; val *= b;
                up += val >> 32;  down += val & 0XFFFFFFFF;
                //aR^2 + bR

                down += c;
                up += down >> 32; down &= 0XFFFFFFFF;

                //down 低32位,  up 高64位, 合计96位.  2进制表示.
                do {
                        //把 2^64拆分为商和余数.
                        uint32_t temp = up >> 32;
                        iQuotient += temp * iMultiple;
                        up &= 0XFFFFFFFF;
                        up += iRemainder_up * temp;
                        down += iRemainder_down * temp;

                        up += down >> 32;
                        down &= 0XFFFFFFFF;
                } while (up >> 32);
                up <<= 32; up += down;
                iQuotient += up / divsor;
                return iQuotient;
        }
        return 0;
}
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-11-21 17:34 , Processed in 0.024335 second(s), 16 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表