[fpc-pascal] interested in building a library for functions?

[Angel Montesinos]

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