Does This Method Provide The Exact Same Return As SPINs Square Root Operator

idbruce

Hello Eveyone

While I am not working on my 3D printer, I have been spending some time porting certain sections of the Teacup firmware for 3D printing, from C/C++ to SPIN. One of the functions that I am porting obtains the square root of an integer, and I am wondering, if all this code is really necessary, since SPIN already has a square root operator. At the beginning of the code, there is a comment about the return of the function and it has me concerned about the return value, and if it would be identical to using the SPIN operator:

\return sqrt(a - 1) < returnvalue <= sqrt(a)

EDIT: In other words... On the surface, it appears that it would return the same, but would SPIN have a similar expression, or is there any possibility they could offer different results.

The Propeller Manual states:

Square Root ‘^^’
The Square Root operator returns the square root of a value. Square Root can be used in both variable and constant expressions. When used with variable expressions or integer constant expressions, Square Root returns the 32-bit truncated integer result. When used with floating-point constant expressions, Square Root returns the 32-bit single-precision floating-point result. Example:

X := ^^Y

Square Root becomes an assignment operator when it is the sole operator to the left of a variable on a line by itself. For example:

^^Y

This would store the square root of the value of Y back into Y.

CON

  _CLKMODE            = XTAL1 + PLL16X
  _XINFREQ            = 5_000_000

  'Full Duplex Serial Pin Assignments
  TX                  = 30
  RX                  = 31

  'Full Duplex Serial Driver Constants  
  COMM_MODE           = 0
  COMM_BAUD_RATE      = 115_200

VAR

  'Global variables for the IntegerSquareRoot method
  BYTE bIntSqRt_C
  BYTE bIntSqRt_Z
  BYTE bIntSqRt_J
  BYTE bIntSqRt_Y2
  
  WORD wIntSqRt_B
  WORD wIntSqRt_X
  WORD wIntSqRt_I
  WORD wIntSqRt_Y2

  LONG lIntSqRt_Y2

OBJ

  Terminal:     "FullDuplexSerial"
    
PUB Main | SquareRoot

  'Start FullDuplexSerial
  Terminal.Start(RX, TX, COMM_MODE, COMM_BAUD_RATE)
  WAITCNT(CNT + (1 * CLKFREQ))

  SquareRoot := IntegerSquareRoot(603729)
  Terminal.Dec(SquareRoot)

PUB IntegerSquareRoot(IntSqRt_A)

  {
  \param a find square root of this number
  \return sqrt(a - 1) < returnvalue <= sqrt(a)
  }

  wIntSqRt_B := IntSqRt_A >> 16
  bIntSqRt_C := wIntSqRt_B >> 8
  wIntSqRt_X := 0
  bIntSqRt_Z := 0
  
  bIntSqRt_J := $8

  REPEAT WHILE bIntSqRt_J

    bIntSqRt_Z |= bIntSqRt_J
    bIntSqRt_Y2 := bIntSqRt_Z * bIntSqRt_Z
    
    IF bIntSqRt_Y2 > bIntSqRt_C
    
      bIntSqRt_Z ^= bIntSqRt_J

    bIntSqRt_J >>= 1  
  
  wIntSqRt_X := bIntSqRt_Z << 4   

  wIntSqRt_I := $8

  REPEAT WHILE wIntSqRt_I

    wIntSqRt_X |= wIntSqRt_I
    wIntSqRt_Y2 := wIntSqRt_X * wIntSqRt_X
    
    IF wIntSqRt_Y2 > wIntSqRt_B
    
      wIntSqRt_X ^= wIntSqRt_I

    wIntSqRt_I >>= 1   
  
  wIntSqRt_X <<= 8
  
  wIntSqRt_I := $80

  REPEAT WHILE wIntSqRt_I

    wIntSqRt_X |= wIntSqRt_I
    lIntSqRt_Y2 := wIntSqRt_X * wIntSqRt_X
    
    IF lIntSqRt_Y2 > IntSqRt_A
    
      wIntSqRt_X ^= wIntSqRt_I

    wIntSqRt_I >>= 1  

  RETURN wIntSqRt_X

kwinn

