4 byte numerator divided by fixed denominator

ttessier

Hello,

I have been going through Tracy Allen's excellent description of how to do double precision math on the Basic Stamp ( http://www.emesystems.com/BS2math6.htm·). Working from the division section ("Method 1: divisor is a constant<65536"), I was able to get division by 1000 working. Unfortunately, my implementation·seems to produce an·answer which was sometimes one higher in value than what I would get in the Windows hex calculator.

For example, dividing the test value 0x88884444 by 1000 I get 0x0022 FC38, Vs. 0x0022 FC37 on the Windows calculator in hex (I am putting·the spaces between each·16 bit number to clearly delineate the words).

When I divide the number 0x8888444488884444 by 1000 I get 0x0022 FC38 BF9F E15A, Vs. 0x0022 FC37 BF9F E159 on the Windows calculator in hex. It seems that the lowest nibble is sometimes one higher than it should be on every second byte. In my application, I would not care if it was just the lowest value character in·24 bits (meaning an error of 1/16,777,216) but it is also happening on nibble 0 of byte 2, so it's actually an error of 1/40,000 after the first run through this value. I actually need to divide by 10^6, which means running the value through the "divide by 1000" algorithm twice thus amplifying the error (dividing by 10,000 will cause the values to explode beyond·8 bits, so I am using a divide-by-1000 twice method).

I'm not sure what I can do to increase the accuracy unless I try going to one of the more complex division algorithms. Maybe some sort of round-off detection test is needed? I thought I would post to the forum first, though, to see if anyone had a suggestion on what I have already done before I embark on a new tack for this problem.·In the code below,·I check the answer by looking at the debug window output, so q0 is overwritten in this test code with each new answer to save on variable code space.

Here is my test code, dividing a 4 byte number 8888444488884444 by 1000:
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' {$STAMP BS2sx}
' {$PBASIC 2.5}
z1 VAR Word· ' two words z1:z0 to make double precision input to start
z0 VAR Word
q1 VAR Word· ' output, quotient, z/divisor
q0 VAR Word
rm VAR Word· ' remainder from the division, z//divisor

z1 = $8888
z0 = $4444
' calculate z/1000, were z=z1:z0, two words
q1=z1/1000
q0=z1//1000
rm=((q0**536*536//1000)+(q0*536//1000)+(z0//1000))
q0=(q0*65)+(q0**35127)+(z0/1000)+(rm/1000)
rm=rm//1000

DEBUG "q1 ", HEX4 q1,CR· 'getting value 0x0022 hex, first byte of answer
DEBUG "q0 ", HEX4 q0,CR· 'getting value 0xF3C8 hex, second byte of answer
DEBUG "rm ", HEX4 rm,CR· 'remainder 0x02EC, passed to next level of division

'''''''''''''continue to next term $8888 hex
'make rm the equivalent of q0, with the old z0 now implied to be the
'new z1, already //1000 from rm in the last section
q0 = rm 'make the transition to next byte - rm is now q0
z0 = $8888
rm=((q0**536*536//1000)+(q0*536//1000)+(z0//1000))
q0=(q0*65)+(q0**35127)+(z0/1000)+(rm/1000)
rm=rm//1000

DEBUG "q0 ", HEX4 q0,CR· 'getting value 0xBF9F hex, third byte of answer
DEBUG "rm ", HEX4 rm,CR· 'remainder 0x0370, passed to next level of division

'''''''''''''continue to next term $4444 hex
'make rm the equivalent of q0, with the old z0 now implied to be the
'new z1, already //1000 from rm in the last section
q0 = rm 'make the transition to next byte - rm is now q0
z0 = $4444
rm=((q0**536*536//1000)+(q0*536//1000)+(z0//1000))
q0=(q0*65)+(q0**35127)+(z0/1000)+(rm/1000)
rm=rm//1000

DEBUG "q0 ", HEX4 q0,CR· 'getting value 0xE15A hex, fourth and final byte of answer
DEBUG "rm ", HEX4 rm,CR· 'final remainder 0x009C
END

Tracy Allen

Hi

I think you may have discovered a typo in the formula for division by 1000, in the placement of parentheses. I think after looking at it again that the formula for remainder rm should be,

rm=((q0**536*536//1000)+(q0*536)//1000+(z0//1000)) ' <-- correct
instead of
rm=((q0**536*536//1000)+(q0*536//1000)+(z0//1000)) ' <--- not

That is, the //1000 is the modulus of the entire term involving q0, not just the modulus of the second term in q0. I'll look at it more closely, but I think that is what is going on. As I mentioned on that web page (thanks for looking!), the formulae for fixed divisors greater than 256 have to be analysed on a case by case basis. Fixed divisors up to 257 are allowed without exception.

You mentioned the number 0x8888444488884444. That is 64 bits, quad precision. Or is it "just" 0x88884444

Another way to divide by a million would be to first divide out all six factors of 2 with a six bit shift right, and then do two successive divides by 125 for the six factors of 5.

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Tracy Allen
www.emesystems.com

ttessier

Thanks so much for your reply! I made the change, and now it works perfectly!

Given the "order of operations" of stamp math, I would not have thought the position of the parentheses would have changed the outcome, but it did! Now I am dividing 1000 into a 4 byte number and getting all the accuracy I get from the Windows calculator in Hex! Not bad for a little 8-bit processor, eh?

As for the quad precision - yes, I am definitely using that, and now getting the same answer as in the hex calculator. I substituted the remainder rm into the next step of the calculation...the code in the original post has comments that should allow following the logic. If·my terse explanation is too·fuzzy to follow, let me know and I'll expand on it. If anything I wrote is useful, please feel free to add it·to your web page. I found your explanations excellent and helpful - you really did a lot of work putting those explanations together and I would not have made it this far without them.

Thanks for the suggestion on the right shift. I'll look at that too and see what I could save in terms of memory space and clock cycles.

Cheers,

-Tom Tessier