1// Copyright 2020 Google LLC 2// 3// This source code is licensed under the BSD-style license found in the 4// LICENSE file in the root directory of this source tree. 5 6display=hexadecimal !; 7print("***IEEE FP32***"); 8print("log2(e) =", round(log2(exp(1)), single, RN)); 9minus_ln2_hi = round(-log(2), single, RN); 10minus_ln2_lo = round(-log(2) - minus_ln2_hi, single, RN); 11print("-log(2):hi =", minus_ln2_hi); 12print("-log(2):lo =", minus_ln2_lo); 13 14print("log2(e) * 8 =", round(log2(exp(1)) * 8, single, RN)); 15minus_ln2_o8_hi = round(-log(2)/8, single, RN); 16minus_ln2_o8_lo = round(-log(2)/8 - minus_ln2_o8_hi, single, RN); 17print("-log(2):hi / 8 =", minus_ln2_o8_hi); 18print("-log(2):lo / 8 =", minus_ln2_o8_lo); 19 20lb = round(-log(2)/2, single, RN); 21ub = round(log(2)/2, single, RN); 22P = fpminimax(expm1(x), [|1,2,3,4,5|], [|SG...|], [lb; ub], relative, 0); 23print("Degree-5 P[expm1(x)] with 0 constraint on [-log(2)/2, log(2)/2] =", horner(P)); 24print(" -log(2)/2 = ", lb); 25print(" +log(2)/2 =", ub); 26print(" relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU)); 27 28lb = round(-log(2)/2, single, RN); 29ub = round(log(2)/2, single, RN); 30P = fpminimax(expm1(x), [|1,2,3,4,5,6|], [|SG...|], [lb; ub], relative, 0); 31print("Degree-6 P[expm1(x)] with 0 constraint on [-log(2)/2, log(2)/2] =", horner(P)); 32print(" -log(2)/2 = ", lb); 33print(" +log(2)/2 =", ub); 34print(" relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU)); 35 36lb = round(-log(2)/2, single, RN); 37ub = round(log(2)/2, single, RN); 38P = fpminimax(expm1(x), [|2,3,4,5|], [|SG...|], [lb; ub], relative, 0+x); 39print("Degree-5 P[expm1(x)] with 0+x constraint on [-log(2)/2, log(2)/2] =", horner(P)); 40print(" -log(2)/2 = ", lb); 41print(" +log(2)/2 =", ub); 42print(" relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU)); 43 44lb = round(-log(2)/2, single, RN); 45ub = round(log(2)/2, single, RN); 46P = fpminimax(expm1(x), [|2,3,4,5,6|], [|SG...|], [lb; ub], relative, 0+x); 47print("Degree-6 P[expm1(x)] with 0+x constraint on [-log(2)/2, log(2)/2] =", horner(P)); 48print(" -log(2)/2 = ", lb); 49print(" +log(2)/2 =", ub); 50print(" relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU)); 51 52lb = round(-log(2)/8, single, RN); 53ub = round(log(2)/8, single, RN); 54P = fpminimax(expm1(x), [|2,3,4|], [|SG...|], [lb; ub], relative, 0+x); 55print("Degree-4 P[expm1(x)] with 0+x constraint on [-log(2)/8, log(2)/8] =", horner(P)); 56print(" -log(2)/8 =", lb); 57print(" +log(2)/8 =", ub); 58print(" relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU)); 59 60lb = round(-log(2)/16, single, RN); 61ub = round(log(2)/16, single, RN); 62P = fpminimax(expm1(x), [|1,2,3,4|], [|SG...|], [lb; ub], relative, 0); 63print("Degree-4 P[expm1(x)] with 0 constraint on [-log(2)/16, log(2)/16] =", horner(P)); 64print(" -log(2)/16 =", lb); 65print(" +log(2)/16 =", ub); 66print(" relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU)); 67 68lb = round(-log(2)/16, single, RN); 69ub = round(log(2)/16, single, RN); 70P = fpminimax(expm1(x), [|2,3,4|], [|SG...|], [lb; ub], relative, 0+x); 71print("Degree-4 P[expm1(x)] with 0+x constraint on [-log(2)/16, log(2)/16] =", horner(P)); 72print(" -log(2)/16 =", lb); 73print(" +log(2)/16 =", ub); 74print(" relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU)); 75 76lb = round(-log(2)/16, single, RN); 77ub = round(log(2)/16, single, RN); 78P = fpminimax(expm1(x), [|1,2,3|], [|SG...|], [lb; ub], relative, 0); 79print("Degree-3 P[expm1(x)] with 0 constraint on [-log(2)/16, log(2)/16] =", horner(P)); 80print(" -log(2)/16 =", lb); 81print(" +log(2)/16 =", ub); 82print(" relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU)); 83 84lb = round(-log(2)/16, single, RN); 85ub = round(log(2)/16, single, RN); 86P = fpminimax(expm1(x), [|2,3|], [|SG...|], [lb; ub], relative, 0+x); 87print("Degree-3 P[expm1(x)] with 0+x constraint on [-log(2)/16, log(2)/16] =", horner(P)); 88print(" -log(2)/16 =", lb); 89print(" +log(2)/16 =", ub); 90print(" relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU)); 91 92lb = round(-log(2)/32, single, RN); 93ub = round(log(2)/32, single, RN); 94P = fpminimax(expm1(x), [|2,3,4|], [|SG...|], [lb; ub], relative, 0+x); 95print("Degree-4 P[expm1(x)] with 0+x constraint on [-log(2)/32, log(2)/32] =", horner(P)); 96print(" -log(2)/32 =", lb); 97print(" +log(2)/32 =", ub); 98print(" relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU)); 99 100lb = round(-log(2)/32, single, RN); 101ub = round(log(2)/32, single, RN); 102P = fpminimax(expm1(x), [|2,3|], [|SG...|], [lb; ub], relative, 0+x); 103print("Degree-3 P[expm1(x)] with 0+x constraint on [-log(2)/32, log(2)/32] =", horner(P)); 104print(" -log(2)/32 =", lb); 105print(" +log(2)/32 =", ub); 106print(" relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU)); 107 108lb = round(-log(2)/64, single, RN); 109ub = round(log(2)/64, single, RN); 110P = fpminimax(expm1(x), [|2,3|], [|SG...|], [lb; ub], relative, 0+x); 111print("Degree-3 P[expm1(x)] with 0+x constraint on [-log(2)/64, log(2)/64] =", horner(P)); 112print(" -log(2)/64 =", lb); 113print(" +log(2)/64 =", ub); 114print(" relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU)); 115 116lb = round(-log(2)/64, single, RN); 117ub = round(log(2)/64, single, RN); 118P = fpminimax(expm1(x), [|1,2|], [|SG...|], [lb; ub], relative, 0); 119print("Degree-2 P[expm1(x)] with 0 constraint on [-log(2)/64, log(2)/64] =", horner(P)); 120print(" -log(2)/64 =", lb); 121print(" +log(2)/64 =", ub); 122print(" relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU)); 123 124lb = round(-log(2)/128, single, RN); 125ub = round(log(2)/128, single, RN); 126P = fpminimax(expm1(x), [|2,3|], [|SG...|], [lb; ub], relative, 0+x); 127print("Degree-3 P[expm1(x)] with 0+x constraint on [-log(2)/128, log(2)/128] =", horner(P)); 128print(" -log(2)/128 =", lb); 129print(" +log(2)/128 =", ub); 130print(" relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU)); 131 132lb = round(-log(2)/128, single, RN); 133ub = round(log(2)/128, single, RN); 134P = fpminimax(expm1(x), [|1,2|], [|SG...|], [lb; ub], relative, 0); 135print("Degree-2 P[expm1(x)] with 0 constraint on [-log(2)/128, log(2)/128] =", horner(P)); 136print(" -log(2)/128 =", lb); 137print(" +log(2)/128 =", ub); 138print(" relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU)); 139 140lb = round(-log(2)/256, single, RN); 141ub = round(log(2)/256, single, RN); 142P = fpminimax(expm1(x), [|1,2|], [|SG...|], [lb; ub], relative, 0); 143print("Degree-3 P[expm1(x)] with 0 constraint on [-log(2)/256, log(2)/256] =", horner(P)); 144print(" -log(2)/256 =", lb); 145print(" +log(2)/256 =", ub); 146print(" relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU)); 147 148lb = round(-log(2)/512, single, RN); 149ub = round(log(2)/512, single, RN); 150P = fpminimax(expm1(x), [|1,2|], [|SG...|], [lb; ub], relative, 0); 151print("Degree-3 P[expm1(x)] with 0 constraint on [-log(2)/512, log(2)/512] =", horner(P)); 152print(" -log(2)/512 =", lb); 153print(" +log(2)/512 =", ub); 154print(" relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU)); 155 156lb = round(-log(2)/512, single, RN); 157ub = round(log(2)/512, single, RN); 158// fpminimax fails, manually find the minimum in the vicinity of degree-3 polynomial coefficients 159P = x * (1 + x * 0x1.FFFFFAp-2); 160print("Degree-3 P[expm1(x)] with 0+x constraint on [-log(2)/512, log(2)/512] =", horner(P)); 161print(" -log(2)/512 =", lb); 162print(" +log(2)/512 =", ub); 163print(" relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU)); 164 165lb = round(-log(2)/1024, single, RN); 166ub = round(log(2)/1024, single, RN); 167P = fpminimax(expm1(x), [|1,2|], [|SG...|], [lb; ub], relative, 0); 168print("Degree-3 P[expm1(x)] with 0 constraint on [-log(2)/1024, log(2)/1024] =", horner(P)); 169print(" -log(2)/1024 =", lb); 170print(" +log(2)/1024 =", ub); 171print(" relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU)); 172 173print("***IEEE FP16***"); 174print("log2(e) =", round(log2(exp(1)), halfprecision, RN)); 175minus_ln2_hi = round(-log(2), halfprecision, RN); 176minus_ln2_lo = round(-log(2) - minus_ln2_hi, halfprecision, RN); 177print("-log(2):hi =", minus_ln2_hi); 178print("-log(2):lo =", minus_ln2_lo); 179 180lb = round(-log(2)/2, halfprecision, RN); 181ub = round(log(2)/2, halfprecision, RN); 182P = fpminimax(expm1(x), [|1,2,3|], [|HP...|], [lb; ub], relative, 0); 183print("Degree-3 P[expm1(x)] with 0 constraint on [-log(2)/2, log(2)/2] =", horner(P)); 184print(" -log(2)/2 = ", lb); 185print(" +log(2)/2 =", ub); 186print(" relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), halfprecision, RU)); 187 188lb = round(-log(2)/2, halfprecision, RN); 189ub = round(log(2)/2, halfprecision, RN); 190P = fpminimax(expm1(x), [|2,3|], [|HP...|], [lb; ub], relative, 0+x); 191print("Degree-3 P[expm1(x)] with 0+x constraint on [-log(2)/2, log(2)/2] =", horner(P)); 192print(" -log(2)/2 = ", lb); 193print(" +log(2)/2 =", ub); 194print(" relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), halfprecision, RU)); 195 196lb = round(-log(2)/2, single, RN); 197ub = round(log(2)/2, single, RN); 198P = fpminimax(expm1(x), [|1,2,3|], [|SG...|], [lb; ub], relative, 0); 199print("Degree-3 P[expm1(x)] with 0 constraint on [-log(2)/2, log(2)/2] =", horner(P)); 200print(" -log(2)/2 = ", lb); 201print(" +log(2)/2 =", ub); 202print(" relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU)); 203 204lb = round(-log(2)/2, single, RN); 205ub = round(log(2)/2, single, RN); 206P = fpminimax(expm1(x), [|1,2|], [|SG...|], [lb; ub], relative, 0); 207print("Degree-2 P[expm1(x)] with 0 constraint on [-log(2)/2, log(2)/2] =", horner(P)); 208print(" -log(2)/2 = ", lb); 209print(" +log(2)/2 =", ub); 210print(" relative error =", round(dirtyinfnorm(P / expm1(x) - 1, [lb; ub]), single, RU)); 211