The square root of a number should be the same in both cases.The only caveats would be rounding/truncating errors. Those should be very small if the algorithms are implemented correctly.

idbruce

kwinn

Thanks for the response.

Although I do not believe I have ever used it, but I may have once or twice, I believe that C/C++ has operators readily available for square roots. This function must have a special purpose, otherwise I cannot imagine someone taking the time to code it all out.

The only caveats would be rounding/truncating errors.

I am just wondering if rounding/truncating is the basis for writing this function in the first place. So instead of being errors, it could possibly be the desired effect.

Heater.

C/C++ does not have operators for square root and other math functions, log, sin, cos, etc.

Rather those functions are provided by libraries. The standard libraries. Hence you will often see "#include <math.h>" in C sources.

One may well want to write ones own square root and other such functions if those standard libraries are not available. I have no idea what is available in the Arduino environment for example.

Or one may do that to get some extra optimization in some cases. Or perhaps you are using some fixed point number format and need specialized functions that work with that.

You could always do some tests to see what different results they produce, if any. Just run the C and Spin codes against a set of random numbers and see what happens.

Anyway in may applications the odd bit up or down in rounding won't make any noticeable difference.

idbruce

Heater

C/C++ does not have operators for square root and other math functions, log, sin, cos, etc.

Rather those functions are provided by libraries. The standard libraries. Hence you will often see "#include <math.h>" in C sources.

I had forgotten about #include <math.h>, which I have used in the past and was an oversight on my part. Thanks for setting me straight

One may well want to write ones own square root and other such functions if those standard libraries are not available. I have no idea what is available in the Arduino environment for example.

Or one may do that to get some extra optimization in some cases. Or perhaps you are using some fixed point number format and need specialized functions that work with that.

Hmmmm.... You stated that very well Heater! Your first reason makes a lot of sense and seems like the primary reason, however I believe you are correct about the second one also. I believe it is a combination of both of those reasons.

EDIT: Or should I say that since it had to be written, might as well take the time to optimize it.

Electrodude

If the only difference between ^^ and your square root function is rounding/truncation (and the rounding/truncation matters), using ^^ and then correcting it should be significantly faster:

PRI sqrt(x) : y | x2

  y := ^^x
  x2 := y*y

  y += { function of (x - x2) or something }

That way, you can take advantage of the PASM square root routine in the Spin interpreter instead of having the overhead of a big Spin loop.

idbruce

Electrodude

I would have to agree. Thanks for your input.

Dave Hein

The standard C libraries contain a sqrt function that takes a double argument, and returns a double result. This would be quite inefficient if all you need is the square root of an integer value. That's probably why the Teacup software contains it's own integer square root function. Here's a 32-bit integer square routine function that I use in spinsim.

uint32_t sqrt32(uint32_t y)
{
    uint32_t x, t1;

    x = 0;
    t1 = 1 << 30;
    while (t1)
    {
        x |= t1;
        if (x <= y)
        {
            y -= x;
            x += t1;
        }
        else
            x -= t1;
        x >>= 1;
        t1 >>= 2;
    }
    return x;
}

It matches the value that is generated by the Spin ^^ operator. You could write a simple C test program to test this against the function that is used in the Teacup software. Or since it appears that you've already converted the Teacup square root function to Spin, you could just write a simple Spin test program to compare it to the Spin ^^ operator. That would quickly answer the question you asked in your first post in this thread. As others have said, even if it doesn't exactly match, the difference probably wouldn't matter.

lonesock

A quick way to pre-round an integer square root is to add in the sqrt of the value 1st.

y := ^^( x + ^^x )

Jonathan

idbruce

Thanks for the input and code Dave.

Jonathan, thanks for the snippet. By the way, what is your part in the KISSlicer project? I will be using that software to create my GCODE.

Dave Hein

Bruce, I looked at the original Teacup C code for the square root function, and it looks like this

