北京邮电大学数值与符号计算实验报告

������2：y = (√2 − 1)x + 2 − √2
������1：y
=(√2
−
1)x
+
1 4
(√2
+
1)
则可得：������0：y = (√2 − 1)������ + 1⁄8 (9 − 3√2)
即 p = √2 − 1，q = 1⁄8 (9 − 3√2
2
此时 误差函数 β(t) = ������0 − √������ = (√2 − 1)������2 − ������ + 1⁄8 (9 − 3√2)，其中令 t=√������，
取������0 = 2������，2s−1 < ⌊√������⌋ ≤ 2������.
c 为 32 位整数，则0 ≤ s ≤ 16，采用二分法求解 s，最多只需要 5 次迭代。 代码实现 unsigned my_isqrt(unsigned c){
//step1,choose x0=2^s; s is in [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16] unsigned x0; asm( "movl %1,%%ebx;" "movl $0,%%esi;" "movl $16,%%edi;" "start_loop:;" "movl %%esi,%%ecx;" "addl %%edi,%%ecx;" "sarl $1,%%ecx;" "movl $1,%%eax;" "shll %%cl,%%eax;" "addl $1,%%eax;" "mull %%eax;" "cmpl %%eax,%%ebx;" "jl left;" "incl %%ecx;" "movl %%ecx,%%esi;" "jmp start_loop;" "left:;" "movl $1,%%eax;" "decl %%ecx;" "shll %%cl,%%eax;" "addl $1,%%eax;" "mull %%eax;" "cmpl %%eax,%%ebx;" "jge end_loop;" "movl %%ecx,%%edi;"
4
3.2 整数 c 的平方根
其牛顿迭代式：
������0 {������������+1 = ⌊1⁄2 (������������ + ⌊������⁄������������⌋)⌋
其中取������0 > ⌊√������⌋
迭代时若������������+1 < ������������ 继续迭代，否则停止，此时������������ = ⌊√������⌋
∅(x)
=
x
−
������(������) 构造迭代式：
������′(������)
������0
{������������+1
=
1 2
(������������
+
������ ������������)
在计算机表示中，浮点数 c 可以表示为尾数与阶数即：
c=（1.������0������1 … ������52）2������ 其中，取 m 为尾数部分，则 m ∈ [1,2]。
则tϵ[1, √2]。
该函数的最大值和最小值为：5⁄8 √2 − 7⁄8 和 7⁄8 − 5⁄8 √2。
则误差|������0 − √������| = 5⁄8 √2 − 7⁄8 = 0.0088835， 二进制为 0.0000001… 此时二进制误差已达到 6 位，64 位机器中浮点数尾数部分为 52 位，则只需要 4 次迭代即可。 实现分为三个部分： 取出阶数 s 计算尾数 m 的平方根
阶数返回，若 s 为偶数，√������ = 2������⁄2√������；否则√������ = 2s−1⁄2√2√������
代码： double my_sqrt(double c){
//step one: set c in [1,2] and s is the multiplicity of c; int s; double result; double p=0.4144142135623730951,q=0.5946699141100893; asm( "movq %2,%%rcx;" "movq $0x7ff0000000000000,%%rax;" "andq %%rcx, %%rax;" "shrq $52,%%rax;" "subq $1023,%%rax;" "movq %%rax,%%rdx;" "movq $0x800fffffffffffff,%%rax;" "andq %%rax,%%rcx;" "movq $0x3ff0000000000000,%%rax;" "orq %%rax,%%rcx;" :"+c"(c),"=d"(s) :"m"(c)); //step 2 //the process of Newton iteration. 3 times iterations reach to 64 bits precisiong asm( "movsd %3,%%xmm0;" "mulsd %1,%%xmm0;" "addsd %2,%%xmm0;" "movl $5,%%ecx;" "movq %3,%%xmm1;"
3
"movabsq $4611686018427387904,%%rax;" //the first operan is 2. that's the format of double float
"movq %%rax,%%xmm2;" "start:;" "cmpl $0,%%ecx;" "jz endloop;" "decl %%ecx;" "movq %%xmm1,%%rax;" "divsd %%xmm0,%%xmm1;" "addsd %%xmm1,%%xmm0;" "movq %%rax,%%xmm1;" "divsd %%xmm2,%%xmm0;" "jmp start;" "endloop:" "movq %%xmm0,%%rdx;" :"+d"(result) :"m"(p),"m"(q),"m"(c)); //step3 //add the multipicity s/2 asm( "movq %%rbx,%%rdx;" "shrq $1,%%rbx;" "movq $0x800fffffffffffff,%%rcx;" "andq %%rcx,%%rax;" "movq $0x3ff,%%rcx;" "addq %%rbx,%%rcx;" "salq $52,%%rcx;" "orq %%rcx,%%rax;" "andq $1,%%rdx;" "cmpl $0,%%edx;" "jz end3;" "movabsq $4609047870845172685,%%rcx;" "movq %%rax,%%xmm1;" "movq %%rcx,%%xmm0;" "mulsd %%xmm1,%%xmm0;" "movq %%xmm0,%%rax;" "end3:;" "movq %%rax,%%rdx;" :"+d"(result) :"a"(result),"b"(s)); return result; } 代码说明，代码使用的是 linux c++嵌套汇编，浮点运算采用的是 xmm 指令。代 码分为三个部分，一是提取 11 位的阶数，然后开根号，再将阶数返回。
方程求根的牛顿迭代法实验报告
目录
方程求根的牛顿迭代法实验报告................................................................................ 1 1 实验环境及要求：.............................................................................................1 2 实验内容：.........................................................................................................1 3 实验步骤：.........................................................................................................2 3.1 浮点数的平方根......................................................................................2 3.2 整数 c 的平方根......................................................................................5 3.3 求解浮点数倒数......................................................................................6 3.4 奇数 C 的逆 ������−1 mod 223 ...................................................................8 4 实验结果及分析说明.......................................................................................10 4.1 浮点数开方............................................................................................10 4.2 整数开根号............................................................................................10 4.3 浮点数求倒数........................................................................................10 4.4 整数倒数................................................................................................10 5 附录一（源代码）........................................................................................... 11
求整数 C 的平方根 ⌊√������⌋ 求浮点数 C 的倒数 ������−1 求奇数 C 模 223的逆 ������−1 mod 223
1
3 实验步骤：
3.1 浮点数的平方根
原理：求解 f(x) = ������2 − ������； f(x)=0 的正根即为 c 的平方根。
根据牛顿迭代格式
选取 ������0= pc + q ，其中，p、q 为常数
其中 曲线为函数 y=√������，其子 x=1,2 处的两点为 A(1,1)，B(2,√2)，记直线 AB
为������2，与������2平行且与曲线相来自百度文库的直线为������1。与这两条直线平行且在正中间的直线������0 即为估计������0的函数。 可求得：
1 实验环境及要求：
Intel i5 core，Ubuntu 12.04 X64 操作系统，使用 Linux c++ 嵌入 AT&T 汇编编码 实现。 在此 64 位 linux 系统中，浮点数 64 位，1 位符号位+11 位阶码（偏移为 1023） +52 位尾数。
2 实验内容：
求浮点数 C 的平方根 √������
5
"jmp start_loop;" "end_loop:;" "incl %%ecx;" "movl $1,%%eax;" "shll %%cl,%%eax;" "endpp:;" :"+a"(x0) :"m"(c)); //step2 the process of Newton iteration, till don't decrease unsigned result; asm( "movl %1,%%ebx;" "movl %2,%%ecx;"