Large number division

[kenjanzen]

I am trying to divide a number by a variable and am having trouble getting the accuracy that I would like. I have searched the forum and most of the answers are beyond my understanding.

I need the following:

W4 = 7,500,000 / W0 ‘ Where 1,500 < W0 < 9,000

I came up with the following code, but any help would be appreciated:

W1 = 50000 / W0 * 150
W2 = 50000 // W0
IF W2 < 436 THEN aa
IF W2 < 655 THEN bb
IF W2 < 873 THEN cc
IF W2 < 1310 THEN dd
IF W2 < 2184 THEN ee
IF W2 < 2621 THEN ff
IF W2 < 4369 THEN gg
IF W2 < 6553 THEN hh
IF W2 < 13106 THEN ii
IF W2 < 32767 THEN jj
GOTO zz

aa:
W3 = W2 * 150 / W0
GOTO zz
bb:
W3 = W2 * 100 / W0 * 3 / 2
GOTO zz
cc:
W3 = W2 * 75 / W0 * 2
GOTO zz
dd:
W3 = W2 * 50 / W0 * 3
GOTO zz
ee:
W3 = W2 * 30 / W0 * 5
GOTO zz
ff:
W3 = W2 * 25 / W0 * 6
GOTO zz
gg:
W3 = W2 * 15 / W0 * 10
GOTO zz
hh:
W3 = W2 * 10 / W0 * 15
GOTO zz
ii:
W3 = W2 * 5 / W0 * 30
GOTO zz
jj:
W3 = W2 * 2 / W0 * 75


zz:
W4 = W1 + W3



Jeremy Leach’s answer in “Help caculating RPM with Pulsin command on PICAXE-8M” is very close to what I need but it does not work for me.

“Drag Strip Timer” is also similar to what I need, but I still need help.

 

[moxhamj]

How much accuracy do you need?

Maximum number in W0 is 65535. So divide 7,500,000 by 65535 and this gives 114.44. Round up to 115. Divide both sides of the equation by 115. So 7,500,000 becomes 65217. The unknown number is between 1500 and 9000. Test with 9000. Divide by 115 = 78.26 which a picaxe will round to 78. Now divide 65217/78 = 836.

The real number is 7,500,000/9000=833.3 recurring. Divide by 836 to get the error which is 99.68%

Is that ok?

So the code is

w1=w0/115
w4=65217/w1

 

[kenjanzen]

The accuracy is still not good enough for what I need. At 1,609 (the lowest accuracy with in the required range)

output is:
65217/(1609/115)=5016

Actual is:
7500000/1609=4661

4661/5016*100=92.9%

I did a bunch of creative IF THEN commands and have the accuracy that I need, but I think I will look at adding a math coprocessor.

Thanks for your help

 

[premelec]

