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 César 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
mpr_init(mpr * op)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
mpr_clear(mpr * op)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
mpr_set(mpr * rop,mpr * op)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
mpr_seti(mpr * rop,long num,long den)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
mpr_setd(mpr * rop,double d)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
mpr_setstr(mpr * rop,char * str,int base)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
mpr_canonicalize(mpr * op)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
mpr_add(mpr * rop,mpr * op1,mpr * op2)160 mpr_add(mpr *rop, mpr *op1, mpr *op2)
161 {
162 mpr_addsub(rop, op1, op2, 0);
163 }
164
165 void
mpr_addi(mpr * rop,mpr * op1,long op2)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
mpr_sub(mpr * rop,mpr * op1,mpr * op2)178 mpr_sub(mpr *rop, mpr *op1, mpr *op2)
179 {
180 mpr_addsub(rop, op1, op2, 1);
181 }
182
183 void
mpr_subi(mpr * rop,mpr * op1,long op2)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
mpr_addsub(mpr * rop,mpr * op1,mpr * op2,int sub)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
mpr_mul(mpr * rop,mpr * op1,mpr * op2)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
mpr_muli(mpr * rop,mpr * op1,long op2)239 mpr_muli(mpr *rop, mpr *op1, long op2)
240 {
241 mpi_muli(mpr_num(rop), mpr_num(op1), op2);
242 }
243
244 void
mpr_div(mpr * rop,mpr * op1,mpr * op2)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
mpr_divi(mpr * rop,mpr * op1,long op2)266 mpr_divi(mpr *rop, mpr *op1, long op2)
267 {
268 mpi_muli(mpr_den(rop), mpr_den(op1), op2);
269 }
270
271 void
mpr_inv(mpr * rop,mpr * op)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
mpr_neg(mpr * rop,mpr * op)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
mpr_abs(mpr * rop,mpr * op)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
mpr_cmp(mpr * op1,mpr * op2)305 mpr_cmp(mpr *op1, mpr *op2)
306 {
307 return (mpr_docmp(op1, op2, 1));
308 }
309
310 int
mpr_cmpi(mpr * op1,long op2)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
mpr_cmpabs(mpr * op1,mpr * op2)325 mpr_cmpabs(mpr *op1, mpr *op2)
326 {
327 return (mpr_docmp(op1, op2, 0));
328 }
329
330 int
mpr_cmpabsi(mpr * op1,long op2)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
mpr_docmp(mpr * op1,mpr * op2,int sign)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
mpr_swap(mpr * op1,mpr * op2)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
mpr_fiti(mpr * op)407 mpr_fiti(mpr *op)
408 {
409 return (mpi_fiti(mpr_num(op)) && mpi_fiti(mpr_den(op)));
410 }
411
412 double
mpr_getd(mpr * op)413 mpr_getd(mpr *op)
414 {
415 return (mpi_getd(mpr_num(op)) / mpi_getd(mpr_den(op)));
416 }
417
418 char *
mpr_getstr(char * str,mpr * op,int base)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