In C, how to divide a float or double number by 2 raised to the power i?

user_7305567 0
Viewed 148

I have written code in C for Atmel's microcontroller SAM E70 that processes 32-bit wide integer values. For further calculations, I normalize the integer value to 0...1.0 as follows:

#define DIV4294967296 ((double) 1.0) / ((double) 4294967296.0)
.
.
double doubleValue;
doubleValue = ((double) intValue) * DIV4294967296;

I understand I could just subtract 32 from the exponent of doubleValue and thereby avoid the more expensive multiplication. I am aware of ldexp() which allows to multiply the exponent by 2 raised to the power i but I couldn't find anything that allows me to explicitly read, manipulate, and write back the exponent of a double. Performing all these steps may actually not be faster than performing the multiplication, so a direct subtraction of 32 from the exponent would be ideal. How is this generally done in C? More importantly, how is this done best with ARM's Cortex V7 instruction set?

Addendum: to answer Eric's question, here's what Atmel Studio 7 shows me as disassembled code for the use of ldexp, scalbn and the multiplication with 0x1p-32:

uint32_t intV = 123456;
 ldr    r3, [pc, #424]
 str    r3, [r7, #28]
double doubleV0 = ((double) intV) * DIV4096;
 ldr    r3, [r7, #36]        
 vmov   s15, r3      
 vcvt.f64.u32   d7, s15      
 vldr   d6, [pc, #272]       
 vmul.f64   d7, d7, d6       
 vstr   d7, [r7, #24]       
double doubleV1 = ldexp(intV, -32);
 ldr    r3, [r7, #28]
 vmov   s15, r3
 vcvt.f64.u32   d7, s15
 mvn    r0, #31
 vmov.f64   d0, d7 
 ldr    r3, [pc, #408]
 blx    r3
 vstr   d0, [r7, #16]
double doubleV2 = scalbn(intV, -32);
 ldr    r3, [r7, #28]
 vmov   s15, r3
 vcvt.f64.u32   d7, s15
 mvn    r0, #31 
 vmov.f64   d0, d7
 ldr    r3, [pc, #384]
 blx    r3
 vstr   d0, [r7, #8]
double doubleV3 = intV * 0x1p-32;
 ldr    r3, [r7, #28]
 vmov   s15, r3
 vcvt.f64.u32   d7, s15
 vldr   d6, [pc, #164]
 vmul.f64   d7, d7, d6
 vstr   d7, [r7]

It looks like none of these match any ARM instructions (like the C function fabs() is compiled straight to the assembly instruction vabs). ldexp and scalbn are coded the same way. The multiplication with 0x1p-32 is coded the same way as my initial multiplication that let me raise my question.

Addendum 2: to show what code it compiles to based on chqrlie's suggestion:

double doubleV4 = ((double) intV);
 vstr   d7, [r7]    
*(uint64_t *)&doubleV4 -= 32ULL << 52;
 mov    r3, r7       
 ldrd   r2, r3, [r3]         
 mov    r1, r7       
 adds   r4, r2, #0       
 adc    r5, r3, #4261412864      
 strd   r4, r5, [r1]

It looks to me like it's the least expensive implementation.

Final verdict: I like chqrlie's answer as it may be of use for those among us where multiplications are too slow. In my case though, I've run an interrupt based routine and measured the execution time for both my initial code and chqrlie's alternative and they run exactly the same time if best optimization (-O3) is used with GCC 9.3.1.

exponent arm double c
1 Answers

If you can assert that the double is stored using the IEEE 754 double-precision binary floating-point format: binary64, with the same endianness an alignment requirements as the 64-bit integers, and its value is large enough for the result to still be a normal value, you can hack the representation directly with this expression which should compile to 2 or 3 instructions:

*(uint64_t *)&doubleValue -= 32ULL << 52;

Yet this form of type punning may cause trouble with agressive optimizers because it violates the C aliasing rules as the value of type double is accessed via a pointer to a different type that is not a character pointer. A better form of type punning can be used, via a union, that will work correctly with most compilers:

union { double d; uint64_t u; } u = doubleValue;
u.u -= 32ULL << 52;
doubleValue = u.d;

To avoid C aliasing issues completely, you can use memcpy:

uint64_t u;
memcpy(&u, &doubleValue, sizeof u);
u -= 32ULL << 52;
memcpy(&doubleValue, &u, sizeof u);

A good optimising compiler should convert these memcpy calls to single instructions.

user_45932670
Related