计算有效值的朋友有福了：32位整数开平方终于做出来（汇编的，在WINAVR中初步测试通过，还在 (amobbs.com 阿莫电子技术论坛) -

lqzh1872 发表于 2009-11-11 22:39:59

计算有效值的朋友有福了：32位整数开平方终于做出来（汇编的，在WINAVR中初步测试通过，还在

;汇编部分
;#include<avr/io.h>
;#include<avr/interrupt.h>
;#include"MyMath.h"

.section .text

;uint16_t mysqrt_i(uint32_t a)
.global mysqrt_i
mysqrt_i:

mov r20,r24 ;r25:r24-->p
mov r21,r25

;r21:r20:r23:r22 32bit
eor r26,r26

eor r25,r25
eor r24,r24

eor r18,r18 ;r20:r19:r18-->4p
eor r19,r19
eor r27,r27

;1

lsl r22
rol r23
rol r20
rol r21
rol r26

lsl r22
rol r23
rol r20
rol r21
rol r26 ;左移2位

cp r18,r26
breq mysqrt_i_next2_0
;商1
ori r24,1

sbc r26,r18
ori r18,4
mysqrt_i_next2_0:

;2
lsl r24 ;

lsl r22
rol r23
rol r20
rol r21
rol r26

lsl r22
rol r23
rol r20
rol r21
rol r26 ;左移2位
cp r18,r26
brlo mysqrt_i_s21
lsl r18
rjmp mysqrt_i_s2_1
mysqrt_i_s21:
;商1
ori r24,1
sbc r26,r18
lsl r18 ;4p
ori r18,0x4
mysqrt_i_s2_1:
;3
lsl r24 ;
;lsl r18 ;4p

lsl r22
rol r23
rol r20
rol r21
rol r26

lsl r22
rol r23
rol r20
rol r21
rol r26 ;左移2位
cp r18,r26
brlo mysqrt_i_s31
lsl r18 ;4P
rjmp mysqrt_i_s3_2
mysqrt_i_s31:
ori r24,1
sbc r26,r18
lsl r18
ori r18,0x4 ;4P

mysqrt_i_s3_2:

;4
lsl r24 ; p

lsl r22
rol r23
rol r20
rol r21
rol r26

lsl r22
rol r23
rol r20
rol r21
rol r26 ;左移2位
cp r18,r26
brlo mysqrt_i_s41
lsl r18
rjmp mysqrt_i_s4_2
mysqrt_i_s41:
ori r24,1
sbc r26,r18
lsl r18 ;4p
ori r18,4
mysqrt_i_s4_2:

;一个字节已经移完，部分商已经产生了4位，R22已经移空
;接下来我移入R22：R26中
;5
lsl r24

lsl r23
rol r20
rol r21
rol r26
rol r22

lsl r23
rol r20
rol r21
rol r26 ;左移2位
rol r22

;
cp r19, r22
brlo mysqrt_i_s51
cp r18,r26
brlo mysqrt_i_s51
lsl r18
rjmp mysqrt_i_s5_2
mysqrt_i_s51:
ori r24,1
sbc r26,r18
sbci r22,0
lsl r18
ori r18,4
mysqrt_i_s5_2:

;6
lsl r24

lsl r23
rol r20
rol r21
rol r26
rol r22

lsl r23
rol r20
rol r21
rol r26 ;左移2位
rol r22

;
cp r19, r22
brlo mysqrt_i_s61
cp r18,r26
brlo mysqrt_i_s61
lsl r18
rjmp mysqrt_i_s6_2
mysqrt_i_s61:
ori r24,1
sbc r26,r18
sbci r22,0
lsl r18
ori r18,4
mysqrt_i_s6_2:

;7
lsl r24

lsl r23
rol r20
rol r21
rol r26
rol r22

lsl r23
rol r20
rol r21
rol r26 ;左移2位
rol r22

