1 /* lrslong.h (lrs long integer arithmetic library */ 2 /* Copyright: David Avis 2000, avis@cs.mcgill.ca */ 3 /* Version 4.0, February 17, 2000 */ 4 5 /* This program is free software; you can redistribute it and/or modify 6 it under the terms of the GNU General Public License as published by 7 the Free Software Foundation; either version 2 of the License, or 8 (at your option) any later version. 9 10 This program is distributed in the hope that it will be useful, 11 but WITHOUT ANY WARRANTY; without even the implied warranty of 12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 GNU General Public License for more details. 14 15 You should have received a copy of the GNU General Public License 16 along with this program; if not, write to the Free Software 17 Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, USA. 18 */ 19 /******************************************************************************/ 20 /* See http://cgm.cs.mcgill.ca/~avis/C/lrs.html for lrs usage instructions */ 21 /******************************************************************************/ 22 /* This package contains the extended precision routines used by lrs 23 and some other miscellaneous routines. The maximum precision depends on 24 the parameter MAX_DIGITS defined below, with usual default of 255L. This 25 gives a maximum of 1020 decimal digits on 32 bit machines. The procedure 26 lrs_mp_init(dec_digits) may set a smaller number of dec_digits, and this 27 is useful if arrays or matrices will be used. 28 */ 29 30 31 /* 32 #ifdef PLRS 33 #include <string> 34 using namespace std; 35 #endif 36 */ 37 38 /***********/ 39 /* defines */ 40 /***********/ 41 /* 42 this is number of longwords. Increasing this won't cost you that much 43 since only variables other than the A matrix are allocated this size. 44 Changing affects running time in small but not very predictable ways. 45 */ 46 47 #define MAX_DIGITS 255L 48 49 /* 50 this is in decimal digits, you pay in memory if you increase this, 51 unless you override by a line with 52 digits n 53 before the begin line of your file. 54 */ 55 #define DEFAULT_DIGITS 100L 56 57 58 59 /**********MACHINE DEPENDENT CONSTANTS***********/ 60 /* MAXD is 2^(k-1)-1 where k is word size */ 61 /* MAXDm is 2^(k/2-1)-1 where k is word size */ 62 /* MAXDa is 2^(k-2)-1 where k is word size */ 63 /* MAXD must be at least 2*BASE^2 */ 64 /* If BASE is 10^k, use "%k.ku" for FORMAT */ 65 /* INTSIZE is number of bytes for integer */ 66 /* 64/128 bit arithmetic */ 67 /************************************************/ 68 #ifdef B128 69 /* 128 bit machines */ /* compiler does not accept big constants! */ 70 #define MAXD 9223372036854775807L /* should be 2^127 -1 but is 2^63 - 1 */ 71 #define MAXDm 9223372036854775807L /* 2^63 - 1 */ 72 #define MAXDa 9223372036854775807L /* should be 2^126 -1 but is 2^63 - 1 */ 73 /* max power of 10 fitting in signed int64 */ 74 #define P10_INT64 1000000000000000000ULL 75 76 #define BASE 1000000000L 77 #define FORMAT "%9.9u" 78 #define BASE_DIG 9 79 #define INTSIZE 16L 80 #define BIT "128bit" 81 #else 82 /* 64 bit machines */ 83 #define MAXD 9223372036854775807LL /* 2^63 - 1 */ 84 #define MAXDm 2147483647LL /* 2^31 - 1 */ 85 #define MAXDa 4611686018427387904LL /* 2^62 - 1 */ 86 #define BASE 1000000000L 87 #define FORMAT "%9.9u" 88 #define BASE_DIG 9 89 #define INTSIZE 16L 90 #define BIT "64bit" 91 #endif 92 93 #define MAXINPUT 1000 /*max length of any input rational */ 94 95 #define POS 1L 96 #define NEG -1L 97 #ifndef TRUE 98 #define TRUE 1L 99 #endif 100 #ifndef FALSE 101 #define FALSE 0L 102 #endif 103 #define ONE 1L 104 #define TWO 2L 105 #define ZERO 0L 106 107 /**********************************/ 108 /* MACROS */ 109 /* dependent on mp implementation */ 110 /**********************************/ 111 112 #ifdef SAFE 113 /* lazy but fast overflow checking */ 114 115 #define mpsafem(a,b) *(a)>MAXDm||*(b)>MAXDm||*(a)<-MAXDm||*(b)<-MAXDm 116 #define mpsafea(a,b) *(a)>MAXDa||*(b)>MAXDa||*(a)<-MAXDa||*(b)<-MAXDa 117 118 #ifdef DEBUG 119 #define mperrorm(a,b) fprintf(stdout," : max(|a|,|b|) > %ld\n",MAXDa);lrs_overflow(1) 120 #define mperrora(a,b) fprintf(stdout," : max(|a|,|b|) > %ld\n",MAXDa);lrs_overflow(1) 121 #define linint(a, ka, b, kb) if( mpsafem(a,b) ) {fprintf(stdout, "\n*linint ");mperrorm(a,b);} else *(a) = *(a) * ka + *(b) * kb 122 #define mulint(a, b, c) if( mpsafem(a,b) ) {fprintf(stdout, "\n*mulint ");mperrorm(a,b);} else *(c) = *(a) * *(b) 123 #define addint(a, b, c) if( mpsafea(a,b) ) {fprintf(stdout, "\n*addint ");mperrora(a,b);} else *(c) = *(a) + *(b) 124 #define subint(a, b, c) if( mpsafea(a,b) ) {fprintf(stdout, "\n*subint ");mperrora(a,b);} else *(c) = *(a) - *(b) 125 #define decint(a, b) if( mpsafea(a,b) ) {fprintf(stdout, "\n*decint ");mperrora(a,b);} else *(a) = *(a) - *(b) 126 #else 127 #define linint(a, ka, b, kb) if( mpsafem(a,b) ) lrs_overflow(1) ; else *(a) = *(a) * ka + *(b) * kb 128 #define mulint(a, b, c) if( mpsafem(a,b) ) lrs_overflow(1) ; else *(c) = *(a) * *(b) 129 #define addint(a, b, c) if( mpsafea(a,b) ) lrs_overflow(1) ; else *(c) = *(a) + *(b) 130 #define subint(a, b, c) if( mpsafea(a,b) ) lrs_overflow(1) ; else *(c) = *(a) - *(b) 131 #define decint(a, b) if( mpsafea(a,b) ) lrs_overflow(1) ; else *(a) = *(a) - *(b) 132 #endif 133 134 #else 135 /* unprotected routines */ 136 #define addint(a, b, c) *(c) = *(a) + *(b) 137 #define subint(a, b, c) *(c) = *(a) - *(b) 138 #define linint(a, ka, b, kb) *(a) = *(a) * ka + *(b) * kb 139 #define mulint(a, b, c) *(c) = *(a) * *(b) 140 #endif 141 142 #define unchecked_decint(a, b) *(a) = *(a) - *(b) /* only safe if a,b come from mulint */ 143 #define divint(a, b, c) *(c) = *(a) / *(b); *(a) = *(a) % *(b) 144 #define exactdivint(a,b,c) *(c) = *(a) / *(b); 145 146 #define abs128(a) (a>0? a : -1*a) 147 #define changesign(a) (*(a) = - *(a)) 148 #define copy(a, b) ((a)[0] = (b)[0]) 149 #define mp_greater(a, b) (*(a) > *(b) ) 150 #define itomp(in, a) *(a) = in 151 152 153 #define one(a) (*(a) == 1) 154 #define negative(a) (*(a) < 0) 155 #define normalize(a) (void) 0 156 #define positive(a) (*(a) > 0) 157 #define sign(a) (*(a) < 0 ? NEG : POS) 158 #ifndef B128 159 #define storesign(a, sa) (*(a) = labs(*(a)) * sa) 160 #else 161 #define storesign(a, sa) (*(a) = abs128(*(a)) * sa) 162 #endif 163 #define zero(a) (*(a) == 0) 164 165 166 /* 167 * convert between decimal and machine (longword digits). Notice lovely 168 * implementation of ceiling function :-) 169 */ 170 #define DEC2DIG(d) ( (d) % BASE_DIG ? (d)/BASE_DIG+1 : (d)/BASE_DIG) 171 #define DIG2DEC(d) ((d)*BASE_DIG) 172 173 #ifndef OMIT_SIGNALS 174 #include <signal.h> 175 #include <stdlib.h> /* labs */ 176 #include <unistd.h> 177 #define errcheck(s,e) if ((long)(e)==-1L){ perror(s);exit(1);} 178 #endif 179 180 #define CALLOC(n,s) xcalloc(n,s,__LINE__,__FILE__) 181 182 /*************/ 183 /* typedefs */ 184 /*************/ 185 #ifndef B128 186 typedef long long lrs_mp[1]; /* type lrs_mp holds one long integer */ 187 typedef long long *lrs_mp_t; 188 typedef long long **lrs_mp_vector; 189 typedef long long ***lrs_mp_matrix; 190 #else 191 typedef __int128 lrs_mp[1]; /* type lrs_mp holds one 128-bit integer */ 192 typedef __int128 *lrs_mp_t; 193 typedef __int128 **lrs_mp_vector; 194 typedef __int128 ***lrs_mp_matrix; 195 #endif 196 197 /*********************/ 198 /*global variables */ 199 /*********************/ 200 201 extern long lrs_digits; /* max permitted no. of digits */ 202 extern long lrs_record_digits; /* this is the biggest acheived so far. */ 203 204 extern FILE *lrs_ifp; /* input file pointer */ 205 extern FILE *lrs_ofp; /* output file pointer */ 206 207 /*********************************************************/ 208 /* Initialization and allocation procedures - must use! */ 209 /******************************************************* */ 210 211 //void mulint(lrs_mp a, lrs_mp b, lrs_mp c); 212 213 long lrs_mp_init (long dec_digits, FILE * lrs_ifp, FILE * lrs_ofp); /* max number of decimal digits, fps */ 214 215 #define lrs_alloc_mp(a) 216 #define lrs_clear_mp(a) 217 218 lrs_mp_t lrs_alloc_mp_t(); /* dynamic allocation of lrs_mp */ 219 lrs_mp_vector lrs_alloc_mp_vector (long n); /* allocate lrs_mp_vector for n+1 lrs_mp numbers */ 220 lrs_mp_matrix lrs_alloc_mp_matrix (long m, long n); /* allocate lrs_mp_matrix for m+1 x n+1 lrs_mp */ 221 222 void lrs_clear_mp_vector (lrs_mp_vector a, long n); 223 void lrs_clear_mp_matrix (lrs_mp_matrix a, long m, long n); 224 225 #ifndef MA 226 #define suf(func) func 227 #endif 228 229 #ifdef MA 230 #ifdef B128 231 #define suf(func) func##_2 232 #else 233 #define suf(func) func##_1 234 #endif 235 #endif 236 237 #define atoaa suf(atoaa) 238 #define atomp suf(atomp) 239 #define comprod suf(comprod) 240 #define divrat suf(divrat) 241 #define gcd suf(gcd) 242 #define getfactorial suf(getfactorial) 243 #define lcm suf(lcm) 244 #define linrat suf(linrat) 245 #define lrs_alloc_mp_matrix suf(lrs_alloc_mp_matrix) 246 #define lrs_alloc_mp_t suf(lrs_alloc_mp_t) 247 #define lrs_alloc_mp_vector suf(lrs_alloc_mp_vector) 248 #define lrs_clear_mp_matrix suf(lrs_clear_mp_matrix) 249 #define lrs_clear_mp_vector suf(lrs_clear_mp_vector) 250 #define lrs_digits suf(lrs_digits) 251 #define lrs_getdigits suf(lrs_getdigits) 252 #define lrs_overflow suf(lrs_overflow) 253 #define lrs_mp_init suf(lrs_mp_init) 254 #define lrs_record_digits suf(lrs_record_digits) 255 #define mptodouble suf(mptodouble) 256 #define mptoi suf(mptoi) 257 #define mpgetstr10 suf(mpgetstr10) 258 #define mulrat suf(mulrat) 259 #define myrandom suf(myrandom) 260 #define notimpl suf(notimpl) 261 #define pmp suf(pmp) 262 #define prat suf(prat) 263 #define cprat suf(cprat) 264 #define cpmp suf(cpmp) 265 #define rattodouble suf(rattodouble) 266 #define readmp suf(readmp) 267 #define readrat suf(readrat) 268 #define plrs_readrat suf(plrs_readrat) 269 #define reduce suf(reduce) 270 #define reducearray suf(reducearray) 271 #define reduceint suf(reduceint) 272 #define stringcpy suf(stringcpy) 273 #define xcalloc suf(xcalloc) 274 275 /*********************************************************/ 276 /* Core library functions - depend on mp implementation */ 277 /******************************************************* */ 278 279 void atomp (const char s[], lrs_mp a); /* convert string to lrs_mp integer */ 280 long compare (lrs_mp a, lrs_mp b); /* a ? b and returns -1,0,1 for <,=,> */ 281 void gcd (lrs_mp u, lrs_mp v); /* returns u=gcd(u,v) destroying v */ 282 void mptodouble (lrs_mp a, double *x); /* convert lrs_mp to double */ 283 long mptoi (lrs_mp a); /* convert lrs_mp to long integer */ 284 char *mpgetstr10(char *, lrs_mp); /* convert lrs_mp to string */ 285 #ifdef PLRS 286 long plrs_readrat (lrs_mp Na, lrs_mp Da, const char * rat); /* take a rational number and convert to lrs_mp */ 287 #endif 288 char *cprat(const char *name, lrs_mp Nt, lrs_mp Dt); /* mp rat to char */ 289 char *cpmp(const char *name, lrs_mp Nt); /* mp int to char */ 290 void pmp (const char *name, lrs_mp a); /* print the long precision integer a */ 291 void prat (const char *name, lrs_mp Nt, lrs_mp Dt); /* reduce and print Nt/Dt */ 292 void readmp (lrs_mp a); /* read an integer and convert to lrs_mp */ 293 long readrat (lrs_mp Na, lrs_mp Da); /* read a rational or int and convert to lrs_mp */ 294 void reduce (lrs_mp Na, lrs_mp Da); /* reduces Na Da by gcd(Na,Da) */ 295 296 /*********************************************************/ 297 /* Standard arithmetic & misc. functions */ 298 /* should be independent of mp implementation */ 299 /******************************************************* */ 300 301 void atoaa (const char in[], char num[], char den[]); /* convert rational string in to num/den strings */ 302 long atos (char s[]); /* convert s to integer */ 303 long comprod (lrs_mp Na, lrs_mp Nb, lrs_mp Nc, lrs_mp Nd); /* +1 if Na*Nb > Nc*Nd,-1 if Na*Nb > Nc*Nd else 0 */ 304 void divrat (lrs_mp Na, lrs_mp Da, lrs_mp Nb, lrs_mp Db, lrs_mp Nc, lrs_mp Dc); 305 /* computes Nc/Dc = (Na/Da) /( Nb/Db ) and reduce */ 306 void getfactorial (lrs_mp factorial, long k); /* compute k factorial in lrs_mp */ 307 void linrat (lrs_mp Na, lrs_mp Da, long ka, lrs_mp Nb, lrs_mp Db, long kb, lrs_mp Nc, lrs_mp Dc); 308 void lcm (lrs_mp a, lrs_mp b); /* a = least common multiple of a, b; b is saved */ 309 void mulrat (lrs_mp Na, lrs_mp Da, lrs_mp Nb, lrs_mp Db, lrs_mp Nc, lrs_mp Dc); 310 /* computes Nc/Dc=(Na/Da)*(Nb/Db) and reduce */ 311 long myrandom (long num, long nrange); /* return a random number in range 0..nrange-1 */ 312 void notimpl (const char *s); /* bail out - help! */ 313 void rattodouble (lrs_mp a, lrs_mp b, double *x); /* convert lrs_mp rational to double */ 314 void reduceint (lrs_mp Na, lrs_mp Da); /* divide Na by Da and return it */ 315 void reducearray (lrs_mp_vector p, long n); /* find gcd of p[0]..p[n-1] and divide through by */ 316 void scalerat (lrs_mp Na, lrs_mp Da, long ka); /* scales rational by ka */ 317 318 /********************************/ 319 /* Matrix and vector allocation */ 320 /********************************/ 321 void lrs_clear_mp_vector (lrs_mp_vector p, long n); 322 lrs_mp_vector lrs_alloc_mp_vector(long n); 323 void lrs_clear_mp_vector(lrs_mp_vector p, long n); 324 lrs_mp_matrix lrs_alloc_mp_matrix(long m, long n); 325 void lrs_clear_mp_matrix(lrs_mp_matrix p, long m, long n); 326 long lrs_mp_init(long dec_digits, FILE *fpin, FILE *fpout); 327 328 329 /* how big are numbers? */ 330 extern long lrs_digits; /* max permitted no. of digits */ 331 extern long lrs_record_digits; /* this is the biggest achieved so far. */ 332 333 /**********************************/ 334 /* Miscellaneous functions */ 335 /******************************** */ 336 337 void lrs_getdigits (long *a, long *b); /* send digit information to user */ 338 339 void stringcpy (char *s, char *t); /* copy t to s pointer version */ 340 341 void *calloc (); 342 void *malloc (); 343 void *xcalloc (long n, long s, long l, const char *f); 344 345 void lrs_default_digits_overflow (); 346 void lrs_exit(int i); 347 void lrs_overflow(int i); 348 /* end of lrs_long.h (vertex enumeration using lexicographic reverse search) */ 349