Stamp II math note #1

multiplying a variable times a constant fraction

(c) 1998 EME Systems, Berkeley CA U.S.A.
Tracy Allen
<stamp index> <home>

Contents (updated 12/14/2005)


*/ and ** operators for multiplying a variable times a fraction, overview

top

These are probably two of the most useful math operators available on the BASIC Stamp.   These operators allow the BASIC Stamp to do fractions, breaking the barrier of integer math.  These operators will come up time and time again as we come up with solutions to practical calculation problems using the BASIC Stamp and as we move on to "higher" math with the BASIC Stamp

To set the stage:

Suppose you have an approximation to PI=3.1416 to four decimal places or  PI=3.14159265359 to 11 decimal places. The first is an approximation using a denominator of 10000 and the second a denominator of 100,000,000,000
 PI ~= 31416 / 10000 = 3.1416
 PI ~=   314159265359 / 100000000000 = 3.14159265359

Recall from high school math that you can choose other denominators, not necessarily a power of 10.   And 22/7 is a good approximation with a denominator of 7.  
 PI ~= 22 / 7 = 3.142857142857
or
 PI ~= 355 / 113 = 3.141592920354

The last one gets it right to  six decimal places with a relatively small denominator.  The point here is that PI can be approximated by using any number in the denominator, and in general a larger number in the denominator allows you to make better approximations to a given number.   Computers work in powers of two, so it is convenient to think about how the computer makes the approximation.  

Think about making approximations using a denominators  256 or 65536, which happen to be the size of the microcomputer byte or word variables. The principal is exactly the same. The best approximations to PI with 65536 as the denominators is,
 PI ~= 205888 / 65536 = 3 + 9280 / 65536 = 3.1416015625
which is 8.90891e-06 too high, or
 PI ~= 205887/65536 = 3 + 9279/65536 = 3.141586303711
which is 6.349879e-06 too low, therefore a little bit better.  Not as good an approximation as 355/113, but still good to several decimal places.

With the denominator of 256 we can't do quite as well:
PI ~= 804 / 256=3 + 36 / 256 = 3.140625

The approximation with 65536 as the denominator is better than with 256 as the denominator, for the same reason as 31416/100000 is better than 314/100.

Here is how this relates to the BASIC Stamp.   The BS2 */ and ** operators multiply times a constant and implicitly divide by a power of two:

Y = X */ Z1   ' X is a variable, Z1 is a constant
' effectively forming the fraction Z1/256
' the Stamp implicitly performs the division by 256

Y = X ** Z2 ' X is a variable, Z2 is a constant
' effectively forming the fraction Z2/65536
' the Stamp implicitly performs the division by 65536

The way the BS2 manual explains the */ command is technical:

The result of the */ operator, Y *? Z, is the middle 16 bits of the 32-bit product of the two 16-bit input values, X and Z

It helps to recall (or to learn!) that dropping the 8 least significant bits of a binary number is equivalent to dividing the number by 28=256. That is precisely what the */ operator does. It implicitly performs the division by 256

The ** operator is similar, in that it takes the top 16 bits of the product. It drops the 16 least significant bits of the 32 bit product, in effect performing an implicit division by 216=65536. In the Stamp manual, the ** operator is explained in the context of double precision math. That is one use of it. This fractional multiplication is another. If you read through these math pages, ** turns up in all sorts of surprising places.



Cookbook example, how to compute Y=X*1.2207:

For example, how to multiply an input 12-bit count (0--4095) by the constant 1.2207? Multiply 1.2207 x 256 = 312.5. Round off. That's the value that goes on the right:

Y = X */ 313 

For example,

Y = 3890 */ 313 

The BS2 result is Y=4756. Compare that to the exact result, 1.2207*3890=4748.52. If you round down and use 312 instead of 313, the BS2 result is 4740, a bit closer to the correct result, but only good to 8 bits of accuracy.

The accuracy is improved by multiplying the */ factor times 10. That is, use 1.2207 x 256 x 10 = 3125 as the factor on the right of the */:

Y = 3890 */ 3125

The BS2 result is now Y = 47485. Again, compare that to the exact result, 1.2207*3890=4748.52. You can keep the extra precision, or, divide by 10 to round the result off to Y=4748.

But there are some caveats. You have to check that the intermediate product won't overflow 24 bits. Otherwise the "middle 16" rule will chop off the most significant bits of the result, giving rubbish. The maximum input in the example was stated to be X = 2^12 = 4096, from an analog to digital converter, so the */ factor also has to be less than 4096, in order to keep the product less than 2^24. That is true in this case. The */ factor is 3125.