cp r19, r22
brlo mysqrt_i_s71
cp r18,r26
brlo mysqrt_i_s71
lsl r18
rol r19
rjmp mysqrt_i_s72
mysqrt_i_s71:
;商1
ori r24,1
sbc r26,r18
sbc r22,r19
lsl r18 ;4p
rol r19
ori r18,4
mysqrt_i_s72:

;8
lsl r24

lsl r23
rol r20
rol r21
rol r26
rol r22

lsl r23
rol r20
rol r21
rol r26 ;左移2位
rol r22

cp r19,r22
brlo mysqrt_i_s81
brne mysqrt_i_s82_2
cp r18,r26
brlo mysqrt_i_s81
mysqrt_i_s82_2:
lsl r18
rol r19
rjmp mysqrt_i_s82

mysqrt_i_s81:
;商1
ori r24,1
sbc r26,r18
sbc r22,r19
lsl r18
rol r19
ori r18,4

mysqrt_i_s82:

;9 r23已经空了,此时可以移进R23:R22:R26
lsl r24
rol r25

lsl r20
rol r21
rol r26
rol r22
rol r23

lsl r20
rol r21
rol r26 ;左移2位
rol r22
rol r23

cp r27, r23
brlo mysqrt_i_s91
cp r19,r22
brlo mysqrt_i_s91
brne mysqrt_i_s92_2
cp r18,r26
brlo mysqrt_i_s91
mysqrt_i_s92_2:
lsl r18
rol r19
rjmp mysqrt_i_s92
mysqrt_i_s91:
;商1
ori r24,1
sbc r26,r18
sbc r22,r19
sbci r23,0

lsl r18
rol r19
ori r18,4
mysqrt_i_s92:

;10
lsl r24
rol r25

lsl r20
rol r21
rol r26
rol r22
rol r23

lsl r20
rol r21
rol r26 ;左移2位
rol r22
rol r23

cp r27, r23
brlo mysqrt_i_s101
cp r19,r22
brlo mysqrt_i_s101
brne mysqrt_i_s102_2
cp r18,r26
brlo mysqrt_i_s101

mysqrt_i_s102_2:
lsl r18
rol r19
rjmp mysqrt_i_s102
mysqrt_i_s101:
;商1
ori r24,1
sbc r26,r18
sbc r22,r19
sbci r23,0
lsl r18
rol r19
ori r18,4
mysqrt_i_s102:

;11
lsl r24
rol r25

lsl r20
rol r21
rol r26
rol r22
rol r23

lsl r20
rol r21
rol r26 ;左移2位
rol r22
rol r23

cp r27, r23
brlo mysqrt_i_s111
cp r19,r22
brlo mysqrt_i_s111
brne mysqrt_i_s112_2
cp r18,r26
brlo mysqrt_i_s111
mysqrt_i_s112_2:
lsl r18
rol r19
rjmp mysqrt_i_s112
mysqrt_i_s111:
;商1
ori r24,1
sbc r26,r18
sbc r22,r19
sbci r23,0
lsl r18
rol r19
ori r18,4
mysqrt_i_s112:

;12
lsl r24
rol r25

lsl r20
rol r21
rol r26
rol r22
rol r23

lsl r20
rol r21
rol r26 ;左移2位
rol r22
rol r23

cp r27, r23
brlo mysqrt_i_s121
cp r19,r22
brlo mysqrt_i_s121
brne mysqrt_i_s122_2
cp r18,r26
brlo mysqrt_i_s121
mysqrt_i_s122_2:
lsl r18
rol r19
rjmp mysqrt_i_s122
mysqrt_i_s121:
;商1
ori r24,1
sbc r26,r18
sbc r22,r19
sbci r23,0
lsl r18
rol r19
ori r18,4 ;4p
mysqrt_i_s122:

;13
;R20已经空了
lsl r24
rol r25

lsl r21
rol r26
rol r22
rol r23
rol r20

lsl r21
rol r26 ;左移2位
rol r22
rol r23
rol r20

