src/math/expm1.sollya

// Copyright 2020 Google LLC
//
// This source code is licensed under the BSD-style license found in the
// LICENSE file in the root directory of this source tree.

display=hexadecimal !;
print("***IEEE FP32***");
print("log2(e) =", round(log2(exp(1)), single, RN));
minus_ln2_hi = round(-log(2), single, RN);
minus_ln2_lo = round(-log(2) - minus_ln2_hi, single, RN);
print("-log(2):hi =", minus_ln2_hi);
print("-log(2):lo =", minus_ln2_lo);

print("log2(e) * 8 =", round(log2(exp(1)) * 8, single, RN));
minus_ln2_o8_hi = round(-log(2)/8, single, RN);
minus_ln2_o8_lo = round(-log(2)/8 - minus_ln2_o8_hi, single, RN);
print("-log(2):hi / 8 =", minus_ln2_o8_hi);
print("-log(2):lo / 8 =", minus_ln2_o8_lo);

lb = round(-log(2)/2, single, RN);
ub = round(log(2)/2, single, RN);
P = fpminimax(expm1(x), [|1,2,3,4,5|], [|SG...|], [lb; ub], relative, 0);
print("Degree-5 P[expm1(x)] with 0 constraint on [-log(2)/2, log(2)/2] =", horner(P));
print("  -log(2)/2 = ", lb);
print("  +log(2)/2 =", ub);
print("  relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU));

lb = round(-log(2)/2, single, RN);
ub = round(log(2)/2, single, RN);
P = fpminimax(expm1(x), [|1,2,3,4,5,6|], [|SG...|], [lb; ub], relative, 0);
print("Degree-6 P[expm1(x)] with 0 constraint on [-log(2)/2, log(2)/2] =", horner(P));
print("  -log(2)/2 = ", lb);
print("  +log(2)/2 =", ub);
print("  relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU));

lb = round(-log(2)/2, single, RN);
ub = round(log(2)/2, single, RN);
P = fpminimax(expm1(x), [|2,3,4,5|], [|SG...|], [lb; ub], relative, 0+x);
print("Degree-5 P[expm1(x)] with 0+x constraint on [-log(2)/2, log(2)/2] =", horner(P));
print("  -log(2)/2 = ", lb);
print("  +log(2)/2 =", ub);
print("  relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU));

lb = round(-log(2)/2, single, RN);
ub = round(log(2)/2, single, RN);
P = fpminimax(expm1(x), [|2,3,4,5,6|], [|SG...|], [lb; ub], relative, 0+x);
print("Degree-6 P[expm1(x)] with 0+x constraint on [-log(2)/2, log(2)/2] =", horner(P));
print("  -log(2)/2 = ", lb);
print("  +log(2)/2 =", ub);
print("  relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU));

lb = round(-log(2)/8, single, RN);
ub = round(log(2)/8, single, RN);
P = fpminimax(expm1(x), [|2,3,4|], [|SG...|], [lb; ub], relative, 0+x);
print("Degree-4 P[expm1(x)] with 0+x constraint on [-log(2)/8, log(2)/8] =", horner(P));
print("  -log(2)/8 =", lb);
print("  +log(2)/8 =", ub);
print("  relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU));

lb = round(-log(2)/16, single, RN);
ub = round(log(2)/16, single, RN);
P = fpminimax(expm1(x), [|1,2,3,4|], [|SG...|], [lb; ub], relative, 0);
print("Degree-4 P[expm1(x)] with 0 constraint on [-log(2)/16, log(2)/16] =", horner(P));
print("  -log(2)/16 =", lb);
print("  +log(2)/16 =", ub);
print("  relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU));

lb = round(-log(2)/16, single, RN);
ub = round(log(2)/16, single, RN);
P = fpminimax(expm1(x), [|2,3,4|], [|SG...|], [lb; ub], relative, 0+x);
print("Degree-4 P[expm1(x)] with 0+x constraint on [-log(2)/16, log(2)/16] =", horner(P));
print("  -log(2)/16 =", lb);
print("  +log(2)/16 =", ub);
print("  relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU));

lb = round(-log(2)/16, single, RN);
ub = round(log(2)/16, single, RN);
P = fpminimax(expm1(x), [|1,2,3|], [|SG...|], [lb; ub], relative, 0);
print("Degree-3 P[expm1(x)] with 0 constraint on [-log(2)/16, log(2)/16] =", horner(P));
print("  -log(2)/16 =", lb);
print("  +log(2)/16 =", ub);
print("  relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU));

lb = round(-log(2)/16, single, RN);
ub = round(log(2)/16, single, RN);
P = fpminimax(expm1(x), [|2,3|], [|SG...|], [lb; ub], relative, 0+x);
print("Degree-3 P[expm1(x)] with 0+x constraint on [-log(2)/16, log(2)/16] =", horner(P));
print("  -log(2)/16 =", lb);
print("  +log(2)/16 =", ub);
print("  relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU));

lb = round(-log(2)/32, single, RN);
ub = round(log(2)/32, single, RN);
P = fpminimax(expm1(x), [|2,3,4|], [|SG...|], [lb; ub], relative, 0+x);
print("Degree-4 P[expm1(x)] with 0+x constraint on [-log(2)/32, log(2)/32] =", horner(P));
print("  -log(2)/32 =", lb);
print("  +log(2)/32 =", ub);
print("  relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU));

lb = round(-log(2)/32, single, RN);
ub = round(log(2)/32, single, RN);
P = fpminimax(expm1(x), [|2,3|], [|SG...|], [lb; ub], relative, 0+x);
print("Degree-3 P[expm1(x)] with 0+x constraint on [-log(2)/32, log(2)/32] =", horner(P));
print("  -log(2)/32 =", lb);
print("  +log(2)/32 =", ub);
print("  relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU));

lb = round(-log(2)/64, single, RN);
ub = round(log(2)/64, single, RN);
P = fpminimax(expm1(x), [|2,3|], [|SG...|], [lb; ub], relative, 0+x);
print("Degree-3 P[expm1(x)] with 0+x constraint on [-log(2)/64, log(2)/64] =", horner(P));
print("  -log(2)/64 =", lb);
print("  +log(2)/64 =", ub);
print("  relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU));

lb = round(-log(2)/64, single, RN);
ub = round(log(2)/64, single, RN);
P = fpminimax(expm1(x), [|1,2|], [|SG...|], [lb; ub], relative, 0);
print("Degree-2 P[expm1(x)] with 0 constraint on [-log(2)/64, log(2)/64] =", horner(P));
print("  -log(2)/64 =", lb);
print("  +log(2)/64 =", ub);
print("  relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU));

lb = round(-log(2)/128, single, RN);
ub = round(log(2)/128, single, RN);
P = fpminimax(expm1(x), [|2,3|], [|SG...|], [lb; ub], relative, 0+x);
print("Degree-3 P[expm1(x)] with 0+x constraint on [-log(2)/128, log(2)/128] =", horner(P));
print("  -log(2)/128 =", lb);
print("  +log(2)/128 =", ub);
print("  relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU));

lb = round(-log(2)/128, single, RN);
ub = round(log(2)/128, single, RN);
P = fpminimax(expm1(x), [|1,2|], [|SG...|], [lb; ub], relative, 0);
print("Degree-2 P[expm1(x)] with 0 constraint on [-log(2)/128, log(2)/128] =", horner(P));
print("  -log(2)/128 =", lb);
print("  +log(2)/128 =", ub);
print("  relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU));

lb = round(-log(2)/256, single, RN);
ub = round(log(2)/256, single, RN);
P = fpminimax(expm1(x), [|1,2|], [|SG...|], [lb; ub], relative, 0);
print("Degree-3 P[expm1(x)] with 0 constraint on [-log(2)/256, log(2)/256] =", horner(P));
print("  -log(2)/256 =", lb);
print("  +log(2)/256 =", ub);
print("  relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU));

lb = round(-log(2)/512, single, RN);
ub = round(log(2)/512, single, RN);
P = fpminimax(expm1(x), [|1,2|], [|SG...|], [lb; ub], relative, 0);
print("Degree-3 P[expm1(x)] with 0 constraint on [-log(2)/512, log(2)/512] =", horner(P));
print("  -log(2)/512 =", lb);
print("  +log(2)/512 =", ub);
print("  relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU));

lb = round(-log(2)/512, single, RN);
ub = round(log(2)/512, single, RN);
// fpminimax fails, manually find the minimum in the vicinity of degree-3 polynomial coefficients
P = x * (1 + x * 0x1.FFFFFAp-2);
print("Degree-3 P[expm1(x)] with 0+x constraint on [-log(2)/512, log(2)/512] =", horner(P));
print("  -log(2)/512 =", lb);
print("  +log(2)/512 =", ub);
print("  relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU));

lb = round(-log(2)/1024, single, RN);
ub = round(log(2)/1024, single, RN);
P = fpminimax(expm1(x), [|1,2|], [|SG...|], [lb; ub], relative, 0);
print("Degree-3 P[expm1(x)] with 0 constraint on [-log(2)/1024, log(2)/1024] =", horner(P));
print("  -log(2)/1024 =", lb);
print("  +log(2)/1024 =", ub);
print("  relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU));

print("***IEEE FP16***");
print("log2(e) =", round(log2(exp(1)), halfprecision, RN));
minus_ln2_hi = round(-log(2), halfprecision, RN);
minus_ln2_lo = round(-log(2) - minus_ln2_hi, halfprecision, RN);
print("-log(2):hi =", minus_ln2_hi);
print("-log(2):lo =", minus_ln2_lo);

lb = round(-log(2)/2, halfprecision, RN);
ub = round(log(2)/2, halfprecision, RN);
P = fpminimax(expm1(x), [|1,2,3|], [|HP...|], [lb; ub], relative, 0);
print("Degree-3 P[expm1(x)] with 0 constraint on [-log(2)/2, log(2)/2] =", horner(P));
print("  -log(2)/2 = ", lb);
print("  +log(2)/2 =", ub);
print("  relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), halfprecision, RU));

lb = round(-log(2)/2, halfprecision, RN);
ub = round(log(2)/2, halfprecision, RN);
P = fpminimax(expm1(x), [|2,3|], [|HP...|], [lb; ub], relative, 0+x);
print("Degree-3 P[expm1(x)] with 0+x constraint on [-log(2)/2, log(2)/2] =", horner(P));
print("  -log(2)/2 = ", lb);
print("  +log(2)/2 =", ub);
print("  relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), halfprecision, RU));

lb = round(-log(2)/2, single, RN);
ub = round(log(2)/2, single, RN);
P = fpminimax(expm1(x), [|1,2,3|], [|SG...|], [lb; ub], relative, 0);
print("Degree-3 P[expm1(x)] with 0 constraint on [-log(2)/2, log(2)/2] =", horner(P));
print("  -log(2)/2 = ", lb);
print("  +log(2)/2 =", ub);
print("  relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU));

lb = round(-log(2)/2, single, RN);
ub = round(log(2)/2, single, RN);
P = fpminimax(expm1(x), [|1,2|], [|SG...|], [lb; ub], relative, 0);
print("Degree-2 P[expm1(x)] with 0 constraint on [-log(2)/2, log(2)/2] =", horner(P));
print("  -log(2)/2 = ", lb);
print("  +log(2)/2 =", ub);
print("  relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU));