| |
I received the following comments and code from Mike Ash, whose Lugnet
'connection' appears not to be working at present. Big thanks Mike. I've yet
to try it, but it looks OK in theory.
**************************************************
On Sun, 12 Nov 2000, mike mcfarlane wrote:
> Hi
> What math functions are available under legOS? I can get the usual +-/ and *
> to work, but not powers or square roots. Are these math functions available
> and what others e.g. sine, cosine etc?
> I thought it would be standard C stuff like pow and sqrt, but they don't
> appear to work and I can't find a header file with any math functions in.
> Not that I like the stuff or anything ;-)just need to use it.
Most math functions aren't available. They would require so much computing
power, and so few people would use them, that they probably aren't worth
the code space in the kernel. You could probably find some sample
implementations of pow and sqrt floating around and just put them in your
program. Same for trig functions, I guess. I have a library of fixed-point
math routines that includes all trig and inverse-trig functions along with
it. It works great in fixed-point, and could probably be converted to
floating if that's what you're looking for. If you're interested, just let
me know, and I'll be happy to send it along.
Hmm, while we're on the subject, is anybody ELSE interested in a library
of fixed-point math routines? They seem to take up significantly less code
space than the 2k needed for the floating-point routines, and they're
probably faster to boot. I use them all over in my program, and I'd be
more than happy to share.
--
"Say not, 'I have found the truth,' but rather, 'I have found a truth.'
"Say not, 'I have found the path of the soul.' Say rather, 'I have met the
soul walking upon my path.'" -- Khalil Gibran
****************************************************************
I'm not sure if this will make it to the list, but we'll see. I sent out a
message twice and it still hasn't shown up.
Anyway, I have uploaded my fixed-point math routines to my web site. This
includes simple stuff like multiplication and division, but also a full
set of trigonometry routines. Even if you'd rather use floating-point, you
may find the technique for the trig functions valuable, since they're
pretty fast and reasonably accurate, and wouldn't be difficult to convert
to floating.
The URLs are:
http://www.mikeash.com/lego/fixedpoint.h
http://www.mikeash.com/lego/fixedpoint.c
The code is fairly well-explained in the comments (I hope!). I don't have
any additional write-up besides them, though. If you have any questions,
please feel free to e-mail me.
I haven't done any speed tests, but I don't expect any significant speed
differences. 32-bit integer math isn't native to the RCX CPU, and there's
extra bitshifting work that has to be done as well, so my guess is that it
won't be much faster than floating-point. However, it does seem to take up
much less code space, which I personally find very valuable.
Comments, suggestions, criticisms, and offers of money are all welcome. :)
--
"Say not, 'I have found the truth,' but rather, 'I have found a truth.'
"Say not, 'I have found the path of the soul.' Say rather, 'I have met the
soul walking upon my path.'" -- Khalil Gibran
****************************************************************
#include "fixedpoint.h"
/* The following tables are the sine and cosine of angles, starting at 90,
and decreasing by powers of two. So entry 0 is for 90 degrees, entry 1
is for 45, entry 2 is for 22.5, etc.
For the reasoning, see RotateVectorInternal and UnRotateVectorInternal.
These tables are only this long because any further entries into the
sine table will go under 256. Because of how the multiplication
routine is done, multiplying by a number less than 256 (which is actually
0.0039) will result in multiplying by zero, which is probably not a good
thing. At the very least, there's no point in continuing the table further.
*/
const long powSinTable[] =
{65536, 46340, 25079, 12785, 6423, 3215, 1608, 804, 402};
const long powCosTable[] =
{0, 46340, 60547, 64276, 65220, 65457, 65516, 65531, 65534};
void RotateUnitVector(longFix theta, longFix *x, longFix *y);
void RotateVectorInternal(longFix theta, longFix *x, longFix *y);
longFix UnRotateUnitVectorInternal(longFix x, longFix y);
longFix UnRotateVectorInternal(longFix x, longFix y);
longFix longFixMulDirty(longFix n1, longFix n2) {
return (n1 >> 8) * (n2 >> 8);
}
longFix longFixDivDirty(longFix n1, longFix n2) {
return (n1) / (n2 >> 16);
}
longFix longFixMul(longFix n1, longFix n2) {
return (n1 >> 8) * (n2 >> 8);
}
longFix longFixDiv(longFix n1, longFix n2) {
return longFixDivDirty(n1, n2);
if((n1 > 0L && n2 > 0L) || (n1 < 0L && n2 > 0L))
return ((abs(n1)/abs(n2)) << 16) | ((abs(n1)/(abs(n2)>>16)) & 0x0000FFFF);
else
return -(((abs(n1)/abs(n2)) << 16) | ((abs(n1)/(abs(n2)>>16)) &
0x0000FFFF));
/* MENTAL NOTE: if calculations involving division start going screwey,
check here. This routine was not exactly carefully validated.
Basically, I came up with some weird equation that produced the
above, and I don't remember how it works anymore. If anybody
could validate (or invalidate!) the above, I would be highly
grateful.
*/
}
/* Trig functions are basically ratios along a unit circle,
so we just rotate unit vectors and do things with the
results. */
longFix longFixSin(longFix theta) {
longFix x=0, y=0;
RotateUnitVector(theta, &x, &y);
return y;
}
longFix longFixCos(longFix theta) {
longFix x=0, y=0;
RotateUnitVector(theta, &x, &y);
return x;
}
longFix longFixTan(longFix theta) {
longFix x=0, y=0;
RotateUnitVector(theta, &x, &y);
if(x != 0)
return y / x;
else
return 0x7FFFFFFF;
/* If x is zero, the result should be infinity. Unfortunately
we can't represent infinity, so this is the next best thing. */
}
/* This function will take any vector and rotate it so that it's pointing
directly to the right (its angle is zero). It will then return the angle
that it needed to be rotated by. The result will be in [0, twopi].
*/
longFix UnRotateVector(longFix x, longFix y) {
longFix theta;
/* UnRotateVectorInternal returns an angle in [-pi, pi], so
it needs to be moved into the range that's expected. */
theta = UnRotateVectorInternal(x, y);
if(theta < 0)
theta += twopi;
return theta;
}
/* This function takes the parameters and simply normalizes theta to the
range of
[-pi, pi] that the internal function expects. */
void RotateVector(longFix theta, longFix *x, longFix *y) {
theta %= twopi;
if(theta > pi)
theta -= twopi;
/*else if(theta < -pi); <-- DO NOT PROGRAM WHILE ON CRACK */
/* The above bug was in this code for a long time, and it
complicated debugging immensely. Just a warning for all
of you. */
else if(theta < -pi)
theta += twopi;
RotateVectorInternal(theta, x, y);
}
void RotateUnitVector(longFix theta, longFix *x, longFix *y) {
*x = 1; *y = 0;
RotateVector(theta, x, y);
}
/* RotateVector will rotate an arbitrary vector about the origin to an
angle between
-pi and pi. The vector is passed into the function with x and y, and
the result is stored in the same place.
Calling this function with a value of theta outside the specificed range
gives
undefined behavior. With this rev of the function, it will simply rotate
pi or
-pi degrees, depending on whether the angle is positive or negative.
If you're wondering how this routine works, here's a quick summary.
Rotating a point about the origin by a specific number of degrees
is pretty trivial. x = x*cos - y*sin, y = x*sin + y*cos
This, of course, assumes you know the cosine and sine of the desired
angle. But the whole reason I wrote this routine was to find the
sine and cosine of a given angle! (Rotate a unit vector about the
origin by a specified angle, and the sine, cosine, tangent, or what-
have-you can be derived easily.) I could write a sine and cosine
routine using a maclaurin series, but that's slow and painful. I
could make a table, but I only have 32k! So, a clever guy (If Rog
from #macdev ever reads this, thanks again!) showed me the light.
Do successive rotations, back and forth. Store a small table, with
angles decreasing by half, and rotate left or right, hunting back and
forth, until you're finally where you need to be.
So this is what this function does. We have sine and cosine tables
above (notice how they're very short, basically insignificant). Now
we loop through the table entries. If we need to rotate to the left,
rotate the current vector to the left. If to the right, rotate to the
right. Once we either hit it right on, or run out of table, leave
the function.
*/
void RotateVectorInternal(longFix theta, longFix *x, longFix *y) {
longFix curAngle, change, newX, newY;
int i;
curAngle = 0; change = pi/2;
for(i = 0; i < sizeof(powSinTable)/sizeof(long); i++) {
if(curAngle > theta) {
newX = longFixMul(*x, powCosTable[i]) - longFixMul(*y, -powSinTable[i]);
newY = longFixMul(*x, -powSinTable[i]) + longFixMul(*y, powCosTable[i]);
curAngle -= change;
}
else if(curAngle < theta) {
newX = longFixMul(*x, powCosTable[i]) - longFixMul(*y, powSinTable[i]);
newY = longFixMul(*x, powSinTable[i]) + longFixMul(*y, powCosTable[i]);
curAngle += change;
}
break;
*x = newX;
*y = newY;
change /= 2;
}
}
/* UnRotateVectorInternal does the opposite of RotateVector (thus the
name). Rotate
vector takes a vector and an angle, and rotates the vector by that angle.
UnRotateVector takes a vector and rotates it until its y coordinate is zero
(i.e. the angle that the vector forms is zero), and then returns that angle.
This will presumably be useful for inverse trig functions.
Inquiries should generally be sent to UnRotateVector, which will return
a more amiable number.
This routine works in the oppsite manner of the function above. That one
rotates
left or right depending on whether the current angle is left or right of
the desired
angle. This one rotates left or right depending on whether the current y value
is above or below the x axis. Once it hits the x axis, we know what the angle
of rotation is and can return it.
*/
longFix UnRotateVectorInternal(longFix xIn, longFix yIn) {
longFix curAngle, change, newX, newY, x, y;
int i;
x = xIn; y = yIn;
curAngle = 0; change = pi/2;
/* Just a note, I had i < sizeof(powSinTable) here, without the
division. That caused problems (four times too many iterations!).
Beware. */
for(i = 0; i < sizeof(powSinTable)/sizeof(long); i++) {
if(y > 0) {
newX = longFixMul(x, powCosTable[i]) + longFixMul(y, powSinTable[i]);
newY = -longFixMul(x, powSinTable[i]) + longFixMul(y, powCosTable[i]);
curAngle += change;
}
else if(y < 0) {
newX = longFixMul(x, powCosTable[i]) - longFixMul(y, powSinTable[i]);
newY = longFixMul(x, powSinTable[i]) + longFixMul(y, powCosTable[i]);
curAngle -= change;
}
else
break;
x = newX;
y = newY;
change /= 2;
}
return curAngle;
}
****************************************************************
#ifndef FIXEDPOINT_H
#define FIXEDPOINT_H
#define pi ((longFix)205887)
#define twopi ((longFix)411775)
/* 2pi turns out to be odd, so it's slightly more
accurate (not noticeable, most likely) to say
twopi than it is to say 2*pi. */
#define abs(x) ((x) < 0 ? -(x) : (x))
/* myDegrees = longFixMul(myRadians, rad2deg);
myRadians = longFixDiv(myDegrees, rad2deg);
*/
#define rad2deg ((longFix)3754936)
/* longFix is just a normal long with sixteen bits
as the fractional, so simple shifting is all
that's needed for the conversion. */
#define longFix2int(x) (((long)(x)) >> 16)
#define int2longFix(x) (((long)(x)) << 16)
typedef long longFix;
longFix longFixMulDirty(longFix, longFix);
longFix longFixDivDirty(longFix, longFix);
/* These first two are less accurate but faster.
Currently multiplication is the same, and the division
routine will likely cause you trouble, so don't use these.
*/
longFix longFixMul(longFix, longFix);
longFix longFixDiv(longFix, longFix);
/* For those of you who have not used fixed-point numbers:
Addition and subtraction between fixed point numbers works
as usual. MyFixed1 = MyFixed2 - MyFixed3. Multiplication
and division is different; use the routines provided.
If you want to add or subtract with a regular integer,
convert the integer to longFix before adding or subtracting.
If you want to multiply or divide with a regular integer,
just go ahead and do it, it works. TwiceAsMuch = 2*SomeLongFix.
If you're using floating point (lord knows why you'd use both
at once with the limited memory space in the RCX, but you never
know), then you can convert back and forth by multiplying and
dividing by 65536. MyFloat = (float)myFixed/65536.0.
MyFixed = MyFloat*65535.
*/
longFix longFixSin(longFix);
longFix longFixCos(longFix);
longFix longFixTan(longFix);
/* First parameter is an angle, the next two are references
to the x and y coords of a vector. The routine will rotate
the vector by the angle. */
void RotateVector(longFix, longFix *, longFix *);
/* The parameters are the x and y coords of a vector. The routine
will rotate the vector until it lies on the positive x axis
and then return how far it had to be rotated. This will basically
give you the angle of the vector. This routine can typically take
the place of most inverse trig functions, and can be used to make
real inverse trig functions if you need them. */
longFix UnRotateVector(longFix, longFix);
#endif
|
|
Message is in Reply To:
 | | Re: Math
|
| Mike, You may have to look at writing your own. I know there's a libmint (integer math) library, but I'm currently at work, and can't check what else is there. Maybe need a libm to put all math functins in?? ROSCO mike mcfarlane (...) (24 years ago, 12-Nov-00, to lugnet.robotics.rcx.legos)
|
4 Messages in This Thread:
  
- Entire Thread on One Page:
- Nested:
All | Brief | Compact | Dots
Linear:
All | Brief | Compact
|
|
|
|
