Path: csiph.com!eternal-september.org!feeder.eternal-september.org!news.eternal-september.org!.POSTED!not-for-mail
From: Paul Rubin <no.email@nospam.invalid>
Newsgroups: comp.lang.forth
Subject: Re: embedded forth
Date: Mon, 30 Jan 2017 14:09:07 -0800
Organization: A noiseless patient Spider
Lines: 43
Message-ID: <87d1f4yzcs.fsf@nightsong.com>
References: <49b3566a-3842-4f9f-a47f-a66590625fef@googlegroups.com> <92d07db1-18b7-410e-a77c-d7ab4f6606d8@googlegroups.com> <87lgtu20wb.fsf@nightsong.com> <6415800d-9ed0-4eef-afcf-6291de9a1ef3@googlegroups.com> <o6lcft$em9$1@cherry.spenarnc.xs4all.nl> <e90a34e0-376e-4910-bd6a-96d8f16bf0bf@googlegroups.com> <o6mla6$f0l$1@dont-email.me> <e8b3a38a-d6b7-41b3-a760-fe07881bdf17@googlegroups.com> <o6nl76$o8v$1@dont-email.me> <d126a340-9abd-4145-99f6-2d01915ebe64@googlegroups.com>
Mime-Version: 1.0
Content-Type: text/plain
Injection-Info: mx02.eternal-september.org; posting-host="97a7bf2542d909943fb48cc1086f163d"; logging-data="6239"; mail-complaints-to="abuse@eternal-september.org";	posting-account="U2FsdGVkX18yR37HD2yrLJe+mAOflqp/"
User-Agent: Gnus/5.13 (Gnus v5.13) Emacs/24.3 (gnu/linux)
Cancel-Lock: sha1:Uosk2+fmdlRy3WfOxMOIWoPFfP0= sha1:uLA9DZXobq25VMb+h+8EPCIkFpg=
Xref: csiph.com comp.lang.forth:54556

Ilya Tarasov <ilya74.tarasov@gmail.com> writes:
> From the first look, we may list degrees, minutes and seconds, but is
> it optimal?

It's convenient, so the question then is whether it's close enough.

Optimal involves non-uniform placement of the interpolation points,
e.g. at the so-called Chebychev nodes according to Wikipedia.  I'd
expect that with non-optimal placement, you can still get the same
accuracy by using more points.  That seems ok to me if you have enough
memory for the bigger tables.

Using naive deg/min/sec interpolation (see below) I get a max error of
about 1.56e-4 which is well above 1 ULP in single precision.  I think
that can be improved by just moving the interpolation point slightly.
There's too stuff going on here at the moment for me to figure that out
right now, but this is a start.

================================================================

: dms.1 ( f:x -- f:int f:frac ) fdup f>s s>f ftuck f- ;
: f>dms ( F:degrees -- F:degrees F:minutes F:seconds )
    dms.1 60e f* dms.1 60e f* dms.1 f+ ;

: dms>f ( F:d F:m F:s -- F:dd )    60e f/ f+ 60e f/ f+ ;
: deg>rad ( F:d -- F:rad ) pi f* 180e f/ ;
: rad>deg ( F:rad -- F:deg ) 180e f* pi f/ ;
: min>rad ( F:min -- F:rad ) 60e f/ deg>rad ;
: sec>rad ( F:min -- F:rad ) 60e f/ min>rad ;

\ approximate sin x using d/m/s interpolation
: xsin ( F:deg -- F:approximate sine )
    f>dms  sec>rad  { F: s } min>rad { F: m } deg>rad { F: d }
    d fcos m fsin f*  
    d fsin m fcos f* f+
    fdup fcos s f*
    f+ ;

: fclearstack  fdepth 0 ?do fdrop loop ;
: cs clearstack fclearstack ;

: err ( deg -- e ) fdup xsin fswap deg>rad fsin f/ 1e f- fabs ;
: maxerr ( -- e ) cs 0e 89999 1 do i s>f 1e-3 f*  err  fmax loop ;