To achieve the best precision, it is always best to use the largest possible */ factor , consistent with no integer overflow. You have to understand the limits of the variables. For example, a 12 bit analog to digital converter is limited to values from 0 to 4095, so the factor on the right can safely be as large as 4096 too, without danger of integer overflow.

 

Here is the same thing using the ** operator. We need to multiply the variable X times the constant, 1.2207. Multiply 0.2207 * 65536 = 14464. That is the factor that goes on the right of the ** operator.

Y = X ** 14464 + X

We have to add an additional "X" at the end because we really want to multiply times 1.2207, not 0.2207. Remember, the ** operator only does fractions less than one. The result is more accurate than was the */ result, for example, when X=3890:

Y = 3890 ** 14464 + 3890

This gives Y=4748. Compare that to the exact result, 1.2207*3890=4748.52. Moreover, if the value of X is limited to X<4096 (from an analog to digital converter), we can coax out one more significant figure, by multiplying X times 10 going into the formula:

Y = (X*10) ** 14464 + (X*10)
' gives Y=47485 when X=3890

This provides enough precision for accurate rounding off the the result, based on the final digit:

Y = (X*10) ** 14464 + (X*10) + 5 / 10
' gives Y=4749 when X=3890

See how it is now rounding up correctly? All the integer math operators round down. So you have to carry forward an extra digit of precision if you want to round up and down in the standard manner with 0.5 in the middle.

One last note. How would you multiply times 1.2207 on the BASIC Stamp without using the */ or ** operators? It is difficult, because the fraction is not a ratio of small whole numbers. The fraction is 12207/10000. You can look at it as:

1.2207 = 1+(2/10)+(2/100)+(7/10000)

And the formula Y=X*1.2207 can be translated to BS2 code as:

Y = X + (2*X/10) + (2*X/100) + (7*X/10000)

Which comes up with the answer Y=4747 when X=3890. Compare to the correct answer, Y=4748.52. There are several steps in the computation where integer round-off can take its toll on the accuracy. The */ or ** operators do really provide the best approximation available for the BASIC Stamp, with much less computation.


Say you have a 12-bit convertor and the reference is trimmed to 5120 millivolts, 1.25 millivolts per bit. The converter output has to be multiplied times 1.25 to get from the raw count to millivolts.

Here is the obvious way to do it:

Y = X*5/4   ' when X=4095, X*5=20475
  1. Y=1*5/4 gives Y=1 millivolt compared to correct value 1.25
  2. Y=2*5/4 gives Y=2 millivolt compared to correct value 2.5
  3. Y=3*5/4 gives Y=3 millivolt compared to correct value 3.75
  4. Y=4*5/4 gives Y=5 millivolt which is the correct value

Although, it is possible to improve the rounding (so it rounds off both up and down as appropriate).

Y = X*5+2/4   ' 2 added for rounding
  1. Y=1*5+2/4 gives Y=1 millivolt <--- 1.25 is rounded down to 1
  2. Y=2*5+2/4 gives Y=3 millivolt <--- 2.5 is rounded up to 3
  3. Y=3*5+2/4 gives Y=4 millivolt <--- 3.75 is rounded up to 4
  4. Y=4*5+2/4 gives Y=5 millivolt <--- 5 is exact

In the above, you can't multiply the X by ten initially, because that would give integer overflow when X was at its upper limit of 4095 for this A/D converter (4095*10*5>65535).

Here is how to get one more significant digit by using */.

Y = X */ 3200 
  1. count gives 12 (millivolts*10), compared to 1.25 millivolts exact
  2. counts give 25 (millivolts*10), compared to 2.50 millivolts exact
  3. counts give 37 (millivolts*10), compared to 3.75 millivolts exact
  4. counts give 50 (millivolts*10), compared to 5.00 millivolts exact

The maximum value in the multiplier will be 4095*3200, which satisfies the requirement that that product be less than 224.

It is not possible to get another significant figure. When X=4095 the result will be Y=51187, which is already close to the integer limit, 65535. The exact value is 4095*1.25=5118.75, but the maximum precision in one 16 bit word is 51187 (note that ** always rounds down).

What if you have to have that last significant figure for a display? Note that the last figure is a 5 when X is odd:

