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

Angel Montesinos kirjoitti torstai 24 helmikuu 2011 20:54:48:

Are there any visualisation features or is the output in text format, too? Sounds like your solver competes with Maxima and other math programs.

Juha

In our previous episode, Angel Montesinos said:

Have a look at FPC package symbolic. It sounds like roughly the same kind of soup. (parsing, differentiating, fast eval)

Are there any visualisation features or is the output in text format, too? Sounds like your solver competes with Maxima and other math programs

Ok, I didn't look at the prog before asking.
I booted Windows and tried it.

> Perhaps the prettiest of them is Superficies_En, wich is an English

version.

I think the Superficies version 2 is newer and preferred for porting. The 3-D images are rotated nicely without OpenGL support. Good.

I tried to convert it but there are 3rd party components which are not ported to Lazarus.
Also, there are references to TBitmap.ScanLine which is not portable. The Superficies v.1 seems to have much more of them, so the code has become more portable.

What are the 64-bit problems you faced?

Regards,
Juha Manninen

Marco van de Voort wrote

Have a look at FPC package symbolic. It sounds like roughly the same kind of soup. (parsing, differentiating, fast eval)

Yes, roughly. Because if I have read it correctly, the package evaluates things ultimately by calls to the system functions like sin, arctan, *, etc. This is very different than building on the fly the minimal object code for evaluating the function, and much slower.

Also, for derivatives, Symbolic builds first the expression that is the derivative of the initial expression. When the derivative is called, it evaluates the new function.

Automatic differentiation does not compute the derivative function, but at each step of the evaluation of the RPN it evaluates the function and the first (or the function, the first and the second) derivative of that step. For instance, if the RPN step is "cosine", the Symbolic approach is to calculate first -sine (and for the second derivative, to calculate then _again_ the cosine, etc.). But in my approach, instead of computing the cosine, it is computed the sincos pair by using the x87 instruction of AMD or Intel. The two values are used directly for computing the first and second derivatives: thus at most _one_ call to an expensive FPU instruction.

In my approach, if a call of a RPN step costs n clock cycles, then the call of the first derivative of that function takes a cost c*n clock cycles, where c depends only on the function called (that is, whether it is * or arctan, etc). And the same (with another greater coefficient) with the second derivative. Suppose that c = 3 for the operator * and one has to compute the value of x*y*z*w, where these symbols denote variables or expressions on which the formula depends. Then Symbolic will compute

x'*y*z*w + x*y'*z*w + x*y*z'*w + x*y*z*w',

that is roughly four times the computation of the function. Automatic differentiation will take only three (= c) times more (this should be taken only as a sloppy explanation). Thus, the more complex the formula the better relative performance.

On the other hand, porting my library to other architectures will be possible if they have an FPU based on a stack of perhaps more than six registers, a collection of instructions that move numbers between them and between them and RAM, and a minimal set of instructions for each of the basic mathematical functions of Pascal: +,-,*,/,sqr, sqrt, sin, cos, tan, arctan, abs, exp, log. But I have no idea about those architectures.

Please, take a look to the file OCFunctDDDoubleSafeIntegrals.pas in Superficies_En. Then perhaps you may understand that it is a sort of minicompiler that converts strings to object code: a minicompiler that cannot be based on fpc. Thus, the problem of porting to 64bits (or by the same token to other architectures) is to change the opcodes used in building that object code, and also to change the relation between the function code and the preparation that fpc must make prior to calling that code, and the use that fpc will do of the results of that code.

An added difficulty that I have thumped with is that the debugger of fpc-Lazarus does not move, by pressing F7, from instruction to instruction in the CPU window, so that is difficult to catch where the errors arise. Another one, lesser, is that the assembler of the debugger CPU window is different from that used in AMD or Intel description of their architecture and machine instructions.

I updated yesterday Superficies so that now it may be used when the DEP (Data Execution Prevention) of Windows is enabled for all programs. I think also that by windows emulation it may run in Linux. I am working now in making DEP compatible the other programs in my web page.

I know about ScanLine and components problems, and think that they are not difficult to bypass by using GR32 or other tricks. But this is not my did not use OpenGl mainly because I believe that is does not provide a mechanism for identifying on what point of the drawn surface is the mouse when it moves on the surface. This is a must for many things in my programs.

I downloaded the files from your website to have a look.

My first thought was 'I am looking at a program written in Reverse Polish Notation that is literally translated into processor machinecode.' RPN is the 'native language' of Hewlett-Packard calculators.

My second thought was that in order to be able to maintain this code you must exactly know what it does, hence you must be a master in the art of differential geometry & a master in processor machinecode.

