ppevvmath.h 17 KB

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