Y = X */ 3200
debug dec Y/10,".",dec1 Y," millivolts"
'displays 0, 1.2, 2.5, 3.8, .... , 5118.8

Z = X // 2 * 5
debug dec Y/10,".",dec1 Y,dec1 Z," millivolts"
'displays 0, 1.25, 2.50, 3.75, .... , 5118.75

Y = Y + X//2 ' rounding up and down
debug dec Y/10,".",dec1 Y," millivolts"
'displays 0, 1.3, 2.5, 3.8, .... , 5118.7

See, there are all sorts of tricks you can play with numbers!


Another example, using RCtime command to determine an unknown capacitance:

You want to use the RCTIME command to determine the value of an unknown capacitance. You have a known resistance of 10kohms. From RCTIME you come up with a number "X" in units of 2 microseconds. Then you have to solve for C,

where "K" is a proportionality factor. This is an opportunity to use the BS2 "*/" command. It is the BSII's way of multiplying a variable "X" times a constant fraction. Both the resistance and "K" are fixed constants. From the BS2 version 1.8 manual, the constant K has a value of around 600. The manual describes how to find out the value of K for your exact setup. Suppose the value is exactly K=1/600.

To use the */ command, you now multiply the exact value of K/R times 256 and move the decimal point to the left as much as possible (limited by integer overflow) and round off to an integer:

Now in your BS2 program, after you get the value "X" from the RCTIME command, you calculate and display capacitance as follows:

cap = X */ 427 
debug DEC cap/10000,".",DEC4 cap

About 3 significant figures--4 is stretching it. Note that the debug command first sends out the integer part, then the decimal point, and then the fractional part.

Now, unlike the case with the A/D converter, there is danger of integer overflow. The value of X can take on any value from 0 to 39291. Why 39291? Because integer overflow occurs if 427 * X > 16777216, which is 2 to the 24th power.

The RCTIME command can return values from 0 to 65535. The value 39291 corresponds to a maximum capacitance of C=6.5 microfarads, measurable in this circuit. Of course, the best approach is to choose a fixed value of R that is most appropriate for the value of capacitance being measured.

advanced technique, rationalize:

Here is another tip for "advanced stampers" . The exact value of the multiplier in this example is 1.6666, which is equal to 5/3. Multiply by 256 and you get 426.666667. For this example, we used the approximate value, 427 on the right side of the */. There is a round-off error. In this situation, where the multiplier is a simple fraction, you can get the best results by making the multiplier "exact". Multiply the factor by three: 3*426.6666667 = 1280. Then divide by 3 afterwards.

cap = X */ 1280 / 3   ' to 4 decimal places, or 
cap = X */ 128 / 3 ' to 3 decimal places

In general this approach gives the best results. Rationalize the multiplier for the */, then divide by a constant after. In the final formula, using 128, there is no integer overflow, even for X=65535.


Another example, Fahrenheit to Celsius conversion:

The point here is to illustrate the tradeoffs of doing the calculations in one way or the other. Reference calculation:

C = F * 5 / 9 - 17.7778

Suppose that F is measured in units of 1/10 of a degree F. I.e. as the LM34 temperature sensor produces 10 millivolts per degree F, and the analog to digital converter resolves to 1.00 millivolts per bit over the range of 0-4.096 volts. Suppose that F will get as high as 300 degrees Fahrenheit (148 degrees C). That is 3000 millivolts, or 3000 counts of the a/d converter. You want to convert this to Celsius and keep the resolution as good as possible, better than 1/10 degree C.

Here are different ways you could implement this on the stamp (assume C and mV are word variables). Observe that division does not work on negative numbers, so for temperatures below freezing at least, the subtraction has to follow the division.

   A.  C = mV * 5 / 9 - 178     ' gives 0.1 degree C
B. C = mV * 50 / 9 - 1778 ' 0.01 degree works up to only 55 C!
C. C = mV */ 142 - 178 ' 0.1 degree C
D. C = mV */ 1422 - 1778 ' 0.01 degree works up to 637 C!
E. C = mV ** 36409 - 178 ' to 0.1 degree C
F. C = mV*10**36409-1778 ' to 0.01 degree C, up to 346 C!

