High-precision reciprocal

Lord Steve

Could anyone provide PASM code to calculate 2^32 divided by x?

mpark

Check out post #13 in this thread: http://forums.parallax.com/showthread.php?123663-What-are-the-standard-clkfreq-values-above-80MHz

It's Spin, not PASM, but the translation should be straightforward.

Lord Steve

PERFECT! Thanks, mpark!

Tracy Allen

Here is that algorithm in both spin and pasm.

Michael, you had asked in that other thread how Kye and I came up with the algorithms. In my case I needed it for the PBASIC on the Stamp, where with 16 bits words it is even more critical. I consider it a fundamental integer algorithm, right up there with multiplication and division. From basic principles, after you do a division, there is a quotient and a remainder. The remainder then is what is to the right of the radix point. This algorithm then calculates the binary fraction to whatever number of binary digits is required. It is just like ordinary long division. At each step it asks if the numerator is greater than the denominator. If no, it inserts a zero into the lsb, if yes, it inserts a 1 and subtracts to get a new numerator. Then it doubles the numerator and repeats.

In the Prop, this is how a program calculates from a desired frequency in Hertz to the angular frequency required by frqx. But it is also often the best way to handle any integer calculation involving fractions. On the Stamp with 16 bits it is more important than on the Prop with 32 bits, but it is still important for highest levels of precision. The main caveat on the Prop is to watch out for the sign bit.

spin:

PRI binRatio(y, x, bits) : f            
    repeat bits                   
        y <<= 1
        f <<= 1
        if y => x    '  
            y -= x
            f++
{
 calculates y/x * 2^bits
 enter with {y < x,    x < 2^30,    x <= 2^30,    bits =< 32}
 exit:  f / 2^bits is the closest binary denominated appoximation less than y/x
 caveat:  take care for the sign bit (see entry conditions)
 }

pasm:

' enter with a and b: a < b, a < 2^30, b <= 2^30
' exit with f: f/2^32 =~ a/b  (enter with nbits := 32)
' if different number of bits, then f: f/(2^n) = a/b

ratioxybits            ' enter here with a, b and n
       mov f,#0    ' frequency parameter
:loop
       shl a,#1              ' double the value of a
       shl f,#1               ' double the estimate of f, appoximation f/(2^n) =~ a/b
       cmpsub a,b    wc, wr       ' if a>=b then a:=a-b, carry:=1, else a unchanged, carry:=0
       addx f,#0             ' add carry to f
       djnz n,#:loop    '  all bits
ratioxybits_ret      ret

lonesock

(Quick note, Tracy: you could speed up the PASM version a bit by using "rcl f, #1" after the cmpsub, instead of a shl and addx. And, of course, if N = 32 then there is no need to zero f at the beginning.)

sorry, back to your regularly scheduled programming [8^)

Jonathan

mpark

Thanks for the additional info, Tracy!

Tracy Allen

Good tip Jonathan. Would this work? The carry needs to be zero going in, and
and wc on the initial mov should do that I think.

' enter with a and b: a < b, a < 2^30, b <= 2^30
' exit with f: f/2^32 =~ a/b  (when entered with n = 32, number of bits)
' if different number of bits, then f: f/(2^n) = a/b

ratioxybits            ' enter here with a, b and n
       mov f,#0  wc  ' frequency parameter
:loop
       shl a,#1              ' double the value of a
       rcl f,#1               ' double the estimate of f, approximation f/(2^n) =~ a/b
       cmpsub a,b    wc, wr       ' if a>=b then a:=a-b, carry:=1, else a unchanged, carry:=0
       djnz n,#:loop    '  all bits
ratioxybits_ret      ret

lonesock

Just move the rcl to right after the cmpsub, that way the C flag is always set right before it is "rotated" into the f flag.

Jonathan

Tracy Allen

I see now it should be in that position after the cmpsub anyway. Otherwise the final bit will always be zero.

' enter with a and b: a < b, a < 2^30, b <= 2^30
' exit with f: f/2^32 =~ a/b  (when entered with n = 32, number of bits)
' if different number of bits, then f: f/(2^n) = a/b
ratioxybits            ' enter here with a, b and n
          mov f,#0  wc  ' frequency parameter
:loop
          shl a,#1              ' double the value of a
          cmpsub a,b    wc, wr       ' if a>=b then a:=a-b, carry:=1, else a unchanged, carry:=0
          rcl f,#1               ' next bit & double estimate of f, approximation f/(2^n) =~ a/b
          djnz n,#:loop    '  all bits
ratioxybits_ret          ret

The indentation of the spin version above was off -- I went back up and corrected it. The f++ is inside the IF ...

Tracy Allen

Someone asked me a couple of questions in PM, puzzled by 1) how 2^32 divided by a number could be considered a reciprocal and 2) so what if it gives a binary fraction, how do you get a decimal fraction (for display)?

Say you have 1 / 7531 which is your reciprocal. The algorithm kicks out 32 binary digits:
result = %00000000000010001011001111000001
And that value is stored as an integer in a long. The binary radix is implied at the left end. The decimal value of that integer is 570305. What this means is that the original number 1/7531 is now approximated by a number with a much larger denominator, 2^32. Due to the fact that the denominator is large, the approximation is good. ("good" can be easily quantified) You can check for yourself that 570305 / 4_294_967_296 is very close to 1 / 7531

To transform that to a decimal fraction, use the Prop ** operator.

decimalResult := result ** 100_000_000

That gives decimalResult = 13278. What that has done is to renormalize the denominator from 2^32, to 10^8.
Now the approximation is 1 / 7531 = 13278 / 100_000_000.
To print that, place the decimal point 8 positons to the left: 1 / 7531 = 0.00013278.

I hope that might a help to other spinners.