Finding the Perfect Hyperbola February 15th, 2010
Patrick Stein

For an application that I’m working on, I needed a way to scale quantities that range (theoretically) over the real numbers (though practically are probably between plus and minus three) into positive numbers. I wanted the function to be everywhere increasing, I wanted f(0) = 1, and I wanted control of the derivative at x = 0.

The easy choice is: f(x) = e^{\alpha x}. This is monotonically increasing. f(0) = 1 and f^\prime(0) = \alpha.

I needed to scale three such quantities and mush them together. I thought it’d be spiffy then to have three different functions that satisfy my criteria. The next logical choice was f(x) = 1 + \mathrm{tanh} (\alpha x). It is everywhere positive and increasing. And, it has f^\prime(0) = \alpha.

Now, I needed third function that was always positive, always increasing, had f(0) = 1 and f^\prime(0) = \alpha. One choice was: f(x) = e^{e^{\alpha x} - 1}. But, that seemed like overkill. It also meant that I really had to keep my \alpha tiny if I didn’t want to scale things into the stratosphere.

Playing with hyperbolas

So, I thought… why don’t I make a hyperbola, rotate it, and shift it so that the apex of one side of the hyperbola is at (0,1). And, I can adjust the parameters of the hyperbola so that f'(0) = \alpha. After a variety of false starts where I tried to keep the hyperbola general until the very end (\frac{x^2}{a^2} - \frac{y^2}{b^2} = r^2, rotated by \theta degrees, and shifted by \beta), I quickly got bogged down in six or seven incredibly ugly equations in eight or nine variables.

So, it was time to start trying to make it easy from the beginning. I noticed that if was going to rotate it by an angle \theta in the clockwise direction, then I needed \phi = \frac{\pi}{2} - \theta to be such that \tan \phi = \alpha if my slope was going to work out right in the end. So, I’m looking at the basic triangle on the right then to determine all of my sines and cosines and such.

Based on that triangle, it was also obvious that the asymptote for my starting hyperbola had to have \frac{b}{a} = \frac{1}{\alpha}. I played around a bit then with making a = \alpha and b = 1. In the end, I found things simplified sooner if I started with a = 1 and b = \frac{1}{\alpha}.

I also needed r to be such that the point (-r,0) rotate up so that its y-coordinate was 1. This meant that r \sin \theta = 1 or r = \frac{1}{\sin \theta} = \sqrt{1 + \alpha^2}.

So, my starting hyperbola then was: x^2 - \alpha^2 y^2 = 1 + \alpha^2.

From there, I had to rotate the x and y by \theta in the clockwise direction. This gave me:

(x\cos\theta - y\sin\theta)^2 - \alpha^2 (x\sin\theta + y\cos\theta)^2 = 1 + \alpha^2

A little multiplying out leads to:

(x^2\cos^2\theta - 2xy\sin\theta\cos\theta + y^2\sin^2\theta) \\ - \alpha^2(x^2\sin^2\theta + 2xy\sin\theta\cos\theta + y^2\cos^2\theta) = 1 + \alpha^2

From there, using cos^2\theta = \frac{\alpha^2}{1 + \alpha^2}, sin^2\theta = \frac{1}{1 + \alpha^2}, and \sin\theta\cos\theta = \frac{\alpha}{1 + \alpha^2}, we come to:

\frac{-2\alpha(1 + \alpha^2)xy + ( 1 - \alpha^4 )y^2}{1 + \alpha^2} = -2\alpha{}xy + (1 - \alpha^2)y^2 = 1 + \alpha^2

The only step remaining was to shift it all over so that when y = 1, we end up with x = 0. Plugging in y = 1, we see that -2x\alpha + 1 - \alpha^2 = 1 + \alpha. That equation balances when x = -\alpha. We want it to balance when x = 0, so we’re going to substitute (x - \alpha) in place of the x in the above equation to get:

-2\alpha(x - \alpha)y + (1 - \alpha^2)y^2 = 1 + \alpha^2

We can easily verify that the point (0,1) is on the curve. And, we can implicitly differentiate the above to get:

-2\alpha(x - \alpha)\frac{dy}{dx} - 2\alpha{}y + 2(1 - \alpha^2)y\frac{dy}{dx} = 0

Plopping x = 0 and y = 1 in there, we find that \frac{dy}{dx} = \alpha.

This is pretty good as it goes. The only step is to complete the square to find a nicer expression for y in terms of x. We start by adding \left[ \frac{\alpha}{\sqrt{1 - \alpha^2}} (x - \alpha) \right]^2 to both sides to get:

\left[ \sqrt(1 - \alpha^2)y - \frac{\alpha}{\sqrt{1 - \alpha^2}(x - \alpha)} \right]^2 = 1 + \alpha^2 + \frac{\alpha}{1 - \alpha^2} (x-\alpha)^2

This is easy enough to solve for y by taking the square root of both sides and shuffling some things about:

y = \frac{\alpha}{1 - \alpha^2}(x - \alpha) + \frac{1}{\sqrt{1 - \alpha^2}} \sqrt{1 + \alpha^2 + \frac{\alpha^2}{1 - \alpha^2}(x - \alpha)^2}