The first is already a steep requirement, the second can be lethal when you are talking about implementations on different processors ;)

Regards / Cees

And a few lines down in the stack of messages you say so yourselve '... at each step of the evaluation of the RPN it evaluates the ...'

Is it then not easier to make a proper RPN-stack processor and feed the RPN-code into it?
There are many versions already at http://www.hpcalc.org/hp48/pc/emulators/

Forge these into a proper 'object' written in Free Pascal and problem solved. No need for endless translations into processor machinecode.

The maintenance is then restricted to maintaining 'the object' and keeping an eye on additions to RPN-stack manipulation.

Are there also changes in individual RPN-instructions? I don't keep track and have not touched RPN in twenty years ;)

Regards / Cees

In our previous episode, Angel Montesinos said:

I know. I've played with going a step inbetween (the OS independant evaluator, my current approach, and your approach, code generation), namely a CPU dependant codegenerator.

I never really started it, but one of the things I noticed is that the CPU instructions for certain goniometric functions don't always accept the full domain. How have you dealt with that in your code generation (or am I wrong, it is over 7 years ago that I made this)

I don't fully understand this. In my approach I calculate d(cos(x))/dx = -sin(x) and then d(-sin(x)/dx = -cos(x) and only then use one floating
- point instruction to calculate the cosinus.

Or do you mean for stuff like Taylor where each term is a further differentation?

The idea with symbolic was mainly to quickly draw graphs. Since the evaluation is then done thousands of times, the symbolic manipulation is not really optimized.

The whole project was a combination of needing something to draw graphs of computed values quick (and derivatives were a requirement). I invested a bit more time because I wanted to play a bit with syntax tree manipulation, just as an exercise.

The evaluation speed was fine for my purposes then, considering the chosen technique (OS independant RPN interpreter). But I'd be interested to see how all this pans out.

Actually the problems I have with symbolic is the lack of custom functions and boolean logic, rather than speed. I have some half forgotten CVS with an half finished attempt at adding that.

Does your engine do boolean logic?

I planned to eliminate this using a common subexpression tactic on the final result. But I never got around to it (and only understood later that it is relatively hard). Maybe indeed the derivation routine should be able to do prepatory work for this.

Hardware FPU stack is afaik a x86 oddity. Most other FPUs are register based. But there is no one true solution. I'd take an evaluator as baseline (for unknown architectures) and then plug in better solutions for specific CPus. Be careful with border conditions though.

The evaluator can then also be used in case e.g. the codegeneration has a problem, even on a known arch.

Hello,

Angel Montesinos wrote

An added difficulty that I have thumped with is that the debugger of fpc-Lazarus does not move, by pressing F7

I think this may be solved already in the latest Lazarus version: "Single Instruction step (over/into)"
(http://wiki.lazarus.freepascal.org/Lazarus_0.9.30_release_notes). What Lazarus version are you using? If it is not RC of 0.9.30, then you should try using it. You can find it at:
http://www.hu.freepascal.org/lazarus/ and it is called "Lazarus-0.9.29-*".

Regards
Žilvinas Ledas

Marco van de Voort wrote

I never really started it, but one of the things I noticed is that the CPU instructions for certain goniometric functions don't always accept the full domain. How have you dealt with that in your code generation (or am I wrong, it is over 7 years ago that I made this)

Yes, the range of argument is -2^63 to 2^63. I follow a suggestion of INTEL: first apply a mod pi or mod 2*pi call for the argument (there is an instruction for this, FPREM), and then apply cos, sin or sincos. This should guarantees better performance, but for greater arguments the FPREM becomes simply innacurate and useless: thus it only serves as a dodge for the exception that would throw the FPU.

My rationale is that in serious mathematics, most time it is useless to compute the derivatives without computing also the value of the function. Thus, the call of sincos in a RPN step will give the function value and, by a further multiplication and sign change, that of its derivative. Also, for other expensive functions, the derivative is a simple function of the function value. For instance, tan' = 1 + tan^2. So, once computed tan (the function), we obtain easily also the derivative.

Also, if with your approach you must compute the derivative of the function
x * cos(x),
you will compute
cos(x)- x*sin(x).
And for the second derivative
-sin(x) - sin(x) - x* cos(x),
because the computation of simplifications is an enormous task: therefore you will call two ot three expensive trig functions. My approach will compute only sincos once.

Another thing that is easy in my approach is to call from the object code any external function given in software: I use this for computing functions given by integrals (though this with much less precision).

Or do you mean for stuff like Taylor where each term is a further differentation?

No, I have seen this in your library, but I don'use Taylor.