Cube Root

[mordred]

bit of a side note one trick I learned in Computer programming is when you wish to multiply and devide by 2 a bit shift left is the same as devide by two and BSR obviously is the same as multiply by two. This was handy to know when speed was a concern before the advent of the math co-processor as the PC's before then would multiply by adding x number of times. Of course this is referring to Binary

 

[Bobbias]

Yes, that's true. In fact, most computer use a system called Shift and Add. Basically, it combines bit shifting with addition to achieve the multiplication, and can be faster than actual multiplication between registers, on occasion.

 

[Peter Nachtwey]

This will be a huuuh? moment.

Here's some results for the two algorithms (DB3=Peter's
Algorithm). I used unity as the initial guess.

You can see that L D[AR2,p#0.0]'s method converges with fewer iterations but this doesn't tell us which one is faster. How long does a divide take relative to a multiply? I count clock cycles.

Also, the trick is getting the initial guess right, then my method takes 5 iterations. This is the key.

L D[AR2,P#0.0], you appear to be taking the cube root of 15625 to get 25. First load 15625 into ACC1 and convert to a float if it isn't already.
Next multiply 0.333333333333333
Then +D 709921077 but treat the result a REAL.
This is your first guess.
Yes, you are adding adding a DINT to a REAL and the result is still a REAL but don't let that bother you.
Now execute the rest of the code and post it as you did before. You will see that both methods will take fewer iterations.

So why does the *0.333333333+L#16#709921077 work?

The S7 should allow you do this trick because STL is close to machine code whereas other PLCs will not let one add a DINT to a REAL.

Manglemender, the code you posted is how XPY is computed.

I have spent a lot of time looking at the source code to C libraries, modifying them and then recompiling to make the routines faster.

 

[Bobbias]

You can see that L D[AR2,p#0.0]'s method converges with fewer iterations but this doesn't tell us which one is faster. How long does a divide take relative to a multiply? I count clock cycles.

Also, the trick is getting the initial guess right, then my method takes 5 iterations. This is the key.

L D[AR2,P#0.0], you appear to be taking the cube root of 15625 to get 25. First load 15625 into ACC1 and convert to a float if it isn't already.
Next multiply 0.333333333333333
Then +D 709921077 but treat the result a REAL.
This is your first guess.
Yes, you are adding adding a DINT to a REAL and the result is still a REAL but don't let that bother you.
Now execute the rest of the code and post it as you did before. You will see that both methods will take fewer iterations.

So why does the *0.333333333+L#16#709921077 work?

The S7 should allow you do this trick because STL is close to machine code whereas other PLCs will not let one add a DINT to a REAL.

Manglemender, the code you posted is how XPY is computed.

I have spent a lot of time looking at the source code to C libraries, modifying them and then recompiling to make the routines faster.

That is some impressive stuff. I remember one of my friends figured out how to optimize some assembly code for converting all the letters to capital (or small) without using a single jump, lol. They did that because the book they were reading on optimization said something along the lines of "now this code, you'd be hard pressed to improve any more than it is" so they took on the challenge, and succeeded.

Out of curiocity how much optimization is applied to PLC code when it's fed through the compiler (I guess this is a platform dependant question, but what's more common, optimizing the code, or not)?

 

[TConnolly]

Not much - partly because most PLC programs are not compiled to machine code but to op codes that are interpreted by the PLC firmware. The interpreter may do some optimization, but I suspect that most don't because they must be highly generalized.

 

[Peter Nachtwey]

Bobbias, did you try it?

Alaric, whether the code is compiled down to pseudo code or machine code there many small optimizations that can be made.
r**3 can still be converted to r*r*r
divides by powers of 2 can still be converted to shift rights. Divides by a constant can be converted to multiplies by 1/constant. MOD(number,2**n) can be converted to AND 2**n-1.

Back on topic.
So why does the *0.333333333+L#16#709921077 work?
The multiplying by 0.3333333333 should be obvious adding L#16#709921077 is not.

 

[Bobbias]

Bobbias, did you try it?

There's a copuple reason I haven't. A) I don't actually know what I'd be doing. I've never worked with ST berfore, B) Don't have a PLC available to play with at the moment, nor does this coomputer have one of the emulators.

I'm a co-op student, and nearly everything I know about PLCs comes from this class, personal research (I'm on here as often as possible) and from the tiny bit we've touched on PLCs in my program. It's pathetic, but my program doesn't actually have a PLC course until after your second (out of 3) CO-OPs (mind you, even though this is technically my second co-op, I didn't find a placement last time, so it doesn't count. This is my first actual co-op.)

From a theory standpoint, I'm assuming that adding a DINT to a REAL would just result in a bitwise addition, regardless of the format... Just like how in C you could say, add an INT to a CHAR, or even add an INT to a BOOL if you felt like it (among various other tricks).

 

[JesperMP]

I am trying your suggestion Peter.

I have this code:

      L MD6              // REAL number to be "qube rooted"
      L 3.333333e-001    // load 0.3333333333333
      *R                 // multiply
      L L#709921077      // load magic number
      +D                 // add as DINT (?)
      T  MD10            // this should be the 1st guess (as REAL)

But result is:

MD     6  FLOATING_POINT     15625.0 
MD     6  DEC                L#1182016512 
MD    10  FLOATING_POINT     1.505932e+029
MD    10  DEC                L#1878215647

I gues the 1st guess should be somewhere near 25.

edit: I spotted that I had put in a *D in stead of a +D. But it didnt improve much.
I hav corrected it in th code and the result above.

 

[Peter Nachtwey]

There's a copuple reason I haven't. A) I don't actually know what I'd be doing. I've never worked with ST berfore, B) Don't have a PLC available to play with at the moment, nor does this coomputer have one of the emulators.