tst r20
brne mysqrt_i_s131
tst r23
brne mysqrt_i_s131
cp r19,r22
brlo mysqrt_i_s131
brne mysqrt_i_s132_2
cp r18,r26
brlo mysqrt_i_s131
mysqrt_i_s132_2:
lsl r18
rol r19
rjmp mysqrt_i_s132
mysqrt_i_s131:
;商1
ori r24,1
sec
sbc r26,r18
sbc r22,r19
sbc r23,0
sbc r20,0
lsl r18
rol r19
ori r18,4
mysqrt_i_s132:

;14
lsl r24
rol r25

lsl r21
rol r26
rol r22
rol r23
rol r20

lsl r21
rol r26 ;左移2位
rol r22
rol r23
rol r20

tst r20
brne mysqrt_i_s141
tst r23
brne mysqrt_i_s141
cp r19,r22
brlo mysqrt_i_s141
brne mysqrt_i_s142_2
cp r18,r26
brlo mysqrt_i_s141
mysqrt_i_s142_2:
lsl r18
rol r19

rjmp mysqrt_i_s142
mysqrt_i_s141:
;商1
ori r24,1
sec
sbc r26,r18
sbc r22,r19
sbci r23,0
sbci r20,0
lsl r18
rol r19

ori r18,4;4p
mysqrt_i_s142:

;15

lsl r24
rol r25

lsl r21
rol r26
rol r22
rol r23
rol r20

lsl r21
rol r26 ;左移2位
rol r22
rol r23
rol r20

tst r20
brne mysqrt_i_s151
tst r23
brne mysqrt_i_s151
cp r19,r22
brlo mysqrt_i_s151
brne mysqrt_i_s152_2
cp r18,r26
brlo mysqrt_i_s151
mysqrt_i_s152_2:
lsl r18
rol r19
rol r27
rjmp mysqrt_i_s152
mysqrt_i_s151:
;商1
ori r24,1
sec
sbc r26,r18
sbc r22,r19
sbci r23,0
sbci r20,0
lsl r18
rol r19
rol r27
ori r18,4 ;4p

mysqrt_i_s152:

;16
lsl r24
rol r25

lsl r21
rol r26
rol r22
rol r23
rol r20

lsl r21
rol r26 ;左移2位
rol r22
rol r23
rol r20

tst r20
brne mysqrt_i_s161
cp r27,r23
brlo mysqrt_i_s161
brne mysqrt_i_s162_2
cp r19,r22
brlo mysqrt_i_s161
brne mysqrt_i_s162_2
cp r18,r26
brlo mysqrt_i_s161
mysqrt_i_s162_2:
lsl r18
rol r19
rol r20
rjmp mysqrt_i_s162
mysqrt_i_s161:
;商1
ori r24,1 ;其实后面的可以不用做了
sec
sbc r26,r18
sbc r22,r19
sbc r23,r27
sbci r20,0
lsl r18
rol r19
rol r20
ori r18,4
mysqrt_i_s162:

;17

ret

lqzh1872 发表于 2009-11-11 22:41:56

/*C语言测试部分*/
#include<avr/io.h>
#include<avr/interrupt.h>
#include"MyMath.h"
int main()
{
unsigned int a,b;
a=0xff00;
b=mysqrt_i(0xff00ff01);
return 0;
}

lqzh1872 发表于 2009-11-11 22:43:15

/*MyMath.h的内容：*/

#ifndef MYMATH_H
#define MYMATH_H
#include <stdint.h>
uint16_t mysqrt_i(uint32_t a);
#endif

lqzh1872 发表于 2009-11-11 22:45:16

使用AVR-GDB调试了很久，希望对大家有帮助，来不及做详细注释了，先放上再说!

coody 发表于 2009-11-11 23:18:07

有什么意义？编译器不提供标准的数学计算功能？

xyqdoudou 发表于 2009-11-11 23:42:48

库函数不一定是速度最快最节省资源的你用AT13做个浮点数除法看看还剩多少

这段程序我相信大家很多人都看过只能用神奇来形容

float InvSqrt (float x)
{
float xhalf = 0.5f*x;
int i = *(int*)&x;
i = 0x5f3759df - (i >> 1); // 计算第一个近似根
x = *(float*)&i;
x = x*(1.5f - xhalf*x*x); // 牛顿迭代法
return x;
}
主要是 0x5f3759df 这个数简直匪夷所思

