d/dx in forth

[Ahmed]

Hello,

In the thread:
https://www.novabbs.com/devel/article-flat.php?group=comp.lang.forth&id=16944&first=1&last=25#start

Paul Robin posted:

For example, sin (sine) is a function, whose input is a real and
whose
output is another real. And "derivative" (d/dx) is a function whose
input is a function, and whose output is another function. So you
get
that d/dx (sin) = cos, or some numerical approximation of same. In
Python or Haskell, this is easy, e.g.:


h = 0.0001
d_over_dx f = g where { g x = (f(x+h) - f(x)) / h }
dsin = d_over_dx sin
print (dsin 0.5) -- prints 0.8775585891507287
print (cos 0.5) -- prints 0.8775825618903728


Pretty close approximation. In Forth it is a bit messier because you
want a signature like


: d/dx ( xt -- xt ) .... ;


and you may have to do some advanced or non-portable things to create
new xt's on the fly like that. But there are various ways you can
get
equivalent functionality using OOP or whatever.


minforth replied:

Clumsy but portable:

0.0001e FVALUE h

DEFER ofunc ' FSIN IS ofunc
DEFER dfunc

: D/DX {: xt f: x -- y :}
x h f+ xt execute x xt execute f- h f/ ;

:NONAME ( x -- y ) ['] ofunc D/DX ; IS dfunc


0.5e dfunc f.
0.5e fcos f.

My question is:
Why not do it directly like this:

: d/dx ( xt -- )
create ,
does> @ dup ( f: x -- y) fdup 1e-4 f+ execute fswap 1e-4 f-
execute f- 0.5e4 f* ;

and the dervative of any function can be obtained directly like this:

' func() d/dx der_func()

for example: (in gforth)
' fsin d/dx der_fsin ok
0.5e fcos f. 0.877582561890373 ok
0.5e der_fsin f. 0.877582560427637 ok

another example:
: f() ( f: x -- y) fdup fdup f* f+ ; \ f(x) = x^2+x
: df() ( f: x -- y) 2e f* 1e f+ ; \ df(x) = 2*x + 1 for comparison
' f() d/dx der_f()

5e der_f() f. 10.9999999999744 ok
5e df() f. 11. ok

Best Regards.

[Anton Ertl]

Post by Ahmed
Hello,
https://www.novabbs.com/devel/article-flat.php?group=comp.lang.forth&id=16944&first=1&last=25#start
For example, sin (sine) is a function, whose input is a real and
whose
output is another real. And "derivative" (d/dx) is a function whose
input is a function, and whose output is another function. So you
get
that d/dx (sin) = cos, or some numerical approximation of same. In
h = 0.0001
d_over_dx f = g where { g x = (f(x+h) - f(x)) / h }
dsin = d_over_dx sin
print (dsin 0.5) -- prints 0.8775585891507287
print (cos 0.5) -- prints 0.8775825618903728
Pretty close approximation. In Forth it is a bit messier because you
want a signature like
: d/dx ( xt -- xt ) .... ;

Most probably he meant ( xt1 -- xt2 )

Post by Ahmed
and you may have to do some advanced or non-portable things to create
new xt's on the fly like that.
: d/dx ( xt -- )
create ,
execute f- 0.5e4 f* ;
' func() d/dx der_func()

or, for his example:

' fsin d/dx fdsin

And given that he gives a name to the result anyway, it's obviously
fine that your solution requires a name and produces a named word
rather than an unnamed word. Your solution uses a slightly different
(unbiased) formula, though.

Post by Ahmed
for example: (in gforth)
' fsin d/dx der_fsin ok
0.5e fcos f. 0.877582561890373 ok
0.5e der_fsin f. 0.877582560427637 ok

The smaller difference from the fcos result demonstrates the advantage
of the unbiased formula.

A solution that actually complies with the ( xt1 -- xt2 ) requirement
and is in standard Forth is as follows (using the biased formula in
the Python or Haskell program above):

0.0001e fvalue h

: d/dx1 ( r1 xt -- r2 )
fdup h f+ dup execute fswap execute f- h f/ ;

: d/dx ( xt1 -- xt2 )

Post by Ahmed
r :noname ( r1 -- r2 ) r> postpone literal postpone d/dx1 postpone ; ;

0.5e ' fsin d/dx execute f. \ prints 0.877558589150729

Maybe he considers the use of DOES> or POSTPONE to be "advanced
things". The disadvantage of the ( xt1 -- xt2 ) variant is that it
may seduce the user to produce multiple instances of the same function
(possibly one instance for every iteration of some loop), with every
instance costing memory that is unlikely to be reclaimed before the
end of the session. With a version that produces a named word that is
far less likely.

The advantage of the unnamed variant is that we can do things like

defer fddsin
' fsin d/dx d/dx is fddsin

without naming the intermediate function. Looking at the result, this
seems to work ok:

0.5e fddsin f. \ prints -0.479513295736922
0.5e fsin fnegate f. \ prints -0.479425538604203

- anton

[Ahmed]

Got it.
Thanks for the clarifications.

[Anton Ertl]

Post by Anton Ertl

Post by Ahmed
Hello,
https://www.novabbs.com/devel/article-flat.php?group=comp.lang.forth&id=16944&first=1&last=25#start
For example, sin (sine) is a function, whose input is a real and
whose
output is another real. And "derivative" (d/dx) is a function whose
input is a function, and whose output is another function. So you
get
that d/dx (sin) = cos, or some numerical approximation of same. In
h = 0.0001
d_over_dx f = g where { g x = (f(x+h) - f(x)) / h }
dsin = d_over_dx sin
print (dsin 0.5) -- prints 0.8775585891507287
print (cos 0.5) -- prints 0.8775825618903728
: d/dx ( xt -- )
create ,
execute f- 0.5e4 f* ;

' fsin d/dx fdsin
0.0001e fvalue h
: d/dx1 ( r1 xt -- r2 )
fdup h f+ dup execute fswap execute f- h f/ ;
: d/dx ( xt1 -- xt2 )

Post by Ahmed
r :noname ( r1 -- r2 ) r> postpone literal postpone d/dx1 postpone ; ;

0.5e ' fsin d/dx execute f. \ prints 0.877558589150729

Another variant is to compile xt1 into xt2; here a standard version
(with Paul Rubin's formula):

0.0001e fvalue h

: d/dx ( xt1 -- xt2 )

Post by Anton Ertl
r :noname ( r1 -- r2 )

postpone fdup postpone h postpone f+ r@ compile,
postpone fswap r> compile, postpone f- postpone h postpone f/
postpone ; ;

0.5e ' fsin d/dx execute f. \ prints "0.877558589150729"

If I do

0.5e ' fsin d/dx disasm/f

on VFX64, I get:

( 0050A358 D9C0 ) FLD ST
( 0050A35A DB2DF0FEFFFF ) FLD TBYTE FFFFFEF0 [RIP] @0050A250
( 0050A360 DEC1 ) FADDP ST(1), ST
( 0050A362 E8F9CBFFFF ) CALL 00506F60 FSIN
( 0050A367 D9C9 ) FXCH ST(1)
( 0050A369 E8F2CBFFFF ) CALL 00506F60 FSIN
( 0050A36E DEE9 ) FSUBP ST(1), ST
( 0050A370 DB2DDAFEFFFF ) FLD TBYTE FFFFFEDA [RIP] @0050A250
( 0050A376 DEF9 ) FDIVP ST(1), ST
( 0050A378 C3 ) RET/NEXT
( 33 bytes, 10 instructions )

Now on to the less portable stuff. The many POSTPONE's in the code
above hamper readability. So let's use ]] ... [[ to get rid of them:

: d/dx ( xt1 -- xt2 )

Post by Anton Ertl
r :noname ( r1 -- r2 ) ]]