You don't need that. I doubt there are any PLCs besides the S7 that can use this trick. Get a free C or C++ compiler.

I'm a co-op student, and nearly everything I know about PLCs comes from this class, personal research (I'm on here as often as possible) and from the tiny bit we've touched on PLCs in my program.

Don't get hung up on the PLC part. It is just a tool. PLCs are robust when it comes to I/O but most are crippled when it comes to computational capabilities.

It's pathetic, but my program doesn't actually have a PLC course until after your second (out of 3) CO-OPs (mind you, even though this is technically my second co-op, I didn't find a placement last time, so it doesn't count. This is my first actual co-op.)

That may be a good thing. Take a programming course using C or C++ first so your eyes are wide open and can see what can be done. I have never taken a PLC course. Therefore I have never been brain damaged in thinking in limited terms imposed by PLCs.

From a theory standpoint, I'm assuming that adding a DINT to a REAL would just result in a bitwise addition, regardless of the format... Just like how in C you could say, add an INT to a CHAR, or even add an INT to a BOOL if you felt like it (among various other tricks).

In C you can type cast variables to convert them from one type to another. In C you convert a dword to a real or float like this (float)dwVariable. This is similar to the very wordy IEC conversions like DWORD_TO_REAL. Usually C will convert the variables for you it depends on the initial and desired type.

 

[JesperMP]

Btw. I am impressed with the accuracy of the iterative methods. They seem to converge quite accurately to 25.0.

If you do

      L     MD6
      LN    
      L     0.33333333333
      *R    
      EXP   
      T     MD14

You get 25.4, not 25.0.

 

[Bobbias]

