/*************************************************************************\ * * Package: MyLib * File: num.c * Environment: ANSI C * * Copyright (c) 2002 Pierre L'Ecuyer, DIRO, Université de Montréal. * e-mail: lecuyer@iro.umontreal.ca * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted without a fee for private, research, * academic, or other non-commercial purposes. * Any use of this software in a commercial environment requires a * written licence from the copyright owner. * * Any changes made to this package must be clearly identified as such. * * In scientific publications which used this software, a reference to it * would be appreciated. * * Redistributions of source code must retain this copyright notice * and the following disclaimer. * * THIS PACKAGE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS OR * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE. * \*************************************************************************/ #include "util.h" #include "bitset.h" #include "num.h" #include #include #include #include #define Deux53 9007199254740992.0 /* 2^53 */ #define Deux17 131072.0 /* 2^17 */ #define UnDeux17 7.62939453125E-6 /* 1 / 2^17 */ #define MASK32 0xffffffffUL double num_TwoExp[num_MaxTwoExp + 1] = { 1.0, 2.0, 4.0, 8.0, 1.6e1, 3.2e1, 6.4e1, 1.28e2, 2.56e2, 5.12e2, 1.024e3, 2.048e3, 4.096e3, 8.192e3, 1.6384e4, 3.2768e4, 6.5536e4, 1.31072e5, 2.62144e5, 5.24288e5, 1.048576e6, 2.097152e6, 4.194304e6, 8.388608e6, 1.6777216e7, 3.3554432e7, 6.7108864e7, 1.34217728e8, 2.68435456e8, 5.36870912e8, 1.073741824e9, 2.147483648e9, 4.294967296e9, 8.589934592e9, 1.7179869184e10, 3.4359738368e10, 6.8719476736e10, 1.37438953472e11, 2.74877906944e11, 5.49755813888e11, 1.099511627776e12, 2.199023255552e12, 4.398046511104e12, 8.796093022208e12, 1.7592186044416e13, 3.5184372088832e13, 7.0368744177664e13, 1.40737488355328e14, 2.81474976710656e14, 5.62949953421312e14, 1.125899906842624e15, 2.251799813685248e15, 4.503599627370496e15, 9.007199254740992e15, 1.8014398509481984e16, 3.6028797018963968e16, 7.2057594037927936e16, 1.44115188075855872e17, 2.88230376151711744e17, 5.76460752303423488e17, 1.152921504606846976e18, 2.305843009213693952e18, 4.611686018427387904e18, 9.223372036854775808e18, 1.8446744073709551616e19 }; double num_TENNEGPOW[] = { 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, 1.0e-9, 1.0e-10, 1.0e-11, 1.0e-12, 1.0e-13, 1.0e-14, 1.0e-15, 1.0e-16 }; int num_IsNumber (char S[]) /********************************************************* * Returns TRUE if the string S begin with a number * * (with the possibility of spaces and a + or - sign * * before the number). * * e.g. * * ' + 2' returns TRUE * * '-+ 2' returns FALSE * * '4hello' returns TRUE * * 'hello' returns FALSE * *********************************************************/ { int Max; int i; int Sign; Max = (int) (strlen (S) - 1); Sign = 0; for (i = 0; i < Max; i++) { if (S[i] != ' ') { if (S[i] == '+' || S[i] == '-') { if (Sign) { return 0; } /* We already saw a sign */ Sign = 1; } else if ((unsigned char) S[i] >= '0' && (unsigned char) S[i] <= '9') { return 1; } else { return 0; } } } /* end for */ return 0; /* There's no digit in S */ } /* end IsNumber() */ void num_IntToStrBase (long k, long b, char S[]) { int Sign; /* insert a '-' if TRUE */ long Char0; long i; long total; long uppbound; if (b < 2 || b > 10) { util_Error ("*** Erreur: IntToStrB demande une b entre 2 et 10 ***"); } Char0 = 48; if (k < 0) { Sign = 1; S[0] = '-'; k = -k; } else { if (k == 0) { S[0] = '0'; S[1] = '\0'; return; } Sign = 0; } i = k; total = 0; while (i > 0) { i = (i / b); ++total; } if (Sign) uppbound = total + 1; else uppbound = total; S[uppbound] = '\0'; for (i = 0; i < total - 1; i++) { S[(uppbound - i) - 1] = (char) ((int) fmod ((double) k, (double) b) + Char0); k = (long) (k / b); } } /*=========================================================================*/ void num_Uint2Uchar (unsigned char *output, unsigned int *input, int L) { int i, j; for (i = 0, j = 0; i < L; i++, j += 4) { output[j + 3] = (unsigned char) (input[i] & 0xff); output[j + 2] = (unsigned char) ((input[i] >> 8) & 0xff); output[j + 1] = (unsigned char) ((input[i] >> 16) & 0xff); output[j] = (unsigned char) ((input[i] >> 24) & 0xff); } } /*=========================================================================*/ void num_WriteD (double x, int I, int J, int K) { int PosEntier = 0, /* Le nombre de positions occupees par la partie entiere de x */ EntierSign, /* Le nombre de chiffres significatifs avant le point */ Neg = 0; /* Nombre n'egatif */ char S[100]; char *p; if (x == 0.0) EntierSign = 1; else { EntierSign = PosEntier = floor (log10 (fabs (x)) + 1); if (x < 0.0) Neg = 1; } if (EntierSign <= 0) PosEntier = 1; if ((x == 0.0) || (((EntierSign + J) >= K) && (I >= (PosEntier + J + Neg + 1)))) printf ("%*.*f", I, J, x); else { /* On doit utiliser la notation scientifique. */ sprintf (S, "%*.*e", I, K - 1, x); p = strstr (S, "e+0"); if (NULL == p) p = strstr (S, "e-0"); /* remove the 0 in e-0 and in e+0 */ if (p) { p += 2; while ((*p = *(p + 1))) p++; printf (" "); /* pour utiliser au moins I espaces */ } printf ("%s", S); } } /***************************************************************************/ void num_WriteBits (unsigned long x, int k) { int i, n = CHAR_BIT * sizeof (unsigned long); unsigned long mask = (unsigned long) 1 << (n - 1); int spaces; lebool flag = FALSE; if (k > 0) { spaces = k - n; for (i = 0; i < spaces; i++) printf (" "); } for (i = 0; i < n; i++) { if (x & mask) { printf ("1"); flag = TRUE; } else if (flag) printf ("0"); else printf (" "); mask >>= 1; } if (k < 0) { spaces = -k - n; for (i = 0; i < spaces; i++) printf (" "); } } /***************************************************************************/ #if LONG_MAX == 2147483647L #define H 32768 /* = 2^d used in MultModL. */ #else #define H 2147483648L #endif long num_MultModL (long a, long s, long c, long m) /* Suppose que 0 < a < m et 0 < s < m. Retourne (a*s + c) % m. */ /* Cette procedure est tiree de : */ /* L'Ecuyer, P. et Cote, S., A Random Number Package with */ /* Splitting Facilities, ACM TOMS, 1991. */ /* On coupe les entiers en blocs de d bits. H doit etre egal a 2^d. */ { long a0, a1, q, qh, rh, k, p; if (a < H) { a0 = a; p = 0; } else { a1 = a / H; a0 = a - H * a1; qh = m / H; rh = m - H * qh; if (a1 >= H) { a1 = a1 - H; k = s / qh; p = H * (s - k * qh) - k * rh; if (p < 0) p = (p + 1) % m + m - 1; } else /* p = (A2 * s * h) % m. */ p = 0; if (a1 != 0) { q = m / a1; k = s / q; p -= k * (m - a1 * q); if (p > 0) p -= m; p += a1 * (s - k * q); if (p < 0) p = (p + 1) % m + m - 1; } /* p = ((A2 * h + a1) * s) % m. */ k = p / qh; p = H * (p - k * qh) - k * rh; if (p < 0) p = (p + 1) % m + m - 1; } /* p = ((A2 * h + a1) * h * s) % m */ if (a0 != 0) { q = m / a0; k = s / q; p -= k * (m - a0 * q); if (p > 0) p -= m; p += a0 * (s - k * q); if (p < 0) p = (p + 1) % m + m - 1; } p = (p - m) + c; if (p < 0) p += m; return p; } /*************************************************************************/ double num_MultModD (double a, double s, double c, double m) { double V; long k; V = a * s + c; if (V >= Deux53 || -V >= Deux53) { k = a * UnDeux17; a -= k * Deux17; V = k * s; k = V / m; V -= k * m; V = V * Deux17 + a * s + c; } k = V / m; V -= k * m; if (V < 0) V += m; return V; } /**************************************************************************/ long num_InvEuclid (long M, long x) /* * Compute the inverse of x mod M by the modified Euclide * algorithm (Knuth V2 p. 325). */ { long u1 = 0, u3 = M, v1 = 1, v3 = x; long t1, t3, qq; if (x == 0) return 0; while (v3 != 0) { qq = u3 / v3; t1 = u1 - v1 * qq; t3 = u3 - v3 * qq; u1 = v1; v1 = t1; u3 = v3; v3 = t3; } if (u1 < 0) u1 += M; if (u3 != 1) { /* In this case, the inverse does not exist! */ fprintf (stderr, "ERROR in num_InvEuclid: inverse does not exist: m = %ld, x = %ld\n", M, x); return 0; } else return u1; } /*------------------------------------------------------------------------*/ unsigned long num_InvExpon (int E, unsigned long Z) /* * Compute the inverse of Z modulo M = 2^E by exponentiation */ { int j; unsigned long res = Z; if (Z == 0) return 0; if (!(Z & 1)) { fprintf (stderr, "ERROR in num_InvExpon: inverse does not exist: E = %d, Z = %ld\n", E, Z); return 0; } for (j = 1; j <= E - 3; j++) res = res * res * Z; return res & bitset_MASK[E]; } /*------------------------------------------------------------------------*/ long num_RoundL (double x) { return (x >= 0) ? (long)(x + 0.5) : (long)(x - 0.5); } double num_RoundD (double x) { double z; (x >= 0) ? modf(x + 0.5, &z) : modf(x - 0.5, &z); return z; } /*------------------------------------------------------------------------*/