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