automated FSL tests
Krishna Myneni
Using Anton Ertl's new ttester, I modified the test code in a few of the FSL
modules, to perform automated tests using the style
T{ ... -> ... }T
Since these demonstration modules involve floating point matrices and arrays, I
created some helper words to generate the Hayes'-like tests, element by element,
for matrices and arrays. Thus, two matrices may be compared using a statement
such as,
4 4 CompareMatrices a{{ b{{
where CompareMatrices generates the element by element tests,
t{ a{{ 0 0 }} F@ -> b{{ 0 0 }} F@ }
t{ a{{ 0 1 }} F@ -> b{{ 0 1 }} F@ }
:
:
Several advantages of the automated tests became apparent during and after
implementation:
1) There is no need to visually compare the output with the expected output.
2) The absolute and relative nearness tolerances may be set for the floating
point comparison.
3) The precision of the comparison may be set to a high value, i.e. the nearness
tolerance may be on the order of 1e-14. Imagine trying to compare two 4x4
matrices visually, to 14 or 15 significant digits for each element, and the
advantage becomes obvious.
4) Along with 3), we may specify more exact reference values for the
computations. For instance, we should no longer be sloppy and use "0.08333333e",
when the reference value is actually 1/12 (that is, 1e 12e F/ ).
5) A long standing mistake in the test code for one of the FSL modules, lufact,
was uncovered by the automated tests. The reference value for the lowest
diagonal element of the LU factored 4x4 matrix was given as "0.00357e", when, in
fact, it is "0.000357..." (or 1/2800). Apparently, this was never caught in a
visual comparison.
FSL modules, featuring new test code, were demonstrated under Gforth:
hilbert.fs
lufact.fs
dets.fs
backsub.fs
invm.fs
These are related modules. The utilities for array and matrix comparison,
fsl-test-utils.fs
is also needed. Currently, there is a dependency on strings.fs (providing the
words STRCPY and STRCAT), which may be eliminated. I found these words to be
convenient, but an equivalent rework may be done using words like PLACE and +PLACE.
Finally, there is a wrapper to run tests for all five of the FSL demonstration
modules,
fsl-tester.fs
The files may be found at:
ftp://ccreweb.org/software/gforth/fsl/
Currently, we lack a method for looping a block of test lines, with a parameter
N, typically used to set the size of a matrix. For example, in invm.fs, we have
the following test code:
TESTING INVM
4 value N
t{ & mat{{ N N }}malloc -> }t
t{ malloc-fail? -> 0 }t
t{ mat{{ N hilbert -> }t
t{ & lmat{{ N N }}malloc -> }t
t{ malloc-fail? -> 0 }t
t{ & piv{ N }malloc -> }t
t{ malloc-fail? -> 0 }t
t{ lui lmat{{ piv{ N lumatrix-init -> }t
t{ mat{{ lui lufact -> }t
t{ lui mat{{ invm -> }t
t{ & mat-ref{{ N N }}malloc -> }t
t{ mat-ref{{ N hilbert-inv -> }t
N N CompareMatrices mat{{ mat-ref{{
t{ & mat{{ }}free -> }t
t{ & lmat{{ }}free -> }t
t{ & piv{ }free -> }t
This block tests inverting a 4x4 Hilbert matrix using lufact. However, we should
really conduct tests for N = 2, 3, ..., 6 and higher possibly. The block of test
lines should be looped with an incrementing N.
Below is given the demonstration output of fsl-tester.fs, using Gforth. The same
tests as in the original test code are being performed.
--
krishna@ungee:~/ccreweb/ftp/software/gforth/fsl> gforth fsl-tester.fs
FSL-UTIL V1.19 20 August 2007 EFC, KM
DYNMEM V1.9b 16 July 2004 EFC
HILBERT V1.1 13 October 1994 EFC
TESTING HILBERT HILBERT-INV HILBERT-DET
LUFACT V1.6d 31 Oct 2006 EFC, KM
TESTING LUMATRIX-INIT LUFACT LU-MALLOC LU-FREE
DETS V1.0b 23 August 2007 EFC redefined mat{{
TESTING DET
BACKSUB V1.2b 23 August 2007 EFC redefined a{{ redefined
pivot{ redefined mat{{ redefined lmat{{ redefined piv{
TESTING problem 1
TESTING problem 2
INVM V1.2b 23 August 2007 EFC redefined a{{ redefined
mat{{ redefined lmat{{ redefined piv{
TESTING INVM
redefined N krishna@ungee:~/ccreweb/ftp/software/gforth/fsl>
--
All of the above tests pass. One can see that the tests are working by
introducing an error by hand, e.g. inserting a line to change the matrix element
prior to the line containing CompareMatrices, and so on.
Krishna Myneni
Krishna Myneni
> Using Anton Ertl's new ttester, I modified the test code in a few of the FSL
> modules, to perform automated tests using the style
>
> T{ ... -> ... }T
The kForth version of the same FSL modules have been revised accordingly. Also,
I'm deprecating the old kForth matrix package in favor of FSL-style arrays and
matrices in non-FSL code. See the kForth page (link below) for details.
There hasn't been much of a reaction (yea or nay) to the idea of automating the
FSL test code. I had hoped that some of the more prominent FSL
users/supporters/maintainers would weigh in with their opinions of this
approach. Nevertheless, for my own use, I plan to extend the automated tests to
other FSL modules. The rough program would be to start with the ODE solver
(runge4), then gradually revise the test code of other modules, in order of
relative importance. There appear to be several Fortran-based testers of ODE
solvers. Right now, I'm looking at implementing the Enright and Pryce test
suite. See http://people.scs.fsu.edu/~burkardt/f77_src/toms648_nsd/toms648_nsd.html
Krishna Myneni
--
kForth Page: http://ccreweb.org/software/kforth/
Gerry
> Using Anton Ertl's new ttester, I modified the test code in a few of the FSL
> modules, to perform automated tests using the style
>
> T{ ... -> ... }T
>
> Since these demonstration modules involve floating point matrices and arrays, I
> created some helper words to generate the Hayes'-like tests, element by element,
> for matrices and arrays. Thus, two matrices may be compared using a statement
> such as,
>
> 4 4 CompareMatrices a{{ b{{
>
> where CompareMatrices generates the element by element tests,
>
> t{ a{{ 0 0 }} F@ -> b{{ 0 0 }} F@ }
> t{ a{{ 0 1 }} F@ -> b{{ 0 1 }} F@ }
> :
> :
>
Why do you have to generate test lines as text (if I understand you
correctly), why can't you compile t{ and }t inside, say, a do loop and
iterate through the matrices by executing the compiled word. I've done
this with the original Hayes { and } when testing integer array
operations.
Gerry
Krishna Myneni
Hmm... I didn't think it would work, although I can't say why. I tried the
following, and it failed with a seg fault error:
\ -----------------------------
( snippet from a test version of fsl-test-utils.4th )
0 value Nrows
0 value Ncols
\ Generate element by element tests for two fp matrices.
0 [IF]
: CompareMatrices ( n m <m1> <m2> -- )
to Ncols to Nrows
bl word s1 strcpy bl word s2 strcpy
Nrows 0 DO
Ncols 0 DO
J 0 <# #S #> s" " strcat I 0 <# #S #> strcat s" " strcat s" }}
F@ " strcat
s" t{ " s1 count strcat s" " strcat 2over strcat
s" -> " strcat s2 count strcat s" " strcat
2swap strcat s" r}t" strcat
( type ) evaluate
LOOP
LOOP
;
[THEN]
0 ptr m1{{
0 ptr m2{{
: CompareMatrices ( n m m1 m2 -- )
to m2{{ to m1{{ to Ncols to Nrows
Nrows 0 DO
Ncols 0 DO
t{ m1{{ J I }} F@ -> m2{{ J I }} F@ r}t
LOOP
LOOP
;
\ ---------------------------
Can you post an example of how you used the original test words inside of a DO loop?
Krishna
Gerry
This is taken from a memory allocation test (the first I could find):
\ --------------------------------------------------------------------
Testing ALLOCATE FREE RESIZE
VARIABLE addr
{ 100 ALLOCATE SWAP addr ! -> 0 }
{ addr @ 1 CELLS 1- and -> 0 } \ Test address is aligned
{ addr @ FREE -> 0 }
{ 99 ALLOCATE SWAP addr ! -> 0 }
{ addr @ 1 CELLS 1- AND -> 0 } \ Test address is aligned
{ addr @ FREE -> 0 }
{ 50 ALLOCATE SWAP addr ! -> 0 }
: writemem 0 DO I 1+ OVER C! 1+ LOOP DROP ; ( ad n -- )
: checkmem ( ad n --- )
0
DO
DUP C@ SWAP >R
{ -> R> I 1+ SWAP >R } \ Luckily tester checks can be compiled
R> 1+
LOOP
DROP
;
addr @ 50 writemem addr @ 50 checkmem
{ addr @ 28 RESIZE SWAP addr ! -> 0 }
addr @ 28 checkmem
{ addr @ 200 RESIZE SWAP addr ! -> 0 }
{ addr @ 28 checkmem
\ --------------------------------------------------------------------
er, what's the r}t in your code?
Gerry
Marcel Hendrix
[..]
> 0 ptr m1{{
> 0 ptr m2{{
> : CompareMatrices ( n m m1 m2 -- )
> to m2{{ to m1{{ to Ncols to Nrows
> Nrows 0 DO
> Ncols 0 DO
> t{ m1{{ J I }} F@ -> m2{{ J I }} F@ r}t
> LOOP
> LOOP
> ;
: CompareMatrices ( n m m1 m2 -- ) 2swap * floats compare throw ;
-marcel
Krishna Myneni
The example above works in the original tester, but does not appear to work in
Anton's revised ttester and tester wrapper. Perhaps Anton could comment on this.
Under ttester, it is required to specify the expected types on the stack for
Forth systems which use a unified fp/data stack. The "r}t" indicates that a real
number is expected on the data stack. For systems which have separated data and
fp stacks, one may use "}t". Real number comparisons are different from ordinary
data stack comparisons, since a nearness test is typically needed rather than
exact matching.
Krishna
Krishna Myneni
That's an exact comparison. What we want typically is a nearness comparison.
Krishna
Anton Ertl
>Gerry wrote:
...
>> { -> R> I 1+ SWAP >R } \ Luckily tester checks can be compiled
>The example above works in the original tester, but does not appear to work in
>Anton's revised ttester and tester wrapper. Perhaps Anton could comment on this.
It appears that the reason is that my revised version takes {
seriously as the start of the test, and compares only the stuff that
has been pushed on the stack(s) between { and -> to the stuff that has
been pushed between -> and }.
In contrast, John Hayes' version ignores { and always consumes all of
the stack at ->, and compares that to the stuff on the stacks at }.
The present example was apparently programmed to accomodate this
misfeature of Hayes' tester, but I don't think that this is an
intended use of Hayes' tester (otherwise, why provide { at all?).
It is certainly possible to write an equivalent test for the revised
version of the tester, and probably also a version that works for both
of them, but the code is not clear enough that I feel inspired to do
it.
Anyway, you can use the revised tester in compiled code, but the error
messages (if any) will not be as enlightening as for interpreted
tests. Example:
require test/ttester.fs
: foo
0
100 0 do ( i^2 )
t{ dup -> i i * }t
i 2* + 1+ ( [i+1]^2 )
loop
drop ;
foo \ no error messages expected
>Under ttester, it is required to specify the expected types on the stack for
>Forth systems which use a unified fp/data stack.
At least if you want approximate matching.
>The "r}t" indicates that a real
>number is expected on the data stack.
Actually, it indicates that an FP number is expected on the
appropriate stack (it works for both separate and combined stacks).
- 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
New standard: http://www.forth200x.org/forth200x.html
EuroForth 2007: http://www.complang.tuwien.ac.at/anton/euroforth2007/
Krishna Myneni
> Krishna Myneni <krishn...@bellsouth.net> writes:
>> Gerry wrote:
> ...
>>> { -> R> I 1+ SWAP >R } \ Luckily tester checks can be compiled
> ...
>> The example above works in the original tester, but does not appear to work in
>> Anton's revised ttester and tester wrapper. Perhaps Anton could comment on this.
>
> It appears that the reason is that my revised version takes {
> seriously as the start of the test, and compares only the stuff that
> has been pushed on the stack(s) between { and -> to the stuff that has
> been pushed between -> and }.
>
> In contrast, John Hayes' version ignores { and always consumes all of
> the stack at ->, and compares that to the stuff on the stacks at }.
> The present example was apparently programmed to accomodate this
> misfeature of Hayes' tester, but I don't think that this is an
> intended use of Hayes' tester (otherwise, why provide { at all?).
>
> It is certainly possible to write an equivalent test for the revised
> version of the tester, and probably also a version that works for both
> of them, but the code is not clear enough that I feel inspired to do
> it.
>
For my own use, I can live without tests being compilable into word defs.
> ...
>> Under ttester, it is required to specify the expected types on the stack for
>> Forth systems which use a unified fp/data stack.
>
> At least if you want approximate matching.
>
>> The "r}t" indicates that a real
>> number is expected on the data stack.
>
> Actually, it indicates that an FP number is expected on the
> appropriate stack (it works for both separate and combined stacks).
>
> - anton
Thanks for the above clarifications Anton.
Krishna
Gerry
> Krishna Myneni <krishnamyn...@bellsouth.net> writes:
> >Gerry wrote:
> ...
> >> { -> R> I 1+ SWAP >R } \ Luckily tester checks can be compiled
> ...
> >The example above works in the original tester, but does not appear to work in
> >Anton's revised ttester and tester wrapper. Perhaps Anton could comment on this.
>
> It appears that the reason is that my revised version takes {
> seriously as the start of the test, and compares only the stuff that
> has been pushed on the stack(s) between { and -> to the stuff that has
> been pushed between -> and }.
>
> In contrast, John Hayes' version ignores { and always consumes all of
> the stack at ->, and compares that to the stuff on the stacks at }.
> The present example was apparently programmed to accomodate this
> misfeature of Hayes' tester,
I didn't realise it was a misfeature and don't see why you say it is,
that's how it was designed and used. I agree your version is an
improvement, but clearly not a backwards compatible replacement.
> but I don't think that this is an
> intended use of Hayes' tester (otherwise, why provide { at all?).
>
John Hayes said why in the original tester.fr - syntactic sugar (at
least that's what it says in the copy I have).
Gerry
Anton Ertl
>Anton Ertl wrote:
>> Krishna Myneni <krishn...@bellsouth.net> writes:
>>> Gerry wrote:
>> ...
>>>> { -> R> I 1+ SWAP >R } \ Luckily tester checks can be compiled
>> ...
>>> The example above works in the original tester, but does not appear to work in
>>> Anton's revised ttester and tester wrapper. Perhaps Anton could comment on this.
>>
>> It appears that the reason is that my revised version takes {
>> seriously as the start of the test, and compares only the stuff that
>> has been pushed on the stack(s) between { and -> to the stuff that has
>> been pushed between -> and }.
>>
>> In contrast, John Hayes' version ignores { and always consumes all of
>> the stack at ->, and compares that to the stuff on the stacks at }.
>> The present example was apparently programmed to accomodate this
>> misfeature of Hayes' tester, but I don't think that this is an
>> intended use of Hayes' tester (otherwise, why provide { at all?).
>>
>> It is certainly possible to write an equivalent test for the revised
>> version of the tester, and probably also a version that works for both
>> of them, but the code is not clear enough that I feel inspired to do
>> it.
>>
>
>
>For my own use, I can live without tests being compilable into word defs.
Note that the problem with Gerry's code is not that the test is
compiled into a word definition, but that it uses { at a place where
the stack is not empty. The two tester versions behave differently in
this case, both interpreted and inside a colon definition. E.g.:
1
{ 1 -> 1 } \ John Hayes' tester does not like this
clearstack
1
{ -> 1 } \ the revised tester does not like this
Note that both testers produce one error message for one of the tests
in this code, but they produce it for different tests.
Anton Ertl
>On 22 Sep, 16:46, an...@mips.complang.tuwien.ac.at (Anton Ertl) wrote:
>> It appears that the reason is that my revised version takes {
>> seriously as the start of the test, and compares only the stuff that
>> has been pushed on the stack(s) between { and -> to the stuff that has
>> been pushed between -> and }.
>>
>> In contrast, John Hayes' version ignores { and always consumes all of
>> the stack at ->, and compares that to the stuff on the stacks at }.
>> The present example was apparently programmed to accomodate this
>> misfeature of Hayes' tester,
On second thought, it's probably just a bug. A misfeature implies the
intent that it be used, but I don't think that was the intent (see
below).
>I didn't realise it was a misfeature and don't see why you say it is,
>that's how it was designed and used.
That's how it was implemented, but not used: All of John Hayes' tests
work with the revised version, so he did not rely on this bug; i.e.,
he did not use it.
>I agree your version is an
>improvement, but clearly not a backwards compatible replacement.
It is backwards compatible for all correct uses. If I had intended
bug-for-bug compatibility, I would not have fixed the bug. You can
easily reintroduce the bug by replacing START-DEPTH @ with 0 in
various places.
>> but I don't think that this is an
>> intended use of Hayes' tester (otherwise, why provide { at all?).
>>
>
>John Hayes said why in the original tester.fr - syntactic sugar (at
>least that's what it says in the copy I have).
Yes, but why have the syntactic sugar at all if it was intended to
have no meaning? No, it was supposed to have a meaning: the start of
the test, and that's what's implemented in the revised version.
Krishna Myneni
> ...
> 5) A long standing mistake in the test code for one of the FSL modules, lufact,
> was uncovered by the automated tests. The reference value for the lowest
> diagonal element of the LU factored 4x4 matrix was given as "0.00357e", when, in
> fact, it is "0.000357..." (or 1/2800). Apparently, this was never caught in a
> visual comparison.
>
>
An updated version of the module GAULEG has been added to the automated tests
for the Gforth collection of FSL modules at ccreweb.org:
ftp://ccreweb.org/software/gforth/fsl/
In performing the automated comparison of the computed G-L weights with high
precision reference values, I found to my surprise that the FSL version only
calculated these weights to an accuracy of about 6 significant digits. Without
any comment, the FSL module had quietly set the tolerance (EPS) to about 1e-6,
even though the original routine from Numerical Recipes had EPS set to 3e-11.
It is fairly clear to me, from the cases I have encountered so far, that the FSL
is in need of rigorous tests.
Krishna Myneni
Marcel Hendrix
> Krishna Myneni wrote:
>> ...
>> 5) A long standing mistake in the test code for one of the FSL modules, lufact,
>> was uncovered by the automated tests. The reference value for the lowest
>> diagonal element of the LU factored 4x4 matrix was given as "0.00357e", when, in
>> fact, it is "0.000357..." (or 1/2800). Apparently, this was never caught in a
>> visual comparison.
In hindsight, it is unfortunate that an algorithm's change history
does not clearly state that an algorithm was independently tested.
I hope this will be changed when you notify cgm of the current problems.
BTW, in visually checking a few of the files I noticed lots of formatting
mistakes and inconsistencies (non-adherence to FSL conventions) that could
use some cleanup. I do not remember that it was so bad, but of course, the
FSL is almost 15 years old already...
> An updated version of the module GAULEG has been added to the automated tests
> for the Gforth collection of FSL modules at ccreweb.org:
> ftp://ccreweb.org/software/gforth/fsl/
> In performing the automated comparison of the computed G-L weights with high
> precision reference values, I found to my surprise that the FSL version only
> calculated these weights to an accuracy of about 6 significant digits. Without
> any comment, the FSL module had quietly set the tolerance (EPS) to about 1e-6,
> even though the original routine from Numerical Recipes had EPS set to 3e-11.
More likely this is because the NR (Fortran) that is quoted as a reference (1986)
did everything in single precision (for speed). Even the NR in C second edition,
1992, still did this. In my adaptation for iForth I rewrote everything in double
precision.
-marcel
Krishna Myneni
>> precision reference values, I found to my surprise that the FSL version only
>> calculated these weights to an accuracy of about 6 significant digits. Without
>> any comment, the FSL module had quietly set the tolerance (EPS) to about 1e-6,
>> even though the original routine from Numerical Recipes had EPS set to 3e-11.
>
> More likely this is because the NR (Fortran) that is quoted as a reference (1986)
> did everything in single precision (for speed). Even the NR in C second edition,
> 1992, still did this. In my adaptation for iForth I rewrote everything in double
> precision.
>
> -marcel
>
I went back to NR in C, 2nd ed. to try to find the source of the low accuracy
weights. It was then that I noticed the discrepancy in the setting of EPS. Maybe
the 1st ed. did set the tolerance lower, but I no longer have a copy of the
first ed. to verify this. In any case, it is useful to know the expected
accuracy of the computation, which is what the automated tests help us to determine.
You should check the EPS setting in your iForth version of gauleg. The double
precision is not doing you any good if this setting is 3.0e-6, as in the current
gauleg from the main FSL site.
Krishna
Marcel Hendrix
> Marcel Hendrix wrote:
[..]
> Even the NR in C second edition,
> 1992, still did this. In my adaptation for iForth I rewrote everything in double
> precision.
[..]
Sorry, the NRC 1992 2nd edition already made gauleg double precision.
> You should check the EPS setting in your iForth version of gauleg. The double
> precision is not doing you any good if this setting is 3.0e-6, as in the current
> gauleg from the main FSL site.
I have a system constant eps ( 2.2204e-016 ) that works in most
cases (see below). In the reworking I also reviewed all magic numbers.
In the printout, the X and W values depend on PRECISION, which says nothing about
the actual accuracy of the values. The final value of the integral should of course
be printed out with the maximum number of available digits.
FORTH> 8 gauleg-test
X:
-9.602899e-001 -7.966665e-001 -5.255324e-001 -1.834346e-001 1.834346e-001 5.255324e-001 ...
7.966665e-001 9.602899e-001
W:
1.012285e-001 2.223810e-001 3.137066e-001 3.626838e-001 3.626838e-001 3.137066e-001 ...
2.223810e-001 1.012285e-001
The values for n = 8 are:
X:
-9.602899e-001 -7.966665e-001 -5.255324e-001 -1.834346e-001 1.834346e-001 5.255324e-001 ...
7.966665e-001 9.602899e-001
W:
1.012285e-001 2.223810e-001 3.137066e-001 3.626838e-001 3.626838e-001 3.137066e-001 ...
2.223810e-001 1.012285e-001
The integral from -1.0 to 1.0 is : 6.66666666666668e-001
(the actual value is : 6.66666666666667e-001)
ok
-marcel
==========================
(*
* LANGUAGE : ANS Forth with DFW extensions
* PROJECT : Forth Environments
* DESCRIPTION : Gauss-Legendre Integration
* CATEGORY : FSL
* AUTHOR : Thursday, March 30, 2006, 23:22 PM, Marcel Hendrix
* LAST CHANGE : Friday, March 31, 2006, 7:02 AM, Marcel Hendrix
*)
NEEDS -fsl_util
REVISION -gauleg "ÄÄÄ Gauss-Legendre Int. Version 1.10 ÄÄÄ"
PRIVATES
DOC
(*
Gauss-Legendre Integration
Forth Scientific Library Algorithm #27
Given lower and upper limits of integration X1 and X2 the routine gauleg returns the abscissas
and weights of the Gauss-Legendre N-point quadrature formula.
The integral of f(x) is then calculated numerically as:
z = \int_x1^x2 f(x) dx = \sum_i=0^n-1 w[i] f( x[i] )
the routine )gl-integrate is provided to do this calculation using the previously calculated values
of x and w.
This program implements the function GAULEG as described in
Press, W.H., B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, 1986;
Numerical Recipes, The Art of Scientific Computing, Cambridge University
Press, 818 pages, ISBN 0-521-30811-9
*)
ENDDOC
0e FVALUE z PRIVATE
0e FVALUE p1 PRIVATE
0e FVALUE pp PRIVATE
0 VALUE x{ PRIVATE -- aliases to user arrays
0 VALUE w{ PRIVATE
DEFER f(x) PRIVATE -- pointer to user function ( F: x -- f[x] )
: calc-pp ( n -- ) ( F: -- r ) \ NOTE: changes Z
LOCAL N
BEGIN
1e TO p1
0e TO pp
N 1+ 1 DO
pp p1 p1 TO pp
I 2* 1- S>F F* z F*
FSWAP I 1- NEGATE S>F F* F+
I S>F F/ TO p1
LOOP
z p1 F* pp F- N S>F F*
z FSQR F1- F/ TO pp
p1 pp F/ FNEGATE
z FDUP FROT F+ FDUP TO z
F- FABS eps F<
UNTIL
pp FSQR ; PRIVATE
: gauleg ( x{ w{ -- ) ( F: x1 x2 -- )
DUP DIMS MAX LOCAL N TO w{ TO x{
N 1 < ABORT" bad value of matrix size (must be > 0) for gauleg"
F2DUP F+ F2/ FLOCAL xm
FSWAP F- F2/ FLOCAL xl
N 1+ 2/
0 DO
PI N S>F 0.5e F+ F/
I S>F 0.75e F+ F* FCOS TO z
N calc-pp
z FSQR FNEGATE F1+ F*
xl F2* FSWAP F/
FDUP w{ I } DF! w{ N 1- I - } DF!
xl z F* FDUP FNEGATE xm F+ x{ I } DF!
xm F+ x{ N 1- I - } DF!
LOOP ;
-- do the integration
: gl-integrate ( xt x{ w{ -- ) ( F: -- z )
DUP DIMS MAX >R TO w{ TO x{ [IS] f(x)
R@ 1 < ABORT" bad value of matrix size (must be > 0) for gl-integrate"
0e R>
0 ?DO
x{ I } DF@ f(x)
w{ I } DF@ F*
F+
LOOP ;
NESTING @ 1 = [IF] -- test code ==============================================
DOUBLE DARRAY xx{
DOUBLE DARRAY wgt{
: func ( -- ) ( F: x -- f[x] ) \ function to integrate
FDUP 0.2e F+ F* ;
: ifunc ( -- ) ( F: x -- I[x] ) \ its (indefinite) integral
FDUP 3e F/ 0.1e F+
FOVER F* F* ;
: values-for-8 ( -- )
CR ." The values for n = 8 are:"
CR ." X:"
CR ." -9.602899e-001 -7.966665e-001 -5.255324e-001 -1.834346e-001 1.834346e-001 5.255324e-001 ..."
CR ." 7.966665e-001 9.602899e-001"
CR ." W:"
CR ." 1.012285e-001 2.223810e-001 3.137066e-001 3.626838e-001 3.626838e-001 3.137066e-001 ..."
CR ." 2.223810e-001 1.012285e-001" ;
: gauleg-test ( n -- ) \ if n = 8, comparison values are also given
LOCAL N
xx{ N }malloc
wgt{ N }malloc
xx{ wgt{ -1e 1e gauleg
CR ." X: " xx{ }print
CR ." W: " wgt{ }print
N 8 = IF values-for-8 ENDIF
['] func xx{ wgt{ gl-integrate
CR ." The integral from -1.0 to 1.0 is : " +E.
CR ." (the actual value is : " 1e ifunc -1e ifunc F- +E. ." )" CR
xx{ }free
wgt{ }free ;
[THEN]
(* End of File *)
Gerry
To summarise as I see it:
1. When John Hayes developed the tester he made the design decision
that the word -> would empty the data stack - a perfectly reasonable
decision in my view.
2. The Hayes tester was generally accepted and used in the Forth
community, so much so that the Forth 200x team decided to adopt it as
the standard tester - again a perfectly reasonable decision.
3. You decided to improve the tester to handle floating point tests
(good idea) and change the specification of { and -> to not empty the
data stack (a reasonable idea but showing cavalier disregard for any
code the change in specification might break).
4. When presented with an example of some test code broken by the
change in specification you loftily declare the original
implementation of { and -> to contain a bug and that any code using
the fact that -> empties the data stack to be a misuse of the tester
and by implication bad. I think that a matter of opinion and disagree
with you.
Having said that I'm willing to live with the change as it is a better
implementation and no big deal to alter my code, or I can just ignore
it - but let's hope that nobody else is affected by your change.
Gerry
(my last post on the topic)
Krishna Myneni
> ... The final value of the integral should of course
> be printed out with the maximum number of available digits.
>
The accuracy of the weights will (should, in most cases) set an upper limit on
the accuracy of the computed integral. Checking the accuracy of the weights is
best done with the automated nearness test, as in the Gforth and kForth versions
of GAULEG.
Krishna
Krishna Myneni
> ... The rough program would be to start with the ODE solver
> (runge4), ...
A preliminary version of RUNGE4, with automated tests, has been posted now. It
turns out that the comment below, from the original (v1.2) test code, regarding
the accuracy of the test case solutions is highly misleading:
FVARIABLE _dt \ With this time increment the dampled vibration test
1.0E-2 _dt F! \ will be in error in the 4th decimal place. Smaller
\ values will give more accurate results.
The test cases for the damped vibrations and for the charging capacitor have
solutions which converge to a fixed point in the long time limit. Therefore, the
accuracy of the test cases is dependent not only on the value of _dt, but also
highly dependent on the ending time (the nearness of the computed solution and
actual solution improves with a longer end time, due to the convergence).
In fact, for short times, we barely get 1 significant digit of accuracy for the
damped vibration problem, using the above value of _dt!
The adaptive step-size solver appears to work very well, and the tolerance for
these tests is set at 1e-14.
Krishna Myneni
Marcel Hendrix
> A preliminary version of RUNGE4, with automated tests, has been posted now. It
> turns out that the comment below, from the original (v1.2) test code, regarding
> the accuracy of the test case solutions is highly misleading:
> FVARIABLE _dt \ With this time increment the dampled vibration test
> 1.0E-2 _dt F! \ will be in error in the 4th decimal place. Smaller
> \ values will give more accurate results.
Something went wrong:
\ http://www.taygeta.com/fsl/library/runge4.seq
FVARIABLE _dt \ With this time increment the dampled vibration test
1.0E-2 _dt F! \ will be in error in the 4th decimal place. Smaller
\ values will give more accurate results.
\ ftp://ccreweb.org/software/gforth/fsl/runge4.fs
50.0E-3 _dt F!
1e-1 rel-near F! \ extremely low accuracy with coarse fixed-step integrator!
1e-1 abs-near F!
Something *should* happen to the accuracy when the timestep is made
five times larger than intended :-)
> The test cases for the damped vibrations and for the charging capacitor have
> solutions which converge to a fixed point in the long time limit. Therefore, the
> accuracy of the test cases is dependent not only on the value of _dt, but also
> highly dependent on the ending time (the nearness of the computed solution and
> actual solution improves with a longer end time, due to the convergence).
> In fact, for short times, we barely get 1 significant digit of accuracy for the
> damped vibration problem, using the above value of _dt!
Are you comparing the last point of the computed trajectory to the fixed point?
I know for certain that the original author meant the discrepancy between output
points when computed with different values of _dt. In dvib_test_twice (only on
taygeta) the value of the timestep is halved again (to 5 ms) to make this point
clear.
-marcel
Krishna Myneni
> Krishna Myneni <krishn...@bellsouth.net> writes Re: automated FSL tests
> [..]
>> A preliminary version of RUNGE4, with automated tests, has been posted now. It
>> turns out that the comment below, from the original (v1.2) test code, regarding
>> the accuracy of the test case solutions is highly misleading:
>
>
>> FVARIABLE _dt \ With this time increment the dampled vibration test
>> 1.0E-2 _dt F! \ will be in error in the 4th decimal place. Smaller
>> \ values will give more accurate results.
>
> Something went wrong:
>
> \ http://www.taygeta.com/fsl/library/runge4.seq
> FVARIABLE _dt \ With this time increment the dampled vibration test
> 1.0E-2 _dt F! \ will be in error in the 4th decimal place. Smaller
> \ values will give more accurate results.
>
> \ ftp://ccreweb.org/software/gforth/fsl/runge4.fs
> 50.0E-3 _dt F!
> 1e-1 rel-near F! \ extremely low accuracy with coarse fixed-step integrator!
> 1e-1 abs-near F!
>
> Something *should* happen to the accuracy when the timestep is made
> five times larger than intended :-)
>
Thanks. I thought I had actually made the time step smaller, when in fact I had
increased it (50e-3 > 1e-2). Going back to the original time step of 1e-2, and
choosing an end time of 0.8, the test works.
1.0E-2 _dt F!
1e-4 rel-near F! \ relatively low accuracy with coarse fixed-step integrator!
1e-4 abs-near F!
\ Don't use too many steps for the test, since the solution converges in
\ the long time limit!
TESTING Damped Vibration (Fixed Step)
t{ 80 dvib_test x{ 0 } F@ -> t_final F@ actual r}t
The test fails if the nearness is set to 1e-5. The original point I made about
these test cases is, I believe, still valid. These are not good test cases
because I can integrate with the same time step for a longer time, say 1000
steps rather than 80, which would take me to an end time of 10. Then, I can set
the nearness to 1e-6, and the test will pass. One may conclude, improperly, that
the accuracy is 6 significant digits along the entire trajectory.
>> The test cases for the damped vibrations and for the charging capacitor have
>> solutions which converge to a fixed point in the long time limit. Therefore, the
>> accuracy of the test cases is dependent not only on the value of _dt, but also
>> highly dependent on the ending time (the nearness of the computed solution and
>> actual solution improves with a longer end time, due to the convergence).
>
> ...
> Are you comparing the last point of the computed trajectory to the fixed point?
> ...
No. We are comparing the last point of the computed trajectory with the actual
solution at the last time (t_final). The integrator, dvib_test, was modified to
store the final time in t_final. We especially need to use t_final for the
adaptive step size integrator, since we can't choose the end time, except to
within the max step size.
Krishna
Krishna Myneni
runge4.fs has been revised:
\ 2007-09-27 km -- cleaned up the test code, using a consistent format for
\ both the E&P and non E&P tests; renamed rksolve1 and
\ rksolve2 to )rksolve1 and )rksolve2 for better syntax,
\ and added )rk_fixed1 and )rk_fixed2 wrappers for the
\ fixed step size tests.
Also, to avoid any possible confusion, please note the following regarding the
revised FSL modules at ccreweb.org,
ftp://ccreweb.org/software/gforth/fsl/
The modules in the above directory are ones which I am revising and testing
under Gforth for my own use. They are experimental, and are not sanctioned by
the FSL. They are also not part of the authoritative Gforth distribution, which
currently does not include the FSL.
The official FSL library may be found at
http://www.taygeta.com/fsl/sciforth.html
Having stated this, I hope that through iterative feedback, some of the changes,
particularly those related to the automated test code, will be adopted by other
Forth system providers. Changes which have some common backing can then be fed
back as proposed changes into the standard library.
Krishna Myneni
Krishna Myneni
Versions of the following FSL files, revised for automated testing, are
available at ftp://ccreweb.org/software/gforth/fsl/
regfalsi.fs
expint.fs
horner.fs
polys.fs
Please view the README file concerning the relation of the files in the above
directory to the FSL.
In performing the automated tests for polys.fs, it appears that one of the
reference values given for the generalized Laguerre polynomials is incorrect,
either due to:
a) Mathematica v2.2 gave the wrong answer
or, more likely,
b) There was a transcription error by the author of the original test code.
\ Generalized Laguerre table comparison values generated from Mathematica V2.2
\ for alpha = 0.75
\ 4.41667e 5.90365e ** 7.61446e ** 9.56526e 11.7909e 14.3154e
4.41667e 5.90365e 7.611446e 9.56526e 11.7909e 14.3154e
If anyone has access to Mathematica v2.2, it would be nice to verify the source
of the above error for the entry "7.61446"
The other polynomials (Chebyshev, Hermite, Laguerre, Legendre, Bessel) now have
reference values with accuracy close to the machine precision.
Krishna Myneni
Marcel Hendrix
> In performing the automated tests for polys.fs, it appears that one of the
> reference values given for the generalized Laguerre polynomials is incorrect,
> either due to:
> a) Mathematica v2.2 gave the wrong answer
Yeah, sure :-)
> or, more likely,
> b) There was a transcription error by the author of the original test code.
Mathematica 3.0 ...
ToString[TableForm[Table[FortranForm[N[ChebyshevT[j,-0.25+0.25*i],20]] ,{i,0,4},{j,3,8}]]]
0.6875 0.53125 -0.953125 -0.0546875 0.98046875 -0.435546875
-1.836909530733565e-16 1. 3.061515884555943e-16 -1. -4.28612223837832e-16 1.
-0.6875 0.53125 0.953125 -0.0546875 -0.98046875 -0.435546875
-1. -0.5 0.5 1. 0.5 -0.5
-0.5625 -0.96875 -0.890625 -0.3671875 0.33984375 0.876953125
ToString[TableForm[Table[PaddedForm[N[HermiteH[j,-0.25+0.25*i],20],{14,10}],{i,0,4},{j,3,8}]]]
2.8750000000 9.0625000000 -27.5312500000 -76.8593750000 368.8046875000 891.6289062500
0.0000000000 12.0000000000 0.0000000000 -120.0000000000 0.0000000000 1680.0000000000
-2.8750000000 9.0625000000 27.5312500000 -76.8593750000 -368.8046875000 891.6289062500
-5.0000000000 1.0000000000 41.0000000000 31.0000000000 -461.0000000000 -895.0000000000
-5.6250000000 -9.9375000000 30.0937500000 144.5156250000 -144.3515625000 -2239.7460937500
ToString[TableForm[Table[PaddedForm[N[LaguerreL[j,-0.25+0.25*i],20] ,{13,9}],{i,0,4},{j,3,8}]]]
1.846354167 2.198079427 2.589363607 3.023323907 3.503265732 4.032691883
1.000000000 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000
0.341145833 0.177246094 0.037263997 -0.080940416 -0.179367756 -0.259886435
-0.145833333 -0.330729167 -0.445572917 -0.504144965 -0.518339224 -0.498362998
-0.476562500 -0.580566406 -0.576684570 -0.501364136 -0.383086177 -0.243678635
ToString[TableForm[Table[PaddedForm[N[LegendreP[j,-0.25+0.25*i],20] ,{13,9}],{i,0,4},{j,3,8}]]]
0.335937500 0.157714844 -0.339721680 0.024276733 0.279918671 -0.152454019
0.000000000 0.375000000 0.000000000 -0.312500000 0.000000000 0.273437500
-0.335937500 0.157714844 0.339721680 0.024276733 -0.279918671 -0.152454019
-0.437500000 -0.289062500 0.089843750 0.323242187 0.223144531 -0.073638916
-0.070312500 -0.350097656 -0.416381836 -0.280776978 -0.034183502 0.197609305
ToString[TableForm[Table[PaddedForm[N[LaguerreL[j,0.75,-0.25+0.25*i],20] ,{13,9}],{i,0,4},{j,3,8}]]]
4.416666667 5.903645833 7.611458333 9.565256076 11.790891617 14.315441410
3.007812500 3.571777344 4.107543945 4.620986938 5.116092682 5.595726371
1.833333333 1.796875000 1.672395833 1.483420139 1.248214286 0.981351919
0.877604167 0.506673177 0.103621419 -0.291271634 -0.651256246 -0.959287825
0.125000000 -0.367187500 -0.779687500 -1.077539062 -1.249302455 -1.298576137
ToString[TableForm[Table[FortranForm[N[BesselJ[j,-0.25+0.25*i],20] ],{i,0,4},{j,3,6}]]]
-0.0003242512526759081 0.00001014077825911821 -2.536516158747241e-7 5.286375870751723e-9
0. 0. 0. 0.
0.0003242512526759081 0.00001014077825911821 2.536516158747241e-7 5.286375870751723e-9
0.002563729994587244 0.0001607364763642876 8.05362724135748e-6 3.360684628618848e-7
0.00848438342327411 0.000801070086542314 0.00006036416651057644 3.785466932038245e-6
With these reference values (accurate to 10^-16):
FORTH> test-cheby test-hermite test-laguerre test-legendre test-glaguerre test-bessel
RMS error: 2.11796688046160e-017
RMS error: 1.63418795292645e-014
RMS error: 4.86136936013886e-011
RMS error: 4.69596297060642e-011
RMS error: 5.01989984144475e-011
RMS error: 1.66345245030888e+003 ok
(all absolute errors)
Obviously the Bessel function is suspect. I tried Mathematica's BesselJ, BesselK, BesselK
and BesselY but no luck?
The errors of Laguerre, Legendre and GLaguerre seem rather high, but they are OK for
single precision applications.
-marcel
Marcel Hendrix
> FORTH> test-cheby test-hermite test-laguerre test-legendre test-glaguerre test-bessel
> RMS error: 2.11796688046160e-017
> RMS error: 1.63418795292645e-014
> RMS error: 4.86136936013886e-011
> RMS error: 4.69596297060642e-011
> RMS error: 5.01989984144475e-011
> RMS error: 1.66345245030888e+003 ok
> The errors of Laguerre, Legendre and GLaguerre seem rather high, but they are OK for
> single precision applications.
Correction; by printing out the Mathematica values with a few more digits:
FORTH> test-cheby test-hermite test-laguerre test-legendre test-glaguerre test-bessel
RMS error: 2.11796688046160e-017
RMS error: 1.63418795292645e-014
RMS error: 4.33185611774133e-017
RMS error: 3.29611662544038e-017
RMS error: 2.05078980308681e-016
RMS error: 1.66345245030888e+003 ok
So Hermite might be in for improvement and Bessel is a mystery. The other functions
are close to machine precision.
-marcel
Krishna Myneni
I'm not entirely sure what you are comparing, but, in the Bessel case, note the
difference between Bessel polynomials and "reverse" Bessel polynomials, which is
what the FSL polys.x computes.
Krishna
Marcel Hendrix
[..]
>> So Hermite might be in for improvement and Bessel is a mystery. The other functions
>> are close to machine precision.
> I'm not entirely sure what you are comparing, but, in the Bessel case, note the
> difference between Bessel polynomials and "reverse" Bessel polynomials, which is
> what the FSL polys.x computes.
Aha! IMHO, that would have deserved a comment.
I renamed bessel to rbessel to help my memory that this is the reverse Bessel polynomial.
ToString[TableForm[Table[FortranForm[N[Sum[(2 n - k)!/(2^(n-k) k! (n-k)!) (-0.25+0.25 i)^k,{k,0,n}],20] ],{i,0,4},{n,3,6}]]]
11.60937499999999 81.41015625 733.4169921875001 8072.67504882813
indeterminate indeterminate indeterminate indeterminate
19.140625 134.22265625 1209.2001953125 13309.59106445312
24.125 170.0625 1536.59375 16945.046875
30.046875 213.59765625 1939.2802734375 21452.23168945312
We now have the problem that for i=1 the Be_n needs to evaluate 0.0^k, which is not defined for k=0.
With some luck, iForth evaluates it to the correct limiting value, but I guess that is not
guaranteed to be so for other Forths. (So the testing may give interesting results here.)
ToString[TableForm[Table[FortranForm[N[Sum[(2 n - k)!/(2^(n-k) k! (n-k)!)(-0.25+0.25 i+10^-16)^k,{k,0,n}], 20] ],{i,0,4},{n,3,6}]]]
15. 105. 945. 10395.
With the adapted reference values from Mathematica 3.0 (and the above crude fix for x=0):
FORTH> test-cheby test-hermite test-laguerre test-legendre test-glaguerre test-rbessel
RMS error: 2.11796688046160e-017
RMS error: 1.63418795292645e-014
RMS error: 4.33185611774133e-017
RMS error: 3.29611662544038e-017
RMS error: 2.05078980308681e-016
RMS error: 5.39983240078416e-013 ok
Given the fact that rBessel has values spanning 10 .. 10.000, the absolute error of 5e-13
for rbessel is satisfying.
IMHO, the Hermite question remains to be solved.
-marcel
Marcel Hendrix
[..]
> FORTH> test-cheby test-hermite test-laguerre test-legendre test-glaguerre test-rbessel
> RMS error: 2.11796688046160e-017
> RMS error: 1.63418795292645e-014
> RMS error: 4.33185611774133e-017
> RMS error: 3.29611662544038e-017
> RMS error: 2.05078980308681e-016
> RMS error: 5.39983240078416e-013 ok
> Given the fact that rBessel has values spanning 10 .. 10.000, the absolute error of 5e-13
> for rbessel is satisfying.
> IMHO, the Hermite question remains to be solved.
Another oeps, sorry. I hadn't noticed that the range of Hermite is between 3 and 1680.
So all is fine.
Here is a slightly faster He_n (my penance ..)
-- nth order Hermite Polynomial
-- Note: He_n+1 = 2 * x * He_n - 2 * n * He_n-1; He_0 = 1, He_1 = 2*x
: He_n ( n -- ) ( F: x -- z )
F2* FDUP 1e 1e FLOCALS| ii He_n-1 He_n 2x |
DUP 0= IF DROP 1e EXIT ENDIF
DUP 1 = IF DROP He_n EXIT ENDIF
He_n-1 He_n
1 DO FDUP 2x F*
FROT ii F2* F* F-
1e +TO ii
LOOP FNIP ;
FORTH> true to verbose? test-hermite
RMS error: 1.63418795292645e-014
Computed:
2.875000e+000 9.062500e+000 -2.753125e+001 -7.685938e+001 3.688047e+002 8.916289e+002 ...
0.000000e+000 1.200000e+001 0.000000e+000 -1.200000e+002 0.000000e+000 1.680000e+003 ...
-2.875000e+000 9.062500e+000 2.753125e+001 -7.685938e+001
Comparison values:
2.875000e+000 9.062500e+000 -2.753125e+001 -7.685937e+001 3.688047e+002 8.916289e+002 ...
0.000000e+000 1.200000e+001 0.000000e+000 -1.200000e+002 0.000000e+000 1.680000e+003 ...
-2.875000e+000 9.062500e+000 2.753125e+001 -7.685937e+001 ok
-marcel
Krishna Myneni
> 0.000000e+000 1.200000e+001 0.000000e+000 -1.200000e+002 0.000000e+000 1.680000e+003 ...
> -2.875000e+000 9.062500e+000 2.753125e+001 -7.685937e+001 ok
>
> -marcel
>
Ok.
I'll update the reference values in my test code for the generalized Laguerre
polynomials with the higher precision results you gave earlier from Mathematica
v3.0. There *is* a comment in my test code about the Bessel polynomials being
reverse Bessel polynomials, added when I discovered the discrepancy myself. But,
apparently the reverse Bessel polynomials are much more useful in engineering,
which is probably why Skip Carter implemented those.
Krishna
Marcel Hendrix
[..]
> There *is* a comment in my test code about the Bessel polynomials being
> reverse Bessel polynomials, added when I discovered the discrepancy myself. But,
> apparently the reverse Bessel polynomials are much more useful in engineering,
> which is probably why Skip Carter implemented those.
Whatever, here is the normal version.
\ Be_n Evaluates the nth order Bessel Polynomial,
\ Be_n(x) = \sum_k=0^n d_k (x/2)^k, d_k = (n+k)! / ((n-k)! k!)
\ Ref. http://en.wikipedia.org/wiki/Bessel_polynomials
-- nth order Bessel Polynomial
-- Note: Be_n+1 = (2n-1)*x*Be_n + Be_n-1; Be_0 = 1, Be_1 = 1+x
: Be_n ( n -- ) ( F: x -- z )
1e FOVER F1+ ( x 1 1+x )
DUP 0= IF DROP FDROP FNIP EXIT ENDIF
DUP 1 = IF DROP FNIP FNIP EXIT ENDIF
3e -FROT ( x 2n-1 Be_n-1 Be_n )
1 DO FDUP 3 FPICK F* 4 FPICK F*
FROT F+
FROT 2e F+ -FROT
LOOP FNIP FNIP FNIP ;
\ Mathematica 3.0
\ ToString[TableForm[Table[FortranForm[N[
\ Sum[(n + k)! / (k! (n-k)!) ((-0.25+0.25 i+10^-20)/2)^k,{k,0,n}],20] ],{i,0,4},{n,3,6}]]]
DOUBLE DARRAY bessel-table{
bessel-table{ 5 4 }}FREAD
0.203125 0.08203125 0.0185546875 0.031005859375
1. 1. 1. 1.
3.671875 8.36328125 22.4892578125 70.208740234375
9.625 36.9375 175.84375 1004.078125
20.265625 111.33203125 771.7568359375 6478.325927734376
bessel-table{ #20 }resize
: test-bessel ( -- )
['] Be_n [IS] poly
bessel-table{ 7 test-polys
verbose? 0= ?EXIT
CR ." Computed: " res{ }print
CR ." Comparison values:" bessel-table{ 0 res{ CDIM 1- }print[] ;
\ FORTH> test-bessel
\ RMS error: 9.11411746365820e-014 ok
-marcel
Krishna Myneni
> m...@iae.nl (Marcel Hendrix) writes Re: automated FSL tests
> [..]
>> FORTH> test-cheby test-hermite test-laguerre test-legendre test-glaguerre test-rbessel
>> RMS error: 2.11796688046160e-017
>> RMS error: 1.63418795292645e-014
>> RMS error: 4.33185611774133e-017
>> RMS error: 3.29611662544038e-017
>> RMS error: 2.05078980308681e-016
>> RMS error: 5.39983240078416e-013 ok
>
>> Given the fact that rBessel has values spanning 10 .. 10.000, the absolute error of 5e-13
>> for rbessel is satisfying.
>
This is a good example of the reason why the floating point test code of David
Williams, which has been integrated into Anton Ertl's ttester.fs, has both
relative and absolute nearness tests.
Krishna
Krishna Myneni
> ...
> \ Be_n Evaluates the nth order Bessel Polynomial,
> \ Be_n(x) = \sum_k=0^n d_k (x/2)^k, d_k = (n+k)! / ((n-k)! k!)
> \ Ref. http://en.wikipedia.org/wiki/Bessel_polynomials
>
> -- nth order Bessel Polynomial
> -- Note: Be_n+1 = (2n-1)*x*Be_n + Be_n-1; Be_0 = 1, Be_1 = 1+x
> : Be_n ( n -- ) ( F: x -- z )
> 1e FOVER F1+ ( x 1 1+x )
> DUP 0= IF DROP FDROP FNIP EXIT ENDIF
> DUP 1 = IF DROP FNIP FNIP EXIT ENDIF
> 3e -FROT ( x 2n-1 Be_n-1 Be_n )
> 1 DO FDUP 3 FPICK F* 4 FPICK F*
> FROT F+
> FROT 2e F+ -FROT
> LOOP FNIP FNIP FNIP ;
Your choice of "Be_n" as the name for the above word conflicts with the existing
FSL polys.x. If you want to add the above function to polys.x, I'd recommend
that you use a different name.
Krishna
Krishna Myneni
> ... But,
> apparently the reverse Bessel polynomials are much more useful in engineering,
> which is probably why Skip Carter implemented those.
An interesting fact:
Low-pass filters, such as those used in high-fidelity (non-distorting) audio
amplifier systems, have transfer functions described by the reverse Bessel
polynomials. More at
http://en.wikipedia.org/wiki/Bessel_filter
Krishna
Jerry Avins
Bessel filters have poorer transition regions than other types such as
Butterworth and Chebychev. They are very good for minimum phase shift
within their passbands. Symmetric FIR filters have no phase distortion
no matter how sharp their transitions are made. With all filters, sharp
transitions in frequency lead to ringing in time.
Jerry
--
Engineering is the art of making what you want from things you can get.
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
Marcel Hendrix
[..]
> Bessel filters have poorer transition regions than other types such as
> Butterworth and Chebychev. They are very good for minimum phase shift
> within their passbands.
Don't you mean 'linear phase shift with frequency' (i.e. a constant
time-delay)?
> Symmetric FIR filters have no phase distortion
> no matter how sharp their transitions are made. With all filters, sharp
> transitions in frequency lead to ringing in time.
I guess recording a signal and running it in reverse through a filter,
than reverse the result of that in time and run it through the same filter
should undo all phase errors (while doubling the amplitude damping).
Do tricks exist to remove the amplitude damping while worsening the phase
shift?
-marcel
John Doty
> Jerry Avins <j...@ieee.org> wrote Re: automated FSL tests
> [..]
>> Bessel filters have poorer transition regions than other types such as
>> Butterworth and Chebychev. They are very good for minimum phase shift
>> within their passbands.
>
> Don't you mean 'linear phase shift with frequency' (i.e. a constant
> time-delay)?
Sure, but sometimes you conceptually remove the group delay and then
consider the filter to have (nearly) constant phase.
>
>> Symmetric FIR filters have no phase distortion
>> no matter how sharp their transitions are made. With all filters, sharp
>> transitions in frequency lead to ringing in time.
>
> I guess recording a signal and running it in reverse through a filter,
> than reverse the result of that in time and run it through the same filter
> should undo all phase errors (while doubling the amplitude damping).
>
> Do tricks exist to remove the amplitude damping while worsening the phase
> shift?
http://en.wikipedia.org/wiki/All-pass_filter
--
John Doty, Noqsi Aerospace, Ltd.
http://www.noqsi.com/
--
Specialization is for robots.
Jerry Avins
I cam make a pretty good guess about what "amplitude damping" means, but
it's not a term I've seen before. Running a signal through a filter
twice doubles the filter's amplitude effect in dB, which is the same as
squaring it when expressed on a linear scale. As for phase shift, the
second pass doubles it if the signal runs forward and cancels it if the
signal runs backward. The backward run can't be started until the first
pass is finished. The technique can be used with recordings, but not
with a real-time signal such as a PA system. It's more than a matter of
latency.
Filters that shift phase without changing amplitude are called
"all-pass" filters. They are not very difficult to design.
Krishna Myneni
> Marcel Hendrix wrote:
>>
>> I guess recording a signal and running it in reverse through a filter,
>> than reverse the result of that in time and run it through the same
> ... As for phase shift, the
> second pass doubles it if the signal runs forward and cancels it if the
> signal runs backward. The backward run can't be started until the first
> pass is finished. The technique can be used with recordings, but not
> with a real-time signal such as a PA system. It's more than a matter of
> latency.
> ...
Time-reversing a signal is a fairly complex matter, but it doesn't necessarily
require recording the signal and playing it backwards, in the traditional sense.
Time-reversing an optical signal is called optical phase conjugation, and can be
accomplished with a technique called degenerate four-wave mixing. An example of
its use is in correcting wavefront distortion. When an optical wavefront passes
through the atmosphere, it becomes distorted due to phase fluctuations. Here,
the atmosphere, or medium through which the wave propagates, is the "filter".
The optical wave can be reflected from a "phase-conjugate mirror", where it
becomes time-reversed and is reflected back through the same distorting medium,
removing the phase (wavefront) errors [1]. Optical phase conjugation can also be
used to remove dispersive effects in fibers, which cause pulse broadening.
I'm not sure if there is an electronic analog to a phase-conjugate mirror, but
it would have to be a nonlinear circuit.
Krishna
[1] See the short description at
http://en.wikipedia.org/wiki/Nonlinear_optics
or for a fuller treatment, the book "Nonlinear Optics", by R. W. Boyd.
John Doty
> Jerry Avins wrote:
>> Marcel Hendrix wrote:
>>> I guess recording a signal and running it in reverse through a filter,
>>> than reverse the result of that in time and run it through the same
>>> filter should undo all phase errors ...
>
>> ... As for phase shift, the
>> second pass doubles it if the signal runs forward and cancels it if the
>> signal runs backward. The backward run can't be started until the first
>> pass is finished. The technique can be used with recordings, but not
>> with a real-time signal such as a PA system. It's more than a matter of
>> latency.
>> ...
>
>
> Time-reversing a signal is a fairly complex matter, but it doesn't necessarily
> require recording the signal and playing it backwards, in the traditional sense.
> Time-reversing an optical signal is called optical phase conjugation, and can be
> accomplished with a technique called degenerate four-wave mixing. An example of
> its use is in correcting wavefront distortion. When an optical wavefront passes
> through the atmosphere, it becomes distorted due to phase fluctuations. Here,
> the atmosphere, or medium through which the wave propagates, is the "filter".
> The optical wave can be reflected from a "phase-conjugate mirror", where it
> becomes time-reversed and is reflected back through the same distorting medium,
> removing the phase (wavefront) errors [1]. Optical phase conjugation can also be
> used to remove dispersive effects in fibers, which cause pulse broadening.
Note that phase conjugation doesn't really time reverse a signal,
because it also "flips" frequencies about the center frequency. For this
reason the double pass technique you describe also flattens the
amplitude response to first order.
>
> I'm not sure if there is an electronic analog to a phase-conjugate mirror, but
> it would have to be a nonlinear circuit.
Seems to me you could do it with a hybrid transformer and a balanced
mixer with local oscillator ("pump") at twice the center frequency. But
I've never seen an application of this. I'll put it in the backbrain
junkyard, to be triumphantly produced when some application needs it
(another opportunity for a "Where did Doty get *that* idea?" moment). We
can file for the patent jointly...
On the other hand, I've actually done a paper design for a satellite
communication system where the satellite would steer the downlink beam
by measuring the uplink carrier phases on a phased antenna array and
conjugating them to get downlink transmitter phases, thus automatically
steering the downlink beam to the ground station. But there, the phase
conjugation was to be handled digitally, and it would be more
complicated than just a mirror because the downlink data would be
different from the uplink data.
Nonlinear optics generally has electronic analogies: after all,
electronics is just a special case of optics, where the wavelengths are
long, and you're usually working with one-dimensional waveguides whose
modes don't correspond neatly to "wave vector" and "polarization". But
the same basic electromagnetic and mathematical principles apply.
Indeed, most nonlinear optics phenomena were first exploited in radio.
>
> Krishna
>
>
> [1] See the short description at
>
> http://en.wikipedia.org/wiki/Nonlinear_optics
>
> or for a fuller treatment, the book "Nonlinear Optics", by R. W. Boyd.
>
Krishna Myneni
> Krishna Myneni wrote:
>> ... The optical wave can be reflected from a "phase-conjugate mirror",
>> where it
>> becomes time-reversed and is reflected back through the same
>> distorting medium,
>> removing the phase (wavefront) errors [1]. Optical phase conjugation
>> can also be
>> used to remove dispersive effects in fibers, which cause pulse
>> broadening.
>
> Note that phase conjugation doesn't really time reverse a signal,
> because it also "flips" frequencies about the center frequency. For this
> reason the double pass technique you describe also flattens the
> amplitude response to first order.
>
I don't understand the above statement. How do you distinguish a sign change in
the frequency from a sign change in time?
In the simple case, for a monochromatic wave
E(r, t) = E_0(r)*exp( i*omega_0*t )
where r is a vector, and, assuming a real amplitude E_0(r), the conjugate wave is
E_c(r, t) = E_0(r)*exp( -i*omega_0*t )
which is equivalent to the time-reversed wave:
E_c(r, t) = E(r, -t)
More generally, if the field is of the form
E(r, t) = E_0(r, t)*exp( i*phi(r, t) )
where the amplitude and phase, E_0(r, t) and phi(r, t), are real functions, then
( E_0(r, t) )^* = E_0(r, t)
( phi(r, t) )^* = phi(r, t)
and the conjugate field becomes
E_c(r, t) = E_0(r, t)*exp( -i*phi(r, t) )
Now, if the amplitude is even in time, and the phase term is odd in time, i.e.
E_0(r, -t) = E_0(r, t)
phi(r, -t) = -phi(r, t)
then, once again we have
E_c(r, t) = E(r, -t)
which is a time-reversed wave.
>>
>> I'm not sure if there is an electronic analog to a phase-conjugate
>> mirror, but
>> it would have to be a nonlinear circuit.
>
> Seems to me you could do it with a hybrid transformer and a balanced
> mixer with local oscillator ("pump") at twice the center frequency. But
> I've never seen an application of this. I'll put it in the backbrain
> junkyard, to be triumphantly produced when some application needs it
> (another opportunity for a "Where did Doty get *that* idea?" moment). We
> can file for the patent jointly...
>
Sure, ... where do I sign?
Krishna
John Doty
> John Doty wrote:
>> Krishna Myneni wrote:
>
>>> ... The optical wave can be reflected from a "phase-conjugate mirror",
>>> where it
>>> becomes time-reversed and is reflected back through the same
>>> distorting medium,
>>> removing the phase (wavefront) errors [1]. Optical phase conjugation
>>> can also be
>>> used to remove dispersive effects in fibers, which cause pulse
>>> broadening.
>> Note that phase conjugation doesn't really time reverse a signal,
>> because it also "flips" frequencies about the center frequency. For this
>> reason the double pass technique you describe also flattens the
>> amplitude response to first order.
>>
>
> I don't understand the above statement. How do you distinguish a sign change in
> the frequency from a sign change in time?
You can't. Nor, for a real valued signal, can you distinguish phase
conjugation from a sign change in time. Optical phase conjugation does
*both* these things, so it does not time reverse the signal.
If you reflect a modulated optical signal off of a phase conjugating
mirror, you don't get tomorrow's data out today, and today's data out
tomorrow. The output signal is simply delayed by the group delay through
the system.
The interesting thing about optical phase conjugation is, of course,
that the phase is conjugated and the frequency flipped about the
reference provided by the pump, not around zero frequency. Just
conjugating phase in the ordinary Fourier transform and negating
frequency is uninteresting because for a real signal that changes nothing.
>>> I'm not sure if there is an electronic analog to a phase-conjugate
>>> mirror, but
>>> it would have to be a nonlinear circuit.
>> Seems to me you could do it with a hybrid transformer and a balanced
>> mixer with local oscillator ("pump") at twice the center frequency.
And actually, this is like the classic AM superheterodyne radio. With
the local oscillator on the high side (the usual practice), the
intermediate frequency signal is flipped in frequency structure and
conjugated in phase relative to the incoming signal. Thus, the music
doesn't come out backwards!
> But
>> I've never seen an application of this. I'll put it in the backbrain
>> junkyard, to be triumphantly produced when some application needs it
>> (another opportunity for a "Where did Doty get *that* idea?" moment). We
>> can file for the patent jointly...
>>
>
> Sure, ... where do I sign?
Need to figure out what problem it solves, first.
Krishna Myneni
I do have an error in how I presented the field change in optical phase
conjugation (should have checked Boyd's book first; see below). However,
I'm somewhat puzzled by your statement that an ideal phase conjugate
mirror (PCM) does not reverse the wavefront, i.e. read the wavefront out
in reverse time order. Of course we are not able to read tomorrow's
paper today, and there must be an inherent delay in the PCM. But the
output signal is not simply a delayed version of the input signal,
otherwise it would just be an ordinary mirror with a delay.
Following Boyd,
The forward propagating wave is given by
E(r, t) = E_0(r,t)*exp( i*omega*t ) + c.c.
where c.c. is the complex conjugate of the previous term. The phase
conjugate mirror does not conjugate the entire expression (as I
expressed in my previous equations), but produces a field in the reverse
direction which is described by
E_c(r, t) = conj(E_0(r,t))*exp( i*omega*t ) + c.c.
and, it follows that
E_c(r, t) = E(r, -t)
This is, of course, an oversimplification of the phenomenon. But,
formally, it shows how the reflected waveform can be considered to be a
time-reversed version of the incident waveform. The delay is not
manifest in the above equations.
Perhaps we are trying to say the same thing, but the language has
ambiguities, which is why a topic like this can only be properly
discussed with mathematics.
Krishna
John Doty
Yes. You must take frequency flipping into account to understand what
the delay is. And then you'll see why you restore the input. Consider
the following one dimensional problem, a concrete version of your original:
I have a 1 s pulse with a rectangular envelope consisting of 1000 cycles
of a 1 kHz tone. Its Fourier transform is a sinc function with first
zeros at 999 and 1001 kHz. The bulk of the power is between 999.5 and
1000.5 Hz. Now suppose I run it through a passive reciprocal filter with
a group delay of 10 s at 1 kHz, increasing by 1 s/Hz. This will distort
the pulse into a chirp that's no longer rectangular and ~2 s long,
because the group delays at 999.5 and 1000.5 Hz differ by a second.
Can we recover the original with phase conjugation? Yes! Feed the output
of the filter to a phase conjugate reflector operating at 1 kHz. The
phase conjugate reflector is memoryless, so its group delay must be
zero, but it flips frequencies about 1 kHz. So, the slow-moving 1000.5
kHz components become faster 999.5 kHz components, and vice-versa. Once
through the filter the second time, the total delay is equalized (at 20
s) and the pulse is restored.
But note that in a three dimensional problem, a phase conjugate mirror
does not generally restore temporal structure: the reflection of a pulse
from the far side of the mirror is going to take longer to return to
sender than the reflection from the near side, even if its constituent
waves come back parallel to the ones sent.
Krishna Myneni
> Yes. You must take frequency flipping into account to understand what
> the delay is. And then you'll see why you restore the input. Consider
> the following one dimensional problem, a concrete version of your original:
>
> I have a 1 s pulse with a rectangular envelope consisting of 1000 cycles
> of a 1 kHz tone. Its Fourier transform is a sinc function with first
> zeros at 999 and 1001 kHz. The bulk of the power is between 999.5 and
> 1000.5 Hz. Now suppose I run it through a passive reciprocal filter with
> a group delay of 10 s at 1 kHz, increasing by 1 s/Hz. This will distort
> the pulse into a chirp that's no longer rectangular and ~2 s long,
> because the group delays at 999.5 and 1000.5 Hz differ by a second.
>
Fine. You are describing a particular incident pulse and a filter transfer
function. The incident pulse can be written as
E_s(r,t) = R(t - z/c)*exp( i*2*pi*1000*t )
where R is a rectangle function which provides the 1000 cycle envelope, and the
exponential term represents the carrier wave. The form of the filter function is
not relevant to this discussion, at least for an ideal phase conjugation process.
> Can we recover the original with phase conjugation? Yes! Feed the output
> of the filter to a phase conjugate reflector operating at 1 kHz. The
> phase conjugate reflector is memoryless, so its group delay must be
> zero, but it flips frequencies about 1 kHz. So, the slow-moving 1000.5
> kHz components become faster 999.5 kHz components, and vice-versa. Once
> through the filter the second time, the total delay is equalized (at 20
> s) and the pulse is restored.
> ...
We seem to agree on the phenomenon here, i.e. that the original pulse can be
restored after reflection from an ideal PCM, and a second pass through the
filter. However, you seem to take exception to calling the ideal PCM process a
"time-reversal". Formally, it has been shown that if we start with a *real*
wave, and conjugate the relevant portions of the wave, then the conjugation
procedure is equivalent to replacing t with -t in the wave function:
incident wave: E(r, t) = E_0(r,t)*exp( i*omega0*t ) + c.c.
PCM reflected wave: E_c(r, t) = conj(E_0(r,t))*exp( i*omega0*t ) + c.c.
Compare E_c with E; it follows that E_c(r, t) = E(r, -t)
Note that E(r, t) is the pulse after it goes through the filter, and before it
is incident on the PCM, i.e., in the frequency domain,
E(r, omega) = H(omega)*E_s(r, omega)
where H(omega) is the filter transfer function.
I believe the conj(E_0(r, t)), and corresponding term in the c.c., is equivalent
to what you mean by "frequency flipping" about the center frequency (omega0). It
is taken into account by virtue of the conjugation. If one accepts the
conjugation procedure described above as representing an *ideal* PCM, then it is
exactly equivalent to replacing t with -t, in the original wave.
Krishna
John Doty
This last step is wrong, and since you don't display your reasoning, I
can't debug it. But suppose E_0(r,t) is real and that also E(r,t) !=
E(r, -t). This time asymmetry of the input is a reasonable assumption,
true except in special cases. Then E_c(r,t) == E(r,t) != E(r,-t).
A better factoring makes things clearer. Let
a(t) be the amplitude. Real.
p(t) be the phase factor. Complex, magnitude 1.
c(t) be the "carrier" wave, exp(i*omega0*t)
Then E(t) == Re(a(t)*p(t)*c(t))
As a(t) is real, we can pull it out: E(t) == a(t)*Re(p(t)*c(t))
Since Re(p(t)*c(t)) <= 1, while a(t) can be any real value, phase
changes cannot time reverse a signal in general.
A phase-conjugate reflector replaces p(t) with conj(p(t)). The
interesting properties emerge from that. But you cannot make a time
machine this way.
Krishna Myneni
> Krishna Myneni wrote:
>> *real*
>> wave, and conjugate the relevant portions of the wave, then the
>> conjugation
>> procedure is equivalent to replacing t with -t in the wave function:
>>
>> incident wave: E(r, t) = E_0(r,t)*exp( i*omega0*t ) + c.c.
>>
>> PCM reflected wave: E_c(r, t) = conj(E_0(r,t))*exp( i*omega0*t ) +
>> c.c.
>>
>> Compare E_c with E; it follows that E_c(r, t) = E(r, -t)
>
> This last step is wrong, and since you don't display your reasoning, I
> can't debug it. But suppose E_0(r,t) is real and that also E(r,t) !=
> E(r, -t). This time asymmetry of the input is a reasonable assumption,
> true except in special cases. Then E_c(r,t) == E(r,t) != E(r,-t).
>
Your example is certainly correct. We must require the real envelope of the
pulse to be symmetric in time for E_c(r, t) = E(r, -t). In my explanation, I
had written the input field in too general a form. But, it should be noted that
it is not necessary for E(r, t) = E(r, -t) in order for E_c(r, t) = E(r, -t),
i.e. we do not have to require E_0(r, t) to be real!
> A better factoring makes things clearer. Let
>
> a(t) be the amplitude. Real.
> p(t) be the phase factor. Complex, magnitude 1.
> c(t) be the "carrier" wave, exp(i*omega0*t)
>
> Then E(t) == Re(a(t)*p(t)*c(t))
>
> As a(t) is real, we can pull it out: E(t) == a(t)*Re(p(t)*c(t))
>
> Since Re(p(t)*c(t)) <= 1, while a(t) can be any real value, phase
> changes cannot time reverse a signal in general.
>
> A phase-conjugate reflector replaces p(t) with conj(p(t)). The
> interesting properties emerge from that. ...
I agree that an *arbitrary* pulse cannot be time reversed. The time reversal
argument only applies if the *real envelope* is time symmetric. I did not make
this distinction earlier! In general, we cannot admit any arbitray complex
function E_0(r, t). The restriction on E_0(r, t), that it's amplitude be time
symmetric, can be seen with the aid of your factoring. I have written it in
slightly modified form, but the equivalence is easy to see.
E(r, t) = E_0(r, t)*exp( i*omega0*t ) + c.c. (1)
Now, let us factor the complex term E_0(r, t) as follows:
E_0(r, t) = A(r, t)*exp( i*phi(r, t) ) (2)
where, without loss of generality, we may take A(r) and phi(r, t) to be real
functions. We then have the following equivalence with your notation:
a(r, t) = A(r, t)
p(r, t) = exp( i*phi(r, t) ) (3)
where we have generalized the functions a() and p() to be functions of both
position and time. The phase conjugated E(r, t) is given by:
E_c(r, t) = A(r, t)*{ conj( p(r, t) )*exp( i*omega0*t ) + c.c. } (4)
Now, when is E_c(r, t) = E(r, -t)? Since the original field, written in the
above notation, is
E(r, t) = A(r, t)*{ p(r, t)*exp( i*omega0*t ) + c.c. } (5)
the equivalency of the conjugate field with the time reversed field is only true
when
A(r, t) = A(r, -t) (6)
Thus, only the real amplitude of the wave needs to be time symmetric. The input
pulse may still be asymmetric, i.e. E(r, t) != E(r, -t), e.g. a pulse with a
symmetric Gaussian envelope, but with chirp.
I think we are now in agreement John.
Krishna
John Doty
Yes, but there's another condition. You also must have phi(r,t) ==
-phi(r,-t). Otherwise, conj(p(r,t)) != p(r,-t).
>
> Thus, only the real amplitude of the wave needs to be time symmetric. The input
> pulse may still be asymmetric, i.e. E(r, t) != E(r, -t), e.g. a pulse with a
> symmetric Gaussian envelope, but with chirp.
A chirp centered on omega0 works, yes. And that's the famous dispersion
correction case. But you still can't recover tomorrow's FM broadcast today!
>
>
> I think we are now in agreement John.
Well, maybe now ;-)
Krishna Myneni
There is a condition on phi(r, t), which I failed to observe. But the condition
does not appear to be the one you state. Note that first terms in braces for
eqns 4 and 5 are multiplied by exp( i*omega0*t ).
I believe the condition required for E_c(r, t) = E(r, -t) is
phi(r, -t) = phi(r, t) (7)
Then,
E(r, -t) = A(r, -t)*{ p(r, -t)*exp( -i*omega0*t ) + c.c. }
= A(r, t)*{ p(r, t)*exp( -i*omega0*t ) + c.c. } (8)
The above expression now equates to E_c(r, t) in eq 4 (note, the first term in
braces in eq. 8 corresponds to the c.c. term in eq 4).
The new requirement leads to the complex amplitude having the property of time
reversibility:
E_0(r, t) = E_0(r, -t) (9)
Note that equation (9) still does not imply that the incident wave is time
reversible, i.e. E(r, t) != E(r, -t).
>>
>> Thus, only the real amplitude of the wave needs to be time symmetric.
>> The input
>> pulse may still be asymmetric, i.e. E(r, t) != E(r, -t), e.g. a pulse
>> with a
>> symmetric Gaussian envelope, but with chirp.
>
> A chirp centered on omega0 works, yes. And that's the famous dispersion
> correction case. But you still can't recover tomorrow's FM broadcast today!
>
>>
>>
>> I think we are now in agreement John.
>
> Well, maybe now ;-)
>
How about now?
Krishna
Krishna Myneni
for the Gforth collection of FSL modules at ccreweb.org:
ftp://ccreweb.org/software/gforth/fsl/
Please see the README file in the above directory.
The automated tests, using high-precision reference values, reveal that the
gamma function calculation only provides 6 or 7 significant digits. The original
module made no mention of the accuracy in its comments section, although the
source reference probably states the accuracy.
While the above accuracy may be adequate for many applications, there exist
relatively simple algorithms for computing the gamma function to a greater
number of significant digits. An alternate implementation may be a worthwhile
addition to the FSL.
Krishna Myneni
John Doty
Yes. And how many times have I explained that when the pulsar is
accelerating at a constant rate the phase residuals are parabolic? You
caught me in an elementary error, confusing phase and frequency.
Ohwa tagu siam!
>
> Then,
>
> E(r, -t) = A(r, -t)*{ p(r, -t)*exp( -i*omega0*t ) + c.c. }
>
> = A(r, t)*{ p(r, t)*exp( -i*omega0*t ) + c.c. } (8)
>
> The above expression now equates to E_c(r, t) in eq 4 (note, the first term in
> braces in eq. 8 corresponds to the c.c. term in eq 4).
>
> The new requirement leads to the complex amplitude having the property of time
> reversibility:
>
> E_0(r, t) = E_0(r, -t) (9)
>
> Note that equation (9) still does not imply that the incident wave is time
> reversible, i.e. E(r, t) != E(r, -t).
>
>
>
>>> Thus, only the real amplitude of the wave needs to be time symmetric.
>>> The input
>>> pulse may still be asymmetric, i.e. E(r, t) != E(r, -t), e.g. a pulse
>>> with a
>>> symmetric Gaussian envelope, but with chirp.
>> A chirp centered on omega0 works, yes. And that's the famous dispersion
>> correction case. But you still can't recover tomorrow's FM broadcast today!
>>
>>>
>>> I think we are now in agreement John.
>> Well, maybe now ;-)
>>
>
> How about now?
Yes.
Krishna Myneni
of FSL library modules. Over 30% of the modules now have test code revised to
perform automated testing. Recent additions, since the last update on 10/21/07, are:
logistic.fs
crc.fs
permcomb.fs
adaptint.fs
elip.fs
gauss.fs
sph_bes.fs
cubic.fs
factorl.fs
shanks.fs
pcylfun.fs
gaussj.fs
The revised test code varies widely in rigor, accuracy, and completeness among
the different modules, but it strives, at a minimum, to provide the same tests
as in the original modules. A single program, fsl-tester.fs, runs the tests for
all of the available modules. The demonstration modules may be found in
ftp://ccreweb.org/software/gforth/fsl/
Please view the README file to find out more about this effort. All of the
revised modules were tested under Gforth, but should require few modifications
for any other 32/64-bit ANS Forth systems which use a separate fp stack.
Corresponding versions for kForth (and other systems with a unified stack, and
having 1 FLOATS = 2 CELLS) may be found, as usual, on the kForth Programming
Examples Page.
Although I haven't explicitly asked for any help previously, it goes without
saying that anyone interested in contributing to this effort, either by working
on remaining modules or by extending/strengthening tests on already revised
modules, or other aspect is most welcome to participate, particularly original
authors. The current code probably makes a lot of implicit assumptions, and
little thought has been given to portability to systems which vary in floating
point precision, e.g. 16-bit systems. A general method of setting the fp
comparison tolerances, based on the system's floating point precision, and on
the module's expected accuracy, is sorely lacking.
Krishna Myneni