fdup h f+ [[ r@ compile, ]] fswap [[ r> compile, ]] f- h f/ ;
[[ ;

This example may also make it clear why the brackets are oriented the
way they are. Compiles on Gforth and VFX64 (and to the same code as
above on VFX64).

Gforth offers an additional nicety. The following code generates the
same code as the previous D/DX.

: d/dx {: xt: xt1 -- xt2 :}
:noname ( r1 -- r2 ) ]] fdup h f+ xt1 fswap xt1 f- h f/ ; [[ ;

The XT: results in XT1 being a defer-flavoured local, so mentioning
the name executes the xt, and POSTPONEing it (e.g., within ]]..[[)
COMPILE,s it.

Another variant is to use Gforth's flat closures:

: d/dx ( xt1 -- xt2 )
[{: xt: xt1 :}d fdup h f+ xt1 fswap xt1 f- h f/ ;] ;

The functionality is the same as my previous D/DX implementations, and
the source code looks similar to the previous one, but the
implementation is quite different:

' fsin d/dx

uses 3 cells in the dictionary (and none in the native code) with the
flat-closure version, while the three definitions before that (which
use POSTPONE in some form or other) take 16 cells of dictionary and
148 bytes of native code. Generating a closure is also much faster.
However, the POSTPONE-based variants will run faster.

Stack purists will dislike the local in the closure. Gforth also
provides a pure-stack variant, and using it here will look as follows:

: d/dx ( xt1 -- xt2 )
[n:d fdup h f+ dup execute fswap execute f- h f/ ;] ;

This variant makes the difference to the POSTPONE-based approaches
more obvious by calling EXECUTE explicitly.

Ahmed Melahis CREATE...DOES> variant can be seen as the classical
Forth syntax for the pure-stack flat-closure approach; he uses a
different formula and inlines H, which obscures this a bit, but using
the same formula and H, a CREATE-DOES variant looks as follows:

: d/dx-named ( xt1 "name" -- xt2 )
create ,
does> ( r1 -- r2 )
@ fdup h f+ dup execute fswap execute f- h f/ ;

where the code from FDUP to F/ is the same as in the pure-stack
closure version above. If you use D/DX-NAMED with Gforth's NONAME,
calling

' fsin noname d/dx-named

consumes three cells, like the flat-closure variants.

Some may conclude that this shows that flat closures are unnecessary
(and some may extend this to other newer features like defer-flavoured
locals and ]]...[[). However, note that flat closures can not just be
allocated in the dictionary, but also on the locals stack or on the
heap, unlike CREATE..DOES>.

Of course, if we judge expressive power by Turing completeness, one
can reject all new features (or leave away many existing ones),
because they do not change Turing completeness, but there is a reason
why we implement more than just a minimal Turing-complete language, so
we obviously use other criteria for judging the desirability of
features.

- anton

[Ruvim]

Post by Anton Ertl

Post by Ahmed
In Forth it is a bit messier because you
want a signature like
: d/dx ( xt -- xt ) .... ;

Most probably he meant ( xt1 -- xt2 )

Do you mean that ( xt -- xt ) is an incorrect stack diagram for this case?

I always thought that such use is correct. Because formally "xt"
(without an index) only denotes the data type, not a data value.

And "( xt -- xt )" does not mean that the data values before execution
and after execution are the same, since they may still be different.

In a stack diagram, the same data type symbols with the same indices
usually denote the same data value. Otherwise, if data values are not
mentioned (like "0", "-1", etc), it is unknown from the stack diagram
where the values are the same. Since different symbols in the does not
mean different data values. For example "abs ( n -- +n )" — in some
cases the data values "before" and "after" are the same.


--
Ruvim

[dxf]

Post by Ruvim

Post by Anton Ertl

  In Forth it is a bit messier because you
  want a signature like
   : d/dx ( xt -- xt ) .... ;

Most probably he meant ( xt1 -- xt2 )

Do you mean that ( xt -- xt ) is an incorrect stack diagram for this case?
I always thought that such use is correct. Because formally "xt" (without an index) only denotes the data type, not a data value.

"xt" is a number on the stack no different from "u" "n" "char" "flag" other
than in what it denotes.

Why differentiate the others when they differ in value but not "xt" ?

[Ruvim]

Post by dxf

Post by Ruvim

Post by Anton Ertl

  In Forth it is a bit messier because you
  want a signature like
   : d/dx ( xt -- xt ) .... ;

Most probably he meant ( xt1 -- xt2 )

Do you mean that ( xt -- xt ) is an incorrect stack diagram for this case?
I always thought that such use is correct. Because formally "xt" (without an index) only denotes the data type, not a data value.

"xt" is a number on the stack no different from "u" "n" "char" "flag" other
than in what it denotes.
Why differentiate the others when they differ in value but not "xt" ?

I don't mean that "xt" is special among the data types. It's just an
example. My point is about the stack notation in general (as a system of
symbols) [1].

In the stack notation, the diagrams "( foo -- foo )" and "( foo1 -- foo2
)", where "foo" is a data type symbol [2], are semantically equivalent.



[1] Forth-2012, 2.2.2 Stack notation
<https://forth-standard.org/standard/notation#subsection.2.2.2>
[2] Forth-2012, 3.1 Data types
<https://forth-standard.org/standard/usage#usage:data>


--
Ruvim

[Anton Ertl]

Post by Ruvim

Post by Anton Ertl

Post by Ahmed
In Forth it is a bit messier because you
want a signature like
: d/dx ( xt -- xt ) .... ;

Most probably he meant ( xt1 -- xt2 )

Do you mean that ( xt -- xt ) is an incorrect stack diagram for this case?

To me (and to vmgen) ( xt -- xt ) describes a word that expects an xt
on the data stack but does not change the data stack at all.

If two items in the stack effect comment have the same name (such as
in this case), they have the same value. That's also the usage in the
standard:

| 6.1.1290 DUP dupe CORE
| ( x -- x x )
|
| R@ r-fetch CORE
| ...
| Execution:
| ( -- x ) ( R: x -- x )

Post by Ruvim
I always thought that such use is correct. Because formally "xt"
(without an index) only denotes the data type, not a data value.

Says who?

My usage and the usage in the standard is that if there are two
different items with the same type, they get two different names,
e.g., by adding a digit. And if there is only one item of that type
(as in R@), there is no need to distinguish it from itself by adding a
digit.

Post by Ruvim
In a stack diagram, the same data type symbols with the same indices
usually denote the same data value. Otherwise, if data values are not
mentioned (like "0", "-1", etc), it is unknown from the stack diagram
where the values are the same. Since different symbols in the does not
mean different data values.

No, but the same name means that they are the same value.

Post by Ruvim
For example "abs ( n -- +n )" — in some
cases the data values "before" and "after" are the same.

The standard has:

| ABS abs CORE
| ( n -- u )

So we don't know whether the standard considers "+n" to be the same
name or a different name from "n".

- anton

[Ruvim]

Post by Anton Ertl

Post by Ruvim

Post by Anton Ertl

Post by Ahmed
In Forth it is a bit messier because you
want a signature like
: d/dx ( xt -- xt ) .... ;

Most probably he meant ( xt1 -- xt2 )

Do you mean that ( xt -- xt ) is an incorrect stack diagram for this case?

To me (and to vmgen) ( xt -- xt ) describes a word that expects an xt
on the data stack but does not change the data stack at all.

Well, vmgen may introduce it's own stack notation.

Post by Anton Ertl
If two items in the stack effect comment have the same name (such as
in this case), they have the same value. That's also the usage in the
| 6.1.1290 DUP dupe CORE
| ( x -- x x )

But the diagram "( x -- 2*x )" is equivalent to it, isn't it?



[...]

Post by Anton Ertl

Post by Ruvim
I always thought that such use is correct. Because formally "xt"
(without an index) only denotes the data type, not a data value.

Says who?

2.2.2 Stack notation says:

| /before/ represents the stack-parameter _data types_
| before execution of the definition and /after/ represents
| them after execution. The symbols used in /before/ and
| /after/ are shown in table 3.1.

The symbols in table 3.1 denote data types, not data objects (values).
In turn, a data type identifies a *set* of data objects.

Thus, a data type symbol in a stack diagram denotes a *set* of data
objects to which the corresponding stack-parameter belongs. (And this
correspondence is established by the position in the /before/ or /after/
tuple and in the stack).

Additional limitations on a stack-parameter can be *only* given in the
corresponding text description (which is defined in 2.2.4.2 Glossary
semantic description).

In a text description, a particular stack-parameter is referred by its
position on the stack, or by its corresponding element from the stack
diagram. Indices of stack diagram elements are only formally needed to
avoid ambiguity in such references. Any other meaning for indices is
only *informal*.


An example of referring stack parameters by position is the glossary
entry 6.2.2300 "TUCK":
| ( x1 x2 -- x2 x1 x2 )
| Copy the first (top) stack item below the second stack item.

NB: the requirement that the certain stack parameters are identical is
given in the text description.

The stack diagram ( 2*x -- 3*x ) is also formally correct for TUCK.
And the diagram ( x x -- x x x ) is equivalent to it.



So a rule of thumb for the /standard/ stack notation is this. If you
want to give an *informal* hint that some stack parameters are
identical, use the same indices (with the same data type symbols, of
course) for the corresponding elements of the stack diagram. If you want
to give an informal hint that some stack parameters (in the same data
type) *can* be not identical, use the different indices. Data type
symbols without indices don't give any hint.





[...]

Post by Anton Ertl

Post by Ruvim
For example "abs ( n -- +n )" — in some
cases the data values "before" and "after" are the same.

| ABS abs CORE
| ( n -- u )
So we don't know whether the standard considers "+n" to be the same
name or a different name from "n".

Sure, we know. "+n" and "n" are symbols for different data types, and +n
is a sub-type of n. See 3.1.1 Data-type relationships

+n ⇒ n ⇒ x;
+n ⇒ u ⇒ x;


In my example I use "+n" rather then "u" as this is also correct for
"ABS", and puts more restrictions on the stack-parameter.


--
Ruvim

[dxf]

Post by Ruvim
...

Post by Anton Ertl
If two items in the stack effect comment have the same name (such as
in this case), they have the same value.  That's also the usage in the
| 6.1.1290 DUP dupe CORE
| ( x -- x x )

But the diagram "( x -- 2*x )" is equivalent to it, isn't it?

Is it normative. AFAIK * is used for special circumstances. You can argue
it's legal but that's just language lawyering. If it confuses ppl what good
is it. It's setting a double standard. How is that new in Forth - you ask.
You got me there :)

