Approximation of Sine function in Integer Forth

[Martin Hepperle]

Recently I wanted to demonstrate a MS-DOS driver for plotting HPGL on a VGA screen.

To show how this driver could be used from any programming language I selected Forth as one example. The classical Forth implementations, like figForth did not support floating point calculations. Thus, everything has to be done with Integers, in my case in the range +/-32767.

HPGL works nicely with integers as it defines its device units as 0.25 mm per unit.

In order to generate some interesting plots with circles and spirals, I wanted to have an approximation of the trigonometric sine and cosine functions.
As I had to stick with integers, I decided to scale the sine by a factor of 1000.
So the approximation produces values in [-1000...+1000].
The input should be the angle in degrees, so that the angular resolution would be 1 degree.

The restriction to signed 16-bit integers generally limits polynomial terms to powers of two and all integer divisions truncate. The often used power series around zero needs at least powers of 3 to produce reasonable accuracy over the range [0 to 90] degrees.

First I tried a simple parabolic curve covering the range [0 to 180] degrees, which worked, but circles looked more like well rounded rectangles.
In a second attempt I generated series expansions at the center point of four intervals over the range [0 to 90] degrees. While this requires tests for the interval, it is still reasonable fast.
Some manual tweaking of the equations to the restriction of integer arithmetic produced a quite acceptable approximation.

The attached PDF document shows my current set of equations and their deviation from the sine function as well as their implementation in Forth and RPL.

Note: I am not looking just for a circle generator (there are a few integer algorithms for digital raster devices), but for the trigonometric functions for general applications.

Maybe there is another more elegant or more accurate approach?

Martin

[edit: programs adapted to match equations below graph]
[edit: added Garth's solution to the graph and page 2]

[Thomas Klemm]

(03-12-2022 04:25 PM)Martin Hepperle Wrote:  Maybe there is another more elegant or more accurate approach?

You could use Bhaskara's Sine and Cosine Approximations.
They are accurate to about 1‰.

Functions in Python

Code:

def sin(x):
    u = (180 - x) * x
    return 4000 * u // (40500 - u)

def cos(x):
    u = x**2
    return 4000 * (8100 - u) // (32400 + u)

Plot the Difference

Code:

import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(0, 90, 91)

Sine

Code:

plt.plot(x, sin(x) -  1000 * np.sin(np.radians(x)), 'b')

Cosine

Code:

plt.plot(x, cos(x) -  1000 * np.cos(np.radians(x)), 'r')

Cheers
Thomas

[Garth Wilson]

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

[Martin Hepperle]

Thank you for your comments!

For starting, I have implemented Garth's Forth code and added the result to the deviation plot.

The overall deviation is much lower than my piecewise approximation but due to exploiting the full range for scaling, Garth's result is also much smoother.

[edit: the attachment in my initial post was updated accordingly]

[the only change I made to Garth's code was to replace -1 SHIFT by 1 RSHIFT in the word U*/_RND as my Forth only has RSHIFT and LSHIFT]

Martin

[Thomas Klemm]

(03-13-2022 10:00 AM)Martin Hepperle Wrote:  the attachment in my initial post was updated accordingly

Not sure if I interpret the new red line correctly but it looks as if you used the sin45 function.
At least a difference of more than 5 for 90° suggests this.

(03-12-2022 09:46 PM)Garth Wilson Wrote:  Even a little above 45°, the accuracy is lost very quickly!

I recommend to use the sin90 function instead.

[Garth Wilson]

(03-13-2022 12:56 PM)Thomas Klemm Wrote:  

(03-12-2022 09:46 PM)Garth Wilson Wrote:  Even a little above 45°, the accuracy is lost very quickly!

I recommend to use the sin90 function instead.

Note that my sin90 and cos90 use the cos45 and sin45 (respectively) to get the desired accuracy above 45 degrees. Then the SIN and COS functions are usually accurate to all 16 bits, and the few times they're off, it's only by one lsb, and even then, it's always in the same direction, ie, high.

[Paul Dale]

Considered a table lookup? Range reduce to forty-five degrees, look up the table, possibly do a linear interpolation between table values.

[Martin Hepperle]

For my purpose, Garth's solution is very usable and "good enough" (famous last words ;-)).
It also fits nicely on a single Forth screen, has only a small memory footprint and takes about the same time as my initial, poor approximation.

I don't need "absolute" accuracy - in this case I would probably consider something like Cordic and a floating point number representation.

In the first update of my PDF document I had Garth's solution overlaid with wrong scaling on the horizontal axis, mapping 0...360 deg to my scale of 0...90 deg. I noticed my mistake and had corrected the document half an hour later, so maybe Thomas downloaded the incorrect, intermediate version?

Martin

[Thomas Klemm]

(03-14-2022 01:16 PM)Martin Hepperle Wrote:  (…) maybe Thomas downloaded the incorrect, intermediate version?

It looks much better now.

[pavel nemec cz]

(03-12-2022 04:25 PM)Martin Hepperle Wrote:  Maybe there is another more elegant or more accurate approach?

Martin

There is an example of transcendential functions calculation in the Lesson 12 of Dr. C.H. Ting collection with the refference to original article in Forth Dimension Volume 04 Issue 1. Many more to learn from Dr. Ting about Forth if you are interested including nice Win32Forth.

KR
Pavel

[Thomas Klemm]

(03-14-2022 05:58 PM)pavel nemec cz Wrote:  There is an example of transcendental functions calculation …

This was an interesting read. Thanks for sharing.

(03-14-2022 12:15 AM)Paul Dale Wrote:  Considered a table lookup?

In the linked document Forth Dimensions we can also find an article on page 6:
Fixed-Point Trig by Table Lookup

[Martin Hepperle]

Thank you for the reference to Forth Dimensions - these old issues contain a lot of useful information and tips and tricks. I had forgotten to look through them again.

The issue V4N1 is indeed interesting as it also shows a pure Forth implementation of Floating Point numbers, which demonstrates how large such an implementatio0n becomes when you use no assembler or coprocessor.


BTW: This web page has an index of articles in V1 to V6 of Forth Dimensions and direct links to the issues in PDF.

Martin

[John Keith]

(03-15-2022 01:15 PM)Martin Hepperle Wrote:  BTW: This web page has an index of articles in V1 to V6 of Forth Dimensions and direct links to the issues in PDF.

Martin

Thanks, great reference!

[deetee]

Some time ago I had to draw a simple ascii-circle for some kind of "clock project" for an Arduino. The primary goal (on Arduinos) is to write codes with small footprint and sufficient precision. I used some kind of lookup table with 16 bytes only to calculate rough values for sine (and cosine).

Maybe this code snippet can give you another idea?
You can save the division (/ 200) by expanding (and/or recalculating) the byte values to int values.

Code:


  const byte sine[] = {0, 21, 42, 62, 81, 100, 118, 134, 149, 162, 173, 183, 190, 196, 199, 200};

  static void printcircle(byte X, byte Y, byte r) { // Print circle at (X|Y) with radius r
  for (byte i = 0; i < 16; i++) {
    byte x = r * sine[15 - i] / 200, y = r * sine[i] / 200;
    byte a = X + x, b = X - x, c = Y + y, d = Y - y;
    printpixel(a, c); printpixel(a, d); printpixel(b, d); printpixel(b, c);
  }

Regards
deetee