ugly math

Nov 12, 2006 42 Replies

store quotient, pop remainder, store, pop remainder,

That's easy enough - have two counters, one which counts down in binary, and the other which counts up in decimal. I did that for an electronics project at school many years ago.

The trick I know from the 8051 days is to mimic the shifting of you binary value in "BCD space". So imagine you shift your 64 bit source into a register one bit at a time from the lsb side. This means the same as: multiply by two and add the new bit. Which also means: add with carry the variable to itself with the new bit in the carry.

The neat thing of the 8051 was that this process could make use of some DA (decimal adjust) opcode.

You mentioned the 68332 has some BCD support (ABCD instruction?). Perhaps you can can do a similar thing. Just loop 64 times over your

20 digit BCD variable doing ABCD's.

Joop

There is *nothing* neat about an 8051!

See above; there's an easier way.

John

Hmmm, at, say, 30 us per loop, you could do the full 64-bit conversion in about 17 million years. Yup, that might win.

John

On that we have agreement ! Had to program the damn things once. Nasty devices - much preferred the Z8

John Larkin wrote:

store quotient, pop remainder, store, pop remainder,

No problem, except for the amount of code. Sorry about that!

Feel free to use it as long as I'm attributed as the author.

Terje

/* ascii.cpp * (C) Terje Mathisen */ #include #include #include

// Portable C version of 32-bit binary to ascii char *dtoa_c(char *buf, unsigned long n) { unsigned _int64 l64, h64; unsigned long l, h, i;

l64 = n;

l64 *= 2814749767; l64 += n >> 1; h = (unsigned long) (l64 >> 48); l = n - h * 100000;

l64 = l; h64 = h;

h64 *= 429497; l64 *= 429497;

buf[0] = (char) (h64 >> 32) + '0'; buf[5] = (char) (l64 >> 32) + '0'; h = (unsigned long) h64; l = (unsigned long) l64;

h = (h + 15) >> 4; l = (l + 15) >> 4;

for (i = 1; i < 5; i++) { h *= 10; l *= 10; buf[i] = (char) (h >> 28) + '0'; buf[i+5] = (char) (l >> 28) + '0'; h &= 0x0fffffff; l &= 0x0fffffff; } buf[10] = '\\0';

return buf; }

// x86 asm to convert 32-bit bin to ascii char *dtoa(char *buf, unsigned long n) { __asm { mov ebx,[n] mov eax,2814749767 mul ebx shr ebx,1 xor ecx,ecx mov edi,[buf] add eax,ebx adc ecx,edx mov eax,100000 shr ecx,16 // ECX = high part mov ebx,[n] imul eax,ecx // High part * 100k sub ebx,eax // Low part mov eax,429497 mul ecx mov ecx,eax add dl,'0' mov eax,429497 mov [edi],dl mul ebx mov ebx,eax add ecx,7 shr ecx,3 add dl,'0' mov [edi+5],dl add ebx,7 shr ebx,3 lea ecx,[ecx+ecx*4] mov edx,ecx and ecx,0fffffffh shr edx,28 lea ebx,[ebx+ebx*4] add dl,'0' mov eax,ebx shr eax,28 mov [edi+1],dl and ebx,0fffffffh add al,'0' mov [edi+6],al lea ecx,[ecx+ecx*4] lea ebx,[ebx+ebx*4] mov edx,ecx mov eax,ebx and ecx,07ffffffh shr edx,27 and ebx,07ffffffh shr eax,27 add dl,'0' add al,'0' mov [edi+2],dl mov [edi+7],al lea ecx,[ecx+ecx*4] lea ebx,[ebx+ebx*4] mov edx,ecx mov eax,ebx and ecx,03ffffffh shr edx,26 and ebx,03ffffffh shr eax,26 add dl,'0' add al,'0' mov [edi+3],dl mov [edi+8],al lea ecx,[ecx+ecx*4] shr ecx,25 lea ebx,[ebx+ebx*4] shr ebx,25 add cl,'0' add bl,'0' mov [edi+10],ah mov [edi+4],cl mov [edi+9],bl } return buf; }