Here are all three curves with \alpha = \frac{1}{4}. The exponential is in black, the hyperbolic tangent is in red, and the hyperbola is in blue:

The first image on the page here was made with Zach’s Vecto library with some post-processing in the GIMP. (Here is the source file: hyperbola.lisp.) The second image was made entirely within the GIMP. And, the last image was made using Foo Plot and the GIMP.

Casting to Integers Considered Harmful August 6th, 2009
Patrick Stein

Background

Many years back, I wrote some ambient music generation code. The basic structure of the code is this: Take one queen and twenty or so drones in a thirty-two dimensional space. Give them each random positions and velocities. Limit the velocity and acceleration of the queen more than you limit the same for the drones. Now, select some point at random for the queen to target. Have the queen accelerate toward that target. Have the drones accelerate toward the queen. Use the average distance from the drones to the queens in the i-th dimension as the volume of the i-th note where the notes are logarithmically spaced across one octave. Clip negative volumes to zero. Every so often, or when the queen gets close to the target, give the queen a new target.

It makes for some interesting ambient noise that sounds a bit like movie space noises where the lumbering enemy battleship is looming in orbit as its center portion spins to create artificial gravity within.

I started working on an iPhone application based on this code. The original code was in C++. The conversion to Objective C was fairly straightforward and fairly painless (as I used the opportunity to try to correct my own faults by breaking things out into separate functions more often).

Visualization troubles

uniform
The original code though chose random positions and velocities from uniform distributions. The iPhone app is going to involve visualization as well as auralization. The picture at the right here is a plot of five thousand points with each coordinate selected from a uniform distribution with range [-20,+20]. Because each axis value is chosen independently, it looks very unnatural.

gauss
What to do? The obvious answer is to use Gaussian random variables instead of uniform ones. The picture at the right here is five thousand points with each coordinate selected from a Gaussian distribution with a standard-deviation of 10. As you can see, this is much more natural looking.

How did I generate the Gaussians?

I have usually used the Box-Muller method of generating two Gaussian-distributed random variables given two uniformly-distributed random variables:

(defun random-gaussian ()
  (let ((u1 (random 1.0))
        (u2 (random 1.0)))
    (let ((mag (sqrt (* -2.0 (log u1))))
          (ang (* 2.0 pi u2)))
      (values (* mag (cos ang))
              (* mag (sin ang))))))

But, I found an article online that shows a more numerically stable version:

(defun random-gaussian ()
  (flet ((pick-in-circle ()
           (loop as u1 = (random 1.0)
                as u2 = (random 1.0)
                as mag-squared = (+ (* u1 u1) (* u2 u2))
                when (< mag-squared 1.0)
                return (values u1 u2 mag-squared))))
    (multiple-value-bind (u1 u2 mag-squared) (pick-in-circle)
      (let ((ww (sqrt (/ (* -2.0 (log mag-squared)) mag-squared))))
        (values (* u1 ww)
                (* u2 ww))))))

For a quick sanity check, I thought, let’s just make sure it looks like a Gaussian. Here, I showed the code in Lisp, but the original code was in Objective-C. I figured, If I just change the function declaration, I can plop this into a short C program, run a few thousand trials into some histogram buckets, and see what I get.

The trouble with zero

So, here comes the problem with zero. I had the following main loop:

#define BUCKET_COUNT 33
#define STDDEV       8.0
#define ITERATIONS   100000
 
  for ( ii=0; ii < ITERATIONS; ++ii ) {
    int bb = val_to_bucket( STDDEV * gaussian() );
    if ( 0 <= bb && bb < BUCKET_COUNT ) {
      ++buckets[ bb ];
    }
  }

I now present you with three different implementations of the val_to_bucket() function.

int val_to_bucket( double _val ) {
  return (int)_val + ( BUCKET_COUNT / 2 );
}
 
int val_to_bucket( double _val ) {
  return (int)( _val + (int)( BUCKET_COUNT / 2 ) );
}
 
int val_to_bucket( double _val ) {
  return (int)( _val + (int)( BUCKET_COUNT / 2 ) + 1 ) - 1;
}

As you can probably guess, after years or reading trick questions, only the last one actually works as far as my main loop is concerned. Why? Every number between -1 and +1 becomes zero when you cast the double to an integer. That’s twice as big a range as any other integer gets. So, for the first implementation, the middle bucket has about twice as many things in it as it should. For the second implementation, the first bucket has more things in it than it should. For the final implementation, the non-existent bucket before the first one is the overloaded bucket. In the end, I used this implementation instead so that I wouldn’t even bias non-existent buckets:

int val_to_bucket( double _val ) {
  return (int)lround(_val) + ( BUCKET_COUNT / 2 );
}

Clifford Algebras for Rotating, Scaling, and Translating Space July 6th, 2009
Patrick Stein

In (very much) earlier articles, I described:

Today, it is time to tackle rotating, translating, and scaling three-dimensional space using Clifford algebras.

Three dimensions now instead of two

