xref: /freebsd/contrib/arm-optimized-routines/pl/math/tools/v_erfc.sollya (revision 072a4ba82a01476eaee33781ccd241033eefcf0b)
1*072a4ba8SAndrew Turner// polynomial for approximating erfc(x)*exp(x*x)
2*072a4ba8SAndrew Turner//
3*072a4ba8SAndrew Turner// Copyright (c) 2022-2023, Arm Limited.
4*072a4ba8SAndrew Turner// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
5*072a4ba8SAndrew Turner
6*072a4ba8SAndrew Turnerdeg = 12; // poly degree
7*072a4ba8SAndrew Turner
8*072a4ba8SAndrew Turneritv = parse(__argv[0]);
9*072a4ba8SAndrew Turner
10*072a4ba8SAndrew Turnerbounds = [|3.725290298461914e-9,
11*072a4ba8SAndrew Turner           0.18920711500272103,
12*072a4ba8SAndrew Turner           0.41421356237309515,
13*072a4ba8SAndrew Turner           0.681792830507429,
14*072a4ba8SAndrew Turner           1,
15*072a4ba8SAndrew Turner           1.378414230005442,
16*072a4ba8SAndrew Turner           1.8284271247461903,
17*072a4ba8SAndrew Turner           2.363585661014858,
18*072a4ba8SAndrew Turner           3,
19*072a4ba8SAndrew Turner           3.756828460010884,
20*072a4ba8SAndrew Turner           4.656854249492381,
21*072a4ba8SAndrew Turner           5.727171322029716,
22*072a4ba8SAndrew Turner           7,
23*072a4ba8SAndrew Turner           8.513656920021768,
24*072a4ba8SAndrew Turner           10.313708498984761,
25*072a4ba8SAndrew Turner           12.454342644059432,
26*072a4ba8SAndrew Turner           15,
27*072a4ba8SAndrew Turner           18.027313840043536,
28*072a4ba8SAndrew Turner           21.627416997969522,
29*072a4ba8SAndrew Turner           25.908685288118864,
30*072a4ba8SAndrew Turner           31|];
31*072a4ba8SAndrew Turner
32*072a4ba8SAndrew Turnera = bounds[itv];
33*072a4ba8SAndrew Turnerb = bounds[itv + 1];
34*072a4ba8SAndrew Turner
35*072a4ba8SAndrew Turnerf = proc(y) {
36*072a4ba8SAndrew Turner  t = y + a;
37*072a4ba8SAndrew Turner  return erfc(t) * exp(t*t);
38*072a4ba8SAndrew Turner};
39*072a4ba8SAndrew Turner
40*072a4ba8SAndrew Turnerpoly = fpminimax(f(x), deg, [|double ...|], [0;b-a]);
41*072a4ba8SAndrew Turner
42*072a4ba8SAndrew Turnerdisplay = hexadecimal;
43*072a4ba8SAndrew Turnerprint("rel error:", accurateinfnorm(1-poly(x)/f(x), [a;b], 30));
44*072a4ba8SAndrew Turnerprint("in [",a,b,"]");
45*072a4ba8SAndrew Turnerprint("coeffs:");
46*072a4ba8SAndrew Turnerfor i from 0 to deg do coeff(poly, i);
47