1 /*************************************************************************\
2 *
3 * Package: MyLib
4 * File: num.c
5 * Environment: ANSI C
6 *
7 * Copyright (c) 2002 Pierre L'Ecuyer, DIRO, Université de Montréal.
8 * e-mail: lecuyer@iro.umontreal.ca
9 * All rights reserved.
10 *
11 * Redistribution and use in source and binary forms, with or without
12 * modification, are permitted without a fee for private, research,
13 * academic, or other non-commercial purposes.
14 * Any use of this software in a commercial environment requires a
15 * written licence from the copyright owner.
16 *
17 * Any changes made to this package must be clearly identified as such.
18 *
19 * In scientific publications which used this software, a reference to it
20 * would be appreciated.
21 *
22 * Redistributions of source code must retain this copyright notice
23 * and the following disclaimer.
24 *
25 * THIS PACKAGE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS OR
26 * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
27 * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
28 *
29 \*************************************************************************/
30
31
32 #include "util.h"
33 #include "bitset.h"
34 #include "num.h"
35 #include <math.h>
36 #include <string.h>
37 #include <stdio.h>
38 #include <limits.h>
39
40
41
42
43 #define Deux53 9007199254740992.0 /* 2^53 */
44 #define Deux17 131072.0 /* 2^17 */
45 #define UnDeux17 7.62939453125E-6 /* 1 / 2^17 */
46 #define MASK32 0xffffffffUL
47
48 double num_TwoExp[num_MaxTwoExp + 1] = {
49 1.0, 2.0, 4.0, 8.0, 1.6e1, 3.2e1,
50 6.4e1, 1.28e2, 2.56e2, 5.12e2, 1.024e3,
51 2.048e3, 4.096e3, 8.192e3, 1.6384e4, 3.2768e4,
52 6.5536e4, 1.31072e5, 2.62144e5, 5.24288e5,
53 1.048576e6, 2.097152e6, 4.194304e6, 8.388608e6,
54 1.6777216e7, 3.3554432e7, 6.7108864e7,
55 1.34217728e8, 2.68435456e8, 5.36870912e8,
56 1.073741824e9, 2.147483648e9, 4.294967296e9,
57 8.589934592e9, 1.7179869184e10, 3.4359738368e10,
58 6.8719476736e10, 1.37438953472e11, 2.74877906944e11,
59 5.49755813888e11, 1.099511627776e12, 2.199023255552e12,
60 4.398046511104e12, 8.796093022208e12,
61 1.7592186044416e13, 3.5184372088832e13,
62 7.0368744177664e13, 1.40737488355328e14,
63 2.81474976710656e14, 5.62949953421312e14,
64 1.125899906842624e15, 2.251799813685248e15,
65 4.503599627370496e15, 9.007199254740992e15,
66 1.8014398509481984e16, 3.6028797018963968e16,
67 7.2057594037927936e16, 1.44115188075855872e17,
68 2.88230376151711744e17, 5.76460752303423488e17,
69 1.152921504606846976e18, 2.305843009213693952e18,
70 4.611686018427387904e18, 9.223372036854775808e18,
71 1.8446744073709551616e19
72 };
73
74
75 double num_TENNEGPOW[] = {
76 1.0, 1.0e-1, 1.0e-2, 1.0e-3, 1.0e-4, 1.0e-5, 1.0e-6, 1.0e-7, 1.0e-8,
77 1.0e-9, 1.0e-10, 1.0e-11, 1.0e-12, 1.0e-13, 1.0e-14, 1.0e-15, 1.0e-16
78 };
79
80
81
num_IsNumber(char S[])82 int num_IsNumber (char S[])
83 /*********************************************************
84 * Returns TRUE if the string S begin with a number *
85 * (with the possibility of spaces and a + or - sign *
86 * before the number). *
87 * e.g. *
88 * ' + 2' returns TRUE *
89 * '-+ 2' returns FALSE *
90 * '4hello' returns TRUE *
91 * 'hello' returns FALSE *
92 *********************************************************/
93 {
94 int Max;
95 int i;
96 int Sign;
97 Max = (int) (strlen (S) - 1);
98 Sign = 0;
99 for (i = 0; i < Max; i++) {
100 if (S[i] != ' ') {
101 if (S[i] == '+' || S[i] == '-') {
102 if (Sign) {
103 return 0;
104 }
105 /* We already saw a sign */
106 Sign = 1;
107 } else if ((unsigned char) S[i] >= '0' &&
108 (unsigned char) S[i] <= '9') {
109 return 1;
110 } else {
111 return 0;
112 }
113 }
114 } /* end for */
115 return 0; /* There's no digit in S */
116 } /* end IsNumber() */
117
118
num_IntToStrBase(long k,long b,char S[])119 void num_IntToStrBase (long k, long b, char S[])
120 {
121 int Sign; /* insert a '-' if TRUE */
122 long Char0;
123 long i;
124 long total;
125 long uppbound;
126 if (b < 2 || b > 10) {
127 util_Error ("*** Erreur: IntToStrB demande une b entre 2 et 10 ***");
128 }
129 Char0 = 48;
130 if (k < 0) {
131 Sign = 1;
132 S[0] = '-';
133 k = -k;
134 } else {
135 if (k == 0) {
136 S[0] = '0';
137 S[1] = '\0';
138 return;
139 }
140 Sign = 0;
141 }
142 i = k;
143 total = 0;
144 while (i > 0) {
145 i = (i / b);
146 ++total;
147 }
148 if (Sign)
149 uppbound = total + 1;
150 else
151 uppbound = total;
152 S[uppbound] = '\0';
153 for (i = 0; i < total - 1; i++) {
154 S[(uppbound - i) - 1] =
155 (char) ((int) fmod ((double) k, (double) b) + Char0);
156 k = (long) (k / b);
157 }
158 }
159
160
161 /*=========================================================================*/
162
num_Uint2Uchar(unsigned char * output,unsigned int * input,int L)163 void num_Uint2Uchar (unsigned char *output, unsigned int *input, int L)
164 {
165 int i, j;
166
167 for (i = 0, j = 0; i < L; i++, j += 4) {
168 output[j + 3] = (unsigned char) (input[i] & 0xff);
169 output[j + 2] = (unsigned char) ((input[i] >> 8) & 0xff);
170 output[j + 1] = (unsigned char) ((input[i] >> 16) & 0xff);
171 output[j] = (unsigned char) ((input[i] >> 24) & 0xff);
172 }
173 }
174
175
176 /*=========================================================================*/
177
num_WriteD(double x,int I,int J,int K)178 void num_WriteD (double x, int I, int J, int K)
179 {
180 int PosEntier = 0, /* Le nombre de positions occupees par la
181 partie entiere de x */
182 EntierSign, /* Le nombre de chiffres significatifs
183 avant le point */
184 Neg = 0; /* Nombre n'egatif */
185 char S[100];
186 char *p;
187
188 if (x == 0.0)
189 EntierSign = 1;
190 else {
191 EntierSign = PosEntier = floor (log10 (fabs (x)) + 1);
192 if (x < 0.0)
193 Neg = 1;
194 }
195 if (EntierSign <= 0)
196 PosEntier = 1;
197
198 if ((x == 0.0) ||
199 (((EntierSign + J) >= K) && (I >= (PosEntier + J + Neg + 1))))
200 printf ("%*.*f", I, J, x);
201
202 else { /* On doit utiliser la notation
203 scientifique. */
204 sprintf (S, "%*.*e", I, K - 1, x);
205 p = strstr (S, "e+0");
206 if (NULL == p)
207 p = strstr (S, "e-0");
208
209 /* remove the 0 in e-0 and in e+0 */
210 if (p) {
211 p += 2;
212 while ((*p = *(p + 1)))
213 p++;
214 printf (" "); /* pour utiliser au moins I espaces */
215 }
216 printf ("%s", S);
217 }
218 }
219
220
221 /***************************************************************************/
222
num_WriteBits(unsigned long x,int k)223 void num_WriteBits (unsigned long x, int k)
224 {
225 int i, n = CHAR_BIT * sizeof (unsigned long);
226 unsigned long mask = (unsigned long) 1 << (n - 1);
227 int spaces;
228 lebool flag = FALSE;
229
230 if (k > 0) {
231 spaces = k - n;
232 for (i = 0; i < spaces; i++)
233 printf (" ");
234 }
235 for (i = 0; i < n; i++) {
236 if (x & mask) {
237 printf ("1");
238 flag = TRUE;
239 } else if (flag)
240 printf ("0");
241 else
242 printf (" ");
243 mask >>= 1;
244 }
245 if (k < 0) {
246 spaces = -k - n;
247 for (i = 0; i < spaces; i++)
248 printf (" ");
249 }
250 }
251
252
253 /***************************************************************************/
254
255 #if LONG_MAX == 2147483647L
256 #define H 32768 /* = 2^d used in MultModL. */
257 #else
258 #define H 2147483648L
259 #endif
260
num_MultModL(long a,long s,long c,long m)261 long num_MultModL (long a, long s, long c, long m)
262 /* Suppose que 0 < a < m et 0 < s < m. Retourne (a*s + c) % m. */
263 /* Cette procedure est tiree de : */
264 /* L'Ecuyer, P. et Cote, S., A Random Number Package with */
265 /* Splitting Facilities, ACM TOMS, 1991. */
266 /* On coupe les entiers en blocs de d bits. H doit etre egal a 2^d. */
267 {
268 long a0, a1, q, qh, rh, k, p;
269 if (a < H) {
270 a0 = a;
271 p = 0;
272 } else {
273 a1 = a / H;
274 a0 = a - H * a1;
275 qh = m / H;
276 rh = m - H * qh;
277 if (a1 >= H) {
278 a1 = a1 - H;
279 k = s / qh;
280 p = H * (s - k * qh) - k * rh;
281 if (p < 0)
282 p = (p + 1) % m + m - 1;
283 } else /* p = (A2 * s * h) % m. */
284 p = 0;
285 if (a1 != 0) {
286 q = m / a1;
287 k = s / q;
288 p -= k * (m - a1 * q);
289 if (p > 0)
290 p -= m;
291 p += a1 * (s - k * q);
292 if (p < 0)
293 p = (p + 1) % m + m - 1;
294 } /* p = ((A2 * h + a1) * s) % m. */
295 k = p / qh;
296 p = H * (p - k * qh) - k * rh;
297 if (p < 0)
298 p = (p + 1) % m + m - 1;
299 } /* p = ((A2 * h + a1) * h * s) % m */
300 if (a0 != 0) {
301 q = m / a0;
302 k = s / q;
303 p -= k * (m - a0 * q);
304 if (p > 0)
305 p -= m;
306 p += a0 * (s - k * q);
307 if (p < 0)
308 p = (p + 1) % m + m - 1;
309 }
310 p = (p - m) + c;
311 if (p < 0)
312 p += m;
313 return p;
314 }
315
316 /*************************************************************************/
317
num_MultModD(double a,double s,double c,double m)318 double num_MultModD (double a, double s, double c, double m)
319 {
320 double V;
321 long k;
322 V = a * s + c;
323 if (V >= Deux53 || -V >= Deux53) {
324 k = a * UnDeux17;
325 a -= k * Deux17;
326 V = k * s;
327 k = V / m;
328 V -= k * m;
329 V = V * Deux17 + a * s + c;
330 }
331 k = V / m;
332 V -= k * m;
333 if (V < 0)
334 V += m;
335 return V;
336 }
337
338
339 /**************************************************************************/
340
num_InvEuclid(long M,long x)341 long num_InvEuclid (long M, long x)
342 /*
343 * Compute the inverse of x mod M by the modified Euclide
344 * algorithm (Knuth V2 p. 325).
345 */
346 {
347 long u1 = 0, u3 = M, v1 = 1, v3 = x;
348 long t1, t3, qq;
349 if (x == 0) return 0;
350
351 while (v3 != 0) {
352 qq = u3 / v3;
353 t1 = u1 - v1 * qq;
354 t3 = u3 - v3 * qq;
355 u1 = v1;
356 v1 = t1;
357 u3 = v3;
358 v3 = t3;
359 }
360 if (u1 < 0)
361 u1 += M;
362
363 if (u3 != 1) { /* In this case, the inverse does not exist! */
364 fprintf (stderr,
365 "ERROR in num_InvEuclid: inverse does not exist: m = %ld, x = %ld\n",
366 M, x);
367 return 0;
368 } else
369 return u1;
370 }
371
372
373 /*------------------------------------------------------------------------*/
374
num_InvExpon(int E,unsigned long Z)375 unsigned long num_InvExpon (int E, unsigned long Z)
376 /*
377 * Compute the inverse of Z modulo M = 2^E by exponentiation
378 */
379 {
380 int j;
381 unsigned long res = Z;
382
383 if (Z == 0) return 0;
384 if (!(Z & 1)) {
385 fprintf (stderr,
386 "ERROR in num_InvExpon: inverse does not exist: E = %d, Z = %ld\n",
387 E, Z);
388 return 0;
389 }
390 for (j = 1; j <= E - 3; j++)
391 res = res * res * Z;
392 return res & bitset_MASK[E];
393 }
394
395
396 /*------------------------------------------------------------------------*/
397
num_RoundL(double x)398 long num_RoundL (double x)
399 {
400 return (x >= 0) ? (long)(x + 0.5) : (long)(x - 0.5);
401 }
402
403
num_RoundD(double x)404 double num_RoundD (double x)
405 {
406 double z;
407 (x >= 0) ? modf(x + 0.5, &z) : modf(x - 0.5, &z);
408 return z;
409 }
410
411
412 /*------------------------------------------------------------------------*/
413