I already said that the temperature would be limited to the range of 0 to 300 degrees F, or the voltage signal from 0 to 3000 mV. But first let's consider the possible input range for each formula, limited by integer overflow of the multiplication. The LM34 sensor won't operate at some of these temperatures, but a thermocouple could.

  1. -17.8 up to 710 degrees C {limited by mV*5<65536}
  2. -17.78 up to only 55 degrees C {limited by mV*50<65536}
  3. -17.8 up to 3623 degrees C. {limited by mV*142<16777216}
  4. -17.78 up to 637 degrees C.{limited by mV*1422<16777216}
  5. -17.8 up to 3623 degrees C {limited by mV<65536}
  6. -17.78 up to 346.3 degrees C {limited by mV*10<65536}

Now look at the accuracy at 100 degrees F: The correct Celsius conversion is 37.78 degrees.

  deg F    mV    deg C     A     B        C     D      E      F
99.9 999 37.72 377 3772 376 3771 377 3772
100.0 1000 37.78 377 3777 376 3776 377 3777
100.1 1001 37.83 378 3783 377 3782 378 3783

Column B has greater precision, but less range. Columns D and F have both good range and good precision. The message is that if you want both high accuracy and a wide range, use the */ or ** forms. If the range and accuracy requirements are not so demanding, then the straight multiply and divide will work just fine.

Another note describes how to display negative values with the decimal point.


Comparison of the */ and the ** operators for fractional multiply

Summary

*/ -- A variable, call it X, is to be multiplied times a constant fraction, A/B. The approximation used for */ is:

    X * A/B = X * F/256  = X */ F
--on a calculator-- ~BS2~~

** -- A variable, call it X, is to be multiplied times a constant fraction, A/B. The approximation used for ** is:

    X * A/B = X * F/65536  = X ** F
--on a calculator-- ~BS2~~


Example: A well tempered musical scale

The octave is the ratio of 2:1 between two sound frequencies. Other tones on the pentatonic scale fall at simple ratios: 5/3 for the sixth, 2:3 for the fifth, 3:4 for the fourth, and 4:5 for the third.

In the well-tempered scale, the notes of the chromatic scale are all powers of the number Z=1.059463. That is because the twelfth power of that number, Z12 is equal to 2, the ratio for the octave. The number Z is the twelfth root of two. If you start at one frequency, say at the high A: 1760 hertz, and multiply that times 1.059463 twelve times, then you come out one octave higher, at 3520 hertz. Or if you divide 1760 twelve times by 1.059463 you come out at 880 hertz, the note A one octave lower in pitch. Here is the ascending chromatic scale of 12 notes:

.
1.00000 Z0 the root note
1.059463 Z1 ...
1.122462 Z2 ...
1.189287 Z3 ...
1.259921 Z4 the ascending major third ~5/4
1.334839 Z5 the ascending major fourth ~4/3
1.414213 Z6 ...
1.498306 Z7 the ascending major fifth ~3/2
1.587400 Z8 ...
1.681791 Z9 the ascending major sixth ~5/3
1.781796 Z10 ...
1.887747 Z11 ...
2.000000 Z12 the octave 2/1

The scale is "well-tempered" because it relies on the predictable and repeatable powers of the twelfth root of two. To the keen ear the notes may sound a little sharp or flat compared to the "perfect" ratios. What does the ear hear as most pleasing? To the intellect of a musician like JS Bach there must have been something very attractive about the mathematical clarity of the well-tempered scale.

Here is the descending scale is produced by the powers of 0.943874 = 1/1.059463. The twelfth power of 0.943874 is 0.5, one octave down. The simple descending minor third, fourth etc. come close to the fractions 4/5, 3/4, 2/3 and 3/5.

.
1.00000 Z0 the root note
0.943874 Z1 ...
0.890898 Z2 ...
0.840896 Z3 ...
0.793699 Z4 the descending minor third ~4/5
0.749152 Z5 the descending minor fourth ~3/4
0.707105 Z6 ...
0.667418 Z7 the descending minor fifth ~2/3
0.629959 Z8 ...
0.594602 Z9 the descending minor sixth ~3/5
0.561229 Z10 ...
0.529723 Z11 ...
0.500000 Z12 the octave 1/2

One way to produce a scale on the stamp is by means of a lookup table. Another is by use of the ** operator to approximate the ratios 1.059463 or 0.943874. Let's cut to the chase and write down a program to generate a chromatic scale. The best approximation to 0.943874 is the ratio, 61858/65536. And the best approximation to 1.059463 is 1+ (3898/65536).

i    var  byte
x var word

top:
x=1760 ' start on A-1760, down two octaves to A-440

for i=1 to 24 ' descending scale
freqout 0,500,x
x=x ** 61858 ' see note below
next
' now back up to A-1760
for i=1 to 24 ' ascending scale
freqout 0,500,x
x=x ** 3898 + x ' see note below
next

