The purpose here is to describe a method of computing trignometric functions by computer. These functions are sine, cosine, arc tangent, etc. There are many ways to solve this problem, but almost every implementation takes advantage of symmetry, so we will do the same and generally limit the discussion to a small convenient range such as angles (in radians) of zero to π/2 or π/4. When it is pertinent, the range will be described precisely and, for completeness, the necessary transforms to produce the full range of angles will also be considered, eventually. With that in mind, the simplest way to generate trigonometric functions is to use a polynomial function that so closely approximates the needed function that the error can hardly, if at all, be detected since the computer can only represent numbers with a limited amount of precision anyway. The problem, which will be taken on faith for the moment, is that the computation of polynomial involves a number of multiplication and/or division instructions that take too long to perform. Another method, at least for sine and cosine, is to use the Taylor's series. It also requires a number of multiplications, but at least generates both sine and cosine simultaneously. The purpose here, however, is not necessarily to show how to compute trignometric functions the fastest, but to show how to compute them in a novel way. If it turns out to be fast, then all the better. Still, this method requires hardly any multiplication or division instructions, which is good if the target machine doesn't have them. In fact, the method simulates a division instruction, so it sounds like good fun.
Before proceeding further, it is time to establish some dogma. The most commonly available trigonometric functions are sine, cosine and arc tangent. Why? Because none of the others are functions; they either have limited domains (e.g. arc sine and arc cosine) or they have singular points (e.g. tangent or secant) where their value is undefined. On one side of these points, the value goes to positive infinity and on the other the value goes to negative infinity. The way these functions are actually implemented is to cheat a little and say that since these points lie at π times something, and π can not be represented exactly in the computer anyway, that these points don't really exist. If a computer approximation of a singular point is put into these functions, the value returned is a very large (probably the largest representable) positive or negative number. Still, the functions cosecant and cotangent have problems if the something times π is zero. Zero is representable exactly and for these functions is a singular point. This case is probably regarded as division by zero. I don't even want to think about it. So we are left with sine, cosine and arc tangent.
But even arc tangent often has problems. The most obvious use is to find the angle of a line, say from the origin to the point (x, y). It would be easy to get the arc tangent of y/x, unless x is zero. So that must be tested as a special case. Also, if (x, y) is a vector rather than a line, then the range of the angle is a complete circle rather than a half of a circle. The problem of dividing y/x is that pairs of negative signs cancel out and disappear forever. This must be avoided with another special case. In C on UNIX, the problem is solved with a special function called atan2, in which x and y are provided independently to the function and it deals with both special cases.
Another thing to consider is that the independent computation of sine and cosine is probably an artifact of the polynomial approximation technique. In this method, there is no relationship between the two but the Taylor's series produces them as a pair. In fact, there is reason to believe that unless the relationship between sine and cosine is deliberately destroyed, any method will yield them as a pair. Similarly, there is reason to believe that, unless deliberately destroyed, any method of generating the arc tangent of x and y will also generate the square root of x² + y²!
The question becomes: just what do we want from a trigonometric function? I suggest there are exactly two trigonometric functions: rectangular to polar and polar to rectangular. If a rectangular coordinate pair is (x, y) and a polar coordinate pair is (r, a), that is the radius and angle, then I propose to provide functions to convert one to the other, and nothing more. Since these functions are defined everywhere, there are no problems with singular points, division by zero or limited domains. Nor are there problems with very large numbers approaching infinity: r can't be any larger than x + y. I also claim that in most cases, coordinate conversion is the most desirable application, but they just aren't available as such in most computer languages. I suspect that the reason they are not available is tied to the rather quaint notion that functions should return exactly one value (and preferably take exactly one argument). That is a fine idea for plotting functions on graph paper; unfortunately, clinging to that notion is causing a lot of trouble. If the problem can be solved using coordinate conversions, then the trick is to use a different piece of graph paper.
Of course, it is still important to be able to compute sine or arc tangent, for example, but these can be easily derived from the coordinate conversion functions. I claim that the conversion using the method that will be described can be performed as quickly as any sine function implementation; however, to derive sine from the conversion may take a little extra time. For example, to get the sine of an angle a, put (1, a) into the polar to rectangular conversion routine and the needed value comes back as y. The cosine would be x. The tangent, cotangent, secant and cosecant are y/x, x/y, 1/x and 1/y, respectively, having coped with x = 0 or y = 0. The inverse functions are trickier, the simplest are arc tangent and arc cotangent. For arc tangent, if the argument is z, then put (1, z) into the rectangular to polar conversion routine and the needed angle is a. For arc cotangent, put in (z, 1). For arc sine and arc cosine it is necessary to know a value, call it s which is the square root of 1 − z². For arc sine, put in (s, z) and for arc cosine put in (z, s). For arc cosecant and arc secant, the reciprocal of the argument must be computed and then the procedure is the same as that for arc sine and arc cosine respectively.
So the dogma is that only the two functions are necessary: one for rectangular to polar conversion, call it rtop, and one for polar to rectangular, call it ptor. This is not to suggest that the traditional trignometric functions should not exist, but that they can be easily derived from the two and that the problems of limited domain and singular points can be handled at that higher level. For example, a math coprocessor might provide only the two functions and let software perform the derivation.
It is worth noting that rtop and ptor are not completely reversable in the strict sense. Polar coordinates are generally limited to a radius not less than zero and an angle not less than zero and less than 2π. In addition, a zero radius normally implies a zero angle. Rtop will produce values in this range although ptor will take anything. If ptor and rtop are run in sequence, the effect is to "normalize" the polar coordinates. This generally isn't a problem. It is worth noting that the traditional arc tangent function has the same difficulty: theoretically any integral number of π's can be added to an arc tangent value and is still correct. Other inverse trignometric functions have similar problems.
Now the idea behind this method in both these conversions is that an angle can be represented as the sum of smaller angles. The trick is to pick the smaller angles in such a way as to make the computations simple. The method is taken from many sources and originally dates back to someone named Briggs three centuries ago who was trying to compute logarithms and exponentials in base ten. The mathematics have been adapted to trigonometry and the specific equations are pretty much as they appear in: W. H. Specker, "A Class of Algorithms for ln x, exp x, sin x, cos x, tan−1 x, and cot−1 x," IEEE Trans. Electronic Computers, EC-14, 85-86 (1965). The fact that the arc tangent function also yields the r value was my own discovery although I am sure everybody who worked on such problems noticed this before I did. The papers I have seen are generally so mathematical as to be nearly illegible and virtually incomprehensible. Rather than annoy the audience with such gibberish, I have found a reasonably simple graphical interpretation of the algorithm which is at least as informative, if not entirely satifactory as a mathematical proof. The questions of accuracy and precision will not be considered except to note that these papers suggest that the error is rather small.
Representing an angle in the polar coordinate representation as the sum of smaller angles is simple in principle since the angle is explicitly given as one of the polar coordinates. Summing is easily done by computer; it takes virtually no time. The thing to be avoided is multiplication, and we will do this no more than twice in each of the two functions rtop and ptor. The problem with summing angles comes with the rectangular coordinates. In fact, it is virtually begging the question to suggest that an angle can be represented in rectangular coordinates; if that could be done, the problem would be solved essentially. Still, it is possible, although it requires multiplication (four of them), to rotate a rectangular coordinate by the angle of a vector given by another rectangular coordinate. If the rotation vector is cleverly chosen, and its angle is known, then it is possible to get the same effect as expressing the angle of the rectangular coordinate explicitly. It is also possible, again by clever selection of the rotation vector, to avoid the multiplications. So clever, in fact, that the four multiplications can be replaced by two additions and two right shift instructions. A right shift of n bits moves the bits of a value to the right n places. Since the computer is binary the effect is to divide the value by two raised to the power of n. This is all very clever and I do not know how or by whom this angle of rotation was selected or discovered. Nor can I suggest a deductive reason to select such an angle. I suspect that this bit of cleverness was an inductive leap in the mind of a mathematician and perhaps we will never know if the angle is the best, but only that it is good. So I can not explain why the algorithm was chosen; I will explain, graphically, why it works.
Actually, there is more than one angle of rotation. The angle will be represented as the sum of smaller angles, so it is necessary to perform a sequence of rotations. The process, in converting a polar coordinate to a rectangle coordinate, is to start with the rectangular coordinate (r, 0). The interesting thing about this particular coordinate is to note that any coordinate whose first number is positive and second number is zero actually requires no conversion at all: it represents the same point in both polar and rectangular coordinate systems. So, given the rectangular coordinate (r, 0), the process is to apply a series of rotations (about the origin) so that the coordinate comes to the correct location. While the rectangular coordinate is being rotated into position, the polar coordinate (r, a) is being rotated back to (r, 0). It stands to reason that if the rotations are by equal amounts but in opposite directions, then when the polar coordinate reached (r, 0), the rectangular coordinate will reach the correct location, that is, the equivalent of the original polar coordinate and the job will be done. Now it is not precisely true, but accurate enough for the moment, to say that each angle of rotation will be exactly half that of the previous rotation. And let us say that we will only be dealing with angles between zero and one quarter of a circle. Suppose we are dealing with an angle of exactly one quarter of a circle, as the original polar coordinate. If we rotate this coordinate one eighth of a circle toward zero and the rectangle coordinate one eighth of a circle in the opposition direction, then we will have performed one half of the necessary rotation; there is still another eighth to go. Next, we can rotate both coordinates one sixteenth of a circle. Then, a thirty-second, and so one. This process could go on forever but never quite be finished; however, each time the answer becomes twice as close to the correct answer, so it doesn't take long before our approximation is as good as our computer can represent it. At that point, there is no reason to continue so the process stops then, or sooner if the needed accuracy becomes sufficient.
Of course, this illustrates only one case: an original angle of one quarter of a circle. Other angles require a decision: if the rotation would cause the polar coordinate to be rotated beyond zero, then that rotation is not performed. In fact a binary number could be constructed composed of bits each of which would indicate if the rotation by that fraction of a circle was performed. It is not hard to see that all the possible binary numbers represent all the possible angles of the original polar coordinate. In fact if the angle of the polar coordinate was measured in units of whole circles, the angle itself would be identical to the binary number that was constructed. That suggests that the binary number should not be constructed at all; instead the bits of the original angle could be examined one by one and if a bit were set, then the rectangular coordinate could be rotated by the fraction of a circle represented by the position of the bit. But it is not that simple; unfortunately, the angle of rotation can not be exactly half that of the previous rotation without mucking up the clever business that allows the four multiplications to be replaced by two additions and two shifts. Each rotation will actually be slightly larger than one half the previous, as we shall see. In fact, the progression of rotations isn't even linear (so it isn't possible to divide the original angle by some magic number and examine the bits of the result). The necessary process is not that simple, but is simple enough. The angle of each rotation is known and, of course, the original angle of the polar coordinate is also known. It is simple enough to discover if the polar coordinate would be rotated beyond zero by simply comparing the two numbers; if the first is larger, the rotation is not performed. If the first is not larger, the polar coordinate is rotated by simply subtracting the angle of the rotation. That is simple enough. Now if there is a good way to rotate the rectangular coordinate, then the problem is solved, at least to convert from polar to rectangular. The precise angles and the method of rotation of the rectangular coordinate will be described soon enough, the purpose here was to describe the technique.
Now converting from rectangular to polar is very similar. In fact, it is virtually identical. If we start with the polar coordinate (r, 0), and the original rectangular coordinate, we simply rotate the rectangular coordinate back to zero and the polar coordinate in an equal and opposite direction. Simple enough. Of course, it is not so easy to test rotate the rectangular coordinate as it was for the polar angle. It would be necessary to actually carry out the rotation first to see if the new point would go beyond zero. If it did, then that rotation would not be used: so the test rotation could be carried out and stored in a temporary variable until the zero test was performed. If it did not cross, then the rectangular coordinate could be updated to the temporary. In fact, there is a simpler method, discussed below. The zero crossing is easy enough to detect since the y value would go from positive to negative, having limited the discussion to angles less than a quarter circle. And, of course, the rotation of the polar angle is a simple matter of adding to it the rotation angle. The angles, as you may have guessed, are the same clever angles used in the polar to rectangular conversion routine discussed earlier. So now the technique of both conversions has been described and we can get to those clever angles. Or can we?
To test your sobriety and to simplify the discussion, there is a flaw in the description of the rectangular to polar conversion. Perhaps you have found it; perhaps you have solved it as well. The flaw is that we started with the polar coordinate (r, 0). Of course, r is not known; it is one of the outputs. Actually, we don't need r as we did for the other conversion. In the polar to rectangular, it was important to start with (r, 0) as the initial rectangular coordinate since in the final rectangular coordinate, the square root of x² + y² must be equal to r, and through each rotation, this is ensured. But in converting from rectangular to polar, as each rotation is performed, it is really just the polar angle that is being determined; in polar coordinates the distance r is independent of the angle and so is not affected by the rotation. It was assumed the polar coordinate (r, 0) was used initially to simplify the discussion, make it easier to compare with the other conversion and keep the focus on the rotation. But r never actually appears in any calcuation; it's just not in the algorithm, so we can call it anything we want. Still, it can not be ignored, since it is an output. And given that there are only two outputs, it would appear that half the problem remains to be computed. In fact, it has been computed: when the calculation of the angle has been completed, the rectangular coordinate has been rotated to (nearly) zero, that is, the coordinate lies on the positive x-axis. This x value is r. Mathematically, in the initial rectangular coordinate, the square root of x² + y² is equal to r. The rotations do not change this relationship. If, after all the rotations, y is zero (or nearly so), then the square root of x² + y² is x, which is the desired r, or nearly so.
The clever angles may now be divulged. The first is the arc tangent of one. The second is the arc tangent of one half. The third is the arc tangent of one fourth, etc. Graphically, these are the angles of a triangle with a base of one and an altitude of one, one half, etc. These numbers are not convenient to compute. In fact, I will not even say how to compute them. Fortunately, you only have to compute them once, ever. The numbers are to be stored in a table in the source code or hardware for the converter. As far as we are concerned these numbers are constants since they will always be the same sequence of numbers. As the sequence progresses, the numbers become progressively smaller until they become smaller than the precision of the machine. It turns out that as the argument to arc tangent becomes very small, the value of the arc tangent (in radians) becomes very close its argument. Since an argument to arc tangent in the sequence tends to half that of the previous argument, only as many numbers are needed as there are bits of precision in the machine. Of course, for floating point numbers, every number has an equal amount of precision no matter how small the number is; here the precision we are talking about is the precision of the angle in the polar coordinate, once the arc tangent numbers become too small to have any meaning in the representation of the angle, there is no point in continuing. The numbers are used by selectively adding them to the polar angle to either rotate the angle back to zero (ptor) or rotate into position from zero (ptor). All that is needed, then, is to look these numbers up from the table in sequence and add, so this is very efficient. It does not matter how difficult it may be to compute these arc tangents since the computation is not carried out during the conversion.
Of course, it is also necessary to use the same sequence of angles to rotate the rectangular coordinates, that's the real trick. To show that the angles will work as advertised, I must resort to some analytic geometry. The rotation of a rectangular coordinate (x, y) by an angle a (about the origin) is (x cos a − y sin a, x sin a + y cos a). The sequence of angles, described mathematically, is arc tangent of two to the power of minus i, where i is zero, one, two, etc. Symbolically, a = atan 2−i, i = 0 to n, and n is the number of bit of precision. Now if a = atan 2−i, then tan a = 2−i. OK, if we take (x cos a − y sin a, x sin a + y cos a) and factor out cos a, we get cos a (x − y sin a / cos a, x sin a / cos a + y). Now what? Well, sin a / cos a = tan a = 2−i, so the new rotated coordinate can be stated as cos a (x − y 2−i, x 2−i + y). Ignoring the cos a for the moment, we can say that to apply the ith rotation in the sequence, the new x value is the old x value minus the old y times 2−i and the new y value is the old y value plus the old x times 2−i. The magic is that to multiply a number by 2−i is the same (in a binary computer) as shifting the number to the right i places. So all that is necessary to compute the new x is to add once and shift once. And for y, subtract and shift. These are easy operations. But what about the cos a?
It is a mathematical identity that sin² a + cos² a = 1. Dividing the whole thing by cos² a yields: tan² a + 1 = cos−2 a. Solving for cos a gives cos a = (tan² a + 1)−½. Since tan a = 2−i, cos a = (2−2i + 1)−½. Graphically, the reciproal of this is the hypotenuse of the aforementioned triangle. To multiply by cos a is to divide out this hypotenuse so that the rotation performs no scaling. Now notice that in the rotation we are multiplying the whole coordinate (x and y) by cos a, if we were to apply all the rotations of the sequence, we would be adding and shifting each x and y and then multiplying by cos a, for each angle in the sequence but we could equally well perform all the adding and shifting for each angle and then perform all the multiplies at once. Actually we could multiply all the cos a's together first and multiply the result by the added and shifted angle. Since the a's are constants, and we want the product of the cos a's, this number, call it K, could be computed and stored in the source code and never need be computed again. Unfortunately, that only works if all the rotations are applied, which only occurs for one particular angle. For all the others we intend to apply only those rotations that would cause the original angle to rotate back to zero.
It's time for another trick: another mathematical identity is that cos a = cos −a. Our problem is that we plan to apply some rotations and not others. Of course not applying a rotation is the same as rotating by an angle of zero. So for any angle in the sequence, we plan to rotate by an angle of either a or zero. But if we could change the choice to an angle of either a or −a, then it wouldn't matter if we were using cos a or cos −a, the value is the same. Graphically, cos a = cos −a because the triangle we've been discussing has a positive altitude, but if it were reflected over its base (flipped over), the effect would be an angle of −a. As was stated above, cos a is the reciprocal of the hypotenuse of this triangle. Clearly, flipping the triangle does not change the length of the hypotenuse, so cos a = cos −a. Also the product of all the cos a's (or cos −a's) would never change. We could compute this product, K, and multiply all our rectangular coordinates by it and get the correct answer. Is it feasible to switch from the choice of a or zero to the choice of a or −a? Indeed, it's not only possible, but it's better.
All we are really trying to achieve is to rotate the original angle back to zero, it is only important that the angle gets closer to zero. In the original plan, we would apply a rotate if the rotation angle was smaller than the given angle. For the case of rotating a polar angle this would simply require comparing the two angles. That is simple but somewhat redundant, since if the rotation is applied, the rotation angle is subtracted from the original angle. Comparing is really a form of subtraction, so there is no point is subtracting twice. When the computer compares, it really subtracts and throws away the result, but notes if the result is positive or negative. If the result is positive then the first number is greater, otherwise it is lesser. (The case of zero is irrelevant here.) It would make more sense to subtract the two explicitly, save the result in a temporary variable, compare it with zero (an elementary machine instruction) and if the rotation is appropriate, update the angle to the temporary. Of course, this amounts to test-rotating the angle. The problem is worse for a rectangular coordinates, since rotation requires more than just a single subtraction. As we have seen, it requires two additions and two shifts. To do that twice would be terrible, so we have already resigned ourselves to using a temporary variable (actually two) in a true test-rotation.
But it is not necessary to perform a test rotation and then selectively apply it. It would be equally appropriate to unconditionally rotate the angle by ½ a, then compare the result with zero. If this new angle were negative, it could be rotated by −½ a (cancelling out the first rotation), and if positive, by ½ a (resulting in a rotation of a). This may seem a bit redundant, but it achieve two things: most importantly it allowed the a or −a type of rotation and the pre-computation of K, and it also simplifies the test so that the angle is simply compared with zero rather than test-rotated. But the difference is that we are now dealing with ½ angles. This does not turn out to be a problem since for any particular angle in the sequence, the next angle is about ½ the current angle. It is not exactly ½, but slightly larger. That only degrades slightly the rate of convergence of the angle to zero. The most interesting effect is that any particular angle in the sequence acts very much like its predecessor, which is fine if the initial angle in the sequence is examined carefully. But in fact, since for each angle we will now be unconditionally applying a rotation of ½ a (and each angle is twice its successor), the net effect of a rotation by a and all its successors is nearly an unconditional rotation by a, and ultimately an unconditional rotation by the initial angle in the sequence. But if we are dealing with original angles that are always positive (and we shall stipulate this), then the initial rotation by a will always be carried out. So no significant modifications are needed to use the a versus −a choice, in fact, the only modification is the benefit of reducing the test to a comparison with zero.
Still the initial angle in the sequence must be understood to be at least half of the largest original angle we wish to convert. It is necessary to limit any original angle to be less than one-quarter circle, so the initial angle of the sequence can be one-eigth circle. So the tangent of a is one, and we know the a = arc tangent of 2−i. Thus 2−i is one, for the initial angle. So the initial i is zero. To compute the value of K, recall that it is the product of all cos a = (2−2i + 1)−½. So K is the product of (2−2i + 1)−½ for i = 0 to n. Again, n is the number of bits of precision required. If n where zero (which is silly), then K = 2−½ = 0.707106... and as n grows larger, K quickly converges to 0.607252... But the value of K should be computed according to n, not by relying on a massive calculation with a huge value of n. K is used to sort of counteract the effects of the rectangular rotations using the additions and shifts, which actually cause the rectangular coordinate to grow larger. To be precise, the purpose of K is to scale the value back to where it ought to be. K is no fundamental constant of nature.
I have deliberately neglected to say very much about how the polar angle is measured, although I have used radians. Actually, it doesn't much matter how the angle is measured. K is not dependent on an angle measure. It is only necessary that the table of arc tangent values and the polar angle handed to ptor and expected from rtop use the same units. That suggests that a pointer to a table could be supplied as an argument to these routines so that they would be completely independent of an method of angle measurement. This is not quite adequate since we mean to deal with angles beyond the range zero to one quarter of circle, at some point. Still, it would not be hard to get that information from the user supplied table. But this is a feature only a computer scientist could appreciate, and wouldn't you be disappointed if I hadn't brought it up? There is no sense in mucking around with tables and I would recommend sticking with radians. I only avoided talking about radians because they tend to suggest magic. They also inspire magic and it took me a long time reading mathematical descriptions of this algorithm to discover that it's unnecessary. If you examine the algorithm you can see that the polar angle only interacts with the table of arc tangents, and the argument of an arc tangent is in the land of rectangular coordinates, so that issue goes no farther.
The other system of units is that of r, x and y. The algorithm has, of course, no interest in these units since any system of polar and corresponding rectangular coordinates could be scaled up or down with no effect. If you have any doubts, get a piece of each type of graph paper, line them up one on top of the other, hold them up to the light and move your arm back and forth. If any of the points move, consult a physician. Of course, r, x and y must all be used consistently in the same units of distance. But be careful of x and y values near the upper bound of representable values because the corresponding r may be too large. If both x and y are at some upper bound, call it z, then the corresponding r is z times the square root of two which is about 1.4. It is probably simpler to assume this number is 2; in other words, that r will require an extra bit.
It is good idea to take care not to mess with the upper bounds, and so a modification is recommended regarding the multiplication by K. In the case of converting from polar to rectangular, we took the original r and multiplied by K. Since K is smaller than one, r becomes smaller. In fact, K is less than 2−½, so the result no longer requires the extra bit. However, in converting from rectangular to polar, if the multiplication by K is postponed, then x could grow into the extra bit. This should be prevented, at the cost of an extra multiply, by multiplying each of x and y by K at the outset, instead of multiplying r by K at the end.
As an historical note and so as not to appear overly contemptuous of mathematicians, it should be pointed out that the basis of this method is rooted in a technique for the computation of logarithms and exponentials that is very similar, mathematically, to this one. A table of logarithms of 1 + 2−i is stored as constants. The decision is of the original, a versus zero variety, rather than the a versus −a that we eventually adopted. It is interesting to note that the method we used produces an oscillating convergence while the other method causes a simple, asymtotic convergence from one side only and that these strongly resemble the character of trignometric functions, which oscillate, and exponential functions which have a monotonic character, respectively. The logarithm and exponential technique does not have any equivalent to the constant K. The issues of precision and range are also problematic. It is a virtually intrinsic property of logarithms and exponentials that the range of input versus output values is grossly different. Still, these functions share with the trignometric the clever trick of looking up and subtracting values read from a table and multiplying other values by adding and shifting.
The most interesting property, however, is that the trignometric functions can, more or less, be expressed in the same form as the exponential because of Euler's formula in which sines and cosines are related to exponentials in the field of complex numbers. The field of complex numbers, itself, is analogous to the two dimensional space in which we discussed coordinate conversion. Actually the field is identical but since it is clearly not necessary to resort to complex numbers in describing the transform I hesitate to suggest that the invocation of Euler's formula is useful. But since the trignometric can be described in terms of the two-dimensional complex field, while the exponential can be confined to the one-dimensional real numbers, the graphical description of the exponential and logarithm functions can probably be shown using a slide rule as a prop.
The value of the logarithm and exponential functions in the implementation of trignometric functions would be to provide a square root function: if y = x½, then y = exp ln x½ or y = exp ½ ln x. Multiplication by ½ of course, would be implemented as a single right shift. The square root may be computed using a simulation of the long division method without encountering range problems. The method is not the best, but requires no multiplication.
The transformation of angles in the range zero to one quarter circle to the full range is simple. It requires two flags, call them xs (for the sign of x) and ys (for the sign of y). When converting to polar, the rectangular coordinates are compared with zero and the result stored in these flags. The coordinates are then replaced by their absolute values for the interesting part of the algorithm. Upon completion, r is in its final form. The angle is corrected, if x was negative, by subtracting it from π (or the measure of one half circle). This reflects the angle over the y axis. Then if y was negative, the angle is subtracted from 2π (or the measure of one full circle). This reflects it over the x axis.
To convert to rectangular, the original polar angle must be corrected. If it is at least π (one half circle), then ys is set and the angle is subtracted from 2π (one full circle). If the angle is at least π/2 (one quarter circle), then xs is set and the angle is subtracted from π. Upon completion of the conversion, the signs of x and y are made negative if xs and/or ys (respectively) was set.
To complete the job, it is a good idea to allow for the special cases. If an input r is zero, then the rectangular coordinates output are both zero. If an input x and y are zero, the polar result is zero. Finally, input polar coordinates can be corrected for negative r (by adding π to the angle and using a positive r), angles beyond one full circle (positive or negative) (by modulo dividing a by 2π) and negative angles (by adding 2π).
For the square root routine, it is reasonable to return a value with the same sign as the input. The sign is stored in a flag while the operation is performed on the absolute value of the input. The sign is applied to the result.