64位数开根号c语言,平方根的C语言实现(三) ——最终程序实现
了解了浮点数的存储以及手算平方根的原理,我们可以考虑程序实现了。
先实现一个64位整数的平方根,根据之前的手算平方根,程序也不是那么难写了。
#include
uint64_t _sqrt_u64(uint64_t a)
{
int i;
uint64_t res;
uint64_t remain;
//0的平方根是0,特殊处理一下
if(a == 0ull)
return 0ull;
//找到最高位的1,并且产生平方根结果最高位的1
for(i=62;;i-=2)
if(a&(3ull<
res = 1ull;
remain = ((a&(3ull<>i) - 1ull;
i -= 2;
break;
}
//根据手算平方根的原理,依次产生各位结果
for(;i>=0;i-=2) {
//右移动两位,并把a接着的两位并入remain
remain = (remain<<2)|((a&(3ull<>i);
if(((res<<2)|1ull) <= remain) {
//产生新一位的1
remain = remain - ((res<<2)|1ull);
res = (res<<1)|1ull;
} else {
//产生新一位的0
res <<= 1;
}
}
return res;
}
其实,可以合在一起写,代码会短一些,但效率会低那么一点点,而且编译器应该不太容易优化。
#include
uint64_t _sqrt_u64(uint64_t a)
{
int i;
uint64_t res;
uint64_t remain;
res = remain = 0ull;
for(i=62;i>=0;i-=2) {
remain = (remain<<2)|((a&(3ull<>i);
if(((res<<2)|1ull) <= remain) {
remain = remain - ((res<<2)|1ull);
res = (res<<1)|1ull;
} else {
res <<= 1;
}
}
return res;
}
不过,我们不需要这个结果。
为了验证其正确性,我们来写个C语言的main
#include
#include
#include
uint64_t _sqrt_u64(uint64_t a);
int main()
{
uint64_t a, b;
scanf("%" PRIu64, &a);
b = _sqrt_u64(a);
printf("%" PRIu64 "\n",b);
return 0;
}
我们shell程序测试一下,我们当然不可能测试过每一个64bits的数,这个运算量太大,不现实。我们可以用随机取一部分来测试。
#!/bin/bash
#编译gcc -O2 -s sqrt_u64.c main_sqrt_u64.c -o a.out
#随机测试10000个数for((i=0;i<10000;i++));do#随机产生bits0~64,如果是0,代表测试的数就是0
#如果不是0,则代表要产生的数二进制可以有多少位
let bits=RANDOM%65
if [ $bits -eq 0 ]; thenx=0y=0
else#产生一个bits位的二进制数x
x=$({
#最高位1echo -n 1#之后每位随机产生for((j=1;j
echo -n $xdone})
#用bc将x转换成十进制
x=$(echo ‘obase=10;ibase=2;‘"$x" |bc)
#用bc计算x的平方根取整,理论上和我们的C语言计算一致
y=$(echo ‘sqrt(‘"$x"‘)‘ |bc)fi#z是我们的C语言计算结果
z=$(echo $x | ./a.out)
#比较,如果不一致,就报错if [ $y -ne $z ];then
echo$x $y $z error
exit1
fi
done
echo OK
测试结果表明,我们的C语言还是可以得到正确的结果的。
再来回忆下第一节里讲过的浮点数结构,
S(1bits) | N(8bits) | A(23bits)
对于浮点数a*2n,
1<=a<2,n为整数,
如果n是偶数,
那么a*2n的平方根是sqrt(a)*2n/2,也满足1<=sqrt(a)<2,n/2是整数;
如果n为奇数,
那么a*2n的平方根是sqrt(2*a)*2(n-1)/2,也满足1<=sqrt(2*a)<2,(n-1)/2是整数。
所以此处要用a或者2*a来开平方根,
回忆一下浮点数的结构,单精度浮点数的精度是23位。
表示的是科学计数法a*2n的a减去1的部分,那么加上整数1可以用二进制24位表示。
于是,我们就想,一个二进制48位或47位长的数,平方根是二进制24位。那么,我们就可以用一个48位或47位的二进制整数的平方根计算结果的小数部分。
nan/inf/-inf以及负数的平方根都是nan,
0.0的平方根是0.0,
-0.0的平方根是-0.0(可能只是某些库里是这样的),
以上都可以在计算的时候特殊化一下。
规格数(就是用科学计数法表示的浮点数)的平方根也是规格数,
S=0,N=0,A>0代表的是A*2-149,也就是(A*2)*2-150,
我们稍微计算一下,可以明白,所有的此类数的平方根都在规格数表示的范围内。
于是,有了以下的代码。
#include
static uint32_t _sqrt_(uint64_t a)
{
int i;
uint64_t res;
uint64_t remain;
res = remain = 0ull;
//之前整数平方根被直接优化,我们只需要求47位或者48位整数的平方根
for(i=46;i>=0;i-=2) {
remain = (remain<<2)|((a&(3ull<>i);
if(((res<<2)|1ull) <= remain) {
remain = remain - ((res<<2)|1ull);
res = (res<<1)|1ull;
} else {
res <<= 1;
}
}
return (uint32_t)res;
}
float mysqrtf(float f)
{
union {
float f;
uint32_t u;
} n;
uint32_t N,A;
int _N, i;
uint64_t _A;
n.f = f;
if(n.u == 0x80000000 || n.u == 0x00000000) /* 0.0/-0.0 */
return n.f;
N = (n.u&(0xff<<23))>>23;
if(N==0xff||(n.u&0x80000000)) { /* inf/-inf/nan/ f < 0.0*/
n.u = 0x7fc00000; /* nan */
return n.f;
}
if(N!=0x0) { /* 用科学计数法表示的规格数 */
A = (n.u&0x7fffff)|0x800000;
_N = (int)N - 127;
if(N&0x1) {
_A = (uint64_t)A<<23;
} else {
_A = (uint64_t)A<<24;
_N--;
}
} else { //A*2^(-149)这种表示方式的浮点数
//还是需要找最高位
for(i=22;;i--)
if(n.u&((0x1)<
break;
//然后需要移位,要区分奇数和偶数
if(i&0x1) {
_N = i-149;
_A = (uint64_t)n.u << (46-i);
} else {
_N = i-150;
_A = (uint64_t)n.u << (47-i);
}
}
//小数部分
A = _sqrt_(_A);
//指数部分
N = (uint32_t)(_N/2+127);
//得到结果
n.u = (A&0x7fffff)|(N<<23);
return n.f;
}
同样,也写个测试用的程序,对inf/-inf/nan/0.0/-0.0以及负数不测了,这些很简单。
#include
#include
#include
#include
#include
#include
int main(int argc, char **argv)
{
union {
float f;
uint32_t u;
} n;
uint32_t A,N;
float f,f2;
int i;
srand((unsigned)time(NULL));
//随机10000个数据
for(i=0;i<10000;i++) {
N = rand()%256;
if(N==255)
N=254;
A = 0x0;
A |= rand()%256;
A |= (rand()%256)>>8;
A |= (rand()%256)>>16;
n.u = (A&0x7fffff)|(N<<23);
f = sqrtf(n.f);
f2 = mysqrtf(n.f);
printf("%.60f %.60f\n",f,f2);
}
return 0;
}
结果发现,我们的程序和数学库里的sqrtf结果有细微差别。
于是,我们决定再加个小东西,就是四舍五入。之前我们用的是47位或者48位数开平方,为了四舍五入,我们需要多一位,于是就用49位或者50位数开平方。
修改一下mysqrtf,增加两位拿去开平方,_sqrt_也动一下。
#include
static uint32_t _sqrt_(uint64_t a)
{
int i;
uint64_t res;
uint64_t remain;
res = remain = 0ull;
//之前整数平方根被直接优化,我们只需要求49位或者50位整数的平方根
for(i=48;i>=0;i-=2) {//这里之前是46,改成48
remain = (remain<<2)|((a&(3ull<>i);
if(((res<<2)|1ull) <= remain) {
remain = remain - ((res<<2)|1ull);
res = (res<<1)|1ull;
} else {
res <<= 1;
}
}
return (uint32_t)res;
}
float mysqrtf(float f)
{
union {
float f;
uint32_t u;
} n;
uint32_t N,A;
int _N, i;
uint64_t _A;
n.f = f;
if(n.u == 0x80000000 || n.u == 0x00000000) /* 0.0/-0.0 */
return n.f;
N = (n.u&(0xff<<23))>>23;
if(N==0xff||(n.u&0x80000000)) { /* inf/-inf/nan/ f < 0.0*/
n.u = 0x7fc00000; /* nan */
return n.f;
}
if(N!=0x0) { /* 用科学计数法表示的规格数 */
A = (n.u&0x7fffff)|0x800000;
_N = (int)N - 127;
if(N&0x1) {
_A = (uint64_t)A<<25;
} else {
_A = (uint64_t)A<<26;
_N--;
}
} else { //A*2^(-149)这种表示方式的浮点数
//还是需要找最高位
for(i=22;;i--)
if(n.u&((0x1)<
break;
//然后需要移位,要区分奇数和偶数
if(i&0x1) {
_N = i-149;
_A = (uint64_t)n.u << (48-i);
} else {
_N = i-150;
_A = (uint64_t)n.u << (49-i);
}
}
//小数部分
A = _sqrt_(_A);
//四舍五入
A = (A+(A&0x1))>>1;
//指数部分
N = (uint32_t)(_N/2+127);
//得到结果
n.u = (A&0x7fffff)|(N<<23);
return n.f;
}
然后再测,准确无误。于是我们可以完工了。
64位数开根号c语言,平方根的C语言实现(三) ——最终程序实现相关推荐
- 3倍根号x加1分之一c语言,用C语言将一个数开根号后再取倒数的方法
在上学的时候,曾经看过有人写过这样的算法,就是将一个数开根号后再取倒数的算法,我本人也觉得十分巧妙,于是就将它积累了下来,让我们来看看是怎么回事: #include #include float my ...
- c语言指针倒数之和,用C语言将一个数开根号后再取倒数的方法
在上学的时候,曾经看过有人写过这样的算法,就是将一个数开根号后再取倒数的算法,我本人也觉得十分巧妙,于是就将它积累了下来,让我们来看看是怎么回事: #include #include float my ...
- scratch lenet(4): 开根号的C语言实现
文章目录 1. 目的 2 二分法求开根号 2.1 数学原理:单调函数 2.2 代码实现:注意事项 2.3 代码实现: 完整代码 2.4 验证结果 3. 牛顿法 3.1 数学原理:迭代求解 3.2 代码 ...
- 计算机上根号是哪一个,电脑上怎么哪个键是数学中的开根号啊
该楼层疑似违规已被系统折叠 隐藏此楼查看此楼 手工开根号法,只适用于任何一个整数或者有限小数开二次方. 因为网上写不出样式复杂的计算式,所以只能尽量书写,然后通过口述来解释: 假设一个整数145645 ...
- 计算机如何实现开根号?
今天看到一个问题:计算机如何实现开根号? 如何求一个数字的算术平方根(又叫开根号,或者开方)? 大家普遍都是用计算器直接计算的,对于程序员来说,就是调用sqrt()方法.但是其内部又是怎么实现的呢?下 ...
- 开根号的笔算算法图解_一个数的开根号怎么计算
一个数的开根号怎么计算2020-11-08 15:46:47文/钟诗贺 带根号的式子可以直接进行开平方的运算.一些特殊的根号运算有;√2≈1.414.1/2-√3≈0.5-1.732≈-1.232.2 ...
- 计算机上开根号是哪个键,电脑上怎么哪个键是数学中的开根号啊
该楼层疑似违规已被系统折叠 隐藏此楼查看此楼 手工开根号法,只适用于任何一个整数或者有限小数开二次方. 因为网上写不出样式复杂的计算式,所以只能尽量书写,然后通过口述来解释: 假设一个整数145645 ...
- maya python 开根号_maya python
胡泳滨MayaPython简易教程,如需转载,请标明出处地址: http://huyongbin.blogbus.com/c3363976/ 谢谢配合! MayaPython第一篇 - 介绍 大家好, ...
- 【codevs3119】高精度开根号(二分答案)
problem 高精度开根号 输入一个数 求平方根 solution 二分答案,如果mid*mid>原数就去找更小的,反之找更大的. 精度小于二忽略不计? 用到高精加,高精乘,加低精,除低精,比 ...
- android开根号,定点数开根号的性能问题
8种机械键盘轴体对比 本人程序员,要买一个写代码的键盘,请问红轴和茶轴怎么选? 开根号有两种比较常见的方式:牛顿迭代法和二分法. 二分法 public static double SqrtBinary ...
最新文章
- 简单的docker命令ubuntu系统
- elementui表格添加滚动条_如何给PDF文档添加超链接?
- 华为鸿蒙os生态,华为鸿蒙系统终于来了! 首款方舟编译器应用正式上架: 鸿蒙OS可用...
- React 第七章 条件渲染
- 软件测试bug文档模板,国家标准测试计划文档模板
- (十七)Activitivi5之组任务分配
- android 项目 功能 源码 eclipse的
- MATLAB 画图 x轴换成 字符串
- 移动平台开发项目(推箱子小游戏)
- 13.文本文件和二进制文件的区别
- idea打包servlet成war包部署在tomcat
- 全国勘察设计注册暖通空调工程师专业基础考试大纲(送审稿)
- STM32学习笔记——HC05
- [隐匿的学习笔记]JVM(2)运行时数据区
- 无法获得 VMCI 驱动程序的版本: 句柄无效解决方法
- 代理服务器的安全证书有问题 错误代码8,如何修复Internet Explorer 8中的证书错误...
- 爱荷华州立大学计算机学院,爱荷华州立大学最新qs世界排名
- Unity发布项目,记录日志并写入文件。
- PhalAPI学习笔记拓展篇 ———ADM模式中NotORM实现简单CURD
- 手残把下载文件夹位置移动到了D盘根目录,导致了一系列问题的解决方法
热门文章
- 微信网页授权 获取 unionId
- [影视源码]全民影院源码 综合影视HTML源码 无需更新搭建即可用
- 操作 神通数据库_国产神通数据库教程
- solidity教程(一)搭建僵尸工厂
- 给Intel AX200装上个Killer 1650X驱动
- ibm x3850 x5连接存储后,linux操作系统无法正常启动,《七小服公开课》— IBM X3850 X5服务器无法开机故障 处理步骤...
- navicat安装(linux)
- java cximage_CxImage的简单用法
- RestSharp介绍
- java架构师证书_java架构师证书怎么考?做架构师有什么要求?