\$title 'Test correctness of sqrec intrinsic' (FNSQREC2,SEQ=726) \$ontext sqrec(x,S) = 1/x if x >= S 1/S - r/S^2 + r^2/s^3 otherwise, r = x - S for S >= 1e-10, default S = 1e-10. Contributor: Steve, Mar 2017 \$offtext \$include fnset_xy.inc reps = 8e-16; sets T / t1 * t6, t11 * t16 t21 * t26 / torig(t) 'unsmoothed inputs' ; table data(T,V) x y f_ fx_ fxx_ rc_ ec_ t1 1 1e-10 t2 1e-6 1e-10 t3 1e-9 1e-10 t4 1e-20 1e-10 t5 EPS 1e-10 t6 -1e+2 1e-10 t11 1 1e-5 t12 1e-3 1e-5 t13 1e-4 1e-5 t14 1e-99 1e-5 t15 EPS 1e-5 t16 -1e+20 1e-5 t21 4000 2000 t22 2001 2000 t23 2000 2000 t24 1999.999 2000 t25 EPS 2000 t26 -1e20 2000 ; scalars r, ros ; * 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+3e-16) * data(T,'x'); }; torig(t) = [data(T,'x') > data(T,'y')]; loop{torig(T), data(T, 'f_') = 1 / data(T,'x'); data(T, 'fx_') = -data(T,'f_') / data(T,'x'); data(T,'fxx_') = -2 * data(T,'fx_') / data(T,'x'); }; loop{t\$[not torig(t)], r = data(t,'x') - data(t,'y'); ros = r / data(t,'y'); data(t, 'f_') = ( (ros - 1) * ros + 1) / data(t,'y'); data(t, 'fx_') = (2*ros-1) / data(t,'y') / data(t,'y'); data(t,'fxx_') = (2/data(t,'y')) / data(t,'y') / data(t,'y'); }; loop {T, data(T,'fyx_') = data(T,'fxy_'); data(T, 'f') = sqrec.value( data(T,'x'),data(T,'y')); data(T, 'fx') = sqrec.grad(1: data(T,'x'),data(T,'y')); data(T, 'fy') = sqrec.grad(2: data(T,'x'),data(T,'y')); data(T,'fxx') = sqrec.hess(1:1:data(T,'x'),data(T,'y')); data(T,'fxy') = sqrec.hess(1:2:data(T,'x'),data(T,'y')); data(T,'fyx') = sqrec.hess(2:1:data(T,'x'),data(T,'y')); data(T,'fyy') = sqrec.hess(2:2:data(T,'x'),data(T,'y')); data(T, 'rc') = mathlastrc; data(T, 'ec') = mathlastec; }; display data; \$include fntest_xy.inc