[Ruvim]

Post by dxf

Post by Ruvim
...

Post by Anton Ertl
If two items in the stack effect comment have the same name (such as
in this case), they have the same value.  That's also the usage in the
| 6.1.1290 DUP dupe CORE
| ( x -- x x )

But the diagram "( x -- 2*x )" is equivalent to it, isn't it?

Is it normative. AFAIK * is used for special circumstances. You can argue
it's legal but that's just language lawyering. If it confuses ppl what good
is it. It's setting a double standard. How is that new in Forth - you ask.
You got me there :)

Some people tends to treat stack diagram items as arbitrary names for
actual stack parameters (i.e., names for the data objects in the
corresponding stack positions), for documentation purposes. They don't
pay attention to data types at all, or use the standard data type
symbols in the names just as informal hint.

A corollary of this view is that the same name in two stack diagram
items inevitably means that the corresponding stack parameters are
identical.

But it is not a standard stack notation. In the standard stack notation,
the names in a stack diagram are data types symbols with optional
indices/subscripts, not names for actual stack parameters or data objects.

And the same data type symbol in two stack diagram items means that the
corresponding stack parameters belong to the same data type. And nothing
else.

If you use the standard stack notation, and think this case is
confusing, don't use data type symbols without indices if the
corresponding stack parameters are always identical. Anyway, it's an
informal hint.