我测试过很准...我看不太懂汇编希望楼主看一下这个也许对你有些启发
另请教高人 0x5f3759df 这玩意怎么来的。。
通常的牛顿迭代法都是需要编写循环语句的，但是在这里，因为只迭代两次，所以并不需要循环。
i = 0x5f3759df - (i >> 1); 简直...我只知道那是个浮点数而且很准想不通
望高人解答

qianhng 发表于 2009-11-12 07:04:58

顶，希望LZ把浮点数的也做出来。

zhiwei 发表于 2009-11-12 08:36:21

我以前也用汇编写过了。论坛坏了就没有了，

zhiwei 发表于 2009-11-12 08:38:30

重新给大家传一次吧。IAR下的。

//U16 Sqrt(U32 num)
//input : num -> r19:r16
//output: res -> r17:r16
//use volatile registers r9:r0,
//r1:r0 -> result
//r30_r3:r2 -> result left shift temp
//r22:r20 -> num left shift extend
//r23-> loop counter
Sqrt32:
      clr r0    ;clr result
      clr r1
      movwr2, r0;clr result shift
      movwr30,r0;//const - 0
      movwr20,r0
      movwr22,r0
      ldi r23,0x11;//do 17 times
wait:
      rcall shiftb;//shift 2 bit to calculate
      rcall shiftb
      dec r23 ;//check if finished
      cpi r23, 0x1
      breqcal
      cpi r20, 0;//if data is not 0, start calculate
      breqwait //shift until msb is not 0,and to calculate
cal:
      inc r2
      cp r20, r2;//extend - (4 result + 1)
      cpc r21, r3
      cpc r22, r30
      brcchig
      clc
      rjmpcalend
hig:
      sub r20, r2
      sbc r21, r3
      sbc r22, r30
      sec
calend:
      rol r0 ;result reflash
      rol r1
      dec r23
      breqcend;//check if finished
      movwr2, r0;//init 4 * result
      clr r30
      rcall shifta;//result and input number shift
      rcall shifta
      rjmpcal
cend:
      movwr16,r0
   ret

shifta://route shift data
      lsl r2
   rol r3
   rol r30
shiftb: lsl r16
   rol r17
   rol r18
   rol r19
   rol r20
   rol r21
   rol r22
      ret

qianhng 发表于 2009-11-12 08:38:44

你的浮点开平方执行一次要大约多少机器周期？

voidx 发表于 2009-11-12 08:40:56

貌似太长了

AVR FFT里面有有符号32位开方汇编，貌似很短。

zhiwei 发表于 2009-11-12 08:41:55

单片机我还没有用过浮点。双字节乘法、除法什么的都是自己写汇编函数，感觉效率会高一些。

snoopyzz 发表于 2009-11-12 08:48:34

感觉靠不准，需要实际测试和库函数的优劣...

pldjn 发表于 2009-11-12 08:51:32

