Vdsp(bf561)中的浮点运算（7）：float乘法运算

快乐虾

http://blog.csdn.net/lights_joy/

[e-mail]

本文适用于

ADSP-BF561

Visual DSP++ 5.0 (update 6)

Vdsp dual processor simulate

欢迎转载，但请保留作者信息

1.1 Vdsp对float乘法运算的处理

在vdsp下，可以很简单地用：

float mul(float x, float y)

{

float r = x * y;

return r;

}

来完成浮点乘法运算，编译器自动将里面的乘法操作转换为___float32_mul的函数调用，这个函数的调用实现在libdsp/fpmult.asm中，在这个文件的开头说明了这个函数的用法：

/***************************************************************************

Copyright (c) 2000-2008 Analog Devices Inc. All rights reserved.

****************************************************************************

File name : fpmult.asm

This function performs 32 bit floating point multiplication. Implemention

is based on the algorithm mentioned in the reference. Some more conditionS

are added in the present algorithm to take care of various testcases.

Registers used:

Operands in R0 & R1

R0 - X operand,

R1 - Y operand

R2 - R7

Special case:

1) If (x AND y) == 0, Return 0,

2) Overflow : If(Exp(x) + Exp(y) > 254,

Return 0X7F80000 or 0xFF80000

depending upon sign of X and y.

3) Underflow : If(Exp(x) + Exp(y) <= -149, RETURN 0.

Reference : Computer Architecture a Quantitative Approach 2nd Ed.

by Jhon L Hennessy and David Patterson

!!NOTE- Uses non-standard clobber set in compiler:

DefaultClobMinusBIMandLoopRegs

Remember to change the #pragma regs_clobbered in fpmult.c in softfloat if you

change this clobber set

**************************************************************/

1.2 当x和y都为0

看实现代码：

R2 = R0^R1; /* sign of the result of X * Y */

P0 = R2; /* Store sign */

R2 = R0 << 1; /* Remove the sign bit of X */

CC = R2 == 0;

R3 = R1 << 1; /* Remove the sign bit of Y */

CC |= AZ;

IF CC JUMP .HANDLE_ZERO_OP;

………..

.HANDLE_ZERO_OP:

// One operand is zero. If the other is NaN or inf, return NaN, otherwise

// zero

R2 = R2 | R3;

R2 = R2 >> 24;

R3 = MAXBIASEXP+1; // Max biased exponent

CC = R2 == R3;

R0 = 0;

IF !CC JUMP .SIGN_RESULT; // Return zero (signed appropriately)

R0 = -1; // return default NAN

RTS; // no regs to restore on this path

当x和y都为0时，直接返回0，此时的操作需要的CYCLE数为12。

这里一个比较有意思的事情是0是可以有正负号的，也就是说如果计算0.0 * (-0.0)，那么结果的符号位保持为负！

1.3 当x为nan或者inf

看代码：

[--SP] = (R7:4); /* Push registers R4-R7 */

R4 = MAXBIASEXP+1; /* Max biased exponent */

/* Get exponents. */

R2 = R2 >> 24; /* Exponent of X operand */

R3 = R3 >> 24; /* Exponent of Y operand */

CC = R2 == R4;

IF CC JUMP .HANDLE_NAN_INF_X;

……………..

.HANDLE_NAN_INF_X:

// Return signed inf if X=inf, and Y is either inf or a non-zero number

// Otherwise return NaN

R5 = R1 << 1; // If Y = 0

CC = AZ; // Return NaN

IF CC JUMP .RET_DEFAULT_NAN;

CC = R3 < R4; // if y exp is valid

R5 = R1 << 9; // or if Y significand is zero

CC |= AZ;

R5 = R0 << 9; // then Y is either inf or a valid number

CC &= AZ; // And if X is inf, then return inf, otherwise NaN

IF !CC JUMP .RET_DEFAULT_NAN;

.RET_INF:

R0 = 0x7F8 (Z);

R0 <<= 21;

(R7:4) = [SP++]; /* Pop registers R4-R7 */

.SIGN_RESULT:

CC = P0 < 0; /* Extract sign, and set into result */

R0 = ROT R0 BY -1; /* by rotating back with CC */

RTS;

.RET_DEFAULT_NAN:

R0 = -1;

(R7:4) = [SP++]; /* Pop registers R4-R7 */

RTS;

l 当x为inf且y为inf或者非0值

结果为inf。

mul(inf, inf)的结果为inf。

mul(inf, 1.0)的结果为inf。

此时的符号位仍然根据两个操作数的符号位来确定。

此时的操作需要的CYCLE数为53。

l 当x为inf且y为0

结果为NAN。

l 当x为inf且y为-inf

结果为-inf

1.4 当y为nan或者inf

看代码：

CC = R3 == R4;

IF CC JUMP .HANDLE_NAN_INF_Y;

…………..

// Handle certain identities concerning multiplying by NaNs and +/-inf

.HANDLE_NAN_INF_Y:

// Swap operands and exponents

R5 = R0;

R0 = R1;

R1 = R5;

R5 = R2;

R2 = R3;

R3 = R5;

.HANDLE_NAN_INF_X:

………..

简单交换xy，然后和上一个判断条件相同。

1.5 向上溢出

当两个操作数的指数相加大于254时，认为溢出：

// Compute result exponent, and check for overflow

R4 = BIASEXP;

R5 = R2 + R3;

R5 = R5 - R4;

R4 <<= 1; // R4 now 254, max allowed exponent

CC = R4 < R5;

IF CC JUMP .RET_INF;

……………

.RET_INF:

R0 = 0x7F8 (Z);

R0 <<= 21;

(R7:4) = [SP++]; /* Pop registers R4-R7 */

.SIGN_RESULT:

CC = P0 < 0; /* Extract sign, and set into result */

R0 = ROT R0 BY -1; /* by rotating back with CC */

RTS;

理解这段代码有点费劲，因为float能表示的最大指数是127，那么应该是大于127就认为溢出，为什么要和254做比较呢？

判断条件应该是在

(e1 – 127) + (e2 - 127) > 127

时溢出，在上面的代码里相当于等价改为e1 + e2 – 127 > 254

为何不干脆改为e1 + e2 > 381呢？

奇怪！

溢出计算时的CYCLE值为52。

1.6 正常计算

在正常计算一个浮点乘法时，所用的CYCLE值为92，居然比做浮点加减所用的CYCLE少，一个出乎意料的结果！

1.7 向下溢出

在代码说明里给出的下溢条件是

Exp(x) + Exp(y) <= -149

为什么为是这样的条件呢？想不通！看来非得去看看Reference提及的书才行。

2 参考资料

Vdsp(bf561)中的浮点运算（6）：float加减运算(2009-8-13)

Vdsp(bf561)中的浮点运算（5）：float类型表示总结(2009-8-12)

Vdsp(bf561)中的浮点运算（4）：FLT_MAX(2009-8-12)

Vdsp(bf561)中的浮点运算（3）：FLT_MIN(2008-12-19)

Vdsp(bf561)中的浮点运算（2）：float的疑问(2008-12-18)

Vdsp(bf561)中的浮点运算（1）：文档的说法(2008-12-16)