Moreover, if you want to use new data types in your stack diagrams, you
should document these data types, their data type symbols and data type
relationships. Otherwise it will be not a standard stack notation.


--
Ruvim

[minforth]

TEMP-CONVERSION ( F: celsius -- F: fahrenheit )
should be clear enough, right?

[Ruvim]

Post by minforth
TEMP-CONVERSION ( F: celsius -- F: fahrenheit )
should be clear enough, right?

Right, while you have not introduced subtypes of i*r.

So in the standard stack notation it would be:

TEMP-CONVERSION ( F: r.celsius -- r.fahrenheit )

1. A stack-id ("F" in this case) is used once.
2. The data type symbol "r" is used (it's defined in 12.3.1)
3. A subscript separated by a dot from a data type symbol.

The standard stack notation does not define any delimiter for
subscripts, but some delimiter should be used in plain text if the
subscript is not a number. It is convenient to use a dot as a separator.


--
Ruvim

[dxf]

Post by Ruvim

Post by minforth
TEMP-CONVERSION ( F: celsius -- F: fahrenheit )
should be clear enough, right?

Right, while you have not introduced subtypes of i*r.
  TEMP-CONVERSION ( F: r.celsius -- r.fahrenheit )

Unnecessarily complicated.

[Ruvim]

Post by dxf

Post by Ruvim

Post by minforth
TEMP-CONVERSION ( F: celsius -- F: fahrenheit )
should be clear enough, right?

Right, while you have not introduced subtypes of i*r.
  TEMP-CONVERSION ( F: r.celsius -- r.fahrenheit )

Unnecessarily complicated.

The Forth standard uses namely this notation.

Stack diagrams are part of the program (or library) documentation. And
what notations are used for documentation is up to authors.

My initial question was about the standard stack notation.


--
Ruvim

[a***@spenarnc.xs4all.nl]

Post by Ruvim

Post by minforth
TEMP-CONVERSION ( F: celsius -- F: fahrenheit )
should be clear enough, right?

Right, while you have not introduced subtypes of i*r.
TEMP-CONVERSION ( F: r.celsius -- r.fahrenheit )
1. A stack-id ("F" in this case) is used once.
2. The data type symbol "r" is used (it's defined in 12.3.1)
3. A subscript separated by a dot from a data type symbol.
The standard stack notation does not define any delimiter for
subscripts, but some delimiter should be used in plain text if the
subscript is not a number. It is convenient to use a dot as a separator.

I don't standardizing stack notations are worth worrying about.
The standards document has to be precise, so there worries are warranted.

Post by Ruvim
Ruvim

Groetjes Albert

[Hans Bezemer]

Post by Anton Ertl

Post by Ahmed
Hello,
https://www.novabbs.com/devel/article-flat.php?group=comp.lang.forth&id=16944&first=1&last=25#start
For example, sin (sine) is a function, whose input is a real and
whose
output is another real. And "derivative" (d/dx) is a function whose
input is a function, and whose output is another function. So you
get
that d/dx (sin) = cos, or some numerical approximation of same. In
h = 0.0001
d_over_dx f = g where { g x = (f(x+h) - f(x)) / h }
dsin = d_over_dx sin
print (dsin 0.5) -- prints 0.8775585891507287
print (cos 0.5) -- prints 0.8775825618903728
Pretty close approximation. In Forth it is a bit messier because you
want a signature like
: d/dx ( xt -- xt ) .... ;

Most probably he meant ( xt1 -- xt2 )

Post by Ahmed
and you may have to do some advanced or non-portable things to create
new xt's on the fly like that.
: d/dx ( xt -- )
create ,
execute f- 0.5e4 f* ;
' func() d/dx der_func()

' fsin d/dx fdsin
And given that he gives a name to the result anyway, it's obviously
fine that your solution requires a name and produces a named word
rather than an unnamed word. Your solution uses a slightly different
(unbiased) formula, though.

Post by Ahmed
for example: (in gforth)
' fsin d/dx der_fsin ok
0.5e fcos f. 0.877582561890373 ok
0.5e der_fsin f. 0.877582560427637 ok

The smaller difference from the fcos result demonstrates the advantage
of the unbiased formula.
A solution that actually complies with the ( xt1 -- xt2 ) requirement
and is in standard Forth is as follows (using the biased formula in
0.0001e fvalue h
: d/dx1 ( r1 xt -- r2 )
fdup h f+ dup execute fswap execute f- h f/ ;
: d/dx ( xt1 -- xt2 )

Post by Ahmed
r :noname ( r1 -- r2 ) r> postpone literal postpone d/dx1 postpone ; ;

0.5e ' fsin d/dx execute f. \ prints 0.877558589150729
Maybe he considers the use of DOES> or POSTPONE to be "advanced
things". The disadvantage of the ( xt1 -- xt2 ) variant is that it
may seduce the user to produce multiple instances of the same function
(possibly one instance for every iteration of some loop), with every
instance costing memory that is unlikely to be reclaimed before the
end of the session. With a version that produces a named word that is
far less likely.
The advantage of the unnamed variant is that we can do things like
defer fddsin
' fsin d/dx d/dx is fddsin
without naming the intermediate function. Looking at the result, this
0.5e fddsin f. \ prints -0.479513295736922
0.5e fsin fnegate f. \ prints -0.479425538604203
- anton

Well, with a little tinkering (using the preprocessor) it works for me too:

include lib/fp2.4th

float array (h) 1 s>f 10000 s>f f/ latest f! does> f@ ;

:macro d/dx >#> @dup :noname fdup (h) f+ _#_ fswap _#_ f- (h) f/ ` ;` ;
:macro d/dx: : >#> _#_ >#> @dup fdup (h) f+ _#_ fswap _#_ f- (h) f/ ` ;` ;

[DEFINED] ZenFP [IF] \ can be used with ZEN or ANS
include lib/zenfsin.4th
[ELSE]
include lib/fsinfcos.4th
[THEN]

include 4pp/lib/float.4pp \ making life easier
include 4pp/lib/ddivdx.4pp \ now include the payload

fclear
13 set-precision

f% 0.5e d/dx fsin execute f. cr \ prints "0.877558589150729"

d/dx: dsin fsin \ now define a word named DSIN
f% 0.5e dsin f. cr \ and execute it..

Hans Bezemer

BTW..

Post by Anton Ertl
#> = parse next word, put it on the string stack

@dup = duplicate top of the string stack
_#_ = print the top of the string stack with a space
` xx` = do not interpret this string

[Krishna Myneni]

On 5/24/24 04:57, Ahmed wrote:
;;;

 that d/dx (sin) = cos, or some numerical approximation of same.  In
   h = 0.0001
   d_over_dx f = g where { g x = (f(x+h) - f(x)) / h }

...

Side note: you may already be aware that the accuracy of the derivative
improves significantly if you average the forward and backward slopes,

df/dx \approx [(f(x+h) - f(x)) + (f(x) - f(x-h))] / 2h

= [f(x+h) - f(x-h)] / 2h

The error goes as O(h^2) rather than as O(h).

This is the method I implement for computing a numerical derivative in
my XYPLOT-32 program:

https://github.com/mynenik/XYPLOT-32/blob/master/modules/fsl/extras/derivative.4th

--
Krishna

[Ahmed]

Hello,
Thanks for the link.

The forward derivative was posted by Paul Robin.
I just asked why not do it simply and directly like this:

: d/dx ( xt -- ) create , does> @ dup ( f: x -- y) fdup 1e-4 f+ execute
fswap 1e-4 f- execute f- 0.5e4 f* ;

where h = 1e-4 for example.

It is the same as that you proposed.

I have a question concerning your derivative implementation (the link
you've given):
in the word derivative, we find:

dx1 F@ F0= dx2 F@ F0= or
IF 1 unloop EXIT THEN

why you use or in the condition? why not and?

I think you implemented:
the derivative at (x1, f(x1)) as

d/dx f(x1) = (( (f(x2)-f(x1))/(x2-x1))+(( (f(x1)-f(x0))/(x1-x0)))/2,
where (x0,f(x0)), (x1, f(x1)) and (x2, f(x2)) are successive points.

and that is diffrent from
d/dx f(x1) = (f(x2)-f(x0))/(x2-x0) when (x1-x0) != (x2-x1)
the latter can be written d/dx f(x1) = (f(x1+h) - f(x1-h))/(2h) wher h =
x1-x0 = x2-x1
the formula you've given in your reply.

Thanks.

[Krishna Myneni]

On 5/25/24 02:34, Ahmed wrote:
...

Post by Ahmed
I have a question concerning your derivative implementation (the link
    IF 1 unloop EXIT  THEN
why you use or in the condition? why not and?
the derivative at (x1, f(x1)) as

The routine stops the calculation whenever it encounters adjacent points
which have dx = 0. Thus either dx1 being zero OR dx2 being zero results
in an error condition -- the derivative calculation is considered invalid.

--
Krishna

[Ahmed]

Hi,
Thanks for the explanation.
So it is like a jump where the derivative is not defined (not in the
sens of distributions where it is given by Dirac delta distribution (or
generalized function)).

Thank you again.

[Ahmed]

Hi again,

But in that case, one still can calculate the derivative
d/dx (f(x1)) = (f(x2)-f(x0))/(x2-x0)
where (x0,f(x0)), (x1, f(x1)) and (x2, f(x2)) are three consecutive
points with:
- x1 = x0 and f(x1) != f(x0) and x2 != x1.
or
- x2 = x1 and f(x2) != f(x1) and x1 != x0.

The problem is when x2 = x1 = x0, and that's why I asked why not "and"
in place of "or".

Thanks again.

Best Regards.

[Krishna Myneni]

Post by Ahmed
Hi again,
But in that case, one still can calculate the derivative d/dx (f(x1)) =
(f(x2)-f(x0))/(x2-x0) where (x0,f(x0)), (x1, f(x1)) and (x2, f(x2)) are
three consecutive
- x1 = x0 and f(x1) != f(x0) and x2 != x1.
or - x2 = x1 and f(x2) != f(x1) and x1 != x0.
The problem is when x2 = x1 = x0, and that's why I asked why not "and"
in place of "or".

You can still calculate a local derivative if one of the intervals, dx1
or dx2 is not zero; however, my preference in this case is that an error
be flagged for the user, since the data is itself problematic. The user
needs to consider why the data represents a multi-valued function.

--
Krishna

[Ahmed]

Ah, I see.
Thanks for the reply.

[Krishna Myneni]

Post by Ahmed
Ah, I see.
Thanks for the reply.

My derivative routine is basic. For experimental data, it assumes you
have done whatever preprocessing on the raw data to provide a reasonably
smooth approximation of an analytic function to it.

In measurements, one is rarely dealing with data which, in raw form, is
an analytic function. It's possible that the data contains repeated
measurements of Y at the same control parameter X, and Y is a random
variable drawn from some distribution. Then, one could write a
derivative routine which computes the statistics of Y at each X, and
then compute an alternative derivative e.g., a derivative based on the
mean y-values.

I have not implemented such a method for XYPLOT, although one can write
an external Forth module to add that functionality to the program.

Then, of course, there is also the case where X itself may have a
distribution ...

--
Krishna

[minforth]

Even signals disturbed by noise and smoothed by splines
often show considerable deviations in their first derivative.
This can even make them unusable. Simple bandpass filters,
e.g. classic Butterworth, are often better. Finding the right
frequency bands is usually a question of empirical experience.

[Krishna Myneni]

Post by minforth
Even signals disturbed by noise and smoothed by splines
often show considerable deviations in their first derivative.
This can even make them unusable. Simple bandpass filters,
e.g. classic Butterworth, are often better. Finding the right
frequency bands is usually a question of empirical experience.

I'm wary of using filters in connection with computing derivatives.
Filters will exhibit latency and therefore alter some properties of a
signal. One such case might be a repetitive signal consisting of a
short-duration pulse followed by a long-duration pulse. If we want to
obtain a measure of the peak to peak time interval between the short and
long pulse within a cycle, the distortion of the waveform by the filter
has to be considered.

Alternatively, if a trigger pulse is available to mark the start of the
cycle, the repetitive signal can be averaged until the noise is low
enough to use a derivative to obtaining the zero crossings. Curve
fitting is generally a better technique though, when a model function
exists.

--
Krishna

[minforth]

Post by Krishna Myneni

Post by minforth
Even signals disturbed by noise and smoothed by splines
often show considerable deviations in their first derivative.
This can even make them unusable. Simple bandpass filters,
e.g. classic Butterworth, are often better. Finding the right
frequency bands is usually a question of empirical experience.

I'm wary of using filters in connection with computing derivatives.
Filters will exhibit latency and therefore alter some properties of a
signal.

That's why Butterworths are a good compromise. They exhibit only
slight phase distortion. Doing eg FFT with raw signals using ends
in garbage.

Post by Krishna Myneni
One such case might be a repetitive signal consisting of a
short-duration pulse followed by a long-duration pulse. If we want to
obtain a measure of the peak to peak time interval between the short and
long pulse within a cycle, the distortion of the waveform by the filter
has to be considered.

That's a typical application for matched filters.

Post by Krishna Myneni
Alternatively, if a trigger pulse is available to mark the start of the
cycle, the repetitive signal can be averaged until the noise is low
enough to use a derivative to obtaining the zero crossings. Curve
fitting is generally a better technique though, when a model function
exists.

My quasi daily "aha"s were hidden 50 Hz (60 Hz in the US) hums within
signals. Curve fitting happily counts them in, whereas a simple notch
filter can take care of that.

[a***@spenarnc.xs4all.nl]

Post by Krishna Myneni

Post by Ahmed
Ah, I see.
Thanks for the reply.

My derivative routine is basic. For experimental data, it assumes you
have done whatever preprocessing on the raw data to provide a reasonably
smooth approximation of an analytic function to it.
In measurements, one is rarely dealing with data which, in raw form, is
an analytic function.

An analytic function is defined in an infinitly many immeasurable
points where only Descartes and Cantor pretend to know the properties off.
Seriously?
At most an analytic function is a pretend thing, not a property of real
data, not even "rarely".

Post by Krishna Myneni
--
Krishna

Groetjes Albert

[Krishna Myneni]

Post by a***@spenarnc.xs4all.nl

Post by Krishna Myneni

Post by Ahmed
Ah, I see.
Thanks for the reply.

My derivative routine is basic. For experimental data, it assumes you
have done whatever preprocessing on the raw data to provide a reasonably
smooth approximation of an analytic function to it.
In measurements, one is rarely dealing with data which, in raw form, is
an analytic function.

An analytic function is defined in an infinitly many immeasurable
points where only Descartes and Cantor pretend to know the properties off.
Seriously?
At most an analytic function is a pretend thing, not a property of real
data, not even "rarely".

Good grief. We shouldn't even be talking about numerical differentiation
then.

--
Krishna

[a***@spenarnc.xs4all.nl]

Post by Krishna Myneni

Post by a***@spenarnc.xs4all.nl

Post by Krishna Myneni

Post by Ahmed
Ah, I see.
Thanks for the reply.

My derivative routine is basic. For experimental data, it assumes you
have done whatever preprocessing on the raw data to provide a reasonably
smooth approximation of an analytic function to it.
In measurements, one is rarely dealing with data which, in raw form, is
an analytic function.

An analytic function is defined in an infinitly many immeasurable
points where only Descartes and Cantor pretend to know the properties off.
Seriously?
At most an analytic function is a pretend thing, not a property of real
data, not even "rarely".

Good grief. We shouldn't even be talking about numerical differentiation
then.

In practical manners it is okay to talk about numerical differentiation.
For instance implementing a pid regulator.
The derivative calculated numerically serves a great purpose. The thought
that a jurky motion of a piston is or isn't an analytical function, and
can extended to minus infinity in the complex plane is foreign to me.

Post by Krishna Myneni
--
Krishna

Groetjes Albert

[ahmed]

Hi,
I think Mr K. M. meant "analytic formula fitting the acquired (measured)
data".
And analytic functions are a class of functions over the complexe plane
(C) with some properties that give them certain characteristics
(derivation for example).
See the link
https://en.wikipedia.org/wiki/Analytic_function

For automatic control we use numeric derivation for modeling,
identification and control of systems based on sampled time data
(transfer functions in z plane, z-transform).
The approximations used are for example:
Euler forward
Euler backward
Tustin
And the digital PID controller is just an example of discrete systems
with a derivation action
(d/dt(in time domaine)<--->s (Laplace transform)<--->(1-z^(-1))/Ts (in z
domaine, Backwrad Euler)
(d/dt(in time domaine)<--->s (Laplace transform)<--->(z-1)/Ts (in z
domaine, Forward Euler)
(d/dt(in time domaine)<--->s (Laplace transform)<--->(2/Ts)*(z-1)/(z+1)
(in z domaine, Tustin)

Ts is the sampling period.

I used a PID controller with arduino-based real-time simulation and that
worked as required and a comparison with a simulation under
matlab/simulink confirmed the results.

Best regards.

[a***@spenarnc.xs4all.nl]

Post by ahmed
Hi,
I think Mr K. M. meant "analytic formula fitting the acquired (measured)
data".
And analytic functions are a class of functions over the complexe plane
(C) with some properties that give them certain characteristics
(derivation for example).
See the link
https://en.wikipedia.org/wiki/Analytic_function

An analytic function has the properties that if you know a
tiny environment -- mathematically speaking --
you can calculate all derivatives and use these in a series
formula to calculate the function anywere.
So If you know the orbit of the earth for the next 10 mS
you could calculate the position to all eternity
(forgetting about 3 body problems and chaos theory.)
This is all nice and dandy for the Rieman Zeta function but not
for this situation.

(I have a degree a theoretical physics. I'm not worth a
dime in physics, but this mathematics subject I was quite
good at.)

Post by ahmed
Best regards.

[Ahmed]

Post by a***@spenarnc.xs4all.nl

Post by ahmed
Hi,
I think Mr K. M. meant "analytic formula fitting the acquired (measured)
data".
And analytic functions are a class of functions over the complexe plane
(C) with some properties that give them certain characteristics
(derivation for example).
See the link
https://en.wikipedia.org/wiki/Analytic_function

An analytic function has the properties that if you know a
tiny environment -- mathematically speaking --
you can calculate all derivatives and use these in a series
formula to calculate the function anywere.
So If you know the orbit of the earth for the next 10 mS
you could calculate the position to all eternity
(forgetting about 3 body problems and chaos theory.)
This is all nice and dandy for the Rieman Zeta function but not
for this situation.
(I have a degree a theoretical physics. I'm not worth a
dime in physics, but this mathematics subject I was quite
good at.)

Agreed.

Best Regards

Post by a***@spenarnc.xs4all.nl

Post by ahmed
Best regards.