Back when we used Clifford algebras to rotate, translate, and scale the plane, we were using the two-dimesional Clifford algebra. With the two-dimensional Clifford algebra, we represented two-dimensional coordinates (x,y) as xe_1 + ye_2. It shouldn’t surprise you then to find we’re going to represent three-dimensional coordinates (x,y,z) as xe_1 + ye_2 + ze_3.

As before, we will have e_1e_1 = 1 and e_2e_2 = 1. Similarly, we will have e_3e_3 = 1. In the two-dimesional case, we showed that e_1e_2 = -e_2e_1. By the same logic as the two-dimensional case, we also find that e_1e_3 = -e_3e_1 and e_2e_3 = - e_3e_2. We could potentially also end up multiplying e_1, e_2, and e_3 all together. This isn’t going to be equal to any combination of the other things we’ve seen so we’ll just leave it written e_1e_2e_3.

Read the rest of this entry ⇒

Quaternions for Rotating, Scaling, and Translating Space June 11th, 2009
Patrick Stein

In earlier posts, I described how complex numbers can be used to rotate, scale, and translate the plane, how Clifford algebras can be used to rotate, scale, and translate the plane, and why I resorted to an awkward trick for the Clifford algebra rotations of the plane. In this post, I am going to explain what the quaternions are and describe how they can be used to represent a rotation in three-dimensional space.

What are the quaternions

Okay, remember how we got the complex numbers? We needed something that was the square root of negative one.

Now, imagine that you are Sir William Rowan Hamilton. The year is 1843. It is springtime. You know how to use the complex numbers to represent points in the plane. And, you know that when you do that, you can use complex numbers to rotate, scale, and translate the points. That’s all well and good, but you don’t live in a two-dimensional world. How are you going to do the same sort of thing with three-dimensional space? How are you going to multiply triples?

You spend months on this. If only you could say, How about I let there be another number that is different from i (and from -i) that has the same property that its square is negative one? You fight with this for months. You try to represent a point with coordinates (x,y,z) as x + yi + zj. But, nothing you come up with makes any sense.

Your kids are harassing you, Daddy, did you figure out how to multiply triples yet? You have to answer them every morning with a polite, No, not yet.

Then, you’re walking along the Royal Canal in Dublin. It’s mid-October already. My, how the year has flown by. Bam, it hits you. If you add a third number like i and j which is equal to i\cdot j, everything works out. You get so excited, that you carve your equations into a stone bridge over the canal:

i^2 = j^2 = k^2 = ijk = -1

Read the rest of this entry ⇒

What Was Up With That Rotation Trick? June 10th, 2009
Patrick Stein

In my prior post about using Clifford algebras to do plane rotations, I finished with a non-intuitive step at the end. Rather than multiplying on the right by an element representing a rotation of angle \theta, I multiplied on the left by an element representing a rotation of angle \frac{\theta}{2} and multiplied on the right by an element representing a rotation of angle -\frac{\theta}{2}.

Why did I do this? Well, I mentioned it would be awkward for the two-dimensional case, but that it will be important when we get to three or more dimensions. Well, work for a moment with \frac{\theta}{2} being a quarter rotation (ninety degrees, \frac{\pi}{2} radians). This means our total rotation is going to be a half turn (180 degrees, \pi radians).

For that \frac{\theta}{2}, r = e_1e_2 and so \overline{r} = -e_1e_2. Let’s just look at what it does to our unit vectors e_1 and e_2 to multiply on the left by r and on the right by \overline{r}.

For e_1, we get -e_1e_2e_1e_1e_2 = -e_1e_2e_2 = -e_1. Similarly, for e_2, we get -e_1e_2e_2e_1e_2 = -e_2.

So far, we were only working in two dimensions. As such, there wasn’t any e_3 to worry about. But, what if there were? What happens to the z-coordinate of something if you rotate things parallel to the xy-plane? It remains unchanged.

Well, what would happen if we multiplied e_3 on the right by \cos\theta + \sin\theta e_1e_2? We would end up with \cos\theta e_3 + \sin\theta e_3e_1e_2 = \cos\theta e_3 + \sin\theta e_1e_2e_3. We’ve ended up scaling e_3 and adding in a trivector e_1e_2e_3. We’ve made a mess.

Let’s try it instead with our trick. We’re going to start with -e_1e_2e_3e_1e_2. Every time we transpose elements with different subscripts, we flip the sign. Every time we get two elements next to each other with the same subscript, they cancel out. So, switching the e_3 with the second e_1, we get e_1e_2e_1e_3e_2. From there, we can switch the first two elements to get -e_2e_1e_1e_3e_2 which is just -e_2e_3e_2. We can switch the e_3 with the second e_2 to get: e_2e_2e_3 which is just e_3. So, our trick leaves e_3 unchanged.

In the above, there is nothing special about the subscript three. It would work for any subscript except one or two. So, the trick allows us to break the rotation up into two parts that still do what we want with e_1 and e_2 but leave our other directions unchanged (or, maybe it’s easier to think of them as changing them and then changing them right back).

Updates In Email

Email:

RSS Feeds

l