搜索
1方程求根的牛顿迭代法实验报告目录方程求根的牛顿迭代法实验报告................................................................................11实验环境及要求：.............................................................................................12实验内容：.........................................................................................................13实验步骤：.........................................................................................................23.1浮点数的平方根......................................................................................23.2整数c的平方根......................................................................................53.3求解浮点数倒数......................................................................................63.4奇数C的逆𝑐−1mod223...................................................................84实验结果及分析说明.......................................................................................104.1浮点数开方............................................................................................104.2整数开根号............................................................................................104.3浮点数求倒数........................................................................................104.4整数倒数................................................................................................105附录一（源代码）...........................................................................................111实验环境及要求：Inteli5core，Ubuntu12.04X64操作系统，使用Linuxc++嵌入AT&T汇编编码实现。在此64位linux系统中，浮点数64位，1位符号位+11位阶码（偏移为1023）+52位尾数。2实验内容：求浮点数C的平方根√𝐶求整数C的平方根⌊√𝐶⌋求浮点数C的倒数𝐶−1求奇数C模223的逆𝐶−1mod22323实验步骤：3.1浮点数的平方根原理：求解f(x)=𝑥2−𝑐；f(x)=0的正根即为c的平方根。根据牛顿迭代格式∅(x)=x−𝑓(𝑥)𝑓′(𝑥)构造迭代式：{𝑥0𝑥𝑛+1=12(𝑥𝑛+𝑐𝑥𝑛)在计算机表示中，浮点数c可以表示为尾数与阶数即：c=（1.𝑥0𝑥1…𝑥52）2𝑠其中，取m为尾数部分，则m∈[1,2]。选取𝑥0=pc+q，其中，p、q为常数其中曲线为函数y=√𝑥，其子x=1,2处的两点为A(1,1)，B(2,√2)，记直线AB为𝑙2，与𝑙2平行且与曲线相切的直线为𝑙1。与这两条直线平行且在正中间的直线𝑙0即为估计𝑥0的函数。可求得：𝑙2：y=(√2−1)x+2−√2𝑙1：y=(√2−1)x+14(√2+1)则可得：𝑙0：y=(√2−1)𝑥+18⁄(9−3√2)即p=√2−1，q=18⁄(9−3√23此时误差函数β(t)=𝑥0−√𝑐=(√2−1)𝑡2−𝑡+18⁄(9−3√2)，其中令t=√𝑐，则tϵ[1,√2]。该函数的最大值和最小值为：58⁄√2−78⁄和78⁄−58⁄√2。则误差|𝑥0−√𝑐|=58⁄√2−78⁄=0.0088835，二进制为0.0000001…此时二进制误差已达到6位，64位机器中浮点数尾数部分为52位，则只需要4次迭代即可。实现分为三个部分：取出阶数s计算尾数m的平方根阶数返回，若s为偶数，√𝑐=2𝑠2⁄√𝑚；否则√𝑐=2s−12⁄√2√𝑚代码：doublemy_sqrt(doublec){//stepone:setcin[1,2]andsisthemultiplicityofc;ints;doubleresult;doublep=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));//step2//theprocessofNewtoniteration.3timesiterationsreachto64bitsprecisiongasm(movsd%3,%%xmm0;mulsd%1,%%xmm0;addsd%2,%%xmm0;movl$5,%%ecx;movq%3,%%xmm1;4movabsq$4611686018427387904,%%rax;//thefirstoperanis2.that'stheformatofdoublefloatmovq%%rax,%%xmm2;start:;cmpl$0,%%ecx;jzendloop;decl%%ecx;movq%%xmm1,%%rax;divsd%%xmm0,%%xmm1;addsd%%xmm1,%%xmm0;movq%%rax,%%xmm1;divsd%%xmm2,%%xmm0;jmpstart;endloop:movq%%xmm0,%%rdx;:+d(result):m(p),m(q),m(c));//step3//addthemultipicitys/2asm(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;jzend3;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));returnresult;}代码说明，代码使用的是linuxc++嵌套汇编，浮点运算采用的是xmm指令。代码分为三个部分，一是提取11位的阶数，然后开根号，再将阶数返回。53.2整数c的平方根其牛顿迭代式：{𝑥0𝑥𝑛+1=⌊12⁄(𝑥𝑛+⌊𝑐𝑥𝑛⁄⌋)⌋其中取𝑥0⌊√𝑐⌋迭代时若𝑥𝑛+1𝑥𝑛继续迭代，否则停止，此时𝑥𝑛=⌊√𝐶⌋取𝑥0=2𝑠，2s−1⌊√𝑐⌋≤2𝑠.c为32位整数，则0≤s≤16，采用二分法求解s，最多只需要5次迭代。代码实现unsignedmy_isqrt(unsignedc){//step1,choosex0=2^s;sisin[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]unsignedx0;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;jlleft;incl%%ecx;movl%%ecx,%%esi;jmpstart_loop;left:;movl$1,%%eax;decl%%ecx;shll%%cl,%%eax;addl$1,%%eax;mull%%eax;cmpl%%eax,%%ebx;jgeend_loop;movl%%ecx,%%edi;6jmpstart_loop;end_loop:;incl%%ecx;movl$1,%%eax;shll%%cl,%%eax;endpp:;:+a(x0):m(c));//step2theprocessofNewtoniteration,tilldon'tdecreaseunsignedresult;asm(movl%1,%%ebx;movl%2,%%ecx;start_loop_2:;movl%%ebx,%%eax;movl$0,%%edx;divl%%ecx;addl%%ecx,%%eax;shrl$1,%%eax;cmpl%%eax,%%ecx;jleend_loop_2;movl%%eax,%%ecx;jmpstart_loop_2;end_loop_2:movl%%ecx,%%eax;:+a(result):m(c),m(x0));returnresult;}代码说明，代码分为2个部分，分别是step1求解x0，step2牛顿迭代。3.3求解浮点数倒数其函数为：f(x)=1𝑥⁄−𝑐其牛顿迭代形式：{𝑥0𝑥𝑛+1=𝑥𝑛(2−𝑐𝑥𝑛)初值𝑥0=pc+q的选取（p，q为常量），使得7∅=max([𝑥0−1𝑐⁄1𝑐⁄]2)=𝑚𝑎𝑥[𝑐(𝑝𝑐+𝑞)−1]2，c∈[1,2]取得最小值，其中，误差函数为𝑥0−1𝑐⁄1𝑐⁄。上式中令𝑓(𝑐)=[𝑐(𝑝𝑐+𝑞)−1]2=𝑝2𝑐4+2𝑝𝑞𝑐3+(𝑞2−2𝑝)𝑐2−2𝑝𝑐+1其取得最大值，并使得其最大值在p,q为变量时最小，则p=−817⁄，q=2417⁄。此时误差|𝑐𝑥0−1|2−8。则要达到52位（64位浮点数有52位尾数）的精度，只需要4次牛顿迭代。代码实现doublemy_inverse(doublec){//stepone:setcin[1,2]andsisthemultipicityofc;ints;asm(movq%2,%%rcx;movq$0x7ff0000000000000,%%rax;andq%%rcx,%%rax;shrq$52,%%rax;subq$1023,%%rax