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 */
45static void mpr_addsub(mpr *rop, mpr *op1, mpr *op2, int sub);
46
47	/* do the hard work of mpr_cmp and mpr_cmpabs */
48static int mpr_docmp(mpr *op1, mpr *op2, int sign);
49
50/*
51 * Implementation
52 */
53void
54mpr_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
67void
68mpr_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
79void
80mpr_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
88void
89mpr_seti(mpr *rop, long num, long den)
90{
91    mpi_seti(mpr_num(rop), num);
92    mpi_seti(mpr_den(rop), den);
93}
94
95void
96mpr_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
126void
127mpr_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
138void
139mpr_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
159void
160mpr_add(mpr *rop, mpr *op1, mpr *op2)
161{
162    mpr_addsub(rop, op1, op2, 0);
163}
164
165void
166mpr_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
177void
178mpr_sub(mpr *rop, mpr *op1, mpr *op2)
179{
180    mpr_addsub(rop, op1, op2, 1);
181}
182
183void
184mpr_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
195static void
196mpr_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
217void
218mpr_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
238void
239mpr_muli(mpr *rop, mpr *op1, long op2)
240{
241    mpi_muli(mpr_num(rop), mpr_num(op1), op2);
242}
243
244void
245mpr_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
265void
266mpr_divi(mpr *rop, mpr *op1, long op2)
267{
268    mpi_muli(mpr_den(rop), mpr_den(op1), op2);
269}
270
271void
272mpr_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
282void
283mpr_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
289void
290mpr_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
304int
305mpr_cmp(mpr *op1, mpr *op2)
306{
307    return (mpr_docmp(op1, op2, 1));
308}
309
310int
311mpr_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
324int
325mpr_cmpabs(mpr *op1, mpr *op2)
326{
327    return (mpr_docmp(op1, op2, 0));
328}
329
330int
331mpr_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
344static int
345mpr_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
394void
395mpr_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
406int
407mpr_fiti(mpr *op)
408{
409    return (mpi_fiti(mpr_num(op)) && mpi_fiti(mpr_den(op)));
410}
411
412double
413mpr_getd(mpr *op)
414{
415    return (mpi_getd(mpr_num(op)) / mpi_getd(mpr_den(op)));
416}
417
418char *
419mpr_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