< mailto:info@emesystems.com
>
The BASIC Stamp division operator "/" returns the integer quotient, and the "//" operator returns the integer remainder of the division. Dividing a numerator N by a denominator D:
N = Q * D + R Q and R integers, 0<=R<D
with the quotient Q and remainder R. Written another way,
N/D = Q + R/D 0<=R/D<1
Q is what the BASIC Stamp returns from division "/", and R is what the BASIC Stamp returns for the "//" operator.
Suppose the problem is one of finding the ratio of two numbers that might be voltages that come from an analog to digital converter, or time intervals or counts that come from pulsin or count functions. You might end up needing to find the value of the fraction, say, 4567/6789, where in general both the numerator and the denominator are variables. On the stamp, the result you want will be, say, a 4 digit number 6727 representing a decimal fraction times 10000. 4567/6789=6727/10000, that is to say 0.6727. However, his problem is hard to solve on the Stamp. In integer division,
4567/6789=0 and 4567//6789= 4567
no help at all. There are some approximations, but let us cut right to the chase.
It is possible to divide two numbers on the BASIC Stamp using basically the same procedure you would use to perform long division by hand. Remember?! That is, divide the left-most digits of the numerator by the denominator to find the first digit(s) of the quotient, and a remainder. Then multiply the remainder times the base, "bring down" the next digit of the numerator, divide the denominator into that, and so on. The process may be carried out to as many places as is desired, to calculate both an integer quotient and as many places past the decimal point as desired.
Here we show how to do this both in base 10 and in base 2, and then give some examples. For use on the Stamp, I have found the base 2 approach to be the most useful.
In the first example, this is done in base 10. The 16 bit word size limits the size of the integer numerator and denominator to +/- 6553. At each step, the remainder from the previous step has to be multiplied times ten, and to keep that below 65535, the 16 bit word size, the denominator has to be less than 6553. Negative numbers are handled outside the main division loop using the usual rules. The main division loop operates only on a pair of positive numbers.
' fixed point divide routine, four places. decimal algorithm
N var word ' input: numerator ' -6554<N<+6554
D var word ' denominator ' -6554<D<+6554
X var word ' place holder
I var word ' integer part
J var nib ' index
F var word ' fractional part ' .0000 to .9999
main: ' demonstration program
N=355 ' demo numbers
D=113
gosub divide
debug sdec I,46,dec4 F
' show result
end
divide:
I=1-(N.bit15^D.bit15*2) ' sign of result
N=ABS N ' divide absolute values
D=ABS D
I=N/D*I ' integer part
X=1000 ' .1,.01,.001,.0001 place holder
F=0 ' fractional part
for J=1 to 4 ' 4 places
N=N//D*10 ' remainder * 10
F=N/D*X+F ' accumulate next digit
X=X/10 ' next place holder
next
return
This spits out a number that approximates the fractional part to 4 decimal places. That is,
N/D =~ I + F/10000
where "N" and "D" are the original numerator and denominator. The variable "I" holds the integer part of the quotient, and the variable F holds the fractional part. Note that this is an approximation. There is a remainder left over at each step, unless of course the decimal fraction is exact, like .5000, or .2500. But if it is a longer fraction, like .3333___, then the result is an approximation. It could be continued to 5 places to find a result like .33333, but that would not work for.66666, because that would exceed the size of a word, which would be .65535. So, in the above, the approximation is limited to 4 places.
This can be modified to spit out digits to a display or to a serial port. In that case the calculation of the fractional part can be carried out to more than 4 decimal places. There does not have to be an accumulation, and the variables X, I and F are not needed. You can just keep on going as long as you want. 1/3 comes out as 0.3333333333...
J=N.bit15 ^ D.bit15 ' sign of result, 1=negative
N=ABS N ' divide absolute values
D=ABS D
debug 13*J+32,DEC N/D,"." ' sign, integer part and decimal point
for J=0 to 15 ' 16 places
N=N//D*10 ' remainder * 10
debug DEC1 N/D ' display digit
next
The following implements the division algorithm, but using base 2 instead of base 10. The core algorithm, show below in red, is one of the most useful algorithms for doing math on the BASIC Stamp. Now the input range is better, because a remainder can be as high as ±32767. Here at each step, the remainder is multiplied times 2, and 2*32768<65535, the largest value that will fit in one word variable. It takes more steps to accomplish the division in base 2, because there are more binary digits.
' different version of same thing, uses
' binary instead of decimal.
' imporved input range -32768<N,D<+32767
' 2s complement fixed point divide
' Numerator/Denominator
' to four decimal places
' Integer.Fractional
' binary algorithm--good input range
' but slower than decimal algorithm
' better for implementation in machine language
N var word ' input: numerator ' -32767<N<+32767
D var word ' denominator ' -32767<D<+32767
I var word ' integer part
J var nib ' index
F var word ' fractional part .0000 to .9999
main: ' demonstration program
N=29870 ' demo numbers
D=3110
gosub divide
debug sdec I,46,dec4 F
' show result
end
divide2: ' subroutine
I=1-(N.bit15 ^ D.bit15*2) ' sign of result
N=abs N ' divide + numbers
D=abs D
I=N/D*I ' integer part
F=0 ' initialize
'---------binary division loop-----------
for J=15 to 0 ' 16 bits
N=N//D<<1 ' remainder*2
F.bit0(J)=N/D ' next bit
next
'----------------------------------------
F=F**10000 ' normalize F*10000/65536
return
The core of the division algorithm is highlighted in red. It is the part that comes up with a fraction F/65536 that approximates the fraction N/D:
N/D = I + F/65536
That is, the original fraction has been approximated by an integer part plus a fraction that has is a power of two in its denominator. In a decimal fraction, the approximation ends up with an implied power of ten in the denominator.
Another way to look at it is that the binary division loop that calculates F is executed sixteen times, and each time, the remainder is shifted left one binary place, that is, multiplied by two. Multiplying 2 times 2 times 2 sixteen times is 65536. Each step around the loop is computing a binary coefficient for that power of two, a polynomial that represents the number:
F = a15*215 + a14*214 + a13*213 + ... + a2*22 + a1*21 + a0*20
Each one of the coefficients is either zero or one.If all coefficients are one, the F=65535. If all coefficients are zero, then F=0. The value of F/65536 is the closest approximation to the value of N/D.
The final step, F=F**10000 renormalizes the fraction, to convert it from one with an implied denominator of 65536 to one with an implied denominator of 10000. That may require some explanation. In renormaization, the operation F=F**10000, effectively multiplies the fractional part F times 10000, and divides implicitly by 65536. (See the separate section on the operation **). The result is the desired numerator for a denominator is 10000, a proper decimal approximation. The stamp can then print out a decimal point followed by the four decimal digits.
The binary division loop by itself is sometimes very useful, by itself, as part of another algorithm. Just remember that the inner loop highlighted in red comes up with the approximation:
N/D =~ I + F/65536
The following example shows that it can be desireable to keep the binary approximation as part of a longer calculation.
Suppose you have a reading taken from a sensor, and there is a calibration constant that converts the raw reading into proper units. The conversion constant is applied by using the ** or */ operator:
Y = X ** calcon
Remember, ** means that you are approximating the conversion factor by a fraction with a denominator of 65536. (calcon/65536) Now suppose you want to implement an automatic calibration procedure, such that the user will put the sensor in a standard condition (pH, humidity, temperature, whatever), and then press a button to perform a calibration. Or you may have the user enter the calibration value at the keyboard, and the Stamp has to calculate the appropriate multiplier. The Stamp will read the value of X=Xc and then has to compute a new value of calcon and store it in memory. You know the calibration value is Y=Yc. We have to solve for calcon
calcon = Yc / Xc * 65536 This is normal math, not stamp math
The factor of 65536 comes from the implied division by 2^16 that comes along with the ** operator. This is another case where long division is essential. You probably can't solve this problem effectively on the Stamp without using long division.
So by applying the binary long division algorithm to Yc/Xc, it automatically comes up with calcon, including the factor of 65536, without any additional normalization.
For example, here is an example for a light level sensor. The formula for conversion of millivolts to units of watts per meter squared is:
WpM2 = millivolts * 0.231
which you can write on the stamp as
WpM2 = millivolts ** 15139 'because 0.231 = 15139/65536
We already talked about that in the article on multiplication by a constant fraction. The plot thickens when we need to have the user enter a new calibration factor from the keyboard, or in reference to a calibration standard. Of course we could simply provide the new calibration value as say, 15027. But usually it comes in a different form, like a new decimal multiplier, 0.229. The number to put on the right hand side of the ** is then calculated as:
multiplier/65536 = .231 as on a calculator
multiplier = 0.231 * 65536 as on a calculator
= 231/1000 * 65536 ditto
It is quite easy to pass the fraction N/D = 231/1000 to the core of the above algorithm (highlighted in red) and come out with the number F=15027. If you have understood the logic, you know that the factor 65536 is automatically.
Alternatively, suppose that a second light meter is used for calibration, and it reads a value W2. The first light meter puts out a certain millivolt reading, but the first light reading is incorrect because the existing calibration constant is wrong. The calibration constant must be adjusted bso that the answer comes out right:
W2 = millivolts ** newmultiplier ' as done on the stamp
newmultiplier = W2 / millivolts * 65536 As on a calculator
Again, it is easy to feed the values W2/millivolts = N/D to the core long division algorithm (binary version). Out pops the new multiplier.
Here's a discussion of how to calculate the reciprocal of a number using the limited Stamp integer math functions. This is 1/X, or more generally, C/X, where C is a constant, and X is a variable. You need this for things like
Note that the */ and ** operator are of no help here. They are only useful for multiplying a variable in the numerator times a fixed fraction; Y=X*(A/B). Now we want the variable X in the denominator.
In practical problems, the variable often appears in both the
numerator and in the denominator, e.g.,
Y = C*X/(X+K), where both C and K are constants. In problems
of this sort, first find C/(X+K), and then multiply by
X.
There are several ways to go about this, so the discussion will be in several parts.
shortcut #1, simple integer division into the largest possible numerator:
If you write a program with,
Y=1/X
on the stamp, you get Y=0 for all values of X except X=1. Instead, you have to choose a large value for C, such as
Y=10000/X.
It is a matter of scaling. If X represents a frequency measurement in hertz, then 1/X represents the period in seconds, and 10000/X represents the period in tenths of milliseconds, a larger integer. When you display the number, you put the decimal point in the correct place and call out the corresponding units. Say X=150 hertz. Then the period is 10000/X, or 66 tenths of a millisecond, which you would display as 0.0066 seconds, or 6.6 millisecond. (The true answer from the division is 6.666.. milliseconds. The stamp truncates the non-integer remainder.)
For the best resolution, use the largest possible numerator. You can't actually use 100000 in the numerator, because it is larger than the 65535 limit of 16-bith representation. But you can use 50000 as the numerator, and then multiply the result by two.
Y=50000 / X * 2 '(meaningful for X>1)
Again trying the example of X=150, the value of 50000/X is 333, and Y=666. This would be displayed as 0.00666 seconds, or 6.66 milliseconds, or 6660 microseconds. See, there is one more digit added by using the larger number in the numerator. But not a full digit really, because the final digit will always be an even number. It is 5 times better, not ten times.
The best results come at around X=224, which is close to the square root of 50000. At that point, the resolution of changes in the output equals the resolution of changes in the input.
X 50000/X Y
223 224 448
224 223 446
225 222 444
Notice that the resolution is around one part in 200. At small values of X, a single step in X gives a huge step in Y. For example, the step from X=2 to X=3 gives a step in Y from 50000 down to 33332. And at large values of X, X can change a lot before there is any change at all in Y. For example, the value of Y stays equal to 24 for all values of X from 3847 to 4166. So with the numerator=50000, and the denominator around 224, the resolution is about 1/2%.
shortcut #2, extended division:
One way to improve the accuracy is to make use of the remainder after the first division.
I=10000/X ' integer part from BS2 division
R=10000//X ' remainder from BS2 division
for example, if X=150, then the integer division is I=10000/150=66, and the remainder is R=10000//150=100. That is, written out as product plus remainder:
66*150+100 = 10000
You can get one more decimal place in the result by shifting the result of the first integer division over one place and dividing X into ten times the remainder.
Y = 10000/X*10 + (10000//X)*10/X (1)
For example, with X=150, the result is:
Y= (10000/150*10) + (10000//X*10/X) = (66*10)+(1000/150) = 666
This is the same answer as would come from dividing X into 100000. Here are some values of X around 316,which is the square root of 100000.
X 10000/X*10 10000//X*10 10000//X*10/X Y
315 310 2350 7 317
316 310 2040 6 316
317 310 1730 5 315
The result acts as if we had actually divided into 100000!
There are bounds on the allowable values of X. The conditions are:
10000/X*10 < 65536 ==> X > 1.
10000//X*10 < 65536 ==> X < 10001
So the synthetic division is good over the range of 1<X<10001.
We could have multiplied by 100 to get even more accuracy:
Y = 10000/X*100 + (10000//X)*100/X (2)
And here the answer would be
Y = (10000/150 * 100) + (10000//150 * 100 / 324) = 6600 + 100*100/150 = 6666
The constraints on the input range of X are more strict here:
10000/X*100 < 65536 ==> X > 10.
10000//X*100 < 65536 ==> X < 715
The effective numerator is one million, and there are four significant figures in the result, so long as X is within the bounds stated. So the limit on the input range is 10<X<715. The number 715 is not obvious--it comes from analysis of the modular equation. The upper limit is at least 655. But the upper limit can be greater than that, depending on the solution of the modular equation.
Synthetic division is not limited to simple powers of ten in the numerator. Here are a couple of numerical examples to show some other possibilities.
Example: find Y=23567/X, when X=234. The correct decimal answer is 100.714
Y = 23567/X*10 +(23567//X*10/X).
if X=234, then Y=1007.
This has one extra significant figure, compared to the simple integer division,
Y = 23567/234 = 100.
Example: find Y=235679/X, when X=567. Note that the numerator is greater than 65535. The correct decimal answer is 415.6596
Y = 23567/X*10 + (23567//X*10+9/X) ' note adding the ones digit, here=9
if X=567, then Y=415 ' answer to the ones place, remainder is truncated
or even better (the value, X=567 is in a range that allows this, because 567 * 100 < 65536):
The above method can be generalized to kick out more significant figures, an iterative process. You start with a dividend D and a divisor X, and you want to find a quotient Qn and remainders Rn. Initially break the value into the leftmost number of digits still less than 65536, and that was in the above example D=235679, D0=23567, with one more digit remaining, call it r0=9.Y = 23567/X*100 + (23567//567*100+90/X)
if X=567, then Y=4156
Q0 = D0 / X ' 23567 / 567 = 41
D1 = D0 // X * 10 + r0 ' 23567 // 567 * 10 + 9 = 3209
Q1 = D1 / X = 5
D2 = D1 // X * 10 + r1 ' 3209 // 567 * 10 + 0 = 3740 but in this example, no more digits r1, r2 ...
Q2 = D2 / X ' 3740 / 567 = 6
D3 = D2 // X * 10 + r2 ' 3740 // 567 * 10 + 0 = 3380
Q3 = D3/ X ' 3380 / 567 = 5
D4 = D3 // X * 10 + r3 ' 5450
Q4 = D4 / X ' 9
The best with simple integer division would be:
Y = 58919/(X/4) = 58919/141 = 417 ' where 235679 = 58919 *4, redistributing the 4, generally not so good
Note in the above that the final digit adds into the term after
the remainder from the first division is multiplied times the factor
of 10. This is just like long division, where you find the remainder
and then bring down the next digit of the dividend. The
iterative process is exactly what you would do in long division.
Example: You want to measure the velocity of a projectile that will be moving at somewhere between 20 and 80 meters per second (56 to 112 miles per hour). An arrangement of photocells ten centimeters apart triggers a flip flop to provide an input to the PULSIN command on the stamp. The pulse will be somewhere in the range of 1250 to 5000 microseconds. The number returned by the pulsin command (measuring in units of 2 microseconds) will be between 625 and 2500. With a simple integer division this is:
V = 50000/X
V=80 at X=625
V=56 at X=879 (the correct decimal answer is V=56.88)
V=20 at X=2500.
While with synthetic division the answer gains one significant figure:
V = 50000/X*10+(50000//X*10/X)
V=800 at X=625
V=568 at X=879 (the correct decimal answer is V=56.88)
V=200 at X=2500.
Example: You are using the RCtime command with an AD590 temperature sensor, and you need to calibrate the temperature given the value, X, returned by the RCtime command. You find that the RCtime returns a value of X=491 when the temperature is 273 Kelvin (ice bath). The calibration constant is 479*273=134043. The equation to find temperature is, K = 134043/X. Using simple integer division, we would have to lose resolution by dividing both sides by 3 to keep the quantities within the 2^16 bounds. But synthetic division can do better.
To do this on the stamp, you program,
K = 13404/X*10 + (13404//X*10+3/X)
K=297 at X=450
or, to get more resolution (because the expected magnitudes of RCtime are within bounds where tbis is possible <655)
K = 13404/X*100 + (13404//X*100+30/X)
K=2978 at X=450 (297.8 Kelvin, 24.7 degrees Celsius)
Although it should be noted that the RCtime is probably not accurate enough in its operation to justify the additional significant figure. But it could be rounded to 0.5.
Note that the synthetic division is simply a two-step case of long division. Insteead of a loop, the calculation is done in one line in two steps, using the quotient and remainder.
Shortcut #3: using 2^24 as the numerator.
Simple integer division is limited to a numerator of 65535. The number of significant figures can be improved by using a larger number in the numerator. Here is a formula that effectivley divides X into 2^24=16777216. The Basic Stamp cannot handle such numbers directly. This formula works in two steps. Here it is:
X var word ' 256<=X<=65535 input note bounds
Y var word ' 256<=Y<=65535 output
Z var word ' renormalized for numerator of 1000000
J var nib ' helper variable
J=ncd X - 8
Y=65535/X*256+(65535//X<<(8-J)/(X>>J))
Z=Y**10000*/100
debug dec X,32, dec Y,32, dec Z
I'll explain in a minute. Note the limits: both the input and the output fall in the range 256 to 65535. This is most useful for commands like PULSIN and RCtime that can easily be set up to return midrange values, with no chance of hitting the small ones.
The turning point for resolution here, the point where input resolution equals output resolution, is approximately 4096, the square root of the effective numerator of 2^24.
Here is an explanation of the formula:
Explanation 1: This is like the process of synthetic division
described above. The numerator is
2^24-1=16777215. The division is done in base 256, rather
than in base 10 or base 100:
Y = 65535/X*256+(65535//X*256+255/X) <- BUT THIS WON'T WORK!
The trouble is that the remainder, the term 65535//X, can in general be larger than 256. And that times 256 gives integer overflow. So we need another way to interpret the remainder term.
The remainder, 65535//X, is less than or equal to X, so (65535//X/X) is a fraction that lies between zero and 1, and 256 times that fraction is a number that lies between zero and 256. This final term that lies between 0 and 256 is a "correction" to the first term, and it is what gives us the added precision in the result.
Another way to look at it is that we need to scale both the numerator and the denominator in the remainder term in the above equation by some factor that will keep the numerator them it less than 65535.
If X is between 256 and 512, then divide both the remainder and the denominator by 2; if X is between 512 and 1024 then divide by 4; and so on. This brings in the NCD function. The factor to divide by is J=(ncd X - 8). In the formula, the divides are shown as shifts. Divide by 2^J is the same as shift right by J. The <<8 is the same as multiply by 256
J=ncd X - 8
Y=65535/X*256+(65535//X<<(8-J)/(X>>J))
'-----------...............remainder/ power-of-2 factor
'----...........times 256,
'------ then divided by (X/same factor)
Z=Y**10000*/100
'--..........................result from above step
'------....................times 10000/65536
'-----...............times 100/256
The last normalization step in effect multiplies the result Y times the constant fraction in two steps,
1000000/16777216 = 10000/65536 * 100/256.
See the BS2 manual for details about the ** and */ commands.
(This from a question raised on the Stamps listserver.) Using an external circuit, a BS2 pulsin command measures the time in microseconds*2 taken for a vehicle to travel 6 inches. The expected limits of speed should be between 3 to 120 mph. How do you calculate mph from the PULSIN count using integer math?
The conversion factor from inches per second to miles per hour is 0.0568182. Putting it together:
speedMPH = 6 inches * 0.0568182 * 10E6 / (2 * time) = 170454 / time
The expected pulsin time ranges from 56818 down to 1420 (in units of 2 microseconds), to give 3 to 120 mph output.
This formula works fine on a calculator/computer capable of floating point math, but it needs special treatmentto get good results using the Stamp's integer math.
The easiest solution is to scale down the numerator to make it less than 65535. Then scale down time by the same factor to maintain the ratio. For example, divide both numerator and denominator by 5.
'speedMPH = (170454/5)/(time/5) = 34090/(time/5):
time var word
speed var word
measure:
pulsin 1,1,time 'X-OR high-time
time=time/5 ' time in units of 10 microseconds
speed=34090/time
debug "speed = ",DEC speed, "mph" 'print speed etc...
This approach gives 1 mph resolution. Maybe 1 mph resolution is good enough. But information is lost. If you need to split hairs, here is one way to go.
time var word ' 1420<time<56818
speed var word J var nib 'helper variable
measure:
pulsin 1,1,time 'X-OR high-time
' keep time in units of 2 microseconds
J=ncd time - 8
speed=65535/time*256+(65535//time<<(8-J)/(time>>J))*/260
debug "speed =",DEC speed/100,".",DEC2 speed," mph" 'print speed etc...
This approach gives you one or two extra digits to the right of the decimal point.
This program uses 2^24=16777216 as the numerator for the division. Observe that the normalization step using ** 10000 and */100 that was shown in the previous discussion is not used here. Here, instead, there is multiplication by a constant fraction (*/260) to get mph*100. Why */260? Note that 16777216 * 260/256 = 170454, which is exactly the number that has to go in the numerator, (speedMPH = 170454/time).
Another approach to this would be synthetic division.
speed = 17045/time*10+(17045//time*10+4/time)
This would give you one extra digit of resolution, compared to the simple integer division, 34090/(time/5).
< mailto:info@emesystems.com
>