1 /* 2 * Copyright (c) 2002 by The XFree86 Project, Inc. 3 * 4 * Permission is hereby granted, free of charge, to any person obtaining a 5 * copy of this software and associated documentation files (the "Software"), 6 * to deal in the Software without restriction, including without limitation 7 * the rights to use, copy, modify, merge, publish, distribute, sublicense, 8 * and/or sell copies of the Software, and to permit persons to whom the 9 * Software is furnished to do so, subject to the following conditions: 10 * 11 * The above copyright notice and this permission notice shall be included in 12 * all copies or substantial portions of the Software. 13 * 14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 17 * THE XFREE86 PROJECT BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, 18 * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF 19 * OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE 20 * SOFTWARE. 21 * 22 * Except as contained in this notice, the name of the XFree86 Project shall 23 * not be used in advertising or otherwise to promote the sale, use or other 24 * dealings in this Software without prior written authorization from the 25 * XFree86 Project. 26 * 27 * Author: Paulo Csar Pereira de Andrade 28 */ 29 30 /* $XFree86$ */ 31 32 #include "mp.h" 33 34 /* 35 * TODO: 36 * Implement a fast gcd and divexact for integers, so that this code 37 * could be changed to do intermediary calculations faster using smaller 38 * numbers. 39 */ 40 41 /* 42 * Prototypes 43 */ 44 /* do the hard work of mpr_add and mpr_sub */ 45 static void mpr_addsub(mpr *rop, mpr *op1, mpr *op2, int sub); 46 47 /* do the hard work of mpr_cmp and mpr_cmpabs */ 48 static int mpr_docmp(mpr *op1, mpr *op2, int sign); 49 50 /* 51 * Implementation 52 */ 53 void 54 mpr_init(mpr *op) 55 { 56 op->num.digs = mp_malloc(sizeof(BNS)); 57 op->num.sign = 0; 58 op->num.size = op->num.alloc = 1; 59 op->num.digs[0] = 0; 60 61 op->den.digs = mp_malloc(sizeof(BNS)); 62 op->den.sign = 0; 63 op->den.size = op->den.alloc = 1; 64 op->den.digs[0] = 1; 65 } 66 67 void 68 mpr_clear(mpr *op) 69 { 70 op->num.sign = 0; 71 op->num.size = op->num.alloc = 0; 72 mp_free(op->num.digs); 73 74 op->den.sign = 0; 75 op->den.size = op->den.alloc = 0; 76 mp_free(op->den.digs); 77 } 78 79 void 80 mpr_set(mpr *rop, mpr *op) 81 { 82 if (rop != op) { 83 mpi_set(mpr_num(rop), mpr_num(op)); 84 mpi_set(mpr_den(rop), mpr_den(op)); 85 } 86 } 87 88 void 89 mpr_seti(mpr *rop, long num, long den) 90 { 91 mpi_seti(mpr_num(rop), num); 92 mpi_seti(mpr_den(rop), den); 93 } 94 95 void 96 mpr_setd(mpr *rop, double d) 97 { 98 double val, num; 99 int e, sign; 100 101 /* initialize */ 102 if (d < 0) { 103 sign = 1; 104 val = -d; 105 } 106 else { 107 sign = 0; 108 val = d; 109 } 110 111 val = frexp(val, &e); 112 while (modf(val, &num) != 0.0 && val <= DBL_MAX / 2.0) { 113 --e; 114 val *= 2.0; 115 } 116 if (e >= 0) { 117 mpi_setd(mpr_num(rop), d); 118 mpi_seti(mpr_den(rop), 1); 119 } 120 else { 121 mpi_setd(mpr_num(rop), sign ? -num : num); 122 mpi_setd(mpr_den(rop), ldexp(1.0, -e)); 123 } 124 } 125 126 void 127 mpr_setstr(mpr *rop, char *str, int base) 128 { 129 char *slash = strchr(str, '/'); 130 131 mpi_setstr(mpr_num(rop), str, base); 132 if (slash != NULL) 133 mpi_setstr(mpr_den(rop), slash + 1, base); 134 else 135 mpi_seti(mpr_den(rop), 1); 136 } 137 138 void 139 mpr_canonicalize(mpr *op) 140 { 141 mpi gcd; 142 143 memset(&gcd, '\0', sizeof(mpi)); 144 145 mpi_gcd(&gcd, mpr_num(op), mpr_den(op)); 146 if (mpi_cmpabsi(&gcd, 1)) { 147 mpi_div(mpr_num(op), mpr_num(op), &gcd); 148 mpi_div(mpr_den(op), mpr_den(op), &gcd); 149 } 150 151 if (op->den.sign) { 152 op->num.sign = !op->num.sign; 153 op->den.sign = 0; 154 } 155 156 mpi_clear(&gcd); 157 } 158 159 void 160 mpr_add(mpr *rop, mpr *op1, mpr *op2) 161 { 162 mpr_addsub(rop, op1, op2, 0); 163 } 164 165 void 166 mpr_addi(mpr *rop, mpr *op1, long op2) 167 { 168 mpi prod; 169 170 memset(&prod, '\0', sizeof(mpi)); 171 172 mpi_muli(&prod, mpr_den(op1), op2); 173 mpi_add(mpr_num(rop), mpr_num(op1), &prod); 174 mpi_clear(&prod); 175 } 176 177 void 178 mpr_sub(mpr *rop, mpr *op1, mpr *op2) 179 { 180 mpr_addsub(rop, op1, op2, 1); 181 } 182 183 void 184 mpr_subi(mpr *rop, mpr *op1, long op2) 185 { 186 mpi prod; 187 188 memset(&prod, '\0', sizeof(mpi)); 189 190 mpi_muli(&prod, mpr_den(op1), op2); 191 mpi_sub(mpr_num(rop), mpr_num(op1), &prod); 192 mpi_clear(&prod); 193 } 194 195 static void 196 mpr_addsub(mpr *rop, mpr *op1, mpr *op2, int sub) 197 { 198 mpi prod1, prod2; 199 200 memset(&prod1, '\0', sizeof(mpi)); 201 memset(&prod2, '\0', sizeof(mpi)); 202 203 mpi_mul(&prod1, mpr_num(op1), mpr_den(op2)); 204 mpi_mul(&prod2, mpr_num(op2), mpr_den(op1)); 205 206 if (sub) 207 mpi_sub(mpr_num(rop), &prod1, &prod2); 208 else 209 mpi_add(mpr_num(rop), &prod1, &prod2); 210 211 mpi_clear(&prod1); 212 mpi_clear(&prod2); 213 214 mpi_mul(mpr_den(rop), mpr_den(op1), mpr_den(op2)); 215 } 216 217 void 218 mpr_mul(mpr *rop, mpr *op1, mpr *op2) 219 { 220 /* check if temporary storage is required */ 221 if (op1 == op2 && rop == op1) { 222 mpi prod; 223 224 memset(&prod, '\0', sizeof(mpi)); 225 226 mpi_mul(&prod, mpr_num(op1), mpr_num(op2)); 227 mpi_mul(mpr_den(rop), mpr_den(op1), mpr_den(op2)); 228 mpi_set(mpr_num(rop), &prod); 229 230 mpi_clear(&prod); 231 } 232 else { 233 mpi_mul(mpr_num(rop), mpr_num(op1), mpr_num(op2)); 234 mpi_mul(mpr_den(rop), mpr_den(op1), mpr_den(op2)); 235 } 236 } 237 238 void 239 mpr_muli(mpr *rop, mpr *op1, long op2) 240 { 241 mpi_muli(mpr_num(rop), mpr_num(op1), op2); 242 } 243 244 void 245 mpr_div(mpr *rop, mpr *op1, mpr *op2) 246 { 247 /* check if temporary storage is required */ 248 if (op1 == op2 && rop == op1) { 249 mpi prod; 250 251 memset(&prod, '\0', sizeof(mpi)); 252 253 mpi_mul(&prod, mpr_num(op1), mpr_den(op2)); 254 mpi_mul(mpr_den(rop), mpr_num(op2), mpr_den(op1)); 255 mpi_set(mpr_num(rop), &prod); 256 257 mpi_clear(&prod); 258 } 259 else { 260 mpi_mul(mpr_num(rop), mpr_num(op1), mpr_den(op2)); 261 mpi_mul(mpr_den(rop), mpr_num(op2), mpr_den(op1)); 262 } 263 } 264 265 void 266 mpr_divi(mpr *rop, mpr *op1, long op2) 267 { 268 mpi_muli(mpr_den(rop), mpr_den(op1), op2); 269 } 270 271 void 272 mpr_inv(mpr *rop, mpr *op) 273 { 274 if (rop == op) 275 mpi_swap(mpr_num(op), mpr_den(op)); 276 else { 277 mpi_set(mpr_num(rop), mpr_den(op)); 278 mpi_set(mpr_den(rop), mpr_num(op)); 279 } 280 } 281 282 void 283 mpr_neg(mpr *rop, mpr *op) 284 { 285 mpi_neg(mpr_num(rop), mpr_num(op)); 286 mpi_set(mpr_den(rop), mpr_den(op)); 287 } 288 289 void 290 mpr_abs(mpr *rop, mpr *op) 291 { 292 if (mpr_num(op)->sign) 293 mpi_neg(mpr_num(rop), mpr_num(op)); 294 else 295 mpi_set(mpr_num(rop), mpr_num(op)); 296 297 /* op may not be canonicalized */ 298 if (mpr_den(op)->sign) 299 mpi_neg(mpr_den(rop), mpr_den(op)); 300 else 301 mpi_set(mpr_den(rop), mpr_den(op)); 302 } 303 304 int 305 mpr_cmp(mpr *op1, mpr *op2) 306 { 307 return (mpr_docmp(op1, op2, 1)); 308 } 309 310 int 311 mpr_cmpi(mpr *op1, long op2) 312 { 313 int cmp; 314 mpr rat; 315 316 mpr_init(&rat); 317 mpi_seti(mpr_num(&rat), op2); 318 cmp = mpr_docmp(op1, &rat, 1); 319 mpr_clear(&rat); 320 321 return (cmp); 322 } 323 324 int 325 mpr_cmpabs(mpr *op1, mpr *op2) 326 { 327 return (mpr_docmp(op1, op2, 0)); 328 } 329 330 int 331 mpr_cmpabsi(mpr *op1, long op2) 332 { 333 int cmp; 334 mpr rat; 335 336 mpr_init(&rat); 337 mpi_seti(mpr_num(&rat), op2); 338 cmp = mpr_docmp(op1, &rat, 0); 339 mpr_clear(&rat); 340 341 return (cmp); 342 } 343 344 static int 345 mpr_docmp(mpr *op1, mpr *op2, int sign) 346 { 347 int cmp, neg; 348 mpi prod1, prod2; 349 350 neg = 0; 351 if (sign) { 352 /* if op1 is negative */ 353 if (mpr_num(op1)->sign ^ mpr_den(op1)->sign) { 354 /* if op2 is positive */ 355 if (!(mpr_num(op2)->sign ^ mpr_den(op2)->sign)) 356 return (-1); 357 else 358 neg = 1; 359 } 360 /* if op2 is negative */ 361 else if (mpr_num(op2)->sign ^ mpr_den(op2)->sign) 362 return (1); 363 /* else same sign */ 364 } 365 366 /* if denominators are equal, compare numerators */ 367 if (mpi_cmpabs(mpr_den(op1), mpr_den(op2)) == 0) { 368 cmp = mpi_cmpabs(mpr_num(op1), mpr_num(op2)); 369 if (cmp == 0) 370 return (0); 371 if (sign && neg) 372 return (cmp < 0 ? 1 : -1); 373 return (cmp); 374 } 375 376 memset(&prod1, '\0', sizeof(mpi)); 377 memset(&prod2, '\0', sizeof(mpi)); 378 379 /* "divide" op1 by op2 380 * if result is smaller than 1, op1 is smaller than op2 */ 381 mpi_mul(&prod1, mpr_num(op1), mpr_den(op2)); 382 mpi_mul(&prod2, mpr_num(op2), mpr_den(op1)); 383 384 cmp = mpi_cmpabs(&prod1, &prod2); 385 386 mpi_clear(&prod1); 387 mpi_clear(&prod2); 388 389 if (sign && neg) 390 return (cmp < 0 ? 1 : -1); 391 return (cmp); 392 } 393 394 void 395 mpr_swap(mpr *op1, mpr *op2) 396 { 397 if (op1 != op2) { 398 mpr swap; 399 400 memcpy(&swap, op1, sizeof(mpr)); 401 memcpy(op1, op2, sizeof(mpr)); 402 memcpy(op2, &swap, sizeof(mpr)); 403 } 404 } 405 406 int 407 mpr_fiti(mpr *op) 408 { 409 return (mpi_fiti(mpr_num(op)) && mpi_fiti(mpr_den(op))); 410 } 411 412 double 413 mpr_getd(mpr *op) 414 { 415 return (mpi_getd(mpr_num(op)) / mpi_getd(mpr_den(op))); 416 } 417 418 char * 419 mpr_getstr(char *str, mpr *op, int base) 420 { 421 int len; 422 423 if (str == NULL) { 424 len = mpi_getsize(mpr_num(op), base) + mpr_num(op)->sign + 1 + 425 mpi_getsize(mpr_den(op), base) + mpr_den(op)->sign + 1; 426 427 str = mp_malloc(len); 428 } 429 430 (void)mpi_getstr(str, mpr_num(op), base); 431 len = strlen(str); 432 str[len] = '/'; 433 (void)mpi_getstr(str + len + 1, mpr_den(op), base); 434 435 return (str); 436 } 437
Indexes created Wed Jul 01 00:25:26 UTC 2026