FMOD in Forth
Krishna Myneni
After re-reading DPANS94, Floating Point Extensions,
I find there is no equivalent of the C library function
fmod() which returns the remainder (fractional part) of
the result of dividing two fp numbers. Needing this, I
had to resort to the clumsy and inefficient source level
definition
: FMOD ( f1 f2 -- f | remainder of f1/f2 )
FOVER FOVER F/ F>D D>F F* F- ;
FMOD is particularly useful for calculations involving
angles, which are quite common. Unless I've missed something,
I recommend implementation of FMOD in Forth systems which
support floating point. The word and its usage is analagous
to the Forth word MOD.
On a related note, last year after discussions with Guido
Draheim, David Williams and myself, we chose to implement
the word FTRUNC for providing floating point truncation
towards zero. This word complements the existing standard
words FROUND and FLOOR (note that FLOOR truncates towards
minus infinity). The source def of FTRUNC is
: FTRUNC F>D D>F ;
but FTRUNC can be more efficiently implemented with the
floating point unit instructions. I recommend the
implementation of this word --- it exists in recent
versions of PFE and kForth.
Krishna Myneni
andr...@littlepinkcloud.invalid
> After re-reading DPANS94, Floating Point Extensions, I find there is
> no equivalent of the C library function fmod() which returns the
> remainder (fractional part) of the result of dividing two fp
> numbers. Needing this, I had to resort to the clumsy and inefficient
> source level definition
> : FMOD ( f1 f2 -- f | remainder of f1/f2 )
> FOVER FOVER F/ F>D D>F F* F- ;
> FMOD is particularly useful for calculations involving angles, which
> are quite common. Unless I've missed something, I recommend
> implementation of FMOD in Forth systems which support floating
> point.
It's very, very hard to do this accurately. Range reduction for trig
functions requires an approximation to pi that is much longer than the
argument if it is to produce a result that is accurate. There's an
example of how to do this in fdlibm.
Andrew.
Krishna Myneni
> Krishna Myneni <krishn...@compuserve.com> wrote:
>
>> ...
>>
>>: FMOD ( f1 f2 -- f | remainder of f1/f2 )
>> FOVER FOVER F/ F>D D>F F* F- ;
>
>
>>FMOD is particularly useful for calculations involving angles, which
>>are quite common. Unless I've missed something, I recommend
>>implementation of FMOD in Forth systems which support floating
>>point.
>
>
> It's very, very hard to do this accurately. Range reduction for trig
> functions requires an approximation to pi that is much longer than the
> argument if it is to produce a result that is accurate. There's an
> example of how to do this in fdlibm.
>
> Andrew.
Thanks for the warning. I'm not advocating using the above
Forth source definition for implementation. It was a quick
and dirty patch for my use. Your comment implies that
the C math library function fmod() can suffer accuracy
problem when working with angular remainders and these
remainders are used in trig functions. Offhand I
don't see why there is any particular accuracy problem
involved in mod pi or mod 2pi, but it is something to
look at closely when doing an implementation. Thanks
for the pointer to the example.
The ANS fp standard is good but suffers from a few deficiencies
such as lack of words to perform FMOD and FTRUNC functions.
Krishna
andr...@littlepinkcloud.invalid
> andr...@littlepinkcloud.invalid wrote:
>> Krishna Myneni <krishn...@compuserve.com> wrote:
>>
>>> ...
> >>
>>>: FMOD ( f1 f2 -- f | remainder of f1/f2 )
>>> FOVER FOVER F/ F>D D>F F* F- ;
>>
>>
>>>FMOD is particularly useful for calculations involving angles, which
>>>are quite common. Unless I've missed something, I recommend
>>>implementation of FMOD in Forth systems which support floating
>>>point.
>>
>> It's very, very hard to do this accurately. Range reduction for trig
>> functions requires an approximation to pi that is much longer than the
>> argument if it is to produce a result that is accurate. There's an
>> example of how to do this in fdlibm.
> Thanks for the warning. I'm not advocating using the above Forth
> source definition for implementation. It was a quick and dirty patch
> for my use. Your comment implies that the C math library function
> fmod() can suffer accuracy problem when working with angular
> remainders and these remainders are used in trig functions.
Very much so. Any use of fmod is rather suspect!
> Offhand I don't see why there is any particular accuracy problem
> involved in mod pi or mod 2pi, but it is something to look at
> closely when doing an implementation.
Well, if you use fmod for angular reduction you lose one bit of
significance in the result for every bit of magnitude in the argument.
Whether this is so bad depends on how the transcendental functions are
specified, but these days it's common to require 1 ulp, and you don't
get that by using fmod.
Andrew.
Bernd Paysan
> Well, if you use fmod for angular reduction you lose one bit of
> significance in the result for every bit of magnitude in the argument.
> Whether this is so bad depends on how the transcendental functions are
> specified, but these days it's common to require 1 ulp, and you don't
> get that by using fmod.
Important hint: Floating point numbers are not numbers, but intervals (+-1/2
ulp). You do not lose any bit of significance, because you didn't have it
before. The only thing for angular reduction is that pi is known better
than to +-1/2 ulp, and therefore, you could produce a better result than
with fmod, due to a more accurate pi, but it's only 1/2 ulp total that you
can gain (the inaccuracy of pi). If you have real rational angles to
reduce, it's probably better to keep the angle in degrees or a comparable
form where you can reduce it with fmod without losing 1/2 ulp (possible
because the reduction factor is an integer), and then multiply the reduced
angle by pi.
--
Bernd Paysan
"If you want it done right, you have to do it yourself"
http://www.jwdt.com/~paysan/
Jerry Avins
> ... these days it's common to require 1 ulp, ...
I don't know that unit. I've seen it in comic books to imitate the sound
of a gulp, a sort of abbreviated "omygosh". I like the sound. If I knew
what it meant, I might use it often. :-)
Jerry
--
... they proceeded on the sound principle that the magnitude of a lie
always contains a certain factor of credibility, ... and that therefor
... they more easily fall victim to a big lie than to a little one ...
A. H.
———————————————————————————————————————————————————————————————————————
Bernd Paysan
> andr...@littlepinkcloud.invalid wrote:
>
>> ... these days it's common to require 1 ulp, ...
>
> I don't know that unit. I've seen it in comic books to imitate the sound
> of a gulp, a sort of abbreviated "omygosh". I like the sound. If I knew
> what it meant, I might use it often. :-)
Try www.acronymfinder.com: ULP "Unit in the Last Place". There are other
ULPs, too, but this one is the relevant here.
Jerry Avins
> Jerry Avins wrote:
>
>
>>andr...@littlepinkcloud.invalid wrote:
>>
>>
>>> ... these days it's common to require 1 ulp, ...
>>
>>I don't know that unit. I've seen it in comic books to imitate the sound
>>of a gulp, a sort of abbreviated "omygosh". I like the sound. If I knew
>>what it meant, I might use it often. :-)
>
>
> Try www.acronymfinder.com: ULP "Unit in the Last Place". There are other
> ULPs, too, but this one is the relevant here.
>
Thank you. It's a lovely word. Here's one that used to be mine, but I
gave it to my niece: potterphernalia, the odds and ends in a
ceramicist's toolkit.
Coos Haak
> Jerry Avins wrote:
>
>> andr...@littlepinkcloud.invalid wrote:
>>
>>> ... these days it's common to require 1 ulp, ...
>>
>> I don't know that unit. I've seen it in comic books to imitate the sound
>> of a gulp, a sort of abbreviated "omygosh". I like the sound. If I knew
>> what it meant, I might use it often. :-)
>
> Try www.acronymfinder.com: ULP "Unit in the Last Place". There are other
> ULPs, too, but this one is the relevant here.
<OT>
Yes, it's the first answer there.
But I'm still puzzled by YYRRW that Guy Macon used in another thread ;-(
</OT>
--
Coos
Brad Eckert
>
> : FMOD ( f1 f2 -- f | remainder of f1/f2 )
> FOVER FOVER F/ F>D D>F F* F- ;
>
> FMOD is particularly useful for calculations involving
> angles, which are quite common. Unless I've missed something,
> I recommend implementation of FMOD in Forth systems which
> support floating point. The word and its usage is analagous
> to the Forth word MOD.
Integers are the natural way to represent angles. One cell spans 360
degrees, so when it overflows it stays on the unit circle.
You could scale the angle by a constant to convert a floating point
angle to this format. 2^#bits/2pi.
1e0 FATAN 8e0 F* FCONSTANT 2pi
0 1 D>F 2pi F/ FCONSTANT ANGLEINT \ 2^32/2pi
2pi 0 1 D>F F/ FCONSTANT ANGLEFLT \ 2pi/2^32
: FMOD F/ ANGLEINT F* \ angle to cell scale
F>D D>S 0 D>F \ mod the angle
ANGLEFLT F* ; \ restore to radians
The question is, does the x86 FPU have FMOD? If it's not in the Forth
standard I guess not.
Brad
andr...@littlepinkcloud.invalid
> andr...@littlepinkcloud.invalid wrote:
>> Well, if you use fmod for angular reduction you lose one bit of
>> significance in the result for every bit of magnitude in the argument.
>> Whether this is so bad depends on how the transcendental functions are
>> specified, but these days it's common to require 1 ulp, and you don't
>> get that by using fmod.
> Important hint: Floating point numbers are not numbers, but
> intervals (+-1/2 ulp).
That's one way to interpret floating-point numbers, but it's equally
reasonable to regard them as a subset of the rationals.
> You do not lose any bit of significance, because you didn't have it
> before.
Well, that depends on whether you regard floating-point numbers as
intervals or as a subset of the rationals.
> The only thing for angular reduction is that pi is known better than
> to +-1/2 ulp, and therefore, you could produce a better result than
> with fmod, due to a more accurate pi, but it's only 1/2 ulp total
> that you can gain (the inaccuracy of pi).
Not so at all. Even if you regard floating-point numbers as
intervals, it's ideal to return the closest floating-point number to
the sine of the midpoint of the interval.
A good floating-point package will thse days return transcendentals
within 1 ulp, regardless of the magnitude of the arguments.
Andrew.
Julian V. Noble
>
[ deleted ]
>
> The question is, does the x86 FPU have FMOD? If it's not in the Forth
> standard I guess not.
>
> Brad
In fact the Intel and clone fpu's have FPREM and FPREM1 that subtract an
integer multiple of ST(1) from ST(0) until a fractional remainder is found.
That is,
ST(0) <- ST(0) - q * ST(1)
where q is the (integer part of the) quotient of ST(0) / ST(1) .
So the answer is, FMOD can be added to the repertoire with a little
CODE.
--
Julian V. Noble
Professor Emeritus of Physics
j...@lessspamformother.virginia.edu
^^^^^^^^^^^^^^^^^^
http://galileo.phys.virginia.edu/~jvn/
"For there was never yet philosopher that could endure the toothache
patiently." -- Wm. Shakespeare, Much Ado about Nothing. Act v. Sc. 1.
andr...@littlepinkcloud.invalid
> Brad Eckert wrote:
>>
> [ deleted ]
>>
>> The question is, does the x86 FPU have FMOD? If it's not in the Forth
>> standard I guess not.
> In fact the Intel and clone fpu's have FPREM and FPREM1 that subtract an
> integer multiple of ST(1) from ST(0) until a fractional remainder is found.
> That is,
> ST(0) <- ST(0) - q * ST(1)
> where q is the (integer part of the) quotient of ST(0) / ST(1) .
> So the answer is, FMOD can be added to the repertoire with a little
> CODE.
FPREM isn't guaranteed to finish the job: you need to loop until the
C2 flag is cleared. This is needed in order not to have excessive
interrupt latency.
Andrew.
Bill McCarthy
> Brad Eckert wrote:
>>
> [ deleted ]
>>
>> The question is, does the x86 FPU have FMOD? If it's not in the Forth
>> standard I guess not.
>>
>> Brad
> In fact the Intel and clone fpu's have FPREM and FPREM1 that subtract an
> integer multiple of ST(1) from ST(0) until a
> fractional remainder is found.
> That is,
> ST(0) <- ST(0) - q * ST(1)
> where q is the (integer part of the) quotient of ST(0) / ST(1) .
> So the answer is, FMOD can be added to the repertoire with a little
> CODE.
I've written such code ages ago for SwiftForth. For
what it's worth, here's the code:
\ Add a symmetric FMOD
xCODE FMOD ( -- ) ( x y -- mod ) 2 >fs
FXCH
BEGIN
FPREM
FSTSWAX
4 # AH TEST
0=
UNTIL
ST(1) FSTP
f> FNEXT
\ Proper modulo
CODE F_MOD ( -- ) ( x y -- mod ) 2 >fs
FTST
FSTSWAX
EAX 8 # SHL
FXCH
FTST
FSTSWAX
EAX 8 # SHR
$0101 # EAX AND \ isolate signs in AH & AL
EAX 16 # SHL \ move them out of the way
BEGIN
FPREM
FSTSWAX
4 # AH TEST
0=
UNTIL
FTST
FSTSWAX
$40 # AH TEST
0= IF
EAX 16 # SHR \ restore the signs in AH & AL
AH AL XOR
0<> IF
ST(1) ST(0) FADD
THEN
ELSE
FABS \ Get rid of negative zero
THEN
ST(1) FSTP
f> FNEXT
--
Best regards,
Bill
Krishna Myneni
The statement that integers are the natural way to
implement angles depends on the context --- it is
correct to say that integers are the more common way
to represent angles. Physicists often think of angles
in terms of radians. The only reason we associate 360
degrees with the circle is that we have close to 360
days for the orbital period of the earth around the
sun (or the sun around the earth as seen by ancient
people). But pi is a universal ratio for all circles.
Does your implementation of FMOD retain more precision than
my original def?
Krishna
Krishna Myneni
Excellent! FMOD will appear in the next development release
of kForth. Which version corresponds to the C math library
function fmod()? There is some value in retaining consistency
across languages.
Krishna
Krishna Myneni
> ...
> Well, if you use fmod for angular reduction you lose one bit of
> significance in the result for every bit of magnitude in the argument.
> Whether this is so bad depends on how the transcendental functions are
> specified, but these days it's common to require 1 ulp, and you don't
> get that by using fmod.
>
> Andrew.
Ok. There are two separate issues here. One is the implementation
of FMOD in Forth. The second is whether FMOD should be used to
modulo angles. The latter is an accuracy issue which can be dealt
with separately from the issue of implementing FMOD in Forth.
Krishna
Ed Beroset
By the context, I guessed "Yeah, Yeah, Right, Right, Whatever."
Ed
Ed Beroset
had intended to write was that I would very highly recommend "What every
computer scientist should know about floating-point" by David Goldberg,
published in ACM Computing Surveys in 1991. (Also available here
http://cch.loria.fr/documentation/IEEE754/ACM/goldberg.pdf and other
places). It's clear, simple, and extremely useful.
Ed
Marcel Hendrix
[..]
> Ok. There are two separate issues here. One is the implementation
> of FMOD in Forth. The second is whether FMOD should be used to
> modulo angles. The latter is an accuracy issue which can be dealt
> with separately from the issue of implementing FMOD in Forth.
FWIW, iForth has the following high-level names and definitions in
miscutil.frt:
[UNDEFINED] FTRUNC [IF] :INLINE FTRUNC ( F: a -- b ) F>S S>F ; [THEN]
[UNDEFINED] FCIRCULAR [IF] :INLINE FCIRCULAR ( F: a b -- a_mod_b ) F2DUP F/ FLOOR F* F- ; [THEN]
CODE FPREM ( F: r1 r2 -- rem[r1/r2] )
fpop,
@fswap,
BEGIN,
$D9 C, $F8 C, ( fprem, )
$DF9B 16B, $E0 C, ( @fstsw, )
$9E C, ( sahf, )
po,
UNTIL,
@fnip,
$FF C, $E3 C, ( ebx jmp, )
END-CODE
: FP/REM ( F: r1 r2 -- rem[r1/r2] ) ( -- div[r1/r2] )
F2DUP F/ F>S FPREM ;
-marcel
Krishna Myneni
Thanks! You may consider making FMOD an alias of FPREM. While
the name FPREM mirrors the fpu instruction, FMOD strongly
suggests its meaning to Forth programmers who are used to MOD,
as well as to C programmers. BTW, it's refreshing to see a
person who's not afraid to mix machine code with high level
code :)
Krishna
crypto_s...@sympatico.ca
The question really is; are floating point appropriate for the kind of
accuracy you need. You seem to be very specific about it being
accurate. If this is the case, you might want a solution involving
symbolic processing in conjunction with keeping for example 2 32 bit
integers as unresolved divisions for each floating point value instead
of an actual floating point values which is bound to be very imprecise
to a certain extent. The math program MAPLE works in this manner as
far as the symbolic processing goes and the results are of course dead
on the mark. Symbolically speaking that is :).
Of course this requires developping a whole new library of functions
to process unresolved fractions but in the end, it will definately
give you what you want.
Stonelock
Coos Haak
>
> The statement that integers are the natural way to
> implement angles depends on the context --- it is
> correct to say that integers are the more common way
> to represent angles. Physicists often think of angles
> in terms of radians. The only reason we associate 360
> degrees with the circle is that we have close to 360
> days for the orbital period of the earth around the
> sun (or the sun around the earth as seen by ancient
> people). But pi is a universal ratio for all circles.
Why use 360 degrees, surveyors use 400 degrees in a circle.
360 is 6 times 60, the base of Babylonian numbering.
Pi is much younger, a Greek letter and I doubt if they
used it in relation to circles.
In a 16 bit implementation, a circle has 65536 degrees ;-)
--
Coos
Jorge Acereda
> In a 16 bit implementation, a circle has 65536 degrees ;-)
I guess it depends on the application... I normally use the full range
for angles, has lots of advantages, but you can't express "number of
spins" with that.
andr...@littlepinkcloud.invalid
In old Forths we often used 14-bit fractions for angles, so the range
reduction was a trivial AND operation.
Andrew.
Richard Owlett
>
> Why use 360 degrees, surveyors use 400 degrees in a circle.
HUH! I'm not a surveyor, but have worked with some.
Never heard of a 400 degrees in a circle.
A Google search turned up that the French, when setting up metric
system, defined an angular measure such that 1 revolution would be 400
grads ( "grad" became "gon" due to naming conflict in some areas).
Bernd Paysan
> Why use 360 degrees, surveyors use 400 degrees in a circle.
> 360 is 6 times 60, the base of Babylonian numbering.
> Pi is much younger, a Greek letter and I doubt if they
> used it in relation to circles.
The relation between curcumference and diameter of a circle has a long
history. You find it in the bible ("He [Solomon] made the Sea of cast
metal, circular in shape, measuring ten cubits from rim to rim
[diameter = 10] and five cubits high. It took a line of thirty cubits
to measure around it. [circumference = 30]"), but for sure, older
cultures like Babylon, Egypt, India, and China had their own
approximations of pi.
Babylon: 3 1/8 (3.125)
Egypt: (4/3)^4 (3.1605..)
India: sqrt(10) (3.1622..)
China: 355/113 (3.4159292..) (but that was not until 500 AD by Tsu
Ch'ung-Chi, before, they used sqrt(10))
Note that 3 1/7 is more accurate than the first three ones.
Krishna Myneni
> ...
> 1e0 FATAN 8e0 F* FCONSTANT 2pi
>
> 0 1 D>F 2pi F/ FCONSTANT ANGLEINT \ 2^32/2pi
> 2pi 0 1 D>F F/ FCONSTANT ANGLEFLT \ 2pi/2^32
>
> : FMOD F/ ANGLEINT F* \ angle to cell scale
> F>D D>S 0 D>F \ mod the angle
> ANGLEFLT F* ; \ restore to radians
>
Brad,
The above is not a floating point modulus function.
You seem to have intended this word to wrap angles,
but at the moment I can't seem to figure out how
it does this. Anyway, it should not be called FMOD.
E.g.
5e 3e FMOD F.
2 ok ( in analogy with the MOD function )
Your definition gives something different.
Krishna
Krishna Myneni
>
> The question really is; are floating point appropriate for the kind of
> accuracy you need. You seem to be very specific about it being
> accurate. If this is the case, you might want a solution involving
> symbolic processing in conjunction with keeping for example 2 32 bit
> integers as unresolved divisions for each floating point value instead
> of an actual floating point values which is bound to be very imprecise
> to a certain extent. The math program MAPLE works in this manner as
> far as the symbolic processing goes and the results are of course dead
> on the mark. Symbolically speaking that is :).
>
> Of course this requires developping a whole new library of functions
> to process unresolved fractions but in the end, it will definately
> give you what you want.
>
> Stonelock
Actually, the accuracy for my needs only has to be about seven
significant digits so just about any implementation of FMOD is fine, and
I don't require symbolic computation to extend the precision out to an
arbitrary number of decimal places. But since we are talking about
source code implementations, it is worthwhile to compare the precision
of their results. Actual implementation in a Forth system should rely on
the fpu instructions, which will provide the greatest number of
significant digits.
Krishna
Coos Haak
> THX, I didn't know, that 355/113 origins from China. But you made a
> typo, it's not 3.4159292 but 3.141592920, the 1st '1' is missing :-(
>
The mayor of Alkmaar (famous for its cheese market) Adriaan Metius (name
probably derived from surveyor) rediscovered the approximation:
http://www.pimath.de/quadratur/pi_geschichte2.html
> However I'm using it in Forth since long years. AFAIR I saw it 1st in
> Brodie's Starting forth: 355 113 */ works nearly perfect.
>
> 355/133 = 3.14159292
> pi = 3.14159265
> Error = 0.00000027 = 8.6E-8 an excellent approximation.
>
> Bye from Germany
>
> Wolfgang
Bye from the Netherlands
Coos
Bernd Paysan
> THX for that interesting link. Seems it is not all cheese what that
> mayor had in mind :-)
It probably was all cheese. He wanted to know the ratio between
perimeter and diameter of a cheese ;-).
Anton Ertl
>
>After re-reading DPANS94, Floating Point Extensions,
>I find there is no equivalent of the C library function
>fmod() which returns the remainder (fractional part) of
>the result of dividing two fp numbers. Needing this, I
>had to resort to the clumsy and inefficient source level
>definition
>
>: FMOD ( f1 f2 -- f | remainder of f1/f2 )
> FOVER FOVER F/ F>D D>F F* F- ;
It divides flags?-) Floats are indicated by an r prefix.
This is an explicitly round-towards-zero-symmetric operation, whereas
MOD in Forth does not define division rounding and FM/MOD is
explicitly floored. Maybe it should be called FREM in analogy with
SM/REM (OTOH, there is the C name fmod() as precedent). OTOH, there
is a C function drem(), which works with round-to-nearest division.
Is there any reason why you find the round-to-0 variant preferable
over the other two?
>FMOD is particularly useful for calculations involving
>angles, which are quite common.
Hmm, a round-towards-zero 2pi FMOD produces results between 2pi and
-2pi, whereas a round-to-nearest 2pi FMOD produces results between pi
and -pi, and a floored 2pi FMOD produces results between 2pi and 0.
- anton
--
M. Anton Ertl http://www.complang.tuwien.ac.at/anton/home.html
comp.lang.forth FAQs: http://www.complang.tuwien.ac.at/forth/faq/toc.html