This is a call for people willing to relieve me in the task of
maintaining a library for parsing, evaluating and differentiating real
functions with almost optimal speed and precision without recourse to
using multiprecision real types. I have it working for years in
Windows 32 bits, and I think that it could be a nice thing to extend
it in fpc-Lazarus to other CPU and FPU architectures.
The hardest thing has been already done and tested: to conceive it and
to build a working model.
1. SOME HISTORY
For almost thirty years I have been a teacher of differential geometry
in Valencia, Spain. I retired two years ago. Along these years I have
written several programs oriented to teaching and research, all for
Windows 32 bits, and all free. Some of them may be retrieved freely
and with source code in
<http://www.uv.es/montesin>
Perhaps the prettiest of them is Superficies_En, wich is an English
version. Dont't be afraid by the language barrier for the others:
almost all sources are written in English.
After studying the code of a library built years ago by Dmitry
Kaschenko and Vladimir Safin, I wrote first another one with a
different parsing mechanism and with some safety measures. Then, I
added the possibility of parsing functions defined with integrals.
Finally, I added differentiation.
All this done in Delphi with some trials in Kylix and Lazarus.
Now, I feel that I lack the energy and memory to learning the 64 bits
amd-intel x67 architecture, and the more so for other ones. Thus, I
may opt by wainting to Delphi 64 bits and then hope for the best, that
is that by examining the debugging CPU window I become able to
translate the library to 64 bits. Or else, to provide others with the
present code (by the way, it has been open in my web page for years)
and with advice about the underlying theory and tricks. If I find
interested people, my first task would be to write with more detail my
notes about the theory and to document better the present code.
2. OVERVIEW
This is how it works from the viewpoint of the programmer. Integrals
are not included here for brevity:
Supppose the programmer needs, for example, to treat a mathematical
problem in which some real (that is, with result Double or Extended)
function in three independent (resp. Double or Extended) variables
denoted x, y, z is to be used. Suppose that the user has input that
function as the string
fstr:= 'cos(py + sin(3.22x+ qz^3)) - 1/2',
where p and q are considered as parameters.
For example, the user wants the program to draw the surface defined
implicitly by the equation
cos(py + sin(3.22x+ qz^3)) - 1/2 = 0
The programmer should create an instance of TOCFunctDDDoubleSafe
by a call as
theF:= TOCFunctDDDoubleSafe.Create;
Then, the programmer parses fstr by a call to
procedure TOCFunctDDDoubleSafe.ParseNewText(t, vars, pars:
AnsiString);
Here, t is the string defining the function, vars is the string
containing the allowed independent variables, and pars is the string
containing the allowed parameters. In our case, we could call
theF.ParseNewText(fstr, 'xyz', 'pq');
If pars <> '', the programmer must call at least one time the
updater of parameter values. That is for instance
theF.UpdateParameters([3.42, -2.2]);
This means that until another call to UpdateParameters, the parameters
p and q will have the values p = 3.42 and q = -2.2.
Now the programmer may obtain the value of the function by a call as
a:= theF.F([1.1, 0, sqrt(5)]),
where a will be the value of the function when x = 1.1, y = 0 and
z = sqrt(5).
For computing the first derivatives of the function, the programmer
must call
theF.BuildFirstDerivative;
Then, the following function is available
TOCFunctDDDoubleSafe.FD(XV: array of Double): Fd;
where
Fd= record
f, v: Double;
end;
And for computing the second derivatives, the programmer must call
theF.BuildSecondDerivative;
and then the following function is available
TOCFunctDDDoubleSafe.FDD(XVA: array of Double): Fdd;
where
Fdd= record
f, v, a: Double;
end;
An example: the value of the function and its first partial derivative
with respect to 'y' at the point (a,b,c) is obtained by calling:
anFD:= theF.FD([a, b, c, 0, 1, 0]);
So, anFD.f shall be the value of the function at (a,b,c), and
anFD.v shall be the value of that partial derivative at (a,b,c).
In the same manner, the value of the function and its first and second
partial derivatives with respect to 'y' at the point (a,b,c) are
obtained by calling:
anFDD:= theF.FDD([a, b, c, 0, 1, 0, 0, 0, 0]);
So, anFDD.f shall be the value of the function at (a,b,c),
anFDD.v shall be the value of the partial derivative with respect to
'y' at (a,b,c), and anFDD.a the value of the second partial
derivative with respect to 'y' at (a,b,c).
The meaning of the three last slots that seem to be useless, and the
manner for computing crossed partial derivatives requires some
elaboration that I can provide to the interested persons.
The important thing is that the computation of the function and its
derivatives is hard-coded in object code without the use of finite
differences. Thus, one may expect for theF.FD a precision similar to
that one could obtaine by having precompiled with Delphi (or fpc) the
function
function afunction(a, b, c: Double): Fd;
begin
fd.x:= cos(p*y + sin(3.22*x+ q*z*sqr(z))) - 1/2
fd.v:= -p*sin(p*y + sin(3.22*x+ q*z*sqr(z)));
end;
and the same for theF.DD. In general, my library shall be faster than
the precompiled function.
3. MATHEMATICAL BASIS OF DIFFERENTIATION CODE
Essentially, it is the application of the chain rule to the object
code. That is, to extend the technique known as "automatic
differentiation" at the object level of code (instead to the high
level of code). I do not know whether this is new: I haven't found it
the computing literature: at least when I began to write it. After I
wrote it, I haven't had the need to search for it. If you know of
something similar, tell me please.
Googling for "automatic differentiation" could give an initial idea of
the underlying mechanism.
If you download the full distro of Superficies_En, you will see in the
source folder a file with title:
"Object code for the first and second differential of a function
with the instruction set x86.rtf"
It is a roughly written report about the thing, excluding integration.
About integration, I use splines for interpolation: see in that folder
the file "AproxFuncBSplinesEnglish.pdf". I think that my technique
for B-spline interpolation is also new, but I am not sure of it.
The relevant pascal unit in the source folder is
OCFunctDDDoubleSafeIntegrals.pas
Best regards
--
montesin at uv dot es
_______________________________________________
fpc-pascal maillist - fpc-pascal@lists.freepascal.org
http://lists.freepascal.org/mailman/listinfo/fpc-pascal
