Emulating Recursion on the TI-84+

November 30, 2008

Yes… it’s true I think about these things in my free time. It actually happened in the middle of a math class I was taking. We were talking about the Fibonacci sequence and I was thinking that it would be neat to have a recursive implementation on my calculator. Now, I know that the iterative approach is actually faster, but the idea of a true recursive implementation was a fun idea. Now, it may be possible using assembly programming for the calculator, but I don’t really know how to do that.

After a little bit of thinking I decided that programming with the calculator BASIC reminded me a bit of MIPS assembly (I can hear some people groan already). You see the calculator lists can be used as a make shift stack and since all the variables A – Z are global you can simply designate one of those as your “stack pointer”. So I adopted the convention that the S variable is the stack pointer and that the R variable is for returning values. Anything else is fair game. I suppose if I really wanted to put thought into this I would pick some of them as temporary and others for function arguments. However I didn’t want to spend that much time working on this.

The basic idea here is that when you want to call a function (IE another program) at the beginning of your function you simply preserve the state of any variables you wish to change by using the list. I happen to use L5, but you could use any list you wanted. The stack pointer starts at 0 and each time you add to the list to increment the stack pointer by the number of items you added. Then at the end of your function you restore the state of any variables you used and “return” by setting the R variable to be the return value.

Now this has some limitations. First of all the TI-BASIC is quite slow even on my 84+ Silver Edition. The second is that a list can only hold 1000 elements so your stack is essentially size 1000. That being said it is still fun to play around with the idea of recursion or the more broad idea of function calls. For instance I wrote a program to do prime factorization which uses a helper program to tell if a number is prime or not.

Here is the code of the main program FIB which is the recursive Fibonacci number finder:

:ClrList L6
:Prompt A
:1 => S
:Disp R

Now here is the program which does the real work, FIBR:

:If A<=2
:1 => R
:A => L6(S)
:B => L6(S+1)
:T => L6(S+2)
:"Increase stack pointer
:S+3 => S
:A => B
:0 => T
:A-1 => A
:T+R => T
:A-1 => A
:T+R => T
:T => R
:L6(S-3) => A
:L6(S-2) => B
:L6(S-1) => T
:S-3 => S

If you actually went through the trouble of putting this into your calculator you run the program FIB and enter a positive integer. The program will then return the nth fibanocci number. This is a simple program but I think it gets at the idea of recursion. In any case it doesn’t have too practical a purpose it’s just for fun 🙂

Euler Project

November 20, 2008

As pointed out by Wes it has indeed been some time since I posted anything on this blog. About once a week it occurs to me that I would really enjoy writing about this or that, but I just haven’t had the time lately. Part of the reason is how long it actually takes me to get a post up. It typically takes me around 30 to 40 minutes to write something mildly interesting, and then another 10 or 20 to proof read it. It is not as if I couldn’t find that sort of time, but typically I only have it at night, and when I am tried I find it very hard to concentrate on what I am writing 🙂

In any case I want to write about Project Euler today. About a week ago my Professor introduced to to Project Euler. They have a list of computer science and math problems that they ask you to solve, and typically the problem involves you needing to have some sort of mathematical intuition as well as coding skill. One simply makes an account on their site and when you think you have a solution to the problem you submit your answer and the site tells you if you got the right answer or not. So far I have been working on the earlier problems in the list when time permits. For a lot of them I tend to hack together a solution in Python within 5 to 10 minutes. I have heard the later problems get harder, and I await the challenge. Project Euler meets the strange need to write random small computer programs that perform an operation which is “cool” to only a small group of people. I also like the mathematical challenges that some of the problems provide, and it gives me a good test of my problem solving skills, which I think are very important to any computer scientist. The other interesting part to these problems is that Project Euler guarantees that your program should find a solution in under one minute of run time. This adds an extra challenge because you can’t necessarily brute force the answer :).

While writing solutions to these problems I have noticed how much slower Python is than C in certain cases. I have read many different guides to optimizing Python code and applied optimizations everywhere I can, and still the C equivalent to my solutions  runs drastically faster. I am not suggesting that C is necessarily a better language or something. Certainly, one expects a compiled language to outperform a scripting language, but since I had never done any formal study before I found this quite interesting. I do enjoy the speed at which I can throw together the solutions in Python, and I’m glad I taught myself Python a year or two ago because I have used it many times since then.

Anyway, I have a lot more to write about, and have been keeping a list. Hopefully over Thanksgiving break I will be able to put in an entry or two. At the very least when this semester ends I will have plenty of time to write blog entries and plan to be more active during that time.

Intersection of Spherical Circles (Finally!)

July 23, 2008

For the past two weeks I have been pondering on and off how I would calculate the location of the intersection of two spherical circles. I really want Spheriosity to be as complete as possible and having the intersection tool only work on Spherical lines/line segments seems like a real lame idea. So I have finally devised an algorithm which can correctly determine the intersection points. I will give a brief overview of it here. I will leave out details about different parts of the algorithm which I feel are either trivial or would need their own blog post to explain ;).

So we start with two circles… What are the possibilities for their arrangements in terms of intersections:

  1. The most simple is if the circles are not anywhere near each other and don’t intersect.
  2. One circle could be “inside” the other. In this case the circles obviously do not intersect.

    One circle inside the other

    One circle inside the other

  3. The two circles could intersect only in one spot. This could happen in a number of ways:
    • Case 1:

      First Intersection Case

      First Intersection Case

    • Case 2:

      Intersection Case 2

      Intersection Case 2

  4. Then there is the case where the circles intersect twice. This is perhaps the most common case in terms of what users will see, but since the other cases are easy to construct we must also deal with them.

