PandA-2024.02
fixedptc.c
Go to the documentation of this file.
1 #include <stdio.h>
2 #include <string.h>
3 #include <unistd.h>
4 #include <stdint.h>
5 #include "fixedptc.h"
6 
7 
8 char str[25];
9 
10 /* Multiplies two fixedpt numbers, returns the result. */
11 fixedpt
13 {
14  return (((fixedptd)A * (fixedptd)B) >> FIXEDPT_FBITS);
15 }
16 
17 
18 /* Divides two fixedpt numbers, returns the result. */
19 fixedpt
21 {
22  return (((fixedptd)A << FIXEDPT_FBITS) / (fixedptd)B);
23 }
24 
25 /*
26  * Note: adding and substracting fixedpt numbers can be done by using
27  * the regular integer operators + and -.
28  */
29 
39 void __attribute__ ((noinline))
40 fixedpt_str(fixedpt A, char *str, int max_dec)
41 {
42  int ndec = 0, slen = 0;
43  char tmp[12] = {0};
44  fixedptud fr, ip;
45  const fixedptud one = (fixedptud)1 << FIXEDPT_BITS;
46  const fixedptud mask = one - 1;
47 
48  if (max_dec == -1)
49 #if FIXEDPT_BITS == 32
50  max_dec = 2;
51 #elif FIXEDPT_BITS == 64
52  max_dec = 10;
53 #else
54 #error Invalid width
55 #endif
56  else if (max_dec == -2)
57  max_dec = 15;
58 
59  if (A < 0) {
60  str[slen++] = '-';
61  A *= -1;
62  }
63 
64  ip = fixedpt_toint(A);
65  do {
66  tmp[ndec++] = '0' + ip % 10;
67  ip /= 10;
68  } while (ip != 0);
69 
70  while (ndec > 0)
71  str[slen++] = tmp[--ndec];
72  str[slen++] = '.';
73 
74  // printf("before shift left %x\n", fr);
75  fr = (fixedpt_fracpart(A) << FIXEDPT_WBITS) & mask;
76 // printf("before iter %x\n", fr);
77  do {
78  fr = (fr & mask) * 10;
79  // printf("iter %x\n", fr);
80  str[slen++] = '0' + (fr >> FIXEDPT_BITS) % 10;
81  // printf("char %c\n", '0' + (fr >> FIXEDPT_BITS) % 10);
82  ndec++;
83  } while (fr != 0 && ndec < max_dec);
84 
85  if (ndec > 1 && str[slen-1] == '0')
86  str[slen-1] = '\0'; /* cut off trailing 0 */
87  else
88  str[slen] = '\0';
89 }
90 
91 /* Converts the given fixedpt number into a string, using a static
92  * (non-threadsafe) string buffer */
93 char*
94 fixedpt_cstr(const fixedpt A, const int max_dec)
95 {
96 
97  fixedpt_str(A, str, max_dec);
98  return (str);
99 }
100 
101 
102 /* Returns the square root of the given number, or -1 in case of error */
103 fixedpt __attribute__ ((noinline))
105 {
106  int invert = 0;
107  int iter = FIXEDPT_FBITS;
108  int l, i;
109 
110  if (A < 0)
111  return (-1);
112  if (A == 0 || A == FIXEDPT_ONE)
113  return (A);
114  if (A < FIXEDPT_ONE && A > 6) {
115  invert = 1;
116  A = fixedpt_div(FIXEDPT_ONE, A);
117  }
118  if (A > FIXEDPT_ONE) {
119  int s = A;
120 
121  iter = 0;
122  while (s > 0) {
123  s >>= 2;
124  iter++;
125  }
126  }
127 
128  /* Newton's iterations */
129  l = (A >> 1) + 1;
130  for (i = 0; i < iter; i++)
131  l = (l + fixedpt_div(A, l)) >> 1;
132  if (invert)
133  return (fixedpt_div(FIXEDPT_ONE, l));
134  return (l);
135 }
136 
137 
138 /* Returns the sine of the given fixedpt number.
139  * Note: the loss of precision is extraordinary! */
140 fixedpt
142 {
143  int sign = 1;
144  fixedpt sqr, result;
145  const fixedpt SK[2] = {
146  fixedpt_rconst(7.61e-03),
147  fixedpt_rconst(1.6605e-01)
148  };
149 
150  fp %= 2 * FIXEDPT_PI;
151  if (fp < 0)
152  fp = FIXEDPT_PI * 2 + fp;
153  if ((fp > FIXEDPT_HALF_PI) && (fp <= FIXEDPT_PI))
154  fp = FIXEDPT_PI - fp;
155  else if ((fp > FIXEDPT_PI) && (fp <= (FIXEDPT_PI + FIXEDPT_HALF_PI))) {
156  fp = fp - FIXEDPT_PI;
157  sign = -1;
158  } else if (fp > (FIXEDPT_PI + FIXEDPT_HALF_PI)) {
159  fp = (FIXEDPT_PI << 1) - fp;
160  sign = -1;
161  }
162  sqr = fixedpt_mul(fp, fp);
163  result = SK[0];
164  result = fixedpt_mul(result, sqr);
165  result -= SK[1];
166  result = fixedpt_mul(result, sqr);
167  result += FIXEDPT_ONE;
168  result = fixedpt_mul(result, fp);
169  return sign * result;
170 }
171 
172 
173 /* Returns the cosine of the given fixedpt number */
174 fixedpt
176 {
177  return (fixedpt_sin(FIXEDPT_HALF_PI - A));
178 }
179 
180 
181 /* Returns the tangens of the given fixedpt number */
182 fixedpt
184 {
185  return fixedpt_div(fixedpt_sin(A), fixedpt_cos(A));
186 }
187 
188 
189 /* Returns the value exp(x), i.e. e^x of the given fixedpt number. */
190 fixedpt
192 {
193  fixedpt xabs, k, z, R, xp, ans;
194  const fixedpt LN2 = fixedpt_rconst(0.69314718055994530942);
195  const fixedpt LN2_INV = fixedpt_rconst(1.4426950408889634074);
196  const fixedpt EXP_P[5] = {
197  fixedpt_rconst(1.66666666666666019037e-01),
198  fixedpt_rconst(-2.77777777770155933842e-03),
199  fixedpt_rconst(6.61375632143793436117e-05),
200  fixedpt_rconst(-1.65339022054652515390e-06),
201  fixedpt_rconst(4.13813679705723846039e-08),
202  };
203 
204  if (fp == 0){
205 // printf("excuted fp = 0 ");
206  return (FIXEDPT_ONE);
207  }
208 
209 
210  xabs = fixedpt_abs(fp);
211 // printf("xabs = fixedpt_abs(fp); = ");
212 // fixedpt_print(xabs);
213 
214  k = fixedpt_mul(xabs, LN2_INV);
215  //switch the exponents to base 2
216  //now calculating 2^k
217 // printf("k = fixedpt_mul(xabs, LN2_INV); = ");
218 // fixedpt_print(k);
219 
220  k += FIXEDPT_ONE_HALF;
221 // printf("k += FIXEDPT_ONE_HALF; = ");
222 // fixedpt_print(k);
223 
224  k &= ~FIXEDPT_FMASK;
225 // printf("k &= ~FIXEDPT_FMASK; = ");
226 // fixedpt_print(k);
227 
228  if (fp < 0){
229  k = -k;
230 // printf("k = -k; = ");
231 // fixedpt_print(k);
232  }
233 
234  fp -= fixedpt_mul(k, LN2);
235 // printf("fp -= fixedpt_mul(k, LN2);");
236 // fixedpt_print(fp);
237 
238  z = fixedpt_mul(fp, fp);
239 // printf("z = fixedpt_mul(fp, fp); = ");
240 // fixedpt_print(z);
241 
242 
243  /* Taylor */
244  R = FIXEDPT_TWO +
245  fixedpt_mul(z, EXP_P[0] + fixedpt_mul(z, EXP_P[1] +
246  fixedpt_mul(z, EXP_P[2] + fixedpt_mul(z, EXP_P[3] +
247  fixedpt_mul(z, EXP_P[4])))));
248 // printf("R = FIXEDPT_TWO +"
249 // "fixedpt_mul(z, EXP_P[0] + fixedpt_mul(z, EXP_P[1] +"
250 // "fixedpt_mul(z, EXP_P[2] + fixedpt_mul(z, EXP_P[3] +"
251 // "fixedpt_mul(z, EXP_P[4])))));= ");
252 // fixedpt_print(R);
253 
254 
255  xp = FIXEDPT_ONE + fixedpt_div(fixedpt_mul(fp, FIXEDPT_TWO), R - fp);
256 // printf("xp = FIXEDPT_ONE + fixedpt_div(fixedpt_mul(fp, FIXEDPT_TWO), R - fp); ");
257 // fixedpt_print(xp);
258 
259 
260  if (k < 0){
261  k = FIXEDPT_ONE >> (-k >> FIXEDPT_FBITS);
262 // printf("k = FIXEDPT_ONE >> (-k >> FIXEDPT_FBITS);");
263 // fixedpt_print(k);
264  }
265  else{
266  k = FIXEDPT_ONE << (k >> FIXEDPT_FBITS);
267 // printf("k = FIXEDPT_ONE << (k >> FIXEDPT_FBITS);");
268 // fixedpt_print(k);
269  }
270  ans = fixedpt_mul(k, xp);
271 // printf("ans = fixedpt_mul(k, xp)");
272 // fixedpt_print(ans);
273  return (ans);
274 }
275 
276 
277 /* Returns the natural logarithm of the given fixedpt number. */
278 fixedpt
280 {
281  fixedpt log2, xi;
282  fixedpt f, s, z, w, R;
283  const fixedpt LN2 = fixedpt_rconst(0.69314718055994530942);
284  const fixedpt LG[7] = {
285  fixedpt_rconst(6.666666666666735130e-01),
286  fixedpt_rconst(3.999999999940941908e-01),
287  fixedpt_rconst(2.857142874366239149e-01),
288  fixedpt_rconst(2.222219843214978396e-01),
289  fixedpt_rconst(1.818357216161805012e-01),
290  fixedpt_rconst(1.531383769920937332e-01),
291  fixedpt_rconst(1.479819860511658591e-01)
292  };
293 
294  if (x < 0)
295  return (0);
296  if (x == 0)
297  return 0xffffffff;
298 
299  log2 = 0;
300  xi = x;
301  while (xi > FIXEDPT_TWO) {
302  xi >>= 1;
303  log2++;
304  }
305  f = xi - FIXEDPT_ONE;
306  s = fixedpt_div(f, FIXEDPT_TWO + f);
307  z = fixedpt_mul(s, s);
308  w = fixedpt_mul(z, z);
309  R = fixedpt_mul(w, LG[1] + fixedpt_mul(w, LG[3]
310  + fixedpt_mul(w, LG[5]))) + fixedpt_mul(z, LG[0]
311  + fixedpt_mul(w, LG[2] + fixedpt_mul(w, LG[4]
312  + fixedpt_mul(w, LG[6]))));
313  return (fixedpt_mul(LN2, (log2 << FIXEDPT_FBITS)) + f
314  - fixedpt_mul(s, f - R));
315 }
316 
317 
318 /* Returns the logarithm of the given base of the given fixedpt number */
319 fixedpt
321 {
322  return (fixedpt_div(fixedpt_ln(x), fixedpt_ln(base)));
323 }
324 
325 
326 /* Return the power value (n^exp) of the given fixedpt numbers */
327 fixedpt
329 {
330  if (exp == 0)
331  return (FIXEDPT_ONE);
332  if (n < 0)
333  return 0;
334  return (fixedpt_exp(fixedpt_mul(fixedpt_ln(n), exp)));
335 }
336 /*
337 void
338 fixedpt_print(fixedpt A)
339 {
340  char num[30];
341 
342  fixedpt_str(A, num, -2);
343  puts(num);
344 }*/
345 /*
346 void
347 fixedpt_print_file(fixedpt A)
348 {
349  char num[30];
350  int j;
351  FILE *output;
352  output = fopen ( "MonteCarlo.txt", "a" );
353 
354  if ( !output )
355  {
356  fprintf ( stderr, "\n" );
357  fprintf ( stderr, "fixedpt_print_file - Fatal error!\n" );
358  fprintf ( stderr, " Could not open the output file.\n" );
359  return;
360  }
361  //c printf("PRINT_FUNCTON **********************************************\n");
362 // for ( j = 0; j <= n; j++ )
363 // {
364  fixedpt_str(A, num, -2);
365  printf("Blahhhh");
366  puts(num);
367  fputs(num, output);
368  fprintf(output, "\n");
369  //fprintf ( output, " %24.16g\n", num );
370  // }
371 }*/
372 
fixedpt fixedpt_sqrt(fixedpt A)
Definition: fixedptc.c:96
#define FIXEDPT_FBITS
Definition: fixedptc.h:104
char base
This version is stamped on May 10, 2016.
Definition: nussinov.c:24
void fixedpt_str(fixedpt A, char *str, int max_dec)
Convert the given fixedpt number to a decimal string.
Definition: fixedptc.c:35
#define R
Definition: mips.c:40
fixedpt fixedpt_log(fixedpt x, fixedpt base)
Definition: fixedptc.c:298
#define FIXEDPT_ONE_HALF
Definition: fixedptc.h:119
fixedpt fixedpt_pow(fixedpt n, fixedpt exp)
Definition: fixedptc.c:304
#define FIXEDPT_BITS
Definition: fixedptc.h:75
fixedpt fixedpt_ln(fixedpt x)
Definition: fixedptc.c:259
void __attribute__((noinline))
Convert the given fixedpt number to a decimal string.
Definition: fixedptc.c:39
#define fixedpt_toint(F)
Definition: fixedptc.h:109
#define A
Definition: generate.c:13
#define FIXEDPT_TWO
Definition: fixedptc.h:120
#define FIXEDPT_PI
Definition: fixedptc.h:121
#define FIXEDPT_FMASK
Definition: fixedptc.h:105
static const uint32_t k[]
Definition: sha-256.c:22
fixedpt fixedpt_tan(fixedpt A)
Definition: fixedptc.c:168
fixedpt fixedpt_exp(fixedpt fp)
Definition: fixedptc.c:174
#define fixedpt_fracpart(A)
Definition: fixedptc.h:116
fixedpt fixedpt_mul(fixedpt A, fixedpt B)
Definition: fixedptc.c:11
fixedpt fixedpt_sin(fixedpt fp)
Definition: fixedptc.c:130
fixedpt fixedpt_div(fixedpt A, fixedpt B)
Definition: fixedptc.c:17
int result[SIZE]
Definition: adpcm.c:800
char * fixedpt_cstr(const fixedpt A, const int max_dec)
Definition: fixedptc.c:88
fixedpt fixedpt_cos(fixedpt A)
Definition: fixedptc.c:162
int64_t fixedptd
Definition: fixedptc.h:80
char str[25]
Definition: fixedptc.c:8
uint64_t fixedptud
Definition: fixedptc.h:82
#define fixedpt_rconst(R)
Definition: fixedptc.h:107
uint32_t sign
x
Return the smallest n such that 2^n >= _x.
#define FIXEDPT_ONE
Definition: fixedptc.h:118
#define fixedpt_abs(A)
Definition: fixedptc.h:126
#define B
Definition: generate.c:14
int32_t fixedpt
Definition: fixedptc.h:79
uint32_t exp
#define FIXEDPT_HALF_PI
Definition: fixedptc.h:123
#define FIXEDPT_WBITS
Definition: fixedptc.h:95

Generated on Mon Feb 12 2024 13:02:49 for PandA-2024.02 by doxygen 1.8.13