Threaded Mode | Linear Mode

Garth Wilson · 03-12-2022, 09:46 PM

Here's some stuff I did many years ago in 16-bit Forth. Hopefully I correctly converted all the special characters from the DOS/ANSI character set.

A 16-bit circle works out nicely, scaling 360° to be exactly 65,536 counts, with 1 lsb being .00549316°, or .32959 arc-minutes, or 19.7754 arc-seconds. The MOD function doesn't even come up, as there's no need for it. 359° plus 2° is still 1°, and 45° minus 90° is still -45° or 315°. You can consider it signed or unsigned, either way. 90° is represented by $4000, 180° is represented by $8000, and 270° or -90° is represented by $C000. A radian is represented by $28BE. When you take the sin or cos, the output range will be ±1, which, for best resolution, we can scale to as high as ±$7FFF or $8000 depending on how you want to work the trapping (since you can't get +$8000 in a 16-bit signed number).

There's a Wikipedia article on binary angular measure (BAM) and "brads" (binary radians) here, and Jack Crenshaw addresses it in chapter 5 of his fine book "Math Toolkit for Real-Time Programming."

Code:

: U*/_RND       ( U1 U2 U3 -- U1*U2/U3 )    \ Fractional part of answer is rounded in, not truncated.

   >R  UM*

   R@  -1 SHIFT  0  D+

   R>  UM/MOD  NIP      ;

\ cos45 below uses the 1st 4 terms of the infinite series, but I changed one of the coefficients to improve the accuracy.

\ The 8002 started as 7FFF.  The way it is below, the max error I've seen is 1 lsb high.  Usually it's right on though.

7FFF  CONSTANT  7FFF

: cos45   ( rad•28BE -- cos•$7FFF )     \ Input (X) is 0 (0°) up to $2000 (45°).

   7FFF 28BE U*/_RND     \ First, scale the input so 1 radian is 7FFF instead of 28BE.  Now max (45°) is $6488.

   DUP  7FFF U*/_RND     \ ^ 7FFF•X²                            TOS range is     0 for X=0 to $4EF6 for X=45°.

   DUP 168 + 2D0 /       \ ^    "    same/6!         (6!=$2D0.) TOS range is     0 for X=0 to   $1C for X=45°.

         555 SWAP -      \ ^    "    7FFF•(1/4!-X²/6!)          TOS range is  $555 for X=0 to  $539 for X=45°.

   OVER 7FFF U*/_RND     \ ^    "    7FFF•X²(1/4!-X²/6!)        TOS range is     0 for X=0 to  $339 for X=45°.

        4000 SWAP -      \ ^    "    7FFF•(1/2!-X²(1/4!-X²/6!)  TOS range is $4000 for X=0 to $3CC7 for X=45°.

        8002 U*/_RND     \ ^ 7FFF•X²(1/2!-X²(1/4!-X²/6!)        TOS range is     0 for X=0 to $257D for X=45°.

        7FFF SWAP -   ;  \ ^ 7FFF•(1-X²(1/2!-X²(1/4!-X²/6!)     TOS range is $7FFF for X=0 to $5A82 for X=45°.

\ sin45 below uses the 1st 3 terms of the infinite series, but I changed a couple of the coefficients by 1 to improve the

\ accuracy.  The 1556 started as 1555, and the final 8000 was 7FFF.  The way it is below, the max error I've seen is 1 lsb

\ high, as long as the input does not exceed 45°.  Even a little above 45°, the accuracy is lost very quickly!

: sin45   ( rad•28BE -- sin•$7FFF )     \ Input (X) is 0 (0°) up to $2000 (45°).

   7FFF 28BE U*/_RND     \ First, scale the input so 1 radian is 7FFF instead of 28BE.  Now max (45°) is $6488.

   DUP                   \ ^ 7FFF•X  7FFF•X                     Keep a copy of unsquared scaled input for use later.

   DUP  7FFF U*/_RND     \ ^    "    7FFF•X²                           TOS range is     0 for X+0 to $4EF6 for X=45°.

   DUP 3C + 78 /         \ ^    "    7FFF•X²  7FFF•X²/5!    (5!=$78.)  TOS range is     0 for X=0 to   $A8 for X=45°.

        1556 SWAP -      \ ^    "    7FFF•X²  7FFF•(1/3!-X²/5!)        TOS range is $1556 for X=0 to $14AE for X=45°.

        7FFF U*/_RND     \ ^    "    7FFF•X²(1/3!-X²/5!)               TOS range is     0 for X=0 to  $CC1 for X=45°.

        7FFF SWAP -      \ ^    "    7FFF•(1-X²(1/3!-X²/5!))           TOS range is $7FFF for X=0 to $733E for X=45°.

        8000 U*/_RND  ;  \ ^ 7FFF•X(1-X²(1/3!-X²/5!)                   TOS range is     0 for X=0 to $5A82 for X=45°.

: sin90   ( rad•28BE -- sin•$7FFF )     \ Input is 0 (0°) up to $4000 (90°).

   DUP  2000  >                         \ Is it over 45°?

   IF   4000 SWAP - cos45  EXIT         \ If so, just take the cosine of the opposite angle.

   THEN sin45   ;                       \ Otherwise just go straight to the sin word above.

: cos90   ( rad•28BE -- cos•$7FFF )     \ Input is 0 (0°) up to $4000 (90°).  I don't remember why it's here, since it's not used

   DUP  2000  >                         \ Is it over 45°?                                                            in COS below.

   IF   4000 SWAP - sin45  EXIT         \ If so, just take the cosine of the opposite angle.

   THEN cos45   ;                       \ Otherwise just go straight to the cos word above.

: SIN      ( rad•28BE -- sin•$7FFF )

   DUP  3FFF AND                        \ ^ RADs HW_FAR_N2_QUADRANT

   OVER 4000 AND                        \ ^  "          "       BACKWARDS?

   IF   4000 SWAP - THEN  sin90         \ ^  "        90°_SIN

   SWAP 0<                              \ ^ 90°_SIN  RADs

   IF NEGATE THEN       ;               \ ^ SIN

: COS    4000 + SIN              ;    ( RAD -- COS )       \ Good for any quadrant.

: TAN       DUP SIN   SWAP COS   ;    ( RAD -- SIN COS )   \ For any quadrant