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 */
fExponential(fInt exponent)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
fNaturalLog(fInt value)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
fDecodeLinearFuse(uint32_t fuse_value,fInt f_min,fInt f_range,uint32_t bitlength)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
fDecodeLogisticFuse(uint32_t fuse_value,fInt f_average,fInt f_range,uint32_t bitlength)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
fDecodeLeakageID(uint32_t leakageID_fuse,fInt ln_max_div_min,fInt f_min,uint32_t bitlength)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
ConvertToFraction(int X)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
fNegate(fInt X)232 static fInt fNegate(fInt X)
233 {
234 fInt CONSTANT_NEGONE = ConvertToFraction(-1);
235 return (fMultiply(X, CONSTANT_NEGONE));
236 }
237
Convert_ULONG_ToFraction(uint32_t X)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
GetScaledFraction(int X,int factor)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 */
fAdd(fInt X,fInt Y)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 */
fSubtract(fInt X,fInt Y)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
Equal(fInt A,fInt B)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
GreaterThan(fInt A,fInt B)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
fMultiply(fInt X,fInt Y)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
fDivide(fInt X,fInt Y)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
ConvertBackToInteger(fInt A)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
fGetSquare(fInt A)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) */
fSqrt(fInt num)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
SolveQuadracticEqn(fInt A,fInt B,fInt C,fInt Roots[])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 */
GetReal(fInt A)494 static int GetReal (fInt A)
495 {
496 return (A.full >> SHIFT_AMOUNT);
497 }
498
Divide(int X,int Y)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
uGetScaledDecimal(fInt A)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
uPow(int base,int power)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
uAbs(int X)534 static int uAbs(int X)
535 {
536 if (X < 0)
537 return (X * -1);
538 else
539 return X;
540 }
541
fRoundUpByStepSize(fInt A,fInt fStepSize,bool error_term)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