#define E9 1000000000

typedef unsigned __int64 uint64; typedef unsigned long uint32;

char *qtoa(char *buf, uint64 q) { char b[32], *bz, t; uint64 m;

if (q

No, no!

Not even close:

A much "better" approach is to use a random number generator to generate

20 random decimal digits, then use the previous idea (but in reverse) to convert that decimal number to binary, and then compare against the original.

When identical, you've found the proper conversion.

Terje

- "almost all programming can be viewed as an exercise in caching"

Well, my board was limited to 12 bits, so it took something like 1/10 second max. Given that the input was by dip switches, the results arrived in realtime, so the algorithm was good enough.

I think that's cheating, because you are doing extra work for no reason.

But you could implement my algorithm with a Turing machine - that should be slow enough for any need.

store quotient, pop remainder, store, pop remainder,

And it works, as written!

John

Dividing by a Billion would get you nearly the maximum 60 bits worth for conversion.

JosephKK Gegen dummheit kampfen die Gotter Selbst, vergebens.   --Schiller

"John Larkin" wrote in message news: snipped-for-privacy@4ax.com...

Here's a little program I did in Microsoft Visual C, (yea I know I should have used gcc), but this should give you the idea:

---------------code starts-------------------------- /* ll2str.c : decompose 64 bit unsigned integer into string of 20 ASCII digits.

Method: Successive subtraction of powers of ten constants of each decimal digit of the output value.

Notes: This method is only efficient on processors where one 64 bit division takes more time than eight compares and one subtract per decimal digit.

*/

#include "stdafx.h"

const unsigned _int64 sar_digits[77] = {10000000000000000000, 5000000000000000000, 4000000000000000000, 2000000000000000000, 1000000000000000000, 500000000000000000, 400000000000000000, 200000000000000000, 100000000000000000, 50000000000000000, 40000000000000000, 20000000000000000, 10000000000000000, 5000000000000000, 4000000000000000, 2000000000000000, 1000000000000000, 5000000000000000, 4000000000000000, 2000000000000000, 1000000000000000, 500000000000000, 400000000000000, 200000000000000, 100000000000000, 50000000000000, 40000000000000, 20000000000000, 10000000000000, 5000000000000, 4000000000000, 2000000000000, 1000000000000, 500000000000, 400000000000, 200000000000, 100000000000, 50000000000, 40000000000, 20000000000, 10000000000, 5000000000, 4000000000, 2000000000, 1000000000, 500000000, 400000000, 200000000, 100000000, 50000000, 40000000, 20000000, 10000000, 5000000, 4000000, 2000000, 1000000, 500000, 400000, 200000, 100000, 50000, 40000, 20000, 10000, 5000, 4000, 2000, 1000, 500, 400, 200, 100, 50, 40, 20, 10};

void ll2str(char * out, unsigned _int64 in) { int i,j;

j=1; *out ='0'; for (i=0; i

Here are the hex values for the decimal constants I used:

0x8AC7230489E80000 0x4563918244F40000 0x3782DACE9D900000 0x1BC16D674EC80000 0x0DE0B6B3A7640000 0x06F05B59D3B20000 0x058D15E176280000 0x02C68AF0BB140000 0x016345785D8A0000 0x00B1A2BC2EC50000 0x008E1BC9BF040000 0x00470DE4DF820000 0x002386F26FC10000 0x0011C37937E08000 0x000E35FA931A0000 0x00071AFD498D0000 0x00038D7EA4C68000 0x0011C37937E08000 0x000E35FA931A0000 0x00071AFD498D0000 0x00038D7EA4C68000 0x0001C6BF52634000 0x00016BCC41E90000 0x0000B5E620F48000 0x00005AF3107A4000 0x00002D79883D2000 0x0000246139CA8000 0x000012309CE54000 0x000009184E72A000 0x0000048C27395000 0x000003A352944000 0x000001D1A94A2000 0x000000E8D4A51000 0x000000746A528800 0x0000005D21DBA000 0x0000002E90EDD000 0x000000174876E800 0x0000000BA43B7400 0x00000009502F9000 0x00000004A817C800 0x00000002540BE400 0x000000012A05F200 0x00000000EE6B2800 0x0000000077359400 0x000000003B9ACA00 0x000000001DCD6500 0x0000000017D78400 0x000000000BEBC200 0x0000000005F5E100 0x0000000002FAF080 0x0000000002625A00 0x0000000001312D00 0x0000000000989680 0x00000000004C4B40 0x00000000003D0900 0x00000000001E8480 0x00000000000F4240 0x000000000007A120 0x0000000000061A80 0x0000000000030D40 0x00000000000186A0 0x000000000000C350 0x0000000000009C40 0x0000000000004E20 0x0000000000002710 0x0000000000001388 0x0000000000000FA0 0x00000000000007D0 0x00000000000003E8 0x00000000000001F4 0x0000000000000190 0x00000000000000C8 0x0000000000000064 0x0000000000000032 0x0000000000000028 0x0000000000000014 0x000000000000000A

Ok, now I just feel stupid.

My previous two posts have errors.

Here's the code that does work.

I know you don't need it and noone but me cares.

/* ll2str.c : decompose 64 bit unsigned integer into string of 20 ASCII digits.

Method: Successive subtraction of powers of ten constants of each decimal digit of the output value.

Notes: This method is only efficient on processors where one 64 bit division takes more time than eight compares and three subtracts per decimal digit.

*/

#include "stdafx.h"

const unsigned _int64 sar_digits[73] = {10000000000000000000, // 0x8AC7230489E80000 5000000000000000000, // 0x4563918244F40000 4000000000000000000, // 0x3782DACE9D900000 2000000000000000000, // 0x1BC16D674EC80000 1000000000000000000, // 0x0DE0B6B3A7640000 500000000000000000, // 0x06F05B59D3B20000 400000000000000000, // 0x058D15E176280000 200000000000000000, // 0x02C68AF0BB140000 100000000000000000, // 0x016345785D8A0000 50000000000000000, // 0x00B1A2BC2EC50000 40000000000000000, // 0x008E1BC9BF040000 20000000000000000, // 0x00470DE4DF820000 10000000000000000, // 0x002386F26FC10000 5000000000000000, // 0x0011C37937E08000 4000000000000000, // 0x000E35FA931A0000 2000000000000000, // 0x00071AFD498D0000 1000000000000000, // 0x00038D7EA4C68000 500000000000000, // 0x0001C6BF52634000 400000000000000, // 0x00016BCC41E90000 200000000000000, // 0x0000B5E620F48000 100000000000000, // 0x00005AF3107A4000 50000000000000, // 0x00002D79883D2000 40000000000000, // 0x0000246139CA8000 20000000000000, // 0x000012309CE54000 10000000000000, // 0x000009184E72A000 5000000000000, // 0x0000048C27395000 4000000000000, // 0x000003A352944000 2000000000000, // 0x000001D1A94A2000 1000000000000, // 0x000000E8D4A51000 500000000000, // 0x000000746A528800 400000000000, // 0x0000005D21DBA000 200000000000, // 0x0000002E90EDD000 100000000000, // 0x000000174876E800 50000000000, // 0x0000000BA43B7400 40000000000, // 0x00000009502F9000 20000000000, // 0x00000004A817C800 10000000000, // 0x00000002540BE400 5000000000, // 0x000000012A05F200 4000000000, // 0x00000000EE6B2800 2000000000, // 0x0000000077359400 1000000000, // 0x000000003B9ACA00 500000000, // 0x000000001DCD6500 400000000, // 0x0000000017D78400 200000000, // 0x000000000BEBC200 100000000, // 0x0000000005F5E100 50000000, // 0x0000000002FAF080 40000000, // 0x0000000002625A00 20000000, // 0x0000000001312D00 10000000, // 0x0000000000989680 5000000, // 0x00000000004C4B40 4000000, // 0x00000000003D0900 2000000, // 0x00000000001E8480 1000000, // 0x00000000000F4240 500000, // 0x000000000007A120 400000, // 0x0000000000061A80 200000, // 0x0000000000030D40 100000, // 0x00000000000186A0 50000, // 0x000000000000C350 40000, // 0x0000000000009C40 20000, // 0x0000000000004E20 10000, // 0x0000000000002710 5000, // 0x0000000000001388 4000, // 0x0000000000000FA0 2000, // 0x00000000000007D0 1000, // 0x00000000000003E8 500, // 0x00000000000001F4 400, // 0x0000000000000190 200, // 0x00000000000000C8 100, // 0x0000000000000064 50, // 0x0000000000000032 40, // 0x0000000000000028 20, // 0x0000000000000014 10};// 0x000000000000000A

void ll2str(char * out, unsigned _int64 in) { int i,j;

j=1; *out ='0'; for (i=0; i

groups:comp.lang.asm.x86 has had a thread about various binary-to-ascii conversion algorithms recently, on anything but a P4 my algorithm which uses reciprocal multiplication to split the input into manageable parts and then run these in parallel wins.

On the P4 a simple digit by digit subtract/test loop is the fastest, due to the horrible cost of both MUL and shift on that cpu.

BTW, using (semi-binary) search to detect each decimal digit is almost certainly _slower_ than the much simpler linear search, due to the cost of all the branch mispredicts: A 0-to-9 loop will have a single branch miss, while a binary search will fail on average about two times per digit, the difference means that testing on average 5 digits is the better idea.

Terje

- "almost all programming can be viewed as an exercise in caching"

Of course, binary to BCD conversion is normally only done for the benefit of a human viewer, so even a lousy algorithm on a P4 will be fast enough to keep up with a typical display. An average P4 system will have plenty of memory as well, so a simple sprintf( asc, "%lld", bin ); should suffice, and you'll get nice formatting options as a bonus.

That used to be my view as well, but then you have the BCD floating point guys, with IBM actually managing to make it part of the new upcoming FP standard.

The standard benchmark in this arena is the telephone carrier billing program, which basically handles a huge amount of incoming ASCII records, while doing some simple math on them.

The canonical assumption is that the amount of calculation/record is so low that is is better to stay in ASCII/BCD the entire time, to save in/out conversion.

Fast BCD/ASCII to binary is of course trivial, on a SIMD machine you can even do multiple digits/iteration.

Having a comparably fast bin->BCD function pretty much takes away the official reason for the BCD math HW. :-)

Terje

- "almost all programming can be viewed as an exercise in caching"

I was able to find the *following* posting by you...

GeneticDr wrote:

[snip]
[snip]

The code above will (as noted) fail if the input value is >= 1e9, you need one more digit subtraction loop to handle those numbers as well.

Re. performance: Your loop will use at least 6 cycles/iteration/digit, so for a random number like 123,456,789 it would use at least 350-400 cycles. This is the same or more than a naive division-based algorithm:

; EAX has input, EDI -> output buffer mov ebx,10 ; Divisor mov edx,30303030h ; '0000' mov [edi],edx ; Prefill buffer with zeroes! mov [edi+4],edx mov [edi+8],dx add edi,9 ; Point at last digit position next: div ebx ; Remainder in EDX, answer in EAX (40 cycles) add dl,'0' mov [edi],dl dec edi test eax,eax jnz next done:

This version uses 43 cycles/significant digit, so for all numbers with

7 or less digits, it will be faster than your subtraction loop.

If you really want to subtract, then I suggest an unrolled binary search version:

mov ebx,01020400h next_digit: mov cl,0

REPT 4 mov esi,[ebp] ; ...,8000,4000,2000,1000,800,400,200,100,80,40, add ebp,4

add cl,bl ; Trial addition sub eax,esi ; Try to subtract table value

sbb edx,edx ; -1 if subtract resulted in carry

and esi,edx ; Subtracted value if carry, zero otherwise and bl,dl ; Ditto

sub cl,bl add eax,esi

shr ebx,8 ; Next digit ENDM

add cl,'0' mov ebx,01020408h

mov [edi],cl inc edi

dec ch jnz next_digit

This version should use about 250 cycles, but quite a bit of table space (4*10 dwords = 160 bytes.

The code I posted early this week uses no tables, and runs in 40-80 cycles, depending upon the cpu type. :-)

However, although the posting was dated either August 31, or September

1, 1997, and so this narrows down the date of the *target* posting, I wasn't able to find it.

John Savard

I can't be bothered to google for it, so here it is again:

The following function converts a 32-bit unsigned number into 10 ascii digits, with leading zeroes if the input is less than 1e9.

It uses 4 integer MULs, but no DIV/MOD operations at all, running time on anything (except P4) from a PentiumII to a Core 2 Quad is about 30-40 cycles. I have a C version of the same algorithm which is portable to any architecture which can do a 32x32->64 multiplication.

Terje

char *dtoa(char *buf, unsigned long n) { __asm { mov ebx,[n] mov eax,2814749767 mul ebx shr ebx,1 xor ecx,ecx mov edi,[buf] add eax,ebx adc ecx,edx mov eax,100000 shr ecx,16 // ECX = high part mov ebx,[n] imul eax,ecx // High part * 100k sub ebx,eax // Low part mov eax,429497 mul ecx mov ecx,eax add dl,'0' mov eax,429497 mov [edi],dl mul ebx mov ebx,eax add ecx,7 shr ecx,3 add dl,'0' mov [edi+5],dl add ebx,7 shr ebx,3 lea ecx,[ecx+ecx*4] mov edx,ecx and ecx,0fffffffh shr edx,28 lea ebx,[ebx+ebx*4] add dl,'0' mov eax,ebx shr eax,28 mov [edi+1],dl and ebx,0fffffffh add al,'0' mov [edi+6],al lea ecx,[ecx+ecx*4] lea ebx,[ebx+ebx*4] mov edx,ecx mov eax,ebx and ecx,07ffffffh shr edx,27 and ebx,07ffffffh shr eax,27 add dl,'0' add al,'0' mov [edi+2],dl mov [edi+7],al lea ecx,[ecx+ecx*4] lea ebx,[ebx+ebx*4] mov edx,ecx mov eax,ebx and ecx,03ffffffh shr edx,26 and ebx,03ffffffh shr eax,26 add dl,'0' add al,'0' mov [edi+3],dl mov [edi+8],al lea ecx,[ecx+ecx*4] shr ecx,25 lea ebx,[ebx+ebx*4] shr ebx,25 add cl,'0' add bl,'0' mov [edi+10],ah mov [edi+4],cl mov [edi+9],bl } return buf; }

- "almost all programming can be viewed as an exercise in caching"

Yes, since division is much slower, if you can turn your number into a fraction, and then multiply by 10 repeatedly, with a guarantee that the rounding error in the initial "turning the number into a fraction" operation is upwards only, but less than one place in the last bit, it should be a valid and faster algorithm.

When I was thinking of an algorithm that went for speed at all costs, I was thinking of something that multiplied by 100 repeatedly - and used one byte-wide decimal add operation (the x86 architecture having such a thing) for the last step of the conversion.

Basically, one has a seven-bit number from 0 to 99; the last three bits need no conversion, they have a binary value from 0 to 7, which is the same in packed BCD; the first four bits would be masked, and used to index into a table...

x'00', x'08', x'16', x'24', x'32', x'40', x'48', x'56', x'64', x'72', x'80', x'88', x'96'

which is only one *small* table.Perhaps this is too cumbersome to be faster than a single multiply instruction, however.

John Savard

Join the Discussion

Have something to add? Share your thoughts — no account required.

Didn't find your answer?

Ask the community — no account required