Chapter 8

P21Forth Math Routines


P21Forth 1.02 is based on an experimental pre-release version of a new ANS compliant version of eForth. Since eForth has only one math prmitive, UM+ , math operations in eForth tend to be slow. MuP21 only has a multiply step instruction and the slowest instructions are addition and the multiply step.

To inprove performance many other ANS standard math words in P21Forth have been rewritten as CODE words. + 1+ 2+ 1- 2- - 2/ U2/ 2* D2* * 4* 8* 16* 32* were written as CODE words to improve performace on these math operations.

Dr. Michael Montvelishsky, in Saransk, greatly improved the performance of some other math operations in P21Forth by contributing UM/MOD and UM* as CODE words. Like most operations these words perform about 50 times faster on MuP21 as CODE words than they do as COLON definitions.

Dr. Michael Montvelishsky's high performance
Math wordset for MuP21

A math wordset is being constructed for 3D CAD and virtual reality applications. This wordset begins with *. and *_ the fractional multiply routines.

There are two of versions of fractal multiplications:

       *.   where $80000 is 1.0
       *_  where $80000 is 0.5

       *_ ( N1 UF -- N2)
           UF is unsigned fractal from 0 to 0.99
           where $80000 is 0.5

       *.  ( N1 UF -- N2)
           UF is unsigned fractal from 0 to 1.99
           where $80000 is 1.0
These N1 is signed(unsigned) integer(fractional) - or any mix of combinations! But UF is unsigned fractal. So if I'll need signed fractal mult. in future then I'll program some hi-lewel word with sign memory etc. Of course this word will use *. and main calculations will be done inside last one.

There are some small useful hi-level application of Cordic. Of course the main applications is vector and orientation matrix rotation, but this stuff demands more time for explanation. Note

HEX
80000 CONSTANT  PI
PI    CONSTANT -PI
CC    CONSTANT 9B760 \ Cordic constant, fractal image of 0.607253

\ It's OK because angle metric nature is rotational!
\ You may to increment or decrement some angle infinitely
\ without some angle value correction...

CODE CORDIC ( X Y A -- X' Y') ... END-CODE
: SIN*    ( N A -- N*SIN{A} ) 0 SWAP CORDIC NIP  CC *. ;
: COS*    ( N A -- N*COS{A} ) 0 SWAP CORDIC DROP CC *. ;
: TG*     ( N A -- N*TAN{A} ) 0 SWAP CORDIC SWAP */ ;
: CTG*    ( N A -- N*CTAN{A}) 0 SWAP CORDIC */ ;
There is new primitive : POLAR ( x y -- a h') where x,y is cortesian coordinates; a is NEGATIATED angle of polar coordinate and h' is vector length of polar coordinates. To get real value of h' multiply it by Cordic constant. There are several of useful hi-level utils:
: CORT>POLAR ( x y -- h a) POLAR CC *. SWAP NEGATE ;
: HIP ( x y -- h) POLAR NIP *. ;
: ARCCTG ( x y -- a) POLAR DROP NEGATE ;
: ARCTG  ( x y -- a) POLAR DROP PI/2 + ;
In ARCCTG/ARCTG case You may to see at x and y as rational image of some real. For example, to calculate arctg(1.23456) You should type
       123456 100000 ARCTG.
1. ANGLE MEASUREMENT UNIT.

There are several of angle presentation way: integer degrees, fixed and rational radian representation. I think most suitable angle representation is full range integer ( in our case PI = $80000 ). The advances of such way are : maximum of accuracy and calculation speed. It because the main of the angle feature is its repeating a = a +/- n*2PI But integers has same nature too! N + $100000 (for the 0 - $FFFFF range) is N ! And NO difference between sign and unsign! You may to see on the $FFFFF as both -1 and MaxUnsign value. And, of course, if You need to have radians or degree You may to get it by scaling easy and without accuracy loosing. So our PI = -PI = $80000.

2. THE CORDIC CONSTANT CORRECTION

You know all our trigs needs some corrections for its results. This corrections is multiplication by Cordic_Constant. The Cordic_Constant value for 20-bits is 0.607252935009249. There are two way to make this corrections more accurate and more fast.

The most accurate way is

       : CORRECT $242C $3B91 */ ;
The fastest way is
       : CORRECT $9B74 *_ ;
3. THE NAMES.

Constants:

       $80000 CONSTANT PI
       $40000 CONSTANT PI/2
Coordinate -rotating function (one of most important!):
: ROTATE ( x1 y1 a -- x2 y2 ) CORDIC CORRECT SWAP CORRECT SWAP ;
Words to convert between coordinate systems :
: CORT>POLAR ( x y -- h a ) POLAR CORRECT SWAP NEGATE ;
: POLAR>CORT ( h a -- x y ) 0 SWAP CORDIC CORRECT SWAP CORRECT SWAP ;
Multiplications by trigs functions:
SIN*  ( n a -- n*sin a)
COT*  ( n a -- n*cos a)
TG*   ( n a -- n*tan a)
CTG*  ( n a -- n*ctan a)
Arc -functions, with x,y -pair representation (You can to see on it as on rational representation) of input arguments
: ATAN2  ( y x -- a ) POLAR DROP PI/2 + ;
: ACTAN2 ( x y -- a ) POLAR DROP NEGATE ;
: HIP    ( x y -- h=sqrt{x*x+y*y}) POLAR NIP CORRECT ;
4. OVERFLOWING

As you may know the way the cordic algorithms are based on vector rotation. cordic rotates the vector until the angle becomes zero and polar rotates one until Y-coordinates becomes zero. But the vector length grows while this rotation! It becomes 1/Cordic_Constant greater the initial length! There are to cordic correction usage way: before rotation and after one. In my examples I used last one. It's more accurate, and more overflowable ... Anyway we must to remember about the effect.

3D to 2D transformations can be performed using the cordic function. A set of words for 3D project will be built on these math functions. These words will provide a basis or 3D CAD and virtual reality applications.

(Dave Lowry's 3D demo does not use these math functions.)


End of Chapter 8

Return to Table of Contents

Previous Chapter

Next Chapter