Cube Root

[L D[AR2P#0.0]]

I wonder if the OP has followed this thread and has it answered his original question ?

 

[bernie_carlton]

I think we have diverted it into strange areas. The discussion of optimum calculation of the cube root completely side stepped the necessity of the OP for a cube root at all. Then Peter's introduction of Kahan's 'magic number' really got me going. I can't let a 'magic' thing just sit there.

 

[TConnolly]

I suspected that rmb550 probably was on the right track in answering the OP's question with a question.

I know what you mean in being unable to let the magic number just be. I've poked at it, looked at in in hex and binary, overlayed it on an IEEE float, etc. I even went back to my C source code library from my earlier programming days and looked at the source for the pow and exp function (exp uses magic numbers too). Whatever bit twiddle Kahan was up to is still a mystery to me.

 

[JesperMP]

Try 15625.

Result of cube root of 15625:
Regular way: 25.0
Kahans way w. divide: 25.57285
Kahans way w. multiply: 25.57279

 

[Peter Nachtwey]

Terry Wood called this a tangent

I think we have diverted it into strange areas. The discussion of optimum calculation of the cube root completely side stepped the necessity of the OP for a cube root at all. Then Peter's introduction of Kahan's 'magic number' really got me going. I can't let a 'magic' thing just sit there.

Actually I brought up the magic number but L D [AR2,P#0.0] found the source of the info. L D [AR2,p#0.0] brought up using iterations. L D [AR2,#p0.0] just needed a good initial estimator.

Bernie, I don't know how long you stared at the numbers but your question showed insight. No one else asked that question or made that observation. I spent plenty of time staring at the numbers years ago because I had to work this out for the non standard DSP.

Result of cube root of 15625:
Regular way: 25.0
Kahans way w. divide: 25.57285
Kahans way w. multiply: 25.57279

Kahan's method is just an estimation. You should still do iterations as L D [AR2,P#0.0] or I have shown to get the exact value. Combining Kahan's cuberoot estimation and the iteration formula that L D [AR2,P#0.0] showed should result in an exact 32 bit answer in two iterations.

I have other methods of estimating the cube root. They involve simplifying the problem by scaling it do a range of 1 to 8, computing the mantissa and then scaling back. There are polynomials that can be very accurate at computing functions over a small span like from 1 to 8. I hope you guys understand that all the combinations of mantissas will be found from 1 to 8. The rest is just dividing the exponent by 3.

Alaric, did you try finding the minimum worst case magic number? All you need to do is check the numbers from 1 to 8.

I would use a logarithmic scale

for ( i=0 ; i<= 3000 ; i+= 1 )
{
     a=2^(i*0.001);   // a goes from 1 to 8
// put the estimation and comparison code here.

}

I think this is enough for now.

 

[JesperMP]

So, here is the STL source code for a CBR, with Kahans method for 1st guess, and a variable number of subsequent iterations.

FUNCTION "CBR" : VOID
TITLE =
VERSION : 0.1


VAR_INPUT
  rX : REAL ;	
  iNoOfIterations : INT ;	
END_VAR
VAR_OUTPUT
  rXcuberoot : REAL ;	
END_VAR
VAR_TEMP
  rQbrApprox : REAL ;	
  rInv3rX : REAL ;	
  rQbrApprox_cubed : REAL ;	
  iLoopCount : INT ;	
END_VAR
BEGIN
NETWORK
TITLE =

      L     #rX; // Kahans method. With a multiply in stead of a divide.
      DTR   ; // convert REAL to REAL (!? - dont think, just do it)
      L     3.333333e-001; // load 0.33333333333
      *R    ; // multiply by 0.333333333 - same as divide by 3
      RND   ; // convert back to DINT (but it is a REAL, just close your eyes and do it)
      L     L#709921077; // load magic number
      +D    ; // add as DINT
      T     #rQbrApprox; // rQbrKahanApprox should be an approximated (rX)^1/3 (as REAL

      L     3.333333e-001; 
      L     #rX; 
      /R    ; 
      T     #rInv3rX; 


      L     #iNoOfIterations; // to preserve cycle time, no test for counts less than 1 !
NEXT: T     #iLoopCount; 

      L     #rQbrApprox; 
      L     #rQbrApprox; 
      *R    ; 
      L     #rQbrApprox; 
      *R    ; 
      T     #rQbrApprox_cubed; 

      L     #rX; 
      L     4.000000e+000; 
      *R    ; 
      L     #rQbrApprox_cubed; 
      -R    ; 
      L     #rInv3rX; 
      *R    ; 
      L     #rQbrApprox; 
      *R    ; 
      T     #rQbrApprox; 

      L     #iLoopCount; 
      LOOP  NEXT; 

      L     #rQbrApprox; 
      T     #rXcuberoot; 
END_FUNCTION

With rX = 15625
If iNoOfIterations = 1 --> rXcuberoot = 24.97335
If iNoOfIterations = 2 --> rXcuberoot = 24.99994
If iNoOfIterations = 3 --> rXcuberoot = 25.0

I wonder how the cycle time will be on a real CPU.
All the above instructions, including several times through the loop, have to be compared against two loads, one *R, one EXP, one LN, and a transfer.

 

[JesperMP]

You might want to read this:

I have estimated the execution time for the two methods, "simple exp ln" and "Kahans 1st guess + iterations". I dont have a CPU at hand, and PLCSIM is not representative how it will run on a real CPU.
I have calculated the execution times for several CPUs.
Namely 315 and 317 of the previous generation, and 315 and 317 of the current generation.
For simplicitys sake, I have assumed that 2 iterations is enough, as it resulted in very little error in the previous test.

And the result is:

For the previous generation CPU's, the "Kahans 1st guess + iterations" method is significantly faster. Namely 76 µs vs. 681 µs (315) and 18 µs vs. 50 µs (317).

For the current generation CPU's, the "Kahans 1st guess + iterations" method is somewhat slower. Namely 14 µs vs. 7 us (315) and 4.2 µs vs. 2.7 µs (317).

 

[Peter Nachtwey]

That is interesting

There must be a lot of overhead interpreting the STL code relative to the execution time on the newer CPUs. This would penalize methods that use a lot of code vs and exp or ln function that gets compiled.

The other possibility is that the new S7s have built the exp and ln functions into their new silicon like intel has. In otherwords the exp and ln functions are CPU instructions instructions and not funcions that are themseleve written in STL.

You need to try a number besides 15625.
If you just try numbers from 1 to 8 you will get a fair idea.