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