1 /* $NetBSD: ppevvmath.h,v 1.3 2021/12/19 12:21:30 riastradh Exp $ */ 2 3 /* 4 * Copyright 2015 Advanced Micro Devices, Inc. 5 * 6 * Permission is hereby granted, free of charge, to any person obtaining a 7 * copy of this software and associated documentation files (the "Software"), 8 * to deal in the Software without restriction, including without limitation 9 * the rights to use, copy, modify, merge, publish, distribute, sublicense, 10 * and/or sell copies of the Software, and to permit persons to whom the 11 * Software is furnished to do so, subject to the following conditions: 12 * 13 * The above copyright notice and this permission notice shall be included in 14 * all copies or substantial portions of the Software. 15 * 16 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 17 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 18 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 19 * THE COPYRIGHT HOLDER(S) OR AUTHOR(S) BE LIABLE FOR ANY CLAIM, DAMAGES OR 20 * OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 21 * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR 22 * OTHER DEALINGS IN THE SOFTWARE. 23 * 24 */ 25 #include <asm/div64.h> 26 27 #define SHIFT_AMOUNT 16 /* We multiply all original integers with 2^SHIFT_AMOUNT to get the fInt representation */ 28 29 #define PRECISION 5 /* Change this value to change the number of decimal places in the final output - 5 is a good default */ 30 31 #define SHIFTED_2 (2 << SHIFT_AMOUNT) 32 #define MAX_VALUE (1 << (SHIFT_AMOUNT - 1)) - 1 /* 32767 - Might change in the future */ 33 34 /* ------------------------------------------------------------------------------- 35 * NEW TYPE - fINT 36 * ------------------------------------------------------------------------------- 37 * A variable of type fInt can be accessed in 3 ways using the dot (.) operator 38 * fInt A; 39 * A.full => The full number as it is. Generally not easy to read 40 * A.partial.real => Only the integer portion 41 * A.partial.decimal => Only the fractional portion 42 */ 43 typedef union _fInt { 44 int full; 45 struct _partial { 46 unsigned int decimal: SHIFT_AMOUNT; /*Needs to always be unsigned*/ 47 int real: 32 - SHIFT_AMOUNT; 48 } partial; 49 } fInt; 50 51 /* ------------------------------------------------------------------------------- 52 * Function Declarations 53 * ------------------------------------------------------------------------------- 54 */ 55 static fInt ConvertToFraction(int); /* Use this to convert an INT to a FINT */ 56 static fInt Convert_ULONG_ToFraction(uint32_t); /* Use this to convert an uint32_t to a FINT */ 57 static fInt GetScaledFraction(int, int); /* Use this to convert an INT to a FINT after scaling it by a factor */ 58 static int ConvertBackToInteger(fInt); /* Convert a FINT back to an INT that is scaled by 1000 (i.e. last 3 digits are the decimal digits) */ 59 60 static fInt fNegate(fInt); /* Returns -1 * input fInt value */ 61 static fInt fAdd (fInt, fInt); /* Returns the sum of two fInt numbers */ 62 static fInt fSubtract (fInt A, fInt B); /* Returns A-B - Sometimes easier than Adding negative numbers */ 63 static fInt fMultiply (fInt, fInt); /* Returns the product of two fInt numbers */ 64 static fInt fDivide (fInt A, fInt B); /* Returns A/B */ 65 static fInt fGetSquare(fInt); /* Returns the square of a fInt number */ 66 static fInt fSqrt(fInt); /* Returns the Square Root of a fInt number */ 67 68 static int uAbs(int); /* Returns the Absolute value of the Int */ 69 static int uPow(int base, int exponent); /* Returns base^exponent an INT */ 70 71 static void SolveQuadracticEqn(fInt, fInt, fInt, fInt[]); /* Returns the 2 roots via the array */ 72 static bool Equal(fInt, fInt); /* Returns true if two fInts are equal to each other */ 73 static bool GreaterThan(fInt A, fInt B); /* Returns true if A > B */ 74 75 static fInt fExponential(fInt exponent); /* Can be used to calculate e^exponent */ 76 static fInt fNaturalLog(fInt value); /* Can be used to calculate ln(value) */ 77 78 /* Fuse decoding functions 79 * ------------------------------------------------------------------------------------- 80 */ 81 static fInt fDecodeLinearFuse(uint32_t fuse_value, fInt f_min, fInt f_range, uint32_t bitlength); 82 static fInt fDecodeLogisticFuse(uint32_t fuse_value, fInt f_average, fInt f_range, uint32_t bitlength); 83 static fInt fDecodeLeakageID (uint32_t leakageID_fuse, fInt ln_max_div_min, fInt f_min, uint32_t bitlength); 84 85 /* Internal Support Functions - Use these ONLY for testing or adding to internal functions 86 * ------------------------------------------------------------------------------------- 87 * Some of the following functions take two INTs as their input - This is unsafe for a variety of reasons. 88 */ 89 static fInt Divide (int, int); /* Divide two INTs and return result as FINT */ 90 static fInt fNegate(fInt); 91 92 static int uGetScaledDecimal (fInt); /* Internal function */ 93 static int GetReal (fInt A); /* Internal function */ 94 95 /* ------------------------------------------------------------------------------------- 96 * TROUBLESHOOTING INFORMATION 97 * ------------------------------------------------------------------------------------- 98 * 1) ConvertToFraction - InputOutOfRangeException: Only accepts numbers smaller than MAX_VALUE (default: 32767) 99 * 2) fAdd - OutputOutOfRangeException: Output bigger than MAX_VALUE (default: 32767) 100 * 3) fMultiply - OutputOutOfRangeException: 101 * 4) fGetSquare - OutputOutOfRangeException: 102 * 5) fDivide - DivideByZeroException 103 * 6) fSqrt - NegativeSquareRootException: Input cannot be a negative number 104 */ 105 106 /* ------------------------------------------------------------------------------------- 107 * START OF CODE 108 * ------------------------------------------------------------------------------------- 109 */ 110 static fInt fExponential(fInt exponent) /*Can be used to calculate e^exponent*/ 111 { 112 uint32_t i; 113 bool bNegated = false; 114 115 fInt fPositiveOne = ConvertToFraction(1); 116 fInt fZERO = ConvertToFraction(0); 117 118 fInt lower_bound = Divide(78, 10000); 119 fInt solution = fPositiveOne; /*Starting off with baseline of 1 */ 120 fInt error_term; 121 122 static const uint32_t k_array[11] = {55452, 27726, 13863, 6931, 4055, 2231, 1178, 606, 308, 155, 78}; 123 static const uint32_t expk_array[11] = {2560000, 160000, 40000, 20000, 15000, 12500, 11250, 10625, 10313, 10156, 10078}; 124 125 if (GreaterThan(fZERO, exponent)) { 126 exponent = fNegate(exponent); 127 bNegated = true; 128 } 129 130 while (GreaterThan(exponent, lower_bound)) { 131 for (i = 0; i < 11; i++) { 132 if (GreaterThan(exponent, GetScaledFraction(k_array[i], 10000))) { 133 exponent = fSubtract(exponent, GetScaledFraction(k_array[i], 10000)); 134 solution = fMultiply(solution, GetScaledFraction(expk_array[i], 10000)); 135 } 136 } 137 } 138 139 error_term = fAdd(fPositiveOne, exponent); 140 141 solution = fMultiply(solution, error_term); 142 143 if (bNegated) 144 solution = fDivide(fPositiveOne, solution); 145 146 return solution; 147 } 148 149 static fInt fNaturalLog(fInt value) 150 { 151 uint32_t i; 152 fInt upper_bound = Divide(8, 1000); 153 fInt fNegativeOne = ConvertToFraction(-1); 154 fInt solution = ConvertToFraction(0); /*Starting off with baseline of 0 */ 155 fInt error_term; 156 157 static const uint32_t k_array[10] = {160000, 40000, 20000, 15000, 12500, 11250, 10625, 10313, 10156, 10078}; 158 static const uint32_t logk_array[10] = {27726, 13863, 6931, 4055, 2231, 1178, 606, 308, 155, 78}; 159 160 while (GreaterThan(fAdd(value, fNegativeOne), upper_bound)) { 161 for (i = 0; i < 10; i++) { 162 if (GreaterThan(value, GetScaledFraction(k_array[i], 10000))) { 163 value = fDivide(value, GetScaledFraction(k_array[i], 10000)); 164 solution = fAdd(solution, GetScaledFraction(logk_array[i], 10000)); 165 } 166 } 167 } 168 169 error_term = fAdd(fNegativeOne, value); 170 171 return (fAdd(solution, error_term)); 172 } 173 174 static fInt fDecodeLinearFuse(uint32_t fuse_value, fInt f_min, fInt f_range, uint32_t bitlength) 175 { 176 fInt f_fuse_value = Convert_ULONG_ToFraction(fuse_value); 177 fInt f_bit_max_value = Convert_ULONG_ToFraction((uPow(2, bitlength)) - 1); 178 179 fInt f_decoded_value; 180 181 f_decoded_value = fDivide(f_fuse_value, f_bit_max_value); 182 f_decoded_value = fMultiply(f_decoded_value, f_range); 183 f_decoded_value = fAdd(f_decoded_value, f_min); 184 185 return f_decoded_value; 186 } 187 188 189 static fInt fDecodeLogisticFuse(uint32_t fuse_value, fInt f_average, fInt f_range, uint32_t bitlength) 190 { 191 fInt f_fuse_value = Convert_ULONG_ToFraction(fuse_value); 192 fInt f_bit_max_value = Convert_ULONG_ToFraction((uPow(2, bitlength)) - 1); 193 194 fInt f_CONSTANT_NEG13 = ConvertToFraction(-13); 195 fInt f_CONSTANT1 = ConvertToFraction(1); 196 197 fInt f_decoded_value; 198 199 f_decoded_value = fSubtract(fDivide(f_bit_max_value, f_fuse_value), f_CONSTANT1); 200 f_decoded_value = fNaturalLog(f_decoded_value); 201 f_decoded_value = fMultiply(f_decoded_value, fDivide(f_range, f_CONSTANT_NEG13)); 202 f_decoded_value = fAdd(f_decoded_value, f_average); 203 204 return f_decoded_value; 205 } 206 207 static fInt fDecodeLeakageID (uint32_t leakageID_fuse, fInt ln_max_div_min, fInt f_min, uint32_t bitlength) 208 { 209 fInt fLeakage; 210 fInt f_bit_max_value = Convert_ULONG_ToFraction((uPow(2, bitlength)) - 1); 211 212 fLeakage = fMultiply(ln_max_div_min, Convert_ULONG_ToFraction(leakageID_fuse)); 213 fLeakage = fDivide(fLeakage, f_bit_max_value); 214 fLeakage = fExponential(fLeakage); 215 fLeakage = fMultiply(fLeakage, f_min); 216 217 return fLeakage; 218 } 219 220 static fInt ConvertToFraction(int X) /*Add all range checking here. Is it possible to make fInt a private declaration? */ 221 { 222 fInt temp; 223 224 if (X <= MAX_VALUE) 225 temp.full = (X << SHIFT_AMOUNT); 226 else 227 temp.full = 0; 228 229 return temp; 230 } 231 232 static fInt fNegate(fInt X) 233 { 234 fInt CONSTANT_NEGONE = ConvertToFraction(-1); 235 return (fMultiply(X, CONSTANT_NEGONE)); 236 } 237 238 static fInt Convert_ULONG_ToFraction(uint32_t X) 239 { 240 fInt temp; 241 242 if (X <= MAX_VALUE) 243 temp.full = (X << SHIFT_AMOUNT); 244 else 245 temp.full = 0; 246 247 return temp; 248 } 249 250 static fInt GetScaledFraction(int X, int factor) 251 { 252 int times_shifted, factor_shifted; 253 bool bNEGATED; 254 fInt fValue; 255 256 times_shifted = 0; 257 factor_shifted = 0; 258 bNEGATED = false; 259 260 if (X < 0) { 261 X = -1*X; 262 bNEGATED = true; 263 } 264 265 if (factor < 0) { 266 factor = -1*factor; 267 bNEGATED = !bNEGATED; /*If bNEGATED = true due to X < 0, this will cover the case of negative cancelling negative */ 268 } 269 270 if ((X > MAX_VALUE) || factor > MAX_VALUE) { 271 if ((X/factor) <= MAX_VALUE) { 272 while (X > MAX_VALUE) { 273 X = X >> 1; 274 times_shifted++; 275 } 276 277 while (factor > MAX_VALUE) { 278 factor = factor >> 1; 279 factor_shifted++; 280 } 281 } else { 282 fValue.full = 0; 283 return fValue; 284 } 285 } 286 287 if (factor == 1) 288 return ConvertToFraction(X); 289 290 fValue = fDivide(ConvertToFraction(X * uPow(-1, bNEGATED)), ConvertToFraction(factor)); 291 292 fValue.full = fValue.full << times_shifted; 293 fValue.full = fValue.full >> factor_shifted; 294 295 return fValue; 296 } 297 298 /* Addition using two fInts */ 299 static fInt fAdd (fInt X, fInt Y) 300 { 301 fInt Sum; 302 303 Sum.full = X.full + Y.full; 304 305 return Sum; 306 } 307 308 /* Addition using two fInts */ 309 static fInt fSubtract (fInt X, fInt Y) 310 { 311 fInt Difference; 312 313 Difference.full = X.full - Y.full; 314 315 return Difference; 316 } 317 318 static bool Equal(fInt A, fInt B) 319 { 320 if (A.full == B.full) 321 return true; 322 else 323 return false; 324 } 325 326 static bool GreaterThan(fInt A, fInt B) 327 { 328 if (A.full > B.full) 329 return true; 330 else 331 return false; 332 } 333 334 static fInt fMultiply (fInt X, fInt Y) /* Uses 64-bit integers (int64_t) */ 335 { 336 fInt Product; 337 int64_t tempProduct; 338 bool X_LessThanOne __unused, Y_LessThanOne __unused; 339 340 X_LessThanOne = (X.partial.real == 0 && X.partial.decimal != 0 && X.full >= 0); 341 Y_LessThanOne = (Y.partial.real == 0 && Y.partial.decimal != 0 && Y.full >= 0); 342 343 /*The following is for a very specific common case: Non-zero number with ONLY fractional portion*/ 344 /* TEMPORARILY DISABLED - CAN BE USED TO IMPROVE PRECISION 345 346 if (X_LessThanOne && Y_LessThanOne) { 347 Product.full = X.full * Y.full; 348 return Product 349 }*/ 350 351 tempProduct = ((int64_t)X.full) * ((int64_t)Y.full); /*Q(16,16)*Q(16,16) = Q(32, 32) - Might become a negative number! */ 352 tempProduct = tempProduct >> 16; /*Remove lagging 16 bits - Will lose some precision from decimal; */ 353 Product.full = (int)tempProduct; /*The int64_t will lose the leading 16 bits that were part of the integer portion */ 354 355 return Product; 356 } 357 358 static fInt fDivide (fInt X, fInt Y) 359 { 360 fInt fZERO, fQuotient; 361 int64_t longlongX, longlongY; 362 363 fZERO = ConvertToFraction(0); 364 365 if (Equal(Y, fZERO)) 366 return fZERO; 367 368 longlongX = (int64_t)X.full; 369 longlongY = (int64_t)Y.full; 370 371 longlongX = longlongX << 16; /*Q(16,16) -> Q(32,32) */ 372 373 div64_s64(longlongX, longlongY); /*Q(32,32) divided by Q(16,16) = Q(16,16) Back to original format */ 374 375 fQuotient.full = (int)longlongX; 376 return fQuotient; 377 } 378 379 static int ConvertBackToInteger (fInt A) /*THIS is the function that will be used to check with the Golden settings table*/ 380 { 381 fInt fullNumber, scaledDecimal, scaledReal; 382 383 scaledReal.full = GetReal(A) * uPow(10, PRECISION-1); /* DOUBLE CHECK THISSSS!!! */ 384 385 scaledDecimal.full = uGetScaledDecimal(A); 386 387 fullNumber = fAdd(scaledDecimal,scaledReal); 388 389 return fullNumber.full; 390 } 391 392 static fInt fGetSquare(fInt A) 393 { 394 return fMultiply(A,A); 395 } 396 397 /* x_new = x_old - (x_old^2 - C) / (2 * x_old) */ 398 static fInt fSqrt(fInt num) 399 { 400 fInt F_divide_Fprime, Fprime; 401 fInt test; 402 fInt twoShifted; 403 int seed, counter, error; 404 fInt x_new, x_old, C, y; 405 406 fInt fZERO = ConvertToFraction(0); 407 408 /* (0 > num) is the same as (num < 0), i.e., num is negative */ 409 410 if (GreaterThan(fZERO, num) || Equal(fZERO, num)) 411 return fZERO; 412 413 C = num; 414 415 if (num.partial.real > 3000) 416 seed = 60; 417 else if (num.partial.real > 1000) 418 seed = 30; 419 else if (num.partial.real > 100) 420 seed = 10; 421 else 422 seed = 2; 423 424 counter = 0; 425 426 if (Equal(num, fZERO)) /*Square Root of Zero is zero */ 427 return fZERO; 428 429 twoShifted = ConvertToFraction(2); 430 x_new = ConvertToFraction(seed); 431 432 do { 433 counter++; 434 435 x_old.full = x_new.full; 436 437 test = fGetSquare(x_old); /*1.75*1.75 is reverting back to 1 when shifted down */ 438 y = fSubtract(test, C); /*y = f(x) = x^2 - C; */ 439 440 Fprime = fMultiply(twoShifted, x_old); 441 F_divide_Fprime = fDivide(y, Fprime); 442 443 x_new = fSubtract(x_old, F_divide_Fprime); 444 445 error = ConvertBackToInteger(x_new) - ConvertBackToInteger(x_old); 446 447 if (counter > 20) /*20 is already way too many iterations. If we dont have an answer by then, we never will*/ 448 return x_new; 449 450 } while (uAbs(error) > 0); 451 452 return (x_new); 453 } 454 455 static void SolveQuadracticEqn(fInt A, fInt B, fInt C, fInt Roots[]) 456 { 457 fInt *pRoots = &Roots[0]; 458 fInt temp, root_first, root_second; 459 fInt f_CONSTANT10, f_CONSTANT100; 460 461 f_CONSTANT100 = ConvertToFraction(100); 462 f_CONSTANT10 = ConvertToFraction(10); 463 464 while(GreaterThan(A, f_CONSTANT100) || GreaterThan(B, f_CONSTANT100) || GreaterThan(C, f_CONSTANT100)) { 465 A = fDivide(A, f_CONSTANT10); 466 B = fDivide(B, f_CONSTANT10); 467 C = fDivide(C, f_CONSTANT10); 468 } 469 470 temp = fMultiply(ConvertToFraction(4), A); /* root = 4*A */ 471 temp = fMultiply(temp, C); /* root = 4*A*C */ 472 temp = fSubtract(fGetSquare(B), temp); /* root = b^2 - 4AC */ 473 temp = fSqrt(temp); /*root = Sqrt (b^2 - 4AC); */ 474 475 root_first = fSubtract(fNegate(B), temp); /* b - Sqrt(b^2 - 4AC) */ 476 root_second = fAdd(fNegate(B), temp); /* b + Sqrt(b^2 - 4AC) */ 477 478 root_first = fDivide(root_first, ConvertToFraction(2)); /* [b +- Sqrt(b^2 - 4AC)]/[2] */ 479 root_first = fDivide(root_first, A); /*[b +- Sqrt(b^2 - 4AC)]/[2*A] */ 480 481 root_second = fDivide(root_second, ConvertToFraction(2)); /* [b +- Sqrt(b^2 - 4AC)]/[2] */ 482 root_second = fDivide(root_second, A); /*[b +- Sqrt(b^2 - 4AC)]/[2*A] */ 483 484 *(pRoots + 0) = root_first; 485 *(pRoots + 1) = root_second; 486 } 487 488 /* ----------------------------------------------------------------------------- 489 * SUPPORT FUNCTIONS 490 * ----------------------------------------------------------------------------- 491 */ 492 493 /* Conversion Functions */ 494 static int GetReal (fInt A) 495 { 496 return (A.full >> SHIFT_AMOUNT); 497 } 498 499 static fInt Divide (int X, int Y) 500 { 501 fInt A, B, Quotient; 502 503 A.full = X << SHIFT_AMOUNT; 504 B.full = Y << SHIFT_AMOUNT; 505 506 Quotient = fDivide(A, B); 507 508 return Quotient; 509 } 510 511 static int uGetScaledDecimal (fInt A) /*Converts the fractional portion to whole integers - Costly function */ 512 { 513 int dec[PRECISION]; 514 int i, scaledDecimal = 0, tmp = A.partial.decimal; 515 516 for (i = 0; i < PRECISION; i++) { 517 dec[i] = tmp / (1 << SHIFT_AMOUNT); 518 tmp = tmp - ((1 << SHIFT_AMOUNT)*dec[i]); 519 tmp *= 10; 520 scaledDecimal = scaledDecimal + dec[i]*uPow(10, PRECISION - 1 -i); 521 } 522 523 return scaledDecimal; 524 } 525 526 static int uPow(int base, int power) 527 { 528 if (power == 0) 529 return 1; 530 else 531 return (base)*uPow(base, power - 1); 532 } 533 534 static int uAbs(int X) 535 { 536 if (X < 0) 537 return (X * -1); 538 else 539 return X; 540 } 541 542 static fInt fRoundUpByStepSize(fInt A, fInt fStepSize, bool error_term) 543 { 544 fInt solution; 545 546 solution = fDivide(A, fStepSize); 547 solution.partial.decimal = 0; /*All fractional digits changes to 0 */ 548 549 if (error_term) 550 solution.partial.real += 1; /*Error term of 1 added */ 551 552 solution = fMultiply(solution, fStepSize); 553 solution = fAdd(solution, fStepSize); 554 555 return solution; 556 } 557 558