So that loosely describes my problem and the different scenarios that need to be considered. So… on to the solution! Much of this boils down to finding the intersection of two planes. Each circle on the sphere defines it’s own plane which basically slices the sphere, and the intersection of that plane and the sphere could also be used to form the circle. So the intersection of the planes formed by both circles forms a line which, if it goes through the sphere, can be used to find our points of intersection.

Defining that line of intersection turned out to be the largest problem. It was easy enough to get a vector going in the correct direction but I still needed a point to form the line and perform any sort of calculation on it. Here are the steps I take to find the intersection and make sure one actually exists:

  1. Consider two circles AB and CD (defined by center point and radius point).
  2. Find the normal vector to the planes formed by each circle. We will call these vectors a and c
  3. Take a x c if the magnitude squared of the resulting vector is 0 you can stop now because the planes are parallel, and the circles could not possibly intersect (or are coincident).
  4. Form an axis of rotation, r, between points A and C.
  5. For both circles form lines by rotating their respective center points forward and back by the angle formed between their center points, the center of the sphere, and their radius points. (call these lines l and m.
  6. Find the intersection of lines l and m and call this point I
  7. Find the intersection of the line formed by using point I and the cross product of a and c with the sphere.
  8. Consider the following cases:
    • Case 1: Line does not intersect the sphere: no intersection points
    • Case 2: Line intersects once with the sphere: one intersection point
    • Case 3: Line intersects twice with the sphere: two intersection points

There you have it! I know I left some details out, but writing out every single step would make this a really long post ;).

Spheriosty and Parallel Transport on the Sphere

July 11, 2008

Today I was writing some unit tests for Spheriosity and I discovered a flaw with the code that currently handles parallel transported lines. For those of you who are unfamiliar with the concept of parallel transport I will do my best at describing this to you.

Basically start with any two intersecting lines:

Two intersecting lines

Two intersecting lines

Now, slide the first line along the second keeping the angle between the two lines the same. You have just parallel transported your first line! By definition parallel transport merely keeps the angle between those two lines the same as they move.

Parallel transported line

Parallel transported line

Why is that so dang special? But wait! Isn’t that just a fancy way to say “keep the lines parallel”? Well, on the Euclidean plane, yes. However on a sphere there is no such thing as parallel. Two straight lines [known as great circles on a sphere] will always intersect. If you don’t believe me you can always try it out for yourself. Any true mathematician probably doesn’t believe me, as I have not provided a proof 🙂 .

Parallel transport of a great circle

Parallel transport of a great circle

Notice in the above picture that we have transported along the horizontal line. It is obvious that the original line (right line) and the transported line (left line) will intersect.

It is also interesting to note that on a sphere parallel transport is really the same as rotation! (Awesome… I know!) So where to rotate… well you can probably figure it out by simply looking at the sphere, but we want to try and put it in writing.  In order to find the axis to rotate on we have to imagine vectors going from the center of our sphere to the two end points that define our line of transport. Lets call those vectors a and b. Then, we take do: a x b (a cross b) and get some vector c. If we place vector c at the center of the sphere then we get our axis of rotation. Let call this primary axis just so we have name for it

So in the context of Spheriosity I want to take a line of transport, a line to transport and a cursor location and transport the ‘line to transport’ to where the cursor is, but making sure it is with respect to the ‘line of transport’. Translating the cursor location to a line parallel transporting is not a one step process, sadly. Since I don’t know where the cursor will be I had to find a way to map out the cursor location to the line being transported. [A big thanks to my professor who figured this one out]

  1. We use the point provided by the user and the axis of rotation for the line of transport to give us a new axis of rotation. Just to have some order here lets call the user point U, the first axis of rotation primary axis, the line of transport line AB, original line we are transporting line CD, and the new axis of rotation secondary axis
  2. Next we measure the angle between the primary axis and the vector defined from the center of the sphere to point C
  3. Now plot the point which is the intersection of the axis of rotation with the sphere. There are two poential points, so plot whichever one you wish. We will call this point E. Rotate point E with the secondary axis by the angle computed from step 2
  4. Now the first point is in place. We will call this point F. We measure the angle between the vectors defined by going from the center of the sphere to point F and point C.
  5. Using the angle from step 4 we rotate point D using the primary axis
  6. Enjoy the glory of seeing the line parallel transported

There are actually some problems with that method. The first of which is figuring out which direction to rotate point D. The .angle() function of the Vector3f class in Java3D only returns an angle from 0 to PI. The way I figure this out is actually by use of tangent vectors and dot products. To get the tangent vector I actually cheat and only approximate the angle vector by rotating the point C a very small amount (.0001 radians) with the primary axis and use the original point C and the rotated point C to make my vector (vector t). Yes, a crude trick, but it’s quick and easy. Next I define a vector from the center of the sphere to the point the cursor is at (vector c). I take the dot product of t and c,  and if the dot product is negative I know to take the angle from .angle() and rotate that in the opposite direction. Perhaps I will write a separate entry and how the dot product is my hero, but for now you will have to trust that it makes sense 🙂 .

The last issue, which I just discovered today, is that if the line we want to transport has one of it’s endpoints as the axis of rotation the aforementioned method does not work. So I simply special case this, and rotate the end point which isn’t the axis of rotation(point D) by the angle formed by the vectors made from the center of the sphere to the cursor point, and point D. Then I use the same tangent vector trick to figure out which way to rotate.

That’s all for now. Hopefully you at least understand parallel transport a little better :). Please let me know if some of that was unclear and I will do my best to make it sound better.