\$title 'Test correctness of sqexp intrinsic' (FNSQEXP2,SEQ=720) \$ontext sqexp(x,S) = exp(x) if x <= S exp(S) * (1 + r + r^2/2) otherwise, r = x - S for S <= 150, default S = 150. Contributor: Steve, Mar 2017 \$offtext \$include fnset_xy.inc scalar xd; reps = 8e-16; sets T / t1 * t6, t11 * t16 t21 * t26 / tmap(t,t) / t3. (t4*t6) t13.(t14*t16) t23.(t24*t26) / tc(t) ; alias (t,t2); tc(t2) = sum{t, tmap(t,t2)}; table data(T,V) x y f_ fx_ fxx_ rc_ ec_ t1 1 150 t2 149 150 t3 150 150 t4 150.001 150 t5 151 150 t6 1e117 150 t11 1 20 t12 19 20 t13 20 20 t14 20.0001 20 t15 30 20 t16 1e145 20 t21 -50 -20 t22 -21 -20 t23 -20 -20 t24 -19.9999 -20 t25 0 -20 t26 1e148 -20 ; * loop{T\$[data(T,'x') <= data(T,'y')], loop{T\$[not tc(t)], data(T, 'f_') = exp(data(T,'x')); data(T, 'fx_') = data(T,'f_'); data(T,'fxx_') = data(T,'f_'); }; loop{tmap(t,t2), xd = data(t2,'x')-data(t,'x'); data(t2, 'f_') = data(t,'f_') * (1 + xd + 0.5*sqr(xd)); data(t2, 'fx_') = data(t,'f_') * (1 + xd); data(t2,'fxx_') = data(t,'f_'); }; loop {T, data(T,'fyx_') = data(T,'fxy_'); data(T, 'f') = sqexp.value( data(T,'x'),data(T,'y')); data(T, 'fx') = sqexp.grad(1: data(T,'x'),data(T,'y')); data(T, 'fy') = sqexp.grad(2: data(T,'x'),data(T,'y')); data(T,'fxx') = sqexp.hess(1:1:data(T,'x'),data(T,'y')); data(T,'fxy') = sqexp.hess(1:2:data(T,'x'),data(T,'y')); data(T,'fyx') = sqexp.hess(2:1:data(T,'x'),data(T,'y')); data(T,'fyy') = sqexp.hess(2:2:data(T,'x'),data(T,'y')); data(T, 'rc') = mathlastrc; data(T, 'ec') = mathlastec; }; display data; \$include fntest_xy.inc