/*!
integer square root algorithm
\param a find square root of this number
\return sqrt(a - 1) < returnvalue <= sqrt(a)

see http://www.embedded-systems.com/98/9802fe2.htm
*/
// courtesy of http://www.embedded-systems.com/98/9802fe2.htm
uint16_t int_sqrt(uint32_t a) {
	uint32_t rem = 0;
	uint32_t root = 0;
	uint16_t i;

	for (i = 0; i < 16; i++) {
		root <<= 1;
		rem = ((rem << 2) + (a >> 30));
		a <<= 2;
		root++;
		if (root <= rem) {
			rem -= root;
			root++;
		}
		else
			root--;
	}
	return (uint16_t) ((root >> 1) & 0xFFFFL);
}

Your Spin code doesn't seem to be a line-for-line translation of this C code. Anyhow, I wrote a test program to compare my sqrt32 function with Teacup's int_sqrt function, and the results match for all 2 billion+ positive values. It only took a few minutes to test out all the values.

So feel free to use Spin's ^^ operator. It matches exactly.

idbruce

Dave

Thanks for the testing and info, I certainly appreciate it.

I don't know where you got that code, but it is different then the code I have, however I downloaded mine about 9 months ago, so it may have changed, and github won't allow my old browser in for some odd reason, so I cannot get the latest code. Anyhow, this is what I started with:

/*!
  integer square root algorithm
  \param a find square root of this number
  \return sqrt(a - 1) < returnvalue <= sqrt(a)

  This is a binary search but it uses only the minimum required bits for
  each step.
*/
uint16_t int_sqrt(uint32_t a) {
  uint16_t b = a >> 16;
  uint8_t c = b >> 8;
  uint16_t x = 0;
  uint8_t z = 0;
  uint16_t i;
  uint8_t j;

  for (j = 0x8; j; j >>= 1) {
    uint8_t y2;

    z |= j;
    y2 = z * z;
    if (y2 > c)
      z ^= j;
  }
  
  x = z << 4;
  for(i = 0x8; i; i >>= 1) {
    uint16_t y2;

    x |= i;
    y2 = x * x;
    if (y2 > b)
      x ^= i;
  }
  
  x <<= 8;
  for(i = 0x80; i; i >>= 1) {
    uint32_t y2;

    x |= i;
    y2 = (uint32_t)x * x;
    if (y2 > a)
      x ^= i;
  }

  return x;
}

Dave Hein

Oh, you're right. The version they are using is the same as the one you posted. It is located in dda_maths.c. The version I posted is in the file alg.c in the research folder. Anyhow, I tested with the other version, and it matches the ^^ operator also.

idbruce

It is located in dda_maths.c

It surprises me that you have those files at your disposal. Do you have an interest in Teacup?

Dave Hein

I just downloaded the ZIP file from from https://github.com/Traumflug/Teacup_Firmware . I then searched the source tree for "sqrt". It didn't occur to me that there might be more than one sqrt routine in the source tree.

I believe Loopy looked at porting this to the Prop several months ago, and I looked at the code back then. I'm not planning on doing any 3D printing in the near future, but I do have some interest in what other people are doing with 3D printing.

idbruce

For Those That May Be Interested

After reviewing some previous documentation pertaining to Teacup, I rediscovered this summarization, which outlines the processing of key information.

At startup, the code in mendel.c is run first. This initialises all the modules that need it, then starts polling the clock flags and feeding incoming serial characters to the gcode parser. The gcode parser processes each character individually, keeping track via internal state rather than buffering a line and skipping back and forth. The gcode parser converts floating values to integer or fixed-point representations as soon as it encounters a non-numeric character. It calls many module functions directly, but the most interesting part is move creation, where it passes a target position and speed to enqueue()[dda_queue.c] which adds it to the queue, and fires up dda_start()[dda.c] if the queue was empty. dda_start initialises the dda, figures out the stepper directions and first step timeout and a few other bits of housekeeping, then sets the timer for the appropriate timeout. When the timer fires, it calls dda_step()[dda.c] which sends all the step signals then figures out the next step timeout based on acceleration and speed settings. When the last step has been made, the dda "dies" (sets 'live' property to 0) after which queue_step[dda_queue.c] advances the queue read pointer and starts the next dda.

Source: Wikipedia - http://reprap.org/wiki/Teacup_Firmware

