\$title 'Test correctness of sqlog10 intrinsic' (FNSQLOG102,SEQ=724) \$ontext sqlog10(x,S) = log_10(x) if x >= S (1/log 10) (ln(S) + r/S - 1/S^2) otherwise, r = x - S for S >= 1e-150, default S = 1e-150. Contributor: Steve, Mar 2017 \$offtext \$include fnset_xy.inc reps = 8e-16; sets T / t1 * t6, t11 * t16 t21 * t26 / torig(t) 'unsmoothed inputs' ; * doing anything with y = 1e-150 is pretty useless: the 2nd deriv is * so large table data(T,V) x y f_ fx_ fxx_ rc_ ec_ t1 1 1e-140 t2 1e-136 1e-140 t3 1e-139 1e-140 t4 1e-160 1e-140 t5 EPS 1e-140 t6 -1e1 1e-140 t11 1 1e-20 t12 1e-19 1e-20 t13 1e-20 1e-20 t14 1e-120 1e-20 t15 EPS 1e-20 t16 -1e40 1e-20 t21 4000 2000 t22 2001 2000 t23 2000 2000 t24 1999.999 2000 t25 EPS 2000 t26 -1e100 2000 ; scalars r L10 lnS ; * problematic to test exactly at the breakpoint: should 2nd deriv be 0 * or taken from the smooth part? loop{T\$[data(T,'x') = data(T,'y')], data(T,'x') = (1+4e-16) * data(T,'x'); }; torig(t) = [data(T,'x') > data(T,'y')]; L10 = log(10); loop{torig(T), data(T, 'f_') = log(data(T,'x')) / L10; data(T, 'fx_') = 1 / (data(T,'x') * L10); data(T,'fxx_') = (-1 / data(T,'x')) / (data(T,'x')*L10); }; loop{t\$[not torig(t)], lnS = log(data(t,'y')); r = data(t,'x') - data(t,'y'); data(t, 'f_') = (lnS + r/data(t,'y') - 0.5*sqr(r/data(t,'y'))) / L10; data(t, 'fx_') = (1 - r/data(t,'y')) / (data(t,'y')*L10); data(T,'fxx_') = (-1 / data(t,'y')) / (data(t,'y')*L10); }; loop {T, data(T,'fyx_') = data(T,'fxy_'); data(T, 'f') = sqlog10.value( data(T,'x'),data(T,'y')); data(T, 'fx') = sqlog10.grad(1: data(T,'x'),data(T,'y')); data(T, 'fy') = sqlog10.grad(2: data(T,'x'),data(T,'y')); data(T,'fxx') = sqlog10.hess(1:1:data(T,'x'),data(T,'y')); data(T,'fxy') = sqlog10.hess(1:2:data(T,'x'),data(T,'y')); data(T,'fyx') = sqlog10.hess(2:1:data(T,'x'),data(T,'y')); data(T,'fyy') = sqlog10.hess(2:2:data(T,'x'),data(T,'y')); data(T, 'rc') = mathlastrc; data(T, 'ec') = mathlastec; }; display data; \$include fntest_xy.inc