blob: b7429a5278289f35d7a40b76063add4ad916b34f [file] [log] [blame]
Eric Huang770911a2015-11-09 17:34:31 -05001/*
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
25#define SHIFT_AMOUNT 16 /* We multiply all original integers with 2^SHIFT_AMOUNT to get the fInt representation */
26
27#define PRECISION 5 /* Change this value to change the number of decimal places in the final output - 5 is a good default */
28
29#define SHIFTED_2 (2 << SHIFT_AMOUNT)
30#define MAX (1 << (SHIFT_AMOUNT - 1)) - 1 /* 32767 - Might change in the future */
31
32/* -------------------------------------------------------------------------------
33 * NEW TYPE - fINT
34 * -------------------------------------------------------------------------------
35 * A variable of type fInt can be accessed in 3 ways using the dot (.) operator
36 * fInt A;
37 * A.full => The full number as it is. Generally not easy to read
38 * A.partial.real => Only the integer portion
39 * A.partial.decimal => Only the fractional portion
40 */
41typedef union _fInt {
42 int full;
43 struct _partial {
44 unsigned int decimal: SHIFT_AMOUNT; /*Needs to always be unsigned*/
45 int real: 32 - SHIFT_AMOUNT;
46 } partial;
47} fInt;
48
49/* -------------------------------------------------------------------------------
50 * Function Declarations
51 * -------------------------------------------------------------------------------
52 */
53fInt ConvertToFraction(int); /* Use this to convert an INT to a FINT */
54fInt Convert_ULONG_ToFraction(uint32_t); /* Use this to convert an uint32_t to a FINT */
55fInt GetScaledFraction(int, int); /* Use this to convert an INT to a FINT after scaling it by a factor */
56int ConvertBackToInteger(fInt); /* Convert a FINT back to an INT that is scaled by 1000 (i.e. last 3 digits are the decimal digits) */
57
58fInt fNegate(fInt); /* Returns -1 * input fInt value */
59fInt fAdd (fInt, fInt); /* Returns the sum of two fInt numbers */
60fInt fSubtract (fInt A, fInt B); /* Returns A-B - Sometimes easier than Adding negative numbers */
61fInt fMultiply (fInt, fInt); /* Returns the product of two fInt numbers */
62fInt fDivide (fInt A, fInt B); /* Returns A/B */
63fInt fGetSquare(fInt); /* Returns the square of a fInt number */
64fInt fSqrt(fInt); /* Returns the Square Root of a fInt number */
65
66int uAbs(int); /* Returns the Absolute value of the Int */
67fInt fAbs(fInt); /* Returns the Absolute value of the fInt */
68int uPow(int base, int exponent); /* Returns base^exponent an INT */
69
70void SolveQuadracticEqn(fInt, fInt, fInt, fInt[]); /* Returns the 2 roots via the array */
71bool Equal(fInt, fInt); /* Returns true if two fInts are equal to each other */
72bool GreaterThan(fInt A, fInt B); /* Returns true if A > B */
73
74fInt fExponential(fInt exponent); /* Can be used to calculate e^exponent */
75fInt fNaturalLog(fInt value); /* Can be used to calculate ln(value) */
76
77/* Fuse decoding functions
78 * -------------------------------------------------------------------------------------
79 */
80fInt fDecodeLinearFuse(uint32_t fuse_value, fInt f_min, fInt f_range, uint32_t bitlength);
81fInt fDecodeLogisticFuse(uint32_t fuse_value, fInt f_average, fInt f_range, uint32_t bitlength);
82fInt fDecodeLeakageID (uint32_t leakageID_fuse, fInt ln_max_div_min, fInt f_min, uint32_t bitlength);
83
84/* Internal Support Functions - Use these ONLY for testing or adding to internal functions
85 * -------------------------------------------------------------------------------------
86 * Some of the following functions take two INTs as their input - This is unsafe for a variety of reasons.
87 */
88fInt Add (int, int); /* Add two INTs and return Sum as FINT */
89fInt Multiply (int, int); /* Multiply two INTs and return Product as FINT */
90fInt Divide (int, int); /* You get the idea... */
91fInt fNegate(fInt);
92
93int uGetScaledDecimal (fInt); /* Internal function */
94int GetReal (fInt A); /* Internal function */
95
96/* Future Additions and Incomplete Functions
97 * -------------------------------------------------------------------------------------
98 */
99int GetRoundedValue(fInt); /* Incomplete function - Useful only when Precision is lacking */
100 /* Let us say we have 2.126 but can only handle 2 decimal points. We could */
101 /* either chop of 6 and keep 2.12 or use this function to get 2.13, which is more accurate */
102
103/* -------------------------------------------------------------------------------------
104 * TROUBLESHOOTING INFORMATION
105 * -------------------------------------------------------------------------------------
106 * 1) ConvertToFraction - InputOutOfRangeException: Only accepts numbers smaller than MAX (default: 32767)
107 * 2) fAdd - OutputOutOfRangeException: Output bigger than MAX (default: 32767)
108 * 3) fMultiply - OutputOutOfRangeException:
109 * 4) fGetSquare - OutputOutOfRangeException:
110 * 5) fDivide - DivideByZeroException
111 * 6) fSqrt - NegativeSquareRootException: Input cannot be a negative number
112 */
113
114/* -------------------------------------------------------------------------------------
115 * START OF CODE
116 * -------------------------------------------------------------------------------------
117 */
118fInt fExponential(fInt exponent) /*Can be used to calculate e^exponent*/
119{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800120 uint32_t i;
121 bool bNegated = false;
Eric Huang770911a2015-11-09 17:34:31 -0500122
Rex Zhu75ac63d2016-01-06 16:38:48 +0800123 fInt fPositiveOne = ConvertToFraction(1);
124 fInt fZERO = ConvertToFraction(0);
Eric Huang770911a2015-11-09 17:34:31 -0500125
Rex Zhu75ac63d2016-01-06 16:38:48 +0800126 fInt lower_bound = Divide(78, 10000);
127 fInt solution = fPositiveOne; /*Starting off with baseline of 1 */
128 fInt error_term;
Eric Huang770911a2015-11-09 17:34:31 -0500129
Rex Zhu75ac63d2016-01-06 16:38:48 +0800130 uint32_t k_array[11] = {55452, 27726, 13863, 6931, 4055, 2231, 1178, 606, 308, 155, 78};
131 uint32_t expk_array[11] = {2560000, 160000, 40000, 20000, 15000, 12500, 11250, 10625, 10313, 10156, 10078};
Eric Huang770911a2015-11-09 17:34:31 -0500132
Rex Zhu75ac63d2016-01-06 16:38:48 +0800133 if (GreaterThan(fZERO, exponent)) {
134 exponent = fNegate(exponent);
135 bNegated = true;
136 }
Eric Huang770911a2015-11-09 17:34:31 -0500137
Rex Zhu75ac63d2016-01-06 16:38:48 +0800138 while (GreaterThan(exponent, lower_bound)) {
139 for (i = 0; i < 11; i++) {
140 if (GreaterThan(exponent, GetScaledFraction(k_array[i], 10000))) {
141 exponent = fSubtract(exponent, GetScaledFraction(k_array[i], 10000));
142 solution = fMultiply(solution, GetScaledFraction(expk_array[i], 10000));
143 }
144 }
145 }
Eric Huang770911a2015-11-09 17:34:31 -0500146
Rex Zhu75ac63d2016-01-06 16:38:48 +0800147 error_term = fAdd(fPositiveOne, exponent);
Eric Huang770911a2015-11-09 17:34:31 -0500148
Rex Zhu75ac63d2016-01-06 16:38:48 +0800149 solution = fMultiply(solution, error_term);
Eric Huang770911a2015-11-09 17:34:31 -0500150
Rex Zhu75ac63d2016-01-06 16:38:48 +0800151 if (bNegated)
152 solution = fDivide(fPositiveOne, solution);
Eric Huang770911a2015-11-09 17:34:31 -0500153
Rex Zhu75ac63d2016-01-06 16:38:48 +0800154 return solution;
Eric Huang770911a2015-11-09 17:34:31 -0500155}
156
157fInt fNaturalLog(fInt value)
158{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800159 uint32_t i;
160 fInt upper_bound = Divide(8, 1000);
161 fInt fNegativeOne = ConvertToFraction(-1);
162 fInt solution = ConvertToFraction(0); /*Starting off with baseline of 0 */
163 fInt error_term;
Eric Huang770911a2015-11-09 17:34:31 -0500164
Rex Zhu75ac63d2016-01-06 16:38:48 +0800165 uint32_t k_array[10] = {160000, 40000, 20000, 15000, 12500, 11250, 10625, 10313, 10156, 10078};
166 uint32_t logk_array[10] = {27726, 13863, 6931, 4055, 2231, 1178, 606, 308, 155, 78};
Eric Huang770911a2015-11-09 17:34:31 -0500167
Rex Zhu75ac63d2016-01-06 16:38:48 +0800168 while (GreaterThan(fAdd(value, fNegativeOne), upper_bound)) {
169 for (i = 0; i < 10; i++) {
170 if (GreaterThan(value, GetScaledFraction(k_array[i], 10000))) {
171 value = fDivide(value, GetScaledFraction(k_array[i], 10000));
172 solution = fAdd(solution, GetScaledFraction(logk_array[i], 10000));
173 }
174 }
175 }
Eric Huang770911a2015-11-09 17:34:31 -0500176
Rex Zhu75ac63d2016-01-06 16:38:48 +0800177 error_term = fAdd(fNegativeOne, value);
Eric Huang770911a2015-11-09 17:34:31 -0500178
Rex Zhu75ac63d2016-01-06 16:38:48 +0800179 return (fAdd(solution, error_term));
Eric Huang770911a2015-11-09 17:34:31 -0500180}
181
182fInt fDecodeLinearFuse(uint32_t fuse_value, fInt f_min, fInt f_range, uint32_t bitlength)
183{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800184 fInt f_fuse_value = Convert_ULONG_ToFraction(fuse_value);
185 fInt f_bit_max_value = Convert_ULONG_ToFraction((uPow(2, bitlength)) - 1);
Eric Huang770911a2015-11-09 17:34:31 -0500186
Rex Zhu75ac63d2016-01-06 16:38:48 +0800187 fInt f_decoded_value;
Eric Huang770911a2015-11-09 17:34:31 -0500188
Rex Zhu75ac63d2016-01-06 16:38:48 +0800189 f_decoded_value = fDivide(f_fuse_value, f_bit_max_value);
190 f_decoded_value = fMultiply(f_decoded_value, f_range);
191 f_decoded_value = fAdd(f_decoded_value, f_min);
Eric Huang770911a2015-11-09 17:34:31 -0500192
Rex Zhu75ac63d2016-01-06 16:38:48 +0800193 return f_decoded_value;
Eric Huang770911a2015-11-09 17:34:31 -0500194}
195
196
197fInt fDecodeLogisticFuse(uint32_t fuse_value, fInt f_average, fInt f_range, uint32_t bitlength)
198{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800199 fInt f_fuse_value = Convert_ULONG_ToFraction(fuse_value);
200 fInt f_bit_max_value = Convert_ULONG_ToFraction((uPow(2, bitlength)) - 1);
Eric Huang770911a2015-11-09 17:34:31 -0500201
Rex Zhu75ac63d2016-01-06 16:38:48 +0800202 fInt f_CONSTANT_NEG13 = ConvertToFraction(-13);
203 fInt f_CONSTANT1 = ConvertToFraction(1);
Eric Huang770911a2015-11-09 17:34:31 -0500204
Rex Zhu75ac63d2016-01-06 16:38:48 +0800205 fInt f_decoded_value;
Eric Huang770911a2015-11-09 17:34:31 -0500206
Rex Zhu75ac63d2016-01-06 16:38:48 +0800207 f_decoded_value = fSubtract(fDivide(f_bit_max_value, f_fuse_value), f_CONSTANT1);
208 f_decoded_value = fNaturalLog(f_decoded_value);
209 f_decoded_value = fMultiply(f_decoded_value, fDivide(f_range, f_CONSTANT_NEG13));
210 f_decoded_value = fAdd(f_decoded_value, f_average);
Eric Huang770911a2015-11-09 17:34:31 -0500211
Rex Zhu75ac63d2016-01-06 16:38:48 +0800212 return f_decoded_value;
Eric Huang770911a2015-11-09 17:34:31 -0500213}
214
215fInt fDecodeLeakageID (uint32_t leakageID_fuse, fInt ln_max_div_min, fInt f_min, uint32_t bitlength)
216{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800217 fInt fLeakage;
218 fInt f_bit_max_value = Convert_ULONG_ToFraction((uPow(2, bitlength)) - 1);
Eric Huang770911a2015-11-09 17:34:31 -0500219
Rex Zhu75ac63d2016-01-06 16:38:48 +0800220 fLeakage = fMultiply(ln_max_div_min, Convert_ULONG_ToFraction(leakageID_fuse));
221 fLeakage = fDivide(fLeakage, f_bit_max_value);
222 fLeakage = fExponential(fLeakage);
223 fLeakage = fMultiply(fLeakage, f_min);
Eric Huang770911a2015-11-09 17:34:31 -0500224
Rex Zhu75ac63d2016-01-06 16:38:48 +0800225 return fLeakage;
Eric Huang770911a2015-11-09 17:34:31 -0500226}
227
228fInt ConvertToFraction(int X) /*Add all range checking here. Is it possible to make fInt a private declaration? */
229{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800230 fInt temp;
Eric Huang770911a2015-11-09 17:34:31 -0500231
Rex Zhu75ac63d2016-01-06 16:38:48 +0800232 if (X <= MAX)
233 temp.full = (X << SHIFT_AMOUNT);
234 else
235 temp.full = 0;
Eric Huang770911a2015-11-09 17:34:31 -0500236
Rex Zhu75ac63d2016-01-06 16:38:48 +0800237 return temp;
Eric Huang770911a2015-11-09 17:34:31 -0500238}
239
240fInt fNegate(fInt X)
241{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800242 fInt CONSTANT_NEGONE = ConvertToFraction(-1);
243 return (fMultiply(X, CONSTANT_NEGONE));
Eric Huang770911a2015-11-09 17:34:31 -0500244}
245
246fInt Convert_ULONG_ToFraction(uint32_t X)
247{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800248 fInt temp;
Eric Huang770911a2015-11-09 17:34:31 -0500249
Rex Zhu75ac63d2016-01-06 16:38:48 +0800250 if (X <= MAX)
251 temp.full = (X << SHIFT_AMOUNT);
252 else
253 temp.full = 0;
Eric Huang770911a2015-11-09 17:34:31 -0500254
Rex Zhu75ac63d2016-01-06 16:38:48 +0800255 return temp;
Eric Huang770911a2015-11-09 17:34:31 -0500256}
257
258fInt GetScaledFraction(int X, int factor)
259{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800260 int times_shifted, factor_shifted;
261 bool bNEGATED;
262 fInt fValue;
Eric Huang770911a2015-11-09 17:34:31 -0500263
Rex Zhu75ac63d2016-01-06 16:38:48 +0800264 times_shifted = 0;
265 factor_shifted = 0;
266 bNEGATED = false;
Eric Huang770911a2015-11-09 17:34:31 -0500267
Rex Zhu75ac63d2016-01-06 16:38:48 +0800268 if (X < 0) {
269 X = -1*X;
270 bNEGATED = true;
271 }
Eric Huang770911a2015-11-09 17:34:31 -0500272
Rex Zhu75ac63d2016-01-06 16:38:48 +0800273 if (factor < 0) {
274 factor = -1*factor;
275 bNEGATED = !bNEGATED; /*If bNEGATED = true due to X < 0, this will cover the case of negative cancelling negative */
276 }
Eric Huang770911a2015-11-09 17:34:31 -0500277
Rex Zhu75ac63d2016-01-06 16:38:48 +0800278 if ((X > MAX) || factor > MAX) {
279 if ((X/factor) <= MAX) {
280 while (X > MAX) {
281 X = X >> 1;
282 times_shifted++;
283 }
Eric Huang770911a2015-11-09 17:34:31 -0500284
Rex Zhu75ac63d2016-01-06 16:38:48 +0800285 while (factor > MAX) {
286 factor = factor >> 1;
287 factor_shifted++;
288 }
289 } else {
290 fValue.full = 0;
291 return fValue;
292 }
293 }
Eric Huang770911a2015-11-09 17:34:31 -0500294
Rex Zhu75ac63d2016-01-06 16:38:48 +0800295 if (factor == 1)
296 return (ConvertToFraction(X));
Eric Huang770911a2015-11-09 17:34:31 -0500297
Rex Zhu75ac63d2016-01-06 16:38:48 +0800298 fValue = fDivide(ConvertToFraction(X * uPow(-1, bNEGATED)), ConvertToFraction(factor));
Eric Huang770911a2015-11-09 17:34:31 -0500299
Rex Zhu75ac63d2016-01-06 16:38:48 +0800300 fValue.full = fValue.full << times_shifted;
301 fValue.full = fValue.full >> factor_shifted;
Eric Huang770911a2015-11-09 17:34:31 -0500302
Rex Zhu75ac63d2016-01-06 16:38:48 +0800303 return fValue;
Eric Huang770911a2015-11-09 17:34:31 -0500304}
305
306/* Addition using two fInts */
307fInt fAdd (fInt X, fInt Y)
308{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800309 fInt Sum;
Eric Huang770911a2015-11-09 17:34:31 -0500310
Rex Zhu75ac63d2016-01-06 16:38:48 +0800311 Sum.full = X.full + Y.full;
Eric Huang770911a2015-11-09 17:34:31 -0500312
Rex Zhu75ac63d2016-01-06 16:38:48 +0800313 return Sum;
Eric Huang770911a2015-11-09 17:34:31 -0500314}
315
316/* Addition using two fInts */
317fInt fSubtract (fInt X, fInt Y)
318{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800319 fInt Difference;
Eric Huang770911a2015-11-09 17:34:31 -0500320
Rex Zhu75ac63d2016-01-06 16:38:48 +0800321 Difference.full = X.full - Y.full;
Eric Huang770911a2015-11-09 17:34:31 -0500322
Rex Zhu75ac63d2016-01-06 16:38:48 +0800323 return Difference;
Eric Huang770911a2015-11-09 17:34:31 -0500324}
325
326bool Equal(fInt A, fInt B)
327{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800328 if (A.full == B.full)
329 return true;
330 else
331 return false;
Eric Huang770911a2015-11-09 17:34:31 -0500332}
333
334bool GreaterThan(fInt A, fInt B)
335{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800336 if (A.full > B.full)
337 return true;
338 else
339 return false;
Eric Huang770911a2015-11-09 17:34:31 -0500340}
341
342fInt fMultiply (fInt X, fInt Y) /* Uses 64-bit integers (int64_t) */
343{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800344 fInt Product;
345 int64_t tempProduct;
346 bool X_LessThanOne, Y_LessThanOne;
Eric Huang770911a2015-11-09 17:34:31 -0500347
Rex Zhu75ac63d2016-01-06 16:38:48 +0800348 X_LessThanOne = (X.partial.real == 0 && X.partial.decimal != 0 && X.full >= 0);
349 Y_LessThanOne = (Y.partial.real == 0 && Y.partial.decimal != 0 && Y.full >= 0);
Eric Huang770911a2015-11-09 17:34:31 -0500350
Rex Zhu75ac63d2016-01-06 16:38:48 +0800351 /*The following is for a very specific common case: Non-zero number with ONLY fractional portion*/
352 /* TEMPORARILY DISABLED - CAN BE USED TO IMPROVE PRECISION
Eric Huang770911a2015-11-09 17:34:31 -0500353
Rex Zhu75ac63d2016-01-06 16:38:48 +0800354 if (X_LessThanOne && Y_LessThanOne) {
355 Product.full = X.full * Y.full;
356 return Product
357 }*/
Eric Huang770911a2015-11-09 17:34:31 -0500358
Rex Zhu75ac63d2016-01-06 16:38:48 +0800359 tempProduct = ((int64_t)X.full) * ((int64_t)Y.full); /*Q(16,16)*Q(16,16) = Q(32, 32) - Might become a negative number! */
360 tempProduct = tempProduct >> 16; /*Remove lagging 16 bits - Will lose some precision from decimal; */
361 Product.full = (int)tempProduct; /*The int64_t will lose the leading 16 bits that were part of the integer portion */
Eric Huang770911a2015-11-09 17:34:31 -0500362
Rex Zhu75ac63d2016-01-06 16:38:48 +0800363 return Product;
Eric Huang770911a2015-11-09 17:34:31 -0500364}
365
366fInt fDivide (fInt X, fInt Y)
367{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800368 fInt fZERO, fQuotient;
369 int64_t longlongX, longlongY;
Eric Huang770911a2015-11-09 17:34:31 -0500370
Rex Zhu75ac63d2016-01-06 16:38:48 +0800371 fZERO = ConvertToFraction(0);
Eric Huang770911a2015-11-09 17:34:31 -0500372
Rex Zhu75ac63d2016-01-06 16:38:48 +0800373 if (Equal(Y, fZERO))
374 return fZERO;
Eric Huang770911a2015-11-09 17:34:31 -0500375
Rex Zhu75ac63d2016-01-06 16:38:48 +0800376 longlongX = (int64_t)X.full;
377 longlongY = (int64_t)Y.full;
Eric Huang770911a2015-11-09 17:34:31 -0500378
Rex Zhu75ac63d2016-01-06 16:38:48 +0800379 longlongX = longlongX << 16; /*Q(16,16) -> Q(32,32) */
Eric Huang770911a2015-11-09 17:34:31 -0500380
Rex Zhu75ac63d2016-01-06 16:38:48 +0800381 div64_s64(longlongX, longlongY); /*Q(32,32) divided by Q(16,16) = Q(16,16) Back to original format */
Eric Huang770911a2015-11-09 17:34:31 -0500382
Rex Zhu75ac63d2016-01-06 16:38:48 +0800383 fQuotient.full = (int)longlongX;
384 return fQuotient;
Eric Huang770911a2015-11-09 17:34:31 -0500385}
386
387int ConvertBackToInteger (fInt A) /*THIS is the function that will be used to check with the Golden settings table*/
388{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800389 fInt fullNumber, scaledDecimal, scaledReal;
Eric Huang770911a2015-11-09 17:34:31 -0500390
Rex Zhu75ac63d2016-01-06 16:38:48 +0800391 scaledReal.full = GetReal(A) * uPow(10, PRECISION-1); /* DOUBLE CHECK THISSSS!!! */
Eric Huang770911a2015-11-09 17:34:31 -0500392
Rex Zhu75ac63d2016-01-06 16:38:48 +0800393 scaledDecimal.full = uGetScaledDecimal(A);
Eric Huang770911a2015-11-09 17:34:31 -0500394
Rex Zhu75ac63d2016-01-06 16:38:48 +0800395 fullNumber = fAdd(scaledDecimal,scaledReal);
Eric Huang770911a2015-11-09 17:34:31 -0500396
Rex Zhu75ac63d2016-01-06 16:38:48 +0800397 return fullNumber.full;
Eric Huang770911a2015-11-09 17:34:31 -0500398}
399
400fInt fGetSquare(fInt A)
401{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800402 return fMultiply(A,A);
Eric Huang770911a2015-11-09 17:34:31 -0500403}
404
405/* x_new = x_old - (x_old^2 - C) / (2 * x_old) */
406fInt fSqrt(fInt num)
407{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800408 fInt F_divide_Fprime, Fprime;
409 fInt test;
410 fInt twoShifted;
411 int seed, counter, error;
412 fInt x_new, x_old, C, y;
Eric Huang770911a2015-11-09 17:34:31 -0500413
Rex Zhu75ac63d2016-01-06 16:38:48 +0800414 fInt fZERO = ConvertToFraction(0);
Eric Huang770911a2015-11-09 17:34:31 -0500415
Rex Zhu75ac63d2016-01-06 16:38:48 +0800416 /* (0 > num) is the same as (num < 0), i.e., num is negative */
Eric Huang770911a2015-11-09 17:34:31 -0500417
Rex Zhu75ac63d2016-01-06 16:38:48 +0800418 if (GreaterThan(fZERO, num) || Equal(fZERO, num))
419 return fZERO;
Eric Huang770911a2015-11-09 17:34:31 -0500420
Rex Zhu75ac63d2016-01-06 16:38:48 +0800421 C = num;
Eric Huang770911a2015-11-09 17:34:31 -0500422
Rex Zhu75ac63d2016-01-06 16:38:48 +0800423 if (num.partial.real > 3000)
424 seed = 60;
425 else if (num.partial.real > 1000)
426 seed = 30;
427 else if (num.partial.real > 100)
428 seed = 10;
429 else
430 seed = 2;
Eric Huang770911a2015-11-09 17:34:31 -0500431
Rex Zhu75ac63d2016-01-06 16:38:48 +0800432 counter = 0;
Eric Huang770911a2015-11-09 17:34:31 -0500433
Rex Zhu75ac63d2016-01-06 16:38:48 +0800434 if (Equal(num, fZERO)) /*Square Root of Zero is zero */
435 return fZERO;
Eric Huang770911a2015-11-09 17:34:31 -0500436
Rex Zhu75ac63d2016-01-06 16:38:48 +0800437 twoShifted = ConvertToFraction(2);
438 x_new = ConvertToFraction(seed);
Eric Huang770911a2015-11-09 17:34:31 -0500439
Rex Zhu75ac63d2016-01-06 16:38:48 +0800440 do {
441 counter++;
Eric Huang770911a2015-11-09 17:34:31 -0500442
Rex Zhu75ac63d2016-01-06 16:38:48 +0800443 x_old.full = x_new.full;
Eric Huang770911a2015-11-09 17:34:31 -0500444
Rex Zhu75ac63d2016-01-06 16:38:48 +0800445 test = fGetSquare(x_old); /*1.75*1.75 is reverting back to 1 when shifted down */
446 y = fSubtract(test, C); /*y = f(x) = x^2 - C; */
Eric Huang770911a2015-11-09 17:34:31 -0500447
Rex Zhu75ac63d2016-01-06 16:38:48 +0800448 Fprime = fMultiply(twoShifted, x_old);
449 F_divide_Fprime = fDivide(y, Fprime);
Eric Huang770911a2015-11-09 17:34:31 -0500450
Rex Zhu75ac63d2016-01-06 16:38:48 +0800451 x_new = fSubtract(x_old, F_divide_Fprime);
Eric Huang770911a2015-11-09 17:34:31 -0500452
Rex Zhu75ac63d2016-01-06 16:38:48 +0800453 error = ConvertBackToInteger(x_new) - ConvertBackToInteger(x_old);
Eric Huang770911a2015-11-09 17:34:31 -0500454
Rex Zhu75ac63d2016-01-06 16:38:48 +0800455 if (counter > 20) /*20 is already way too many iterations. If we dont have an answer by then, we never will*/
456 return x_new;
457
458 } while (uAbs(error) > 0);
459
460 return (x_new);
Eric Huang770911a2015-11-09 17:34:31 -0500461}
462
463void SolveQuadracticEqn(fInt A, fInt B, fInt C, fInt Roots[])
464{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800465 fInt *pRoots = &Roots[0];
466 fInt temp, root_first, root_second;
467 fInt f_CONSTANT10, f_CONSTANT100;
Eric Huang770911a2015-11-09 17:34:31 -0500468
Rex Zhu75ac63d2016-01-06 16:38:48 +0800469 f_CONSTANT100 = ConvertToFraction(100);
470 f_CONSTANT10 = ConvertToFraction(10);
Eric Huang770911a2015-11-09 17:34:31 -0500471
Rex Zhu75ac63d2016-01-06 16:38:48 +0800472 while(GreaterThan(A, f_CONSTANT100) || GreaterThan(B, f_CONSTANT100) || GreaterThan(C, f_CONSTANT100)) {
473 A = fDivide(A, f_CONSTANT10);
474 B = fDivide(B, f_CONSTANT10);
475 C = fDivide(C, f_CONSTANT10);
476 }
Eric Huang770911a2015-11-09 17:34:31 -0500477
Rex Zhu75ac63d2016-01-06 16:38:48 +0800478 temp = fMultiply(ConvertToFraction(4), A); /* root = 4*A */
479 temp = fMultiply(temp, C); /* root = 4*A*C */
480 temp = fSubtract(fGetSquare(B), temp); /* root = b^2 - 4AC */
481 temp = fSqrt(temp); /*root = Sqrt (b^2 - 4AC); */
Eric Huang770911a2015-11-09 17:34:31 -0500482
Rex Zhu75ac63d2016-01-06 16:38:48 +0800483 root_first = fSubtract(fNegate(B), temp); /* b - Sqrt(b^2 - 4AC) */
484 root_second = fAdd(fNegate(B), temp); /* b + Sqrt(b^2 - 4AC) */
Eric Huang770911a2015-11-09 17:34:31 -0500485
Rex Zhu75ac63d2016-01-06 16:38:48 +0800486 root_first = fDivide(root_first, ConvertToFraction(2)); /* [b +- Sqrt(b^2 - 4AC)]/[2] */
487 root_first = fDivide(root_first, A); /*[b +- Sqrt(b^2 - 4AC)]/[2*A] */
Eric Huang770911a2015-11-09 17:34:31 -0500488
Rex Zhu75ac63d2016-01-06 16:38:48 +0800489 root_second = fDivide(root_second, ConvertToFraction(2)); /* [b +- Sqrt(b^2 - 4AC)]/[2] */
490 root_second = fDivide(root_second, A); /*[b +- Sqrt(b^2 - 4AC)]/[2*A] */
Eric Huang770911a2015-11-09 17:34:31 -0500491
Rex Zhu75ac63d2016-01-06 16:38:48 +0800492 *(pRoots + 0) = root_first;
493 *(pRoots + 1) = root_second;
Eric Huang770911a2015-11-09 17:34:31 -0500494}
495
496/* -----------------------------------------------------------------------------
497 * SUPPORT FUNCTIONS
498 * -----------------------------------------------------------------------------
499 */
500
501/* Addition using two normal ints - Temporary - Use only for testing purposes?. */
502fInt Add (int X, int Y)
503{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800504 fInt A, B, Sum;
Eric Huang770911a2015-11-09 17:34:31 -0500505
Rex Zhu75ac63d2016-01-06 16:38:48 +0800506 A.full = (X << SHIFT_AMOUNT);
507 B.full = (Y << SHIFT_AMOUNT);
Eric Huang770911a2015-11-09 17:34:31 -0500508
Rex Zhu75ac63d2016-01-06 16:38:48 +0800509 Sum.full = A.full + B.full;
Eric Huang770911a2015-11-09 17:34:31 -0500510
Rex Zhu75ac63d2016-01-06 16:38:48 +0800511 return Sum;
Eric Huang770911a2015-11-09 17:34:31 -0500512}
513
514/* Conversion Functions */
515int GetReal (fInt A)
516{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800517 return (A.full >> SHIFT_AMOUNT);
Eric Huang770911a2015-11-09 17:34:31 -0500518}
519
520/* Temporarily Disabled */
521int GetRoundedValue(fInt A) /*For now, round the 3rd decimal place */
522{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800523 /* ROUNDING TEMPORARLY DISABLED
524 int temp = A.full;
525 int decimal_cutoff, decimal_mask = 0x000001FF;
526 decimal_cutoff = temp & decimal_mask;
527 if (decimal_cutoff > 0x147) {
528 temp += 673;
529 }*/
Eric Huang770911a2015-11-09 17:34:31 -0500530
Rex Zhu75ac63d2016-01-06 16:38:48 +0800531 return ConvertBackToInteger(A)/10000; /*Temporary - in case this was used somewhere else */
Eric Huang770911a2015-11-09 17:34:31 -0500532}
533
534fInt Multiply (int X, int Y)
535{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800536 fInt A, B, Product;
Eric Huang770911a2015-11-09 17:34:31 -0500537
Rex Zhu75ac63d2016-01-06 16:38:48 +0800538 A.full = X << SHIFT_AMOUNT;
539 B.full = Y << SHIFT_AMOUNT;
Eric Huang770911a2015-11-09 17:34:31 -0500540
Rex Zhu75ac63d2016-01-06 16:38:48 +0800541 Product = fMultiply(A, B);
Eric Huang770911a2015-11-09 17:34:31 -0500542
Rex Zhu75ac63d2016-01-06 16:38:48 +0800543 return Product;
Eric Huang770911a2015-11-09 17:34:31 -0500544}
Rex Zhu75ac63d2016-01-06 16:38:48 +0800545
Eric Huang770911a2015-11-09 17:34:31 -0500546fInt Divide (int X, int Y)
547{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800548 fInt A, B, Quotient;
Eric Huang770911a2015-11-09 17:34:31 -0500549
Rex Zhu75ac63d2016-01-06 16:38:48 +0800550 A.full = X << SHIFT_AMOUNT;
551 B.full = Y << SHIFT_AMOUNT;
Eric Huang770911a2015-11-09 17:34:31 -0500552
Rex Zhu75ac63d2016-01-06 16:38:48 +0800553 Quotient = fDivide(A, B);
Eric Huang770911a2015-11-09 17:34:31 -0500554
Rex Zhu75ac63d2016-01-06 16:38:48 +0800555 return Quotient;
Eric Huang770911a2015-11-09 17:34:31 -0500556}
557
558int uGetScaledDecimal (fInt A) /*Converts the fractional portion to whole integers - Costly function */
559{
560 int dec[PRECISION];
561 int i, scaledDecimal = 0, tmp = A.partial.decimal;
562
563 for (i = 0; i < PRECISION; i++) {
Rex Zhu75ac63d2016-01-06 16:38:48 +0800564 dec[i] = tmp / (1 << SHIFT_AMOUNT);
565 tmp = tmp - ((1 << SHIFT_AMOUNT)*dec[i]);
566 tmp *= 10;
567 scaledDecimal = scaledDecimal + dec[i]*uPow(10, PRECISION - 1 -i);
568 }
Eric Huang770911a2015-11-09 17:34:31 -0500569
Rex Zhu75ac63d2016-01-06 16:38:48 +0800570 return scaledDecimal;
Eric Huang770911a2015-11-09 17:34:31 -0500571}
572
573int uPow(int base, int power)
574{
575 if (power == 0)
576 return 1;
577 else
578 return (base)*uPow(base, power - 1);
579}
580
581fInt fAbs(fInt A)
582{
583 if (A.partial.real < 0)
584 return (fMultiply(A, ConvertToFraction(-1)));
585 else
586 return A;
587}
588
589int uAbs(int X)
590{
591 if (X < 0)
592 return (X * -1);
593 else
594 return X;
595}
596
597fInt fRoundUpByStepSize(fInt A, fInt fStepSize, bool error_term)
598{
Rex Zhu75ac63d2016-01-06 16:38:48 +0800599 fInt solution;
Eric Huang770911a2015-11-09 17:34:31 -0500600
Rex Zhu75ac63d2016-01-06 16:38:48 +0800601 solution = fDivide(A, fStepSize);
602 solution.partial.decimal = 0; /*All fractional digits changes to 0 */
Eric Huang770911a2015-11-09 17:34:31 -0500603
Rex Zhu75ac63d2016-01-06 16:38:48 +0800604 if (error_term)
605 solution.partial.real += 1; /*Error term of 1 added */
Eric Huang770911a2015-11-09 17:34:31 -0500606
Rex Zhu75ac63d2016-01-06 16:38:48 +0800607 solution = fMultiply(solution, fStepSize);
608 solution = fAdd(solution, fStepSize);
Eric Huang770911a2015-11-09 17:34:31 -0500609
Rex Zhu75ac63d2016-01-06 16:38:48 +0800610 return solution;
Eric Huang770911a2015-11-09 17:34:31 -0500611}
612