Re: [Gretl-devel] fdjac: a better example

On tests: look reasonably well For now I see 2 points: 1) Since 10^3*sqrt($macheps) is ~10^-5 on 32 bit systems there is a problem with, e.g. sqrt(x) near zero. May be, simply make warning in the help text? 2) strictly speaking sqrt(abs(x)) is incorrect and should be replaced by a linear function of abs(x). Still, I don't know what exactly it should be. For quality 0 the following works reasonably well: (for sin(x) for x in (0, 5*10^7) ) if abs(x)<1       h=sqrt($macheps) else      h= (1+(2^(-12))*(abs(x)-1))*sqrt($macheps) endif Also, I think it will be nice to also have  the left-sided version of quality 0: for 1/x with small negative x, for example Oleh 30 грудня 2015, 12:48:19, від "Riccardo (Jack) Lucchetti" < r.lucchetti(a)univpm.it >: On Wed, 30 Dec 2015, oleg_komashko(a)ukr.net wrote: Of course ANY numerical approximation to a derivative is a ratio; therefore, what is crucial to its precision is the relative order of magnitude of numerator and denominator. It's easy to contruct pathological cases and there's no silver bullet that works in all cases. That said, the following settings seem to work reasonably well in a decently-sized array of cases (diff against git HEAD --- not applied yet). <diff> diff --git a/minpack/fdjac2.c b/minpack/fdjac2.c index 9e9676a..54e64ec 100644 --- a/minpack/fdjac2.c +++ b/minpack/fdjac2.c @@ -144,10 +144,10 @@ int fdjac2_(S_fp fcn, int m, int n, int quality, for (j = 0; j < n; j++) { temp = x[j]; if (quality == 2) { - if (fabs(temp) > 1000) { - h = 1000 * sqrt(fabs(temp/1000)) * eps; + if (fabs(temp) > 128) { + h = 128* sqrt(fabs(temp/128)) * eps / sqrt(6); } else { - h = 1000 * eps; + h = 128 * eps / sqrt(6); } } else { h = eps * fabs(temp); </diff> I'm attaching a battery of tests. ------------------------------------------------------- Riccardo (Jack) Lucchetti Dipartimento di Scienze Economiche e Sociali (DiSES) Università Politecnica delle Marche (formerly known as Università di Ancona) r.lucchetti(a)univpm.it http://www2.econ.univpm.it/servizi/hpp/lucchetti -------------------------------------------------------