There is some good documentation, source code, and hyperlinks, within the Teacup firmware, however I still firmly believe that differences in architecture will make a direct porting inefficient, as the example at hand clearly demonstrates, or should I say as Dave has so kindly proven. I think the previous summarization and portions of the Teacup firmware, should be used as a roadmap for creating a new firmware that is unique to the architecture of the Propeller chip.

idbruce

As mentioned elsewhere, there are a lot of files associated with Teacup and a lot of unnecessary code. After some previous reviewing of the source code and my latest efforts, I now have a much better understanding of the key aspects of operation. I have always believed that the method dda_create, within dda.c, was the most critical method for porting to SPIN, and I am now convinced. This method should have been my starting point for a port, but sadly it was not.

I now believe that once this method is ported, it should provide a well rounded overview of what is really necessary for some Propeller firmware, by working within the expansive boundaries of Propeller chip, and utilizing it to it's full potential, for this particular type of application.

Of course, the rest of dda.c will need to be ported, but much of the remaining code will be changed.

Just $0.02 worth.

lonesock

idbruce wrote: »

...
Jonathan, thanks for the snippet. By the way, what is your part in the KISSlicer project? I will be using that software to create my GCODE.

My part of KISSlicer is basically coding in the bugs. And all the other parts. [8^)

There has been a lot of progress recently, a lot of people have been providing useful feedback and bug reports and feature requests, etc. Have you tried the 1.5 Beta 2.xx on the forums?

Sorry to hijack the thread, everybody.

Jonathan

idbruce

My part of KISSlicer is basically coding in the bugs. And all the other parts. [8^)

Have you tried the 1.5 Beta 2.xx on the forums?

No, I am still currently building and testing portions of the printer.

Dave Hein

Bruce, I'm curious why you are converting the Teacup software to Spin, and not just leaving it in C. I noticed that you participated in this thread that attempted to port the C Teacup code to the Propeller. Why not continue from where that thread left off? It seems that converting to Spin code would be error-prone and unnecessary unless the program wouldn't fit in C. Why are you converting it to Spin?

idbruce

Dave

Bruce, I'm curious why you are converting the Teacup software to Spin, and not just leaving it in C. I noticed that you participated in this thread that attempted to port the C Teacup code to the Propeller. Why not continue from where that thread left off? It seems that converting to Spin code would be error-prone and unnecessary unless the program wouldn't fit in C. Why are you converting it to Spin?

In my opinion, there is just too much junk that I do not want or need, which just eats up valuable resourses. Additionally, in the form of raw firmware, considering the complexity of the code, I doubt that I would be able to make it work, with the added factor, that there is just some things that I would like to do differently.

To be perfectly honest, over the last 2 weeks, I have been all over the place with my thinking, first port and then don't port. For the most part, any attempt at a port just gives me a much better understanding of the code. However a couple of days ago, I really started digging deep into the source code, tracking the process and events, and I believe I walked away with a good deal of knowledge (maybe ). Anyhow, for now, at this point in time, I have decided to write my own firmware for the Propeller from scratch, using the Teacup firmware as a guide in the right direction.

I do not know if you noticed, but I recently released a new stepper driver object within the forum. I am currently modifying that stepper driver once again, to be of more use as a 3D printer or CNC driver, so that I may use it in the future firmware, but I am running into problems, because I have forgotten some stuff about cog usage and variables.

Perhaps you might be willing to assist me??? Here is the scenario....

I want to run each driver (axis) within it's own cog, but I also want each of them to start movement simultaneously. I am very close to figuring it all out, but my cogs are dying (after the InitializePins method) and I know I need some sort of repeat loop to keep them alive, meanwhile updating some variables and having a trigger event. Basically I need to update (3) variables within a child object, while the child object is in a loop. Instead of continually checking the values for alteration, I will be using a WAITPEQ to inform the loop that the variables have been changed and start movement.

The attached code is a bit choppy, because it is a work in progress, but I think you will have a good understanding of what I am attempting and what is holding me back. As it stands now, the drivers function and the axes move, but not in parallel as intended. Like I said, I know it needs a repeat loop. Any suggestions???