Not sure just what you are doing... parenthasis does not work so do the division and then the multiplication - 65000/1619*115 [I've forgotten the numbers!] ... sorry if I've misunderstood...

 

[hippy]

I did a bunch of creative IF THEN commands and have the accuracy that I need, but I think I will look at adding a math coprocessor.

That seems to be over-kill and unnecessary cost for a problem which has already been solved in PICAXE software. It would be worthwhile studying the RPM calculating and dragstrip timer code; even though you don't currently understand it you will if you put your mind to it and people here can help with that.

If necessary, someone could probably wrap up the code in a format whereby it can simply be dropped into your own program.

 

[Dippy]

Yes, agreed.
Out of interest, what level of accuracy is acceptable? >98% ? >99.9% ?

 

[westaust55]

Better accuracy but still around 3.6% at worse case

Within the accuracy of the PICAXE integer maths, the best that I can suggest is:

W4 = 63382 / (W0/114)
Which equates to

W4 = 63382 / W0 * 114

you are always going to have the same result for approximately 100 steps so 1500 to 1595 will give the same value and 1596 to 1709 give the same value.

Over the total of 7501 steps,
- the maximum error across the range will be 3.675%
- the average error across the range will be 1.740%
- 111 values from the 7501 steps will be greater than 1.000% error

By comparison, using the factor of 115 as previously suggested while theoretically correct gives error of:
- the maximum error across the range will be 7.610%
- the average error across the range will be 3.211%
so the “WA modified” formula roughly halves the error.

 

[hippy]

By emulating 32-bit math the error will be 0% ( ignoring any fractional parts ). It's just that the code's a little (!) more complex.

 

[westaust55]

32-bit emulation for large numbers

Spurred on by Hippy&#8217;s comment about 32-bit emulation and zero error with a desire to achieve the best results, I had a look at some often quoted past threads which under the titles Drag Strip Timer (in Archive 9) and Help caculating RPM with PULSIN command on PICAXE-8M ( inArchive 19).

The code from the above latter reference thread can be expressed as:

; 32 bit maths emulation for large numbers

SYMBOL lsw   = w0
SYMBOL msw = w1
SYMBOL x     = w2
SYMBOL y     = w3
SYMBOL z     = w4
SYMBOL ans = w5

Init: x = 7500     ; numerator broken into x and y parts
       y = 1000     ; x * y = 7,500,000
       z = 1500     ; divisor  a value from 115 to 65335

Main:

;multiply:  product = x * y
       lsw = x * y      ; low word of product
       msw = x ** y     ; high word of product

;divide ans = product / z

       ans = lsw / z
       ans = - 1 // z * msw/z + ans          ; standard program line
;
;      ans = - 1 // z/5 * msw/z*5 + ans   ; tweaked program line
;
       ans = - 1 / z * msw + ans

END

This was supposed to have zero error but my findings are that there can be considerable error as per the following:

_Divisor____Correct___PICAXE_____Error___PICAXE_____Error
_(= z)_____Answer___Mod code___ (%)____Std code_____(%)
__115______65217_____62494____4.176____65216_____0.002
__150______50000_____49996____0.008____49998_____0.004
__500______15000_____14996____0.027____14998_____0.013
_1500_______5000______4996____0.080_____4955_____0.900
_1609_______4661______4657____0.092_____4578_____1.786
_6000_______1250______1194____4.480_____1150_____8.000
_7500_______1000_______955____4.500______920_____8.000
_9000________833_______831____0.280______803_____3.640
65000________115_______115____1.200______114_____1.200
65535________114_______114____0.368______114_____0.368


Even with tweaking, at 4.5% error (that is the max within the samples I tested and could be worse for some other values) the error is greater than the one off calcultion I proposed earlier.

In conclusion:
- The 32-bit emulation code above is a little more generic/versatile, but
- More complex (as Hipy mentioned)
- Will require more memory space
- Still needs some tweaking and (see the slight change in one line of the code sample above) to try and improve the results.
- Could use more variable space in a general purpose version
- Better results may still be achieved by a &#8220;purpose&#8221; create line of code with tuned constants
- The divisor has a lower limit/value ( ie must be greater than numerator / 65535)

 

[Mycroft2152]

Here's a link to 32 bit division from the PIClist:

http://www.piclist.com/techref/microchip/math/div/div16or32by16to16.htm

Not sure if this will help.

Myc

 

[hippy]

That code is from Jeremy and I must confess to not understanding it as I've not studied it. I think my own PICAXE version used repeated subtraction so was slow but I have better understood 'the proper way to do 32-bit division' ( shifts and subtracts ) so that's much faster and I have that working on a 32-bit processor.

Converting that code to work with a 16-bit PICAXE is the challenge and my mind's on other things at present. The algorithm is as follows ...

lhs := lhs / rhs, where lhs and rhs are 32-bit, acc can be 16-bit or 32-bit, loops is 8-bit

LhsBecomesLhsDividedByRhs:
  Gosub LhsBecomesLhsModuloRhs
  lhs = acc
  Return

LhsBecomesLhsModuloRhs:
  acc   = 0
  loops = 1
  Do While rhs < $80000000
    rhs = rhs * 2
    loops = loops + 1
  Loop
  Do
    acc = acc * 2
    If lhs >= rhs Then
      lhs = lhs - rhs
      acc = acc | 1
    End If
    rhs = rhs / 2
    loops = loops - 1
  Loop Until loops = 0
  Return

 

[sjremington]

Along the iterative line suggested by Hippy, a good place to start would be with W4 = 63382 / W0 * 114

Then increase or decrease W4 using binary chopping to get the closest result.

7,500,000 = $7270E0

Thus:

W4**W0 must give $72 (114 decimal)
then
W4*W0 should be as close as possible to $70E0=28896

You also get the remainder of $7270E0/W0 essentially for free.

Jim

 

[moxhamj]

Hippy, that long division code looks very similar to a Z80 machine code version I looked up last night. Lots of shifts, subtracts and jumps based on carrys. It was 11pm and I must say machine code is one of those things you forget over the years (which is why I like picaxes!), but a generic 32bit by 16bit divide routine would be a useful addition to the code snippets section.

 

[hippy]

Yes, I think that's the 'standard' algorithm, or at least a commonly used one. If someone wants to work out what the PICAXE code would be where lhs and rhs are each formed by two word variables the rest is fairly easy ...

If lhs >= rhs Then

lhs = lhs - rhs

 

[kenjanzen]

I thank you all for your help with this. In part of this project what I am attempting to do is measure the speed of a bullet. We have two chronograph screens (they sense when the bullet passes through them) set up 15 feet apart. I am using the pulsin command on a 28X1 running at 20 MHz. When the bullet passes through the first screen the signal is set high and when the bullet passes through the second screen the signal goes low. The pulsin is used to time how long it takes the bullet to travel the 15 feet between the two screens. I am then displaying the speed in feet/second.

The formula that calculates this is:

Distance between the screens divided by travel time.

The distance is 15 feet

The pulsin command of the 28X1 running at 20 MHz will measure the travel time in 2us units

The formula is then:

15 feet / W0 * 2 us= feet/us

7.5 feet / W0 us * 1000000 us/seconds = feet/second

Or:

7,500,000 / W0

After doing some more research the following is what I have come up with. I’m on vacation and don’t have a picaxe with me to try it, but I believe that it should work. I believe it follows what Hippy mentioned about subtracts and shifts. And yes I believe that it will be slow. But it is the beginning of an understanding for me.

INTERRUPT:
PULSIN 1,1,T_unit ' Measure time between Chronos

L_wrd = 28897 ' Low word of 7,500,000
FOR H_wrd = 0 TO 114 ' High word of 7,500,000
A1:
L_wrd = L_wrd - T_unit
Speed = Speed + 1
IF L_wrd >= T_unit THEN A1
L_wrd = L_wrd - T_unit ' L_wrd will overflow (negative)
L_wrd = L_wrd + 65535 ' L_wrd will overflow (back positive)
Speed = Speed + 1
NEXT ' Carry one from High word

“Not” Everything worth doing is worth doing well.