[转自http://yueweitang.org/blog/posts/inverse-square-root-algorithm-analysis.html]

下面这个求1/\sqrt{x}的函数号称比直接调用sqrt库函数快4倍，来自游戏Quake III的源代码。

float InvSqrt(float x){
float xhalf = 0.5f * x;
int i = *(int*)&x;
i = 0x5f3759df - (i >> 1);
x = *(float*)&i;
x = x * (1.5f - xhalf * x * x);
return x;
}

我们这里分析一下它的原理（指程序的正确性，而不是解释为何快）。

分析程序之前，我们必须解释一下float数据在计算机里的表示方式。一般而言，一个float数据x共32个bit，和int数据一样。其中前23位为有效数字M_x，后面接着一个8位数据E_x表示指数，最后一位表示符号，由于这里被开方的数总是大于0，所以我们暂不考虑最后一个符号位。此时

x=1.M_x 2^{E_x-127}

如果我们把计算机内的浮点数x看做一个整数I_x，那么

I_x = 2^{23}E_x+M_x

现在开始逐步分析函数。这个函数的主体有四个语句，分别的功能是：

int i = *(int*)&x; 这条语句把x转成i=I_x。

i = 0×5f3759df - (i>>1); 这条语句从I_x计算I_{1/\sqrt{x}}。

y = *(float*)&i; 这条语句将I_{1/\sqrt{x}}转换为1/\sqrt{x}。

y = y*(1.5f - xhalf*y*y); 这时候的y是近似解；此步就是经典的牛顿迭代法。迭代次数越多越准确。

关键是第二步 i = 0x5f3759df - (i>>1); 这条语句从I_x计算I_{1/\sqrt{x}}，原理:

令y=1/\sqrt{x}，用x=(1+m_x)2^{e_x}和y=(1+m_y)2^{e_y}带入之后两边取对数，再利用近似表示\log_2(1+z)\sim z+\delta，算一算就得到

I_y = \frac{2}{3}(127-\delta)2^{23}-I_x/2

若取\delta=0.0450465679168701171875，\frac{2}{3}(127-\delta)2^{23}就是程序里所用的常量0x5f3759df。至于为何选择这个\delta，则应该是曲线拟合实验的结果。

qianhng 发表于 2009-11-12 09:03:23

我在avr studio里测试了下楼主的mysqrt_i，执行一次374个机器周期，感觉他把一些循环展开，用空间换时间，所以代码就长了。

cowboy 发表于 2009-11-12 09:53:02

mark

phone 发表于 2009-11-12 09:56:12

代吗确实够长的了。

lqzh1872 发表于 2009-11-12 11:39:18

我做过了测试使用它一次要356个Cycle

lqzh1872 发表于 2009-11-12 11:40:12

16MHz下要22.24uS

zengxy 发表于 2009-11-12 11:52:05

Mark

lqzh1872 发表于 2009-11-12 11:53:46

去掉后面部分还可以再省8个周期

lqzh1872 发表于 2009-11-12 11:59:13

空间换时间，因为我要的速度非常快，保证系统响应时间很重要

coody 发表于 2009-11-12 12:47:49

【13楼】 pldjn

老美一个教授专门花时间研究了这个天才的函数，很是天才。

我比较奇怪，函数作者是写游戏软件的，他怎么就不说明一下呢？

ivws 发表于 2009-11-12 12:59:17

有必要记号一下!

xyqdoudou 发表于 2009-11-12 15:23:17

是写游戏软件的这个函数的作者可能是因为要处理3D画面用到开方吧

qianhng 发表于 2009-11-12 16:51:19

做了个浮点开平方，执行一次要大约537机器周期，可以直接替代winavr里的sqrt()，还在测试看有无错误

qianhng 发表于 2009-11-12 20:59:27

不比不知道，一比吓一跳：
float i,j,k,l,M;
uint32_t i1;
uint16_t i2;
i1=635;
i2=mysqrt_i(i1);
j=3.1;
执行到这里 Cycle counter = 502
i=sqrt(j);
执行到这里 Cycle counter = 4810
k=fsqrt_qianhng(j);
执行到这里 Cycle counter = 5336
M=InvSqrt (j);
执行到这里 Cycle counter = 12742
l=(i-k)*10000000.0;
i2 =(uint16_t) l;
while(1);

通过计算得到函数的机器周期：
sqrt(j)          4308
fsqrt_qianhng(j) 526
InvSqrt (j)    7406

计算结果：
sqrt(j)          0x3f e1 5e 04
fsqrt_qianhng(j) 0x3f e1 5e 04
InvSqrt (j)    0x7f 80 00 00

看样子是InvSqrt (j)出错了^_^

gaolf_2000 发表于 2009-11-12 21:33:57

mark一下

thomasdu 发表于 2009-11-12 21:37:17

这个帖子很有营养哦

qianhng 发表于 2009-11-12 22:28:51

上辛苦工作n小时的劳动成果
;floatting point sqrt
; by qianhng,qianhng@163.com
;2009-11-12 VER 1.3.R
;版权所有、欢迎翻录，请保留所有说明，若程序有问题，麻烦举报于qianhng@163.com
;另：欢迎免费使用，若是对你有所帮助，同样欢迎随意打赏到支付宝帐号qianhng@163.com

使用方法：
http://cache.amobbs.com/bbs_upload782111/files_21/ourdev_502730.jpg
(原文件名:fsqrt_qianhng.jpg)
文件：
点击此处下载 ourdev_502729.rar(文件大小:836字节) (原文件名:fsqrt_qianhng_s.rar)

usingavr 发表于 2009-11-13 08:28:00

InvSqrt原先是在PC上用的，你在AVR上用的时候int是16位的，就不对了
而且InvSqrt算的是开方的倒数，在3D计算上开方的倒数比开方更有用

hollypower 发表于 2009-11-13 09:47:06

记号

hollypower 发表于 2009-11-13 09:55:26

留名

wear778899 发表于 2009-11-13 10:19:24

记号

zhuzi1441 发表于 2009-11-13 10:49:43

mark

nicksean 发表于 2009-11-13 15:51:49

mark

XA144F 发表于 2009-11-13 16:37:42

1=1×1
1+3=2×2
1+3+5=3×3
……
以此类推：1+3+5+……+29+31=16×16

即对于8位（最大255），最多16次减法，16次加法，16次移位和16次比较。

但是扩展到16位，理论上是各项工作完成256次，那就太复杂了……

呵呵，我的程序就是这么做的。

lqzh1872 发表于 2009-11-13 17:22:11

使用int32_t就对了

cjr82123 发表于 2009-11-13 18:41:27

mark

voidx 发表于 2009-11-13 20:49:26

8楼的程序很好，

mark

xml2028 发表于 2009-11-13 21:24:19

mark

fugeone 发表于 2009-11-13 23:09:02

mark

letyoufly 发表于 2009-11-14 12:35:57

phone 发表于 2009-12-24 15:51:11

试用了29楼的程序，发现编译后占的FLASH和SRAM大了很多，SRAM占了240字节。

beiyunhe 发表于 2011-2-20 16:53:40

Mark

rqzzx 发表于 2011-4-1 22:02:27

学习了，谢谢

jackielau 发表于 2011-4-1 22:29:23

mark

yuzr 发表于 2011-4-2 00:48:19

mark

hzn1948 发表于 2011-4-2 08:32:00

好帖！

yangjinke1118 发表于 2011-4-2 10:00:12

mark

tang0571 发表于 2011-4-2 10:04:57

支持这类讨论～～顶～～

anning 发表于 2011-4-2 13:50:46

Mark

xiaobenyi 发表于 2011-4-2 13:59:02

36的方法很有意思哦,不过只能是整数的平方

dkfwy 发表于 2013-3-18 20:24:52

下载试用一下。

sync765 发表于 2013-3-20 15:17:53

很强悍的帖子
学习

xwkm 发表于 2013-3-20 17:20:46

qianhng 发表于 2009-11-12 20:59 static/image/common/back.gif
不比不知道，一比吓一跳：
float i,j,k,l,M;
uint32_t i1;

InvSqrt没错。它求出来的是倒数……

ahong2hao 发表于 2013-4-27 17:38:30

mark ,刚好用到.

kevinhh 发表于 2013-4-28 19:00:58

相当好的东西，MARK

三且寺 发表于 2017-2-10 14:30:11

做信号有效值用的上

凌晨一点 发表于 2017-6-22 19:27:42

xyqdoudou 发表于 2009-11-11 23:42
库函数不一定是速度最快最节省资源的你用AT13做个浮点数除法看看还剩多少

这段程序我相信大家很多人都 ...

不错的算法，验证下先

jackman 发表于 2017-6-23 00:46:59

根据中值定理计算的比较简洁。

页: [1]

amobbs.com 阿莫电子技术论坛's Archiver

计算有效值的朋友有福了：32位整数开平方终于做出来（汇编的，在WINAVR中初步测试通过，还在