Lost Precision With Integer Math. 😕

lardom

4,294,967,295 / 65535 = 65537. There is no remainder. Excellent.
4,294,967,295 * sine 0.7071 = 3,036,971,374.29. This is the return I want.
Using integers, 4,294,967,295 / 65535 = 65537 * 7071 = 463,412,127 / 10,000 = 46,341.2127 Move the decimal four places to the right and then four places to the left.
46,341 * 65535 = 3,036,975,435.
3,036,971,374.29 - 3,036,975,435 = -4,060.71
I want that 0.2127 for precision. I will store the final result in a variable array but I want to know if there is a way to calculate both the left and right side of the decimal point in a single process. Am I missing something?

Duane Degn

I don't have time to look at your problem in detail but one method (TtaMethod) I use regularly was made by Tracy Allen (always a huge help). I added the "MultDivSigned" method to allow negative numbers.

I'm pretty sure TTA is Tracy's initials.

PUB MultDivSigned(N, X, localD) | sign  ' return X*N/D

  sign := 1
  
  if N < 0
    -sign
    -N
  if X < 0
    -sign
    -X
  if localD < 0
    -sign
    -localD
    
  result := TtaMethod(N, X, localD) * sign
   
PUB TtaMethod(N, X, localD)   ' return X*N/D where all numbers and result are positive =<2^31
  return (N / localD * X) + (binNormal(N//localD, localD, 31) ** (X*2))

PUB BinNormal (y, x, b) : f                  ' calculate f = y/x * 2^b
' b is number of bits
' enter with y,x: {x > y, x < 2^31, y <= 2^31}
' exit with f: f/(2^b) =<  y/x =< (f+1) / (2^b)
' that is, f / 2^b is the closest appoximation to the original fraction for that b.
  repeat b
    y <<= 1
    f <<= 1
    if y => x    '
      y -= x
      f++
  if y << 1 => x    ' Round off. In some cases better without.
      f++
  

JonnyMac

The P1 doesn't support unsigned integers directly (the P2 does), so you'll have to write some code to deal with that.

Tracy Allen

Almost like on the BASIC Stamp, you can use the ** operator to maintain precisions when multiplying times a constant factor.

Say the factor is sin(45°) = 0.7071. Pre-multiply that times 2^31 = 2147483648. So 2147483648 * 0.7071 = 1518485687
That is the factor to use with **. As follows:

2147483647 ** 1518485687 * 2 * 2 = 3036971372
That is close to the answer you wanted, which was 3,036,971,374.

The extra factors of 2 come from the business that the P1 treats longs as signed integers, and the trick with ** only works with positive numbers, that is, less than 2^31. But the ** operator normalizes to twice that, 2^32.

If 0.7071 is supposed to be sin(45°) = 0.707106781186548, then the factor is 0.707106781186548 * 2^31 = 1518500250, and the ** estimate is,
2147483647 ** 1518500250 * 2 * 2 = 3037000496
whereas 2^32 * sin(45°) = 3037000499.9760.

lardom

@"Tracy Allen" I was just studying and commenting out different lines of "Sinewave0" with my calculator. I took a closer look at this line just before I decided to login to the forum:

frqb := $8000_0000 - (Fullsine(phsa>>19) << 15)

Magnificent.
I'm going to set aside the thought that I have to design an algorithm based on $FFFF_FFFF and instead concentrate on $8000_0000. I still don't understand what "<<15" does but I realize I first have to figure out how it's connected to "FRQA := 54".
Thank you for explaining "**". I absolutely had no idea how that worked.
I have read your posts for years. Thank you.

Tracy Allen

frqb := $8000_0000 for a cog counter in DUTY mode causes the outb pin to toggle at clkfreq/2, that is 50% high and 50% low. So it comes out at 1/2 scale when filtered by an RC circuit. The rest of the expression for frqb has the effect of modulating up and down from that halfway point, modulating the density of highs vs lows. The sine table returns 16 bit values. Shifting left by 16 bits would give 100% modulation, from zero to 100% density, however, the pulse spacing is quite long near the extremes. Shifting by 15 bits gives modulation from 25% to 75%, and good pulse density throughout. That is why there is <<15 in the above expression. You can observe the effects on depth on your oscilloscope by trying shifts other than 15 bits.

As to the frqa := 54, that is the number that is added to the phase accumulator, phsa, on each clock cycle. The phase accumulator has a length of 2^32 = 4294967296. Notice that 4294967296 / 54 = 79536431.407407407407407...
No coincidence, it is so close to 80_000_000, which is the clkfreq in the demo. Adding 54 on each clock cycle takes very close to 1 second to get all the way through the phase accumulator, and that means 1Hz output from the NCO mode counter, and it also means 1Hz sine output from the coupled DUTY mode DAC. It's actually a tiny bit off from 1 Hz, by about 0.6%.

Tracy Allen

About the method BinNormal(...) that Duane brought up. That is what you can use to calculate the factor to go with ** as a precision fractional multiplier. The way I use that mostly is for calibration. Say you have a 24 bit weighing scale, and you have a standard 1kg weight for calibration. The user presses the "calibrate" button when the weight is on the scale. Say the Prop reads 1003479 from the scale, and the display needs to read 1000.00. The muliplier for subsequent readings is 100000/103479. Plugging those numbers into the BinNormal method for 31 bits effectively solves this equation:
factor / 2^31 = 100000/103479.
and it finds that factor = 2075284501. (rounding up). It stores that factor in DAT eeprom.
Now readings from the scale can subsequently be corrected:
correctedReading := reading * 2 ** 2075284501. (the factor of 2 comes from using 31 bits instead of 32 bits, to avoid negative territory)

Note that you can't do it simply, with correctedReading := reading * 100000 /103479, because that might overflow 32 bits or go negative.

lardom

@"Tracy Allen" You posted this in 2014. I plugged in some numbers and something clicked:💡

'Tracy Allen multiply high
OBJ
   pst : "parallax serial terminal"

VAR
 long Y,X
    
PUB tryDoublePrecision
  pst.start(115200)
  waitcnt(clkfreq/10+cnt)
    
  repeat
    pst.str(string("enter x: ")) ' can be positive or negative
    x := pst.decIn    
    pst.str(string("enter y: "))   
    y := pst.decIn    
    pst.str(string("answer z= $ "))
    pst.hex(x ** y, 8) ' high 32 bits, msb is now the sign
    pst.hex(x * y, 8)  ' append low 32 bits, shown as hex
    pst.newline
    pst.dec(x)
    pst.newline
    pst.dec(y)
    pst.newline

Hexidecimal and binary values displayed as expected! That meant all I had to do was multiply the raw data in the sine table by 65,535. Signed integers were actually a non issue.
@E801 in the sine table is 46340 * 65535 = 3,036,981,900
3,036,981,900 / 4,294,967,295 = 0.70710245071. Time to celebrate.
You also wrote something that made Phil Pilgrim's "UMath" easier to use, I think. I'm trying to find that post again.

lardom

Just to demonstrate how much ground I have to cover to truly understand Tracy Allen's "SineWave0", I'm posting the results of my attempt "Sine Table Practice 8", for comparison. My version does not drop to zero using the same RC filter for both versions. 😕 I'll figure out why. I just wanted to show where I am.