Well, as far as I remember, you didnt need to explicitly cast some of those though, like adding an INT to a CHAR. Could just be the really old C++ book I had (it tought you quite a few things that are considered old fashioned in today's C++ lol).

I've done some C/C++ before. I learned to ptogram in C++, actually, though I've hardly touched it since. I've also programmed a bit in Java, VB, and C#. I've taken programming in highschool (Covered VB, Java, and tought myself Turing for a programming contest in 5 minutes) as well as college (Terrible VB.net class. My teacher was smart, but had NO idea how to actually teach.)

I considered going to college for computer programming, but felt that I really didn't want to spend my day every day sitting at a desk programming. When I saw that there was both the hands on aspect and the programming aspect to PLCs I eally liked that idea.

The company I'm at right now specializes in industrial automation. We come in as contractors instead of an in-house team doing stuff at companies, meaning that we get lots of odd jobs and such. We do tend to concentrate on PLCs, but we don't just do them.

 

[kamenges]

Originally posted by Peter Nachtwey:

I have never taken a PLC course. Therefore I have never been brain damaged in thinking in limited terms imposed by PLCs.

Blessing and a curse, Peter. Blessing and a curse. You also seem to get extremely frustracted when you can't do something the 'easy' C way in a plc, even though there may be a perfectly valid way of doing it.

Remember, it's just a tool. Both a manual and an electric screwdriver will drive screws. One just makes your forearm a bit sore.

As for the +D thing, I'm not sure how much of the total equation is supposed to be encompassed in the constant. Is the constant supposed to be a hex number or a decimal number?

Keith

Keith

 

[TurpoUrpo]

Yeah

L#16#709921077

It can be hex if its DW#16#709921077, or decimal if its L#709921077

 

[TConnolly]

Alaric, whether the code is compiled down to pseudo code or machine code there many small optimizations that can be made.
r**3 can still be converted to r*r*r
divides by powers of 2 can still be converted to shift rights. Divides by a constant can be converted to multiplies by 1/constant. MOD(number,2**n) can be converted to AND 2**n-1.

I don't disagree Peter. I'm don't mean to say that PLC manufacturers cannot optimize the interpreter. What I'm getting at is that I don't think many PLC manufacturers bother to have the interpreter try and optimize the user program much if any. I suspect that most of them create a generalized interpreter that scans the user program in a rigid and literal manner and leave it at that.

Manufacturers that create speedy efficient controllers like the RMC are the exception rather than the rule [/shameless plug]. We shouldn't assume that the PLC will optimize our programs for us, unlike any half way decent C compiler does.

For example, since I'm not privy to the inner working of XPY except what AB publishes about it (2^Y*log2(X)) then I'll write MUL N7:0 N7:0 instead of XPY N7:0 2 because I know how it will process and have no reason to think that the XPY will do anything except 2^2*log2(N7:0) . Likewise, I convert polynomials from Ax³ + Bx² + Cx + D to ((Ax +B)X +C)X + D because I don't depend on the PLC to due anything at all to optimize it.

 

[Peter Nachtwey]

C rules

Blessing and a curse, Peter. Blessing and a curse. You also seem to get extremely frustracted when you can't do something the 'easy' C way in a plc, even though there may be a perfectly valid way of doing it.

Don't get me started.
C isn't always the easy way either. It depends on what I need to do. I can switch between many tools. I use Mathcad, wxMaxima, and Scilab more often now. I use our own ST for example programs.

Remember, it's just a tool. Both a manual and an electric screwdriver will drive screws. One just makes your forearm a bit sore.

I have heard that before some where. Some tools are much better than others.

As for the +D thing, I'm not sure how much of the total equation is supposed to be encompassed in the constant. Is the constant supposed to be a hex number or a decimal number?

A DWORD.

I went through my library and I found and cube root initial estimator that provides 6 bits of accuracy. The C code does two iterations using the method L D[AR2,P#0.0] suggested. Every iteration doubles the precision of the estimate so only two iterations are required for the method L D [AR2,P#0.0] showed.

I found another method that basically calculates two 4th order polynomials and divides one the other. It requires only one iteration.

float cuberoot(float a)
{
	float r = a;
	unsigned int* p = (unsigned int *) &r;
	*p = *p*0.33333333333 + 709921077; // 709921077 is the magic number
	float inv3a = 1.0/(3.0*a);
    printf("\na=%10.7f\n",a);
    printf("\nr=%10.7f\n",r);
	r = r*(4*a-r*r*r)*inv3a;
    printf("\nr=%10.7f\n",r);
	r = r*(4*a-r*r*r)*inv3a;
    printf("\nr=%10.7f\n",r);
	r = r*(4*a-r*r*r)*inv3a;
    printf("\nr=%10.7f\n",r);
	r = r*(4*a-r*r*r)*inv3a;
    printf("\nr=%10.7f\n",r);
	r = r*(4*a-r*r*r)*inv3a;
	return r;
}

int main (int argc, const char * argv[])
{
	float a= 15625.0;
	float r0=pow(a,(1.0/3.0));
	float r = cuberoot(a);
    printf("\npow(a,1/3)=%10.7f,  cuberoot(a)=%10.7f\n",r0,r);
    return 0;
}

I printed the iterations

a=15625.0000000

r=25.5728531  //  You can see the initial estimate is very close with just a multiply and add.  No exp required.

r=24.9733448

r=24.9999447

r=25.0000000

r=25.0000000

pow(a,1/3)=25.0000000,  cuberoot(a)=25.0000000

Bobbias, cut and paste that code and try it. It is pure gold. However, it works only floor floats. A double version has a different magic number that is added. Also, the float must be in IEEE format. On the TI DSP that we use I had to find a different magic formula because the old TI DSP doesn't use a standard floating point format.

It will be interesting to see what L D [AR2,P#0.0] comes up with when he uses the magic cube root estimator. The magic estimator will work with the method he showed too.