goto top ' do it again

This should play a scale two octaves down through A 880 to A 400, and then back up two octaves to the starting note. Almost. There is a small problem. Due to cumulative round-off error, the scale does not go down exactly two octaves. Empirically, the values 61881 and 3940 work better over that two octaves than do the calculated values, 61858 and 3898. The best value depends on the range of notes. Errors can accumulate in a sequence of calculations. It is better to set up your calculation so that it can all be done in just one or a few steps. The above scale has no problem if you set the initial note to x=1760 each time around the loop. The other errors are too small for most ears to hear.


Example: The 1% RN55 resistor series


Here is a similar program, this time to compute the most significant figures of the 1% resistors.  The ratio in this case is 1.024275221382, the 96th root of 10.   Multiply that number by itself 96 times and you get 10.    Round off those numbers to three significant digits, and you get the series of values of 1% resistors.   Here are the values from 100 ohms to 1000 ohms as you would find them in an  electronics parts catalog.  Other decades have the same significant figures, only with the decimal point moved to the left or to the right.  


100 178 316 562
102 182 324 576
105 187 332 590
107 191 340 604
110 196 348 619
113 200 357 634
115 205 365 649
118 210 374 665
121 215 383 681
124 221 392 698
127 226 402 715
130 232 412 732
133 237 422 750
137 243 432 768
140 249 442 787
143 255 453 806
147 261 464 825
150 267 475 845
154 274 487 866
158 280 499 887
162 287 511 909
165 294 523 931
169 301 536 953
174 309 549 976
1000


It is also possible to use the inverse, 0.976300098962 = 1 /  1.024275221382.  Multiply that number by itself 96 times and you get 0.1.   Round those off and you get the same series as above.    Even though the above series of numbers contains only 3 digits, the calculation must be carried out at much higher precision, because the rounding off is often quite close in splitting hairs between rounding up versus rounding down.

In the Stamp, the factor to use with the ** operator is 63983.    That is, 63983/65536 = 0.976303100586, which is close to the correct factor.   But is is precise enough?   Can it be **'d by itself 96 times to generate the above series of numbers?   Almost.   First of all, notice that the program starts with wx=65000,  and then interatively calculates wx = wx** 64983.   The rounding off is done by  dividing by 65, so the final result goes from 1000 down to 100.  The reason for starting at the highest possible number word variable is to acheive the highest possible precision.  The ** operator on the Stamp rounds down after each operation, and the final result after 96 iterations comes out a little too low.    The program below includes a couple of small systematic corrections are the result of empirical tinkering.   The loop index, idx, goes from 1 to 97.   The least significant of idx is used to apply a correction to the central calculation, making it wx = wx ** 63983 + idx.bit0 instead of simply wx = wx ** 63983.  That compensates for the Stamp's tendency to round down everything in integer math, by adding one to every other calculation.   The final result after 96 iterations is wx~6500 as it should be.    The other small adjustment is in the roundoff formula, which also uses the idx.bit0 to alternate the roundoff criterion.

  
'{$STAMP BS2pe}
'{$PBASIC 2.5}
root96 CON 63983 ' 65536/[10^(1/96)], nintysixth root of ten
wx VAR WORD ' for high precision calculation
wy VAR WORD ' resistance values rounded off to 3 digits
idx VAR WORD ' for-next index
idx0 VAR idx.BIT0 ' for trick, empirical for best result

wx=65000
FOR idx=1 TO 97 ' 96 values plus one in the next decade 100-->1000
wy=wx + 32 - idx0 / 65
DEBUG DEC 97-idx, TAB,DEC wx,TAB
DEBUG DEC wy,CR
wx = wx ** root96 + idx0
NEXT
END


If you want to study these values in an Excel spreadsheet,  here are formulae to enter in the columns, and extend it down to 97 cells.
k	10^(k/96)			1% series	
0	=POWER(10,A1/96)*100		=INT(B1+0.5)
1	=POWER(10,A2/96)*100		=INT(B2+0.5)

Note especially numbers like  280.4895377 or 231.5172384, which round off to 280 ohms or 232 ohms respectively, but the precision has to be high to do it properly, instead of 281 or 231 ohms.  It is of course possible and probably more practical to implement the resistor sequence using a lookup table.  But what fun would that be?


<top> <index> <home> logo < mailto:info@emesystems.com>