SecBSD/xenocara
2
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

sync code with last improvements from OpenBSD

Browse source
This commit is contained in:
purplerain 2023-08-28 05:57:34 +00:00
commit 88965415ff
Signed by: purplerain
GPG key ID: F42C07F07E2E35B7
26235 changed files with 29195616 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

822
app/xedit/lisp/mp/mp.c Normal file
View file

@ -0,0 +1,822 @@
/*
* Copyright (c) 2002 by The XFree86 Project, Inc.
*
* Permission is hereby granted, free of charge, to any person obtaining a
* copy of this software and associated documentation files (the "Software"),
* to deal in the Software without restriction, including without limitation
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
* and/or sell copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included in
* all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
* THE XFREE86 PROJECT BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF
* OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
* SOFTWARE.
*
* Except as contained in this notice, the name of the XFree86 Project shall
* not be used in advertising or otherwise to promote the sale, use or other
* dealings in this Software without prior written authorization from the
* XFree86 Project.
*
* Author: Paulo César Pereira de Andrade
*/
/* $XFree86: xc/programs/xedit/lisp/mp/mp.c,v 1.2 2002/11/08 08:01:00 paulo Exp $ */
#include "mp.h"
/*
* TODO:
* o Optimize squaring
* o Write better division code and move from mpi.c to here
* o Make multiplication code don't required memory to be zeroed
* + The first step is easy, just multiply the low word,
* then the high word, that may overlap with the result
* of the first multiply (in case of carry), and then
* just make sure carry is properly propagated in the
* subsequent multiplications.
* + Some code needs also to be rewritten because some
* intermediate addition code in mp_mul, mp_karatsuba_mul,
* and mp_toom_mul is assuming the memory is zeroed.
*/
/*
* Prototypes
*/
/* out of memory handler */
static void mp_outmem(void);
/* memory allocation fallback functions */
static void *_mp_malloc(size_t);
static void *_mp_calloc(size_t, size_t);
static void *_mp_realloc(void*, size_t);
static void _mp_free(void*);
/*
* Initialization
*/
static mp_malloc_fun __mp_malloc = _mp_malloc;
static mp_calloc_fun __mp_calloc = _mp_calloc;
static mp_realloc_fun __mp_realloc = _mp_realloc;
static mp_free_fun __mp_free = _mp_free;
/*
* Implementation
*/
static void
mp_outmem(void)
{
fprintf(stderr, "out of memory in MP library.\n");
exit(1);
}
static void *
_mp_malloc(size_t size)
{
return (malloc(size));
}
void *
mp_malloc(size_t size)
{
void *pointer = (*__mp_malloc)(size);
if (pointer == NULL)
mp_outmem();
return (pointer);
}
mp_malloc_fun
mp_set_malloc(mp_malloc_fun fun)
{
mp_malloc_fun old = __mp_malloc;
__mp_malloc = fun;
return (old);
}
static void *
_mp_calloc(size_t nmemb, size_t size)
{
return (calloc(nmemb, size));
}
void *
mp_calloc(size_t nmemb, size_t size)
{
void *pointer = (*__mp_calloc)(nmemb, size);
if (pointer == NULL)
mp_outmem();
return (pointer);
}
mp_calloc_fun
mp_set_calloc(mp_calloc_fun fun)
{
mp_calloc_fun old = __mp_calloc;
__mp_calloc = fun;
return (old);
}
static void *
_mp_realloc(void *old, size_t size)
{
return (realloc(old, size));
}
void *
mp_realloc(void *old, size_t size)
{
void *pointer = (*__mp_realloc)(old, size);
if (pointer == NULL)
mp_outmem();
return (pointer);
}
mp_realloc_fun
mp_set_realloc(mp_realloc_fun fun)
{
mp_realloc_fun old = __mp_realloc;
__mp_realloc = fun;
return (old);
}
static void
_mp_free(void *pointer)
{
free(pointer);
}
void
mp_free(void *pointer)
{
(*__mp_free)(pointer);
}
mp_free_fun
mp_set_free(mp_free_fun fun)
{
mp_free_fun old = __mp_free;
__mp_free = fun;
return (old);
}
long
mp_add(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2)
{
BNI value; /* intermediate result */
BNS carry; /* carry flag */
long size; /* result size */
if (len1 < len2)
MP_SWAP(op1, op2, len1, len2);
/* unroll start of loop */
value = (BNI)op1[0] + op2[0];
rop[0] = value;
carry = value >> BNSBITS;
/* add op1 and op2 */
for (size = 1; size < len2; size++) {
value = (BNI)op1[size] + op2[size] + carry;
rop[size] = value;
carry = value >> BNSBITS;
}
if (rop != op1) {
for (; size < len1; size++) {
value = (BNI)op1[size] + carry;
rop[size] = value;
carry = value >> BNSBITS;
}
}
else {
/* if rop == op1, than just adjust carry */
for (; carry && size < len1; size++) {
value = (BNI)op1[size] + carry;
rop[size] = value;
carry = value >> BNSBITS;
}
size = len1;
}
if (carry)
rop[size++] = carry;
return (size);
}
long
mp_sub(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2)
{
long svalue; /* intermediate result */
BNS carry; /* carry flag */
long size; /* result size */
/* special case */
if (op1 == op2) {
rop[0] = 0;
return (1);
}
/* unroll start of loop */
svalue = (long)op1[0] - op2[0];
rop[0] = svalue;
carry = svalue < 0;
/* subtracts op2 from op1 */
for (size = 1; size < len2; size++) {
svalue = (long)(op1[size]) - op2[size] - carry;
rop[size] = svalue;
carry = svalue < 0;
}
if (rop != op1) {
for (; size < len1; size++) {
svalue = op1[size] - carry;
rop[size] = svalue;
carry = svalue < 0;
}
}
else {
/* if rop == op1, than just adjust carry */
for (; carry && size < len1; size++) {
svalue = (long)op1[size] - carry;
rop[size] = svalue;
carry = svalue < 0;
}
size = len1;
}
/* calculate result size */
while (size > 1 && rop[size - 1] == 0)
--size;
return (size);
}
long
mp_lshift(BNS *rop, BNS *op, BNI len, long shift)
{
long i, size;
BNI words, bits; /* how many word and bit shifts */
words = shift / BNSBITS;
bits = shift % BNSBITS;
size = len + words;
if (bits) {
BNS hi, lo;
BNI carry;
int adj;
for (i = 1, carry = CARRY >> 1; carry; i++, carry >>= 1)
if (op[len - 1] & carry)
break;
adj = (bits + (BNSBITS - i)) / BNSBITS;
size += adj;
lo = hi = op[0];
rop[words] = lo << bits;
for (i = 1; i < len; i++) {
hi = op[i];
rop[words + i] = hi << bits | (lo >> (BNSBITS - bits));
lo = hi;
}
if (adj)
rop[size - 1] = hi >> (BNSBITS - bits);
}
else
memmove(rop + size - len, op, sizeof(BNS) * len);
if (words)
memset(rop, '\0', sizeof(BNS) * words);
return (size);
}
long
mp_rshift(BNS *rop, BNS *op, BNI len, long shift)
{
int adj = 0;
long i, size;
BNI words, bits; /* how many word and bit shifts */
words = shift / BNSBITS;
bits = shift % BNSBITS;
size = len - words;
if (bits) {
BNS hi, lo;
BNI carry;
for (i = 0, carry = CARRY >> 1; carry; i++, carry >>= 1)
if (op[len - 1] & carry)
break;
adj = (bits + i) / BNSBITS;
if (size - adj == 0) {
rop[0] = 0;
return (1);
}
hi = lo = op[words + size - 1];
rop[size - 1] = hi >> bits;
for (i = size - 2; i >= 0; i--) {
lo = op[words + i];
rop[i] = (lo >> bits) | (hi << (BNSBITS - bits));
hi = lo;
}
if (adj)
rop[0] |= lo << (BNSBITS - bits);
}
else
memmove(rop, op + len - size, size * sizeof(BNS));
return (size - adj);
}
/* rop must be a pointer to len1 + len2 elements
* rop cannot be either op1 or op2
* rop must be all zeros */
long
mp_base_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2)
{
long i, j; /* counters */
BNI value; /* intermediate result */
BNS carry; /* carry value */
long size = len1 + len2;
/* simple optimization: first pass does not need to deference rop[i+j] */
if (op1[0]) {
value = (BNI)(op1[0]) * op2[0];
rop[0] = value;
carry = (BNS)(value >> BNSBITS);
for (j = 1; j < len2; j++) {
value = (BNI)(op1[0]) * op2[j] + carry;
rop[j] = value;
carry = (BNS)(value >> BNSBITS);
}
rop[j] = carry;
}
/* do the multiplication */
for (i = 1; i < len1; i++) {
if (op1[i]) {
/* unrool loop initialization */
value = (BNI)(op1[i]) * op2[0] + rop[i];
rop[i] = value;
carry = (BNS)(value >> BNSBITS);
/* multiply */
for (j = 1; j < len2; j++) {
value = (BNI)(op1[i]) * op2[j] + rop[i + j] + carry;
rop[i + j] = value;
carry = (BNS)(value >> BNSBITS);
}
rop[i + j] = carry;
}
}
if (size > 1 && rop[size - 1] == 0)
--size;
return (size);
}
/* Karatsuba method
* t + ((a0 + a1) (b0 + b1) - t - u) x + ux²
* where t = a0b0 and u = a1b1
*
* Karatsuba method reduces the number of multiplications. Example:
* Square a 40 length number
* instead of a plain 40*40 = 1600 multiplies/adds, it does:
* 20*20+20*20+20*20 = 1200
* but since it is recursive, every 20*20=400 is reduced to
* 10*10+10*10+10*10=300
* and so on.
* The multiplication by x and x² is a just a shift, as it is a
* power of two, and is implemented below by just writting at the
* correct offset */
long
mp_karatsuba_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2)
{
BNI x; /* shift count */
BNI la0, la1, lb0, lb1; /* length of a0, a1, b0, and b1 */
BNS *t; /* temporary memory for t product */
BNS *u; /* temporary memory for u product */
BNS *r; /* pointer to rop */
long xlen, tlen, ulen;
/* calculate value of x, that is 2^(BNSBITS*x) */
if (len1 >= len2)
x = (len1 + 1) >> 1;
else
x = (len2 + 1) >> 1;
/* calculate length of operands */
la0 = x;
la1 = len1 - x;
lb0 = x;
lb1 = len2 - x;
/* allocate buffer for t and (a0 + a1) */
tlen = la0 + lb0;
t = mp_malloc(sizeof(BNS) * tlen);
/* allocate buffer for u and (b0 + b1) */
if (la1 + lb1 < lb0 + lb1 + 1)
ulen = lb0 + lb1 + 1;
else
ulen = la1 + lb1;
u = mp_malloc(sizeof(BNS) * ulen);
/* calculate a0 + a1, store result in t */
tlen = mp_add(t, op1, op1 + x, la0, la1);
/* calculate b0 + b1, store result in u */
ulen = mp_add(u, op2, op2 + x, lb0, lb1);
/* store (a0 + a1) * (b0 + b1) in rop */
r = rop + x; /* multiplied by 2^(BNSBITS*x) */
xlen = mp_mul(r, t, u, tlen, ulen);
/* must zero t and u memory, this is required for mp_mul */
/* calculate t = a0 * b0 */
tlen = la0 + lb0;
memset(t, '\0', sizeof(BNS) * tlen);
tlen = mp_mul(t, op1, op2, la0, lb0);
/* calculate u = a1 * b1 */
ulen = la1 + lb1;
memset(u, '\0', sizeof(BNS) * ulen);
ulen = mp_mul(u, op1 + x, op2 + x, la1, lb1);
/* subtract t from partial result */
xlen = mp_sub(r, r, t, xlen, tlen);
/* subtract u form partial result */
xlen = mp_sub(r, r, u, xlen, ulen);
/* add ux^2 to partial result */
r = rop + (x << 1); /* multiplied by x^2 = 2^(BNSBITS*x*2) */
xlen = len1 + len2;
xlen = mp_add(r, r, u, xlen, ulen);
/* now add t to final result */
xlen = mp_add(rop, rop, t, xlen, tlen);
mp_free(t);
mp_free(u);
if (xlen > 1 && rop[xlen - 1] == 0)
--xlen;
return (xlen);
}
/* Toom method (partially based on GMP documentation)
* Evaluation at k = [ 0 1/2 1 2 oo ]
* U(x) = (U2k + U1)k + U0
* V(x) = (V2k + V1)k + V0
* W(x) = U(x)V(x)
*
* Sample:
* 123 * 456
*
* EVALUATION:
* U(0) = (1*0+2)*0+3 => 3
* U(1) = 1+(2+3*2)*2 => 17
* U(2) = 1+2+3 => 6
* U(3) = (1*2+2)*2+3 => 11
* U(4) = 1+(2+3*0)*0 => 1
*
* V(0) = (4*0+5)*0+6 => 6
* V(1) = 4+(5+6*2)*2 => 38
* V(2) = 4+5+6 => 15
* V(3) = (4*2+5)*2+6 => 32
* V(4) = 4+(5+6*0)*0 => 4
*
* U = [ 3 17 6 11 1 ]
* V = [ 6 38 15 32 4 ]
* W = [ 18 646 90 352 4 ]
*
* After that, we have:
* a = 18 (w0 already known)
* b = 16w0 + 8w1 + 4w2 + 2w3 + w4
* c = w0 + w1 + w2 + w3 + w4
* d = w0 + 2w1 + 4w2 + 8w3 + 16w4
* e = 4 (w4 already known)
*
* INTERPOLATION:
* b = b -16a - e (354)
* c = c - a - e (68)
* d = d - a - 16e (270)
*
* w = (b + d) - 8c = (10w1+8w2+10w3) - (8w1+8w2+8w3) = 2w1+2w3
* w = 2c - w (56)
* b = b/2 = 4w1+w+w3
* b = b-c = 4w1+w+w3 - w1+w2+w3 = 3w1+w2
* c = w/2 (w2 = 28)
* b = b-c = 3w1+c - c = 3w1
* b = b/3 (w1 = 27)
* d = d/2
* d = d-b-w = b+w+4w3 - b-w = 4w3
* d = d/4 (w3 = 13)
*
* RESULT:
* w4*10^4 + w3*10³ + w2*10² + w1*10 + w0
* 40000 + 13000 + 2800 + 270 + 18
* 10 is the base where the calculation was done
*
* This sample uses small numbers, so it does not show the
* advantage of the method. But for example (in base 10), when squaring
* 123456789012345678901234567890
* The normal method would do 30*30=900 multiplications
* Karatsuba method would do 15*15*3=675 multiplications
* Toom method would do 10*10*5=500 multiplications
* Toom method has a larger overhead if compared with Karatsuba method,
* due to evaluation and interpolation, so it should be used for larger
* numbers, so that the computation time of evaluation/interpolation
* would be smaller than the time spent using other methods.
*
* Note that Karatsuba method can be seen as a special case of
* Toom method, i.e:
* U1U0 * V1V0
* with k = [ 0 1 oo ]
* U = [ U0 U1+U0 U1 ]
* V = [ V0 V1+V0 V1 ]
* W = [ U0*V0 (U1+U0)*(V1+V0) (U1+V1) ]
*
* w0 = U0*V0
* w = (U1+U0)*(V1+V0)
* w2 = (U1*V1)
*
* w1 = w - w0 - w2
* w2x² + w1x + w0
*
* See Knuth's Seminumerical Algorithms for a sample implemention
* using 4 stacks and k = [ 0 1 2 3 ... ], based on the size of the
* input.
*/
long
mp_toom_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2)
{
long size, xsize, i;
BNI value; /* used in division */
BNS carry;
BNI x; /* shift count */
BNI l1, l2;
BNI al, bl, cl, dl, el, Ul[3], Vl[3];
BNS *a, *b, *c, *d, *e, *U[3], *V[3];
/* x is the base i.e. 2^(BNSBITS*x) */
x = (len1 + len2 + 4) / 6;
l1 = len1 - (x << 1); /* length of remaining piece of op1 */
l2 = len2 - (x << 1); /* length of remaining piece of op2 */
/* allocate memory for storing U and V */
U[0] = mp_malloc(sizeof(BNS) * (x + 2));
V[0] = mp_malloc(sizeof(BNS) * (x + 2));
U[1] = mp_malloc(sizeof(BNS) * (x + 1));
V[1] = mp_malloc(sizeof(BNS) * (x + 1));
U[2] = mp_malloc(sizeof(BNS) * (x + 2));
V[2] = mp_malloc(sizeof(BNS) * (x + 2));
/* EVALUATE U AND V */
/* Numbers are in the format U2x²+U1x+U0 and V2x²+V1x+V0 */
/* U[0] = U2+U1*2+U0*4 */
/* store U1*2 in U[1], this value is used twice */
Ul[1] = mp_lshift(U[1], op1 + x, x, 1);
/* store U0*4 in U[0] */
Ul[0] = mp_lshift(U[0], op1, x, 2);
/* add U1*2 to U[0] */
Ul[0] = mp_add(U[0], U[0], U[1], Ul[0], Ul[1]);
/* add U2 to U[0] */
Ul[0] = mp_add(U[0], U[0], op1 + x + x, Ul[0], l1);
/* U[2] = U2*4+U1*2+U0 */
/* store U2*4 in U[2] */
Ul[2] = mp_lshift(U[2], op1 + x + x, l1, 2);
/* add U1*2 to U[2] */
Ul[2] = mp_add(U[2], U[2], U[1], Ul[2], Ul[1]);
/* add U0 to U[2] */
Ul[2] = mp_add(U[2], U[2], op1, Ul[2], x);
/* U[1] = U2+U1+U0 */
Ul[1] = mp_add(U[1], op1, op1 + x, x, x);
Ul[1] = mp_add(U[1], U[1], op1 + x + x, Ul[1], l1);
/* Evaluate V[x], same code as U[x] */
Vl[1] = mp_lshift(V[1], op2 + x, x, 1);
Vl[0] = mp_lshift(V[0], op2, x, 2);
Vl[0] = mp_add(V[0], V[0], V[1], Vl[0], Vl[1]);
Vl[0] = mp_add(V[0], V[0], op2 + x + x, Vl[0], l2);
Vl[2] = mp_lshift(V[2], op2 + x + x, l2, 2);
Vl[2] = mp_add(V[2], V[2], V[1], Vl[2], Vl[1]);
Vl[2] = mp_add(V[2], V[2], op2, Vl[2], x);
Vl[1] = mp_add(V[1], op2, op2 + x, x, x);
Vl[1] = mp_add(V[1], V[1], op2 + x + x, Vl[1], l2);
/* MULTIPLY U[] AND V[] */
/* calculate (U2+U1*2+U0*4) * (V2+V1*2+V0*4) */
b = mp_calloc(1, sizeof(BNS) * (Ul[0] * Vl[0]));
bl = mp_mul(b, U[0], V[0], Ul[0], Vl[0]);
mp_free(U[0]);
mp_free(V[0]);
/* calculate (U2+U1+U0) * (V2+V1+V0) */
c = mp_calloc(1, sizeof(BNS) * (Ul[1] * Vl[1]));
cl = mp_mul(c, U[1], V[1], Ul[1], Vl[1]);
mp_free(U[1]);
mp_free(V[1]);
/* calculate (U2*4+U1*2+U0) * (V2*4+V1*2+V0) */
d = mp_calloc(1, sizeof(BNS) * (Ul[2] * Vl[2]));
dl = mp_mul(d, U[2], V[2], Ul[2], Vl[2]);
mp_free(U[2]);
mp_free(V[2]);
/* calculate U0 * V0 */
a = mp_calloc(1, sizeof(BNS) * (x + x));
al = mp_mul(a, op1, op2, x, x);
/* calculate U2 * V2 */
e = mp_calloc(1, sizeof(BNS) * (l1 + l2));
el = mp_mul(e, op1 + x + x, op2 + x + x, l1, l2);
/* INTERPOLATE COEFFICIENTS */
/* b = b - 16a - e */
size = mp_lshift(rop, a, al, 4);
bl = mp_sub(b, b, rop, bl, size);
bl = mp_sub(b, b, e, bl, el);
/* c = c - a - e*/
cl = mp_sub(c, c, a, cl, al);
cl = mp_sub(c, c, e, cl, el);
/* d = d - a - 16e */
dl = mp_sub(d, d, a, dl, al);
size = mp_lshift(rop, e, el, 4);
dl = mp_sub(d, d, rop, dl, size);
/* w = (b + d) - 8c */
size = mp_add(rop, b, d, bl, dl);
xsize = mp_lshift(rop + size, c, cl, 3); /* rop has enough storage */
size = mp_sub(rop, rop, rop + size, size, xsize);
/* w = 2c - w*/
xsize = mp_lshift(rop + size, c, cl, 1);
size = mp_sub(rop, rop + size, rop, xsize, size);
/* b = b/2 */
bl = mp_rshift(b, b, bl, 1);
/* b = b - c */
bl = mp_sub(b, b, c, bl, cl);
/* c = w / 2 */
cl = mp_rshift(c, rop, size, 1);
/* b = b - c */
bl = mp_sub(b, b, c, bl, cl);
/* b = b/3 */
/* maybe the most expensive calculation */
i = bl - 1;
value = b[i];
b[i] = value / 3;
for (--i; i >= 0; i--) {
carry = value % 3;
value = ((BNI)carry << BNSBITS) + b[i];
b[i] = (BNS)(value / 3);
}
/* d = d/2 */
dl = mp_rshift(d, d, dl, 1);
/* d = d - b - w */
dl = mp_sub(d, d, b, dl, bl);
dl = mp_sub(d, d, rop, dl, size);
/* d = d/4 */
dl = mp_rshift(d, d, dl, 2);
/* STORE RESULT IN ROP */
/* first clear memory used as temporary variable w and 8c */
memset(rop, '\0', sizeof(BNS) * (len1 + len2));
i = x * 4;
xsize = (len1 + len2) - i;
size = mp_add(rop + i, rop + i, e, xsize, el) + i;
i = x * 3;
xsize = size - i;
size = mp_add(rop + i, rop + i, d, xsize, dl) + i;
i = x * 2;
xsize = size - i;
size = mp_add(rop + i, rop + i, c, xsize, cl) + i;
i = x;
xsize = size - i;
size = mp_add(rop + i, rop + i, b, xsize, bl) + i;
size = mp_add(rop, rop, a, size, al);
mp_free(e);
mp_free(d);
mp_free(c);
mp_free(b);
mp_free(a);
if (size > 1 && rop[size - 1] == 0)
--size;
return (size);
}
long
mp_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2)
{
if (len1 < len2)
MP_SWAP(op1, op2, len1, len2);
if (len1 < KARATSUBA || len2 < KARATSUBA)
return (mp_base_mul(rop, op1, op2, len1, len2));
else if (len1 < TOOM && len2 < TOOM && len2 > ((len1 + 1) >> 1))
return (mp_karatsuba_mul(rop, op1, op2, len1, len2));
else if (len1 >= TOOM && len2 >= TOOM && (len2 + 2) / 3 == (len1 + 2) / 3)
return (mp_toom_mul(rop, op1, op2, len1, len2));
else {
long xsize, psize, isize;
BNS *ptr;
/* adjust index pointer and estimated size of result */
isize = 0;
xsize = len1 + len2;
mp_mul(rop, op1, op2, len2, len2);
/* adjust pointers */
len1 -= len2;
op1 += len2;
/* allocate buffer for intermediate multiplications */
if (len1 > len2)
ptr = mp_calloc(1, sizeof(BNS) * (len2 + len2));
else
ptr = mp_calloc(1, sizeof(BNS) * (len1 + len2));
/* loop multiplying len2 size operands at a time */
while (len1 >= len2) {
isize += len2;
psize = mp_mul(ptr, op1, op2, len2, len2);
mp_add(rop + isize, rop + isize, ptr, xsize - isize, psize);
len1 -= len2;
op1 += len2;
/* multiplication routines require zeroed memory */
memset(ptr, '\0', sizeof(BNS) * (MIN(len1, len2) + len2));
}
/* len1 was not a multiple of len2 */
if (len1) {
isize += len2;
psize = mp_mul(ptr, op2, op1, len2, len1);
mp_add(rop + isize, rop + isize, ptr, xsize, psize);
}
/* adjust result size */
if (rop[xsize - 1] == 0)
--xsize;
mp_free(ptr);
return (xsize);
}
}

437
app/xedit/lisp/mp/mp.h Normal file
View file

@ -0,0 +1,437 @@
/*
* Copyright (c) 2002 by The XFree86 Project, Inc.
*
* Permission is hereby granted, free of charge, to any person obtaining a
* copy of this software and associated documentation files (the "Software"),
* to deal in the Software without restriction, including without limitation
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
* and/or sell copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included in
* all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
* THE XFREE86 PROJECT BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF
* OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
* SOFTWARE.
*
* Except as contained in this notice, the name of the XFree86 Project shall
* not be used in advertising or otherwise to promote the sale, use or other
* dealings in this Software without prior written authorization from the
* XFree86 Project.
*
* Author: Paulo César Pereira de Andrade
*/
/* $XFree86: xc/programs/xedit/lisp/mp/mp.h,v 1.5tsi Exp $ */
#include <stdio.h>
#include <math.h>
#ifdef sun
#include <ieeefp.h>
#endif
#include <float.h>
#include <stdlib.h>
#include <limits.h>
#include <ctype.h>
#include <string.h>
#ifndef __mp_h_
#define __mp_h_
#ifdef __GNUC__
#define INLINE __inline__
#else
#define INLINE /**/
#endif
/* this normally is better for multiplication and also
* simplify addition loops putting the larger value first */
#define MP_SWAP(op1, op2, len1, len2) { \
BNS *top = op1; \
BNI tlen = len1; \
\
op1 = op2; \
len1 = len2; \
op2 = top; \
len2 = tlen; \
}
/*
* At least this length to use Karatsuba multiplication method
*/
#define KARATSUBA 32
/*
* At least this length to use Toom multiplication method
*/
#define TOOM 128
#if ULONG_MAX > 4294967295UL
/* sizeof(long) == 8 and sizeof(int) == 4 */
# define BNI unsigned long
# define BNS unsigned int
# define MINSLONG 0x8000000000000000UL
# define MAXSLONG 0x7fffffffffffffffUL
# define CARRY 0x100000000
# define LMASK 0xffffffff00000000UL
# define SMASK 0x00000000ffffffffUL
# define BNIBITS 64
# define BNSBITS 32
# ifndef LONG64
# define LONG64
# endif
#else
/* sizeof(long) == 4 and sizeof(short) == 2 */
# define BNI unsigned long
# define BNS unsigned short
# define MINSLONG 0x80000000UL
# define MAXSLONG 0x7fffffffUL
# define CARRY 0x10000
# define LMASK 0xffff0000UL
# define SMASK 0x0000ffffUL
# define BNIBITS 32
# define BNSBITS 16
#endif
#ifdef MAX
#undef MAX
#endif
#define MAX(a, b) ((a) > (b) ? (a) : (b))
#ifdef MIN
#undef MIN
#endif
#define MIN(a, b) ((a) < (b) ? (a) : (b))
/*
* Types
*/
typedef struct _mpi {
unsigned int size : 31;
unsigned int sign : 1;
BNI alloc;
BNS *digs; /* LSF format */
} mpi;
typedef struct _mpr {
mpi num;
mpi den;
} mpr;
typedef void *(*mp_malloc_fun)(size_t);
typedef void *(*mp_calloc_fun)(size_t, size_t);
typedef void *(*mp_realloc_fun)(void*, size_t);
typedef void (*mp_free_fun)(void*);
/*
* Prototypes
*/
/* GENERIC FUNCTIONS */
/* memory allocation wrappers */
void *mp_malloc(size_t size);
void *mp_calloc(size_t nmemb, size_t size);
void *mp_realloc(void *pointer, size_t size);
void mp_free(void *pointer);
mp_malloc_fun mp_set_malloc(mp_malloc_fun);
mp_calloc_fun mp_set_calloc(mp_calloc_fun);
mp_realloc_fun mp_set_realloc(mp_realloc_fun);
mp_free_fun mp_set_free(mp_free_fun);
/* adds op1 and op2, stores result in rop
* rop must pointer to at least len1 + len2 + 1 elements
* rop can be either op1 or op2 */
long mp_add(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2);
/* subtracts op2 from op1, stores result in rop
* rop must pointer to at least len1 + len2 elements
* op1 must be >= op2
* rop can be either op1 or op2 */
long mp_sub(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2);
/* shift op to the left shift bits
* rop must have enough storage for result
* rop can be op */
long mp_lshift(BNS *rop, BNS *op, BNI len, long shift);
/* shift op to the right shift bits
* shift must be positive
* rop can be op */
long mp_rshift(BNS *rop, BNS *op, BNI len, long shift);
/* use simple generic multiplication method
* rop cannot be the same as op1 or op2
* rop must be zeroed
* op1 can be op2 */
long mp_base_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2);
/* use Karatsuba method
* MIN(len1, len2) must be larger than (MAX(len1, len2) + 1) >> 1
* MAX(len1, len2) should be at least 2
* rop cannot be the same as op1 or op2
* rop must be zeroed
* op1 can be op2 */
long mp_karatsuba_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2);
/* use Toom method
* len1 / 3 should be equal to len2 / 3
* len1 / 3 should be at least 1
* rop cannot be the same as op1 or op2
* rop must be zeroed
* op1 can be op2 */
long mp_toom_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2);
/* chooses the available multiplication methods based on it's input
* rop must be a pointer to len1 + len2 elements
* rop cannot be the same as op1 or op2
* rop must be zeroed
* op1 can be op2 */
long mp_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2);
/* INTEGER FUNCTIONS */
/* initialize op and set it to 0 */
void mpi_init(mpi *op);
/* clear memory associated to op */
void mpi_clear(mpi *op);
/* set rop to the value of op */
void mpi_set(mpi *rop, mpi *op);
/* set rop to the value of si */
void mpi_seti(mpi *rop, long si);
/* set rop to the floor(fabs(d)) */
void mpi_setd(mpi *rop, double d);
/* initialize rop to number representation in str in the given base.
* leading zeros are skipped.
* if sign present, it is processed.
* base must be in the range 2 to 36. */
void mpi_setstr(mpi *rop, char *str, int base);
/* adds two mp integers */
void mpi_add(mpi *rop, mpi *op1, mpi *op2);
/* adds op1 and op2 */
void mpi_addi(mpi *rop, mpi *op1, long op2);
/* subtracts two mp integers */
void mpi_sub(mpi *rop, mpi *op1, mpi *op2);
/* subtracts op2 from op1 */
void mpi_subi(mpi *rop, mpi *op1, long op2);
/* multiply two mp integers */
void mpi_mul(mpi *rop, mpi *op1, mpi *op2);
/* multiply op1 by op2 */
void mpi_muli(mpi *rop, mpi *op1, long op2);
/* divides num by den and sets rop to result */
void mpi_div(mpi *rop, mpi *num, mpi *den);
/* divides num by den and sets rop to the remainder */
void mpi_rem(mpi *rop, mpi *num, mpi *den);
/* divides num by den, sets quotient to qrop and remainder to rrop
* qrop is truncated towards zero.
* qrop and rrop are optional
* qrop and rrop cannot be the same variable */
void mpi_divqr(mpi *qrop, mpi *rrop, mpi *num, mpi *den);
/* divides num by then and stores result in rop */
void mpi_divi(mpi *rop, mpi *num, long den);
/* divides num by den and returns remainder */
long mpi_remi(mpi *num, long den);
/* divides num by den
* stores quotient in qrop and returns remainder */
long mpi_divqri(mpi *qrop, mpi *num, long den);
/* sets rop to num modulo den */
void mpi_mod(mpi *rop, mpi *num, mpi *den);
/* returns num modulo den */
long mpi_modi(mpi *num, long den);
/* sets rop to the greatest common divisor of num and den
* result is always positive */
void mpi_gcd(mpi *rop, mpi *num, mpi *den);
/* sets rop to the least common multiple of num and den
* result is always positive */
void mpi_lcm(mpi *rop, mpi *num, mpi *den);
/* sets rop to op raised to exp */
void mpi_pow(mpi *rop, mpi *op, unsigned long exp);
/* sets rop to the integer part of the nth root of op.
* returns 1 if result is exact, 0 otherwise */
int mpi_root(mpi *rop, mpi *op, unsigned long nth);
/* sets rop to the integer part of the square root of op.
* returns 1 if result is exact, 0 otherwise */
int mpi_sqrt(mpi *rop, mpi *op);
/* bit shift, left if shift positive, right if negative
* a fast way to multiply and divide by powers of two */
void mpi_ash(mpi *rop, mpi *op, long shift);
/* sets rop to op1 logand op2 */
void mpi_and(mpi *rop, mpi *op1, mpi *op2);
/* sets rop to op1 logior op2 */
void mpi_ior(mpi *rop, mpi *op1, mpi *op2);
/* sets rop to op1 logxor op2 */
void mpi_xor(mpi *rop, mpi *op1, mpi *op2);
/* sets rop to one's complement of op */
void mpi_com(mpi *rop, mpi *op);
/* sets rop to -op */
void mpi_neg(mpi *rop, mpi *op);
/* sets rop to the absolute value of op */
void mpi_abs(mpi *rop, mpi *op);
/* compares op1 and op2
* returns >0 if op1 > op2, 0 if op1 = op2, and <0 if op1 < op2 */
int mpi_cmp(mpi *op1, mpi *op2);
/* mpi_cmp with a long integer operand */
int mpi_cmpi(mpi *op1, long op2);
/* compares absolute value of op1 and op2
* returns >0 if abs(op1) > abs(op2), 0 if abs(op1) = abs(op2),
* and <0 if abs(op1) < abs(op2) */
int mpi_cmpabs(mpi *op1, mpi *op2);
/* mpi_cmpabs with a long integer operand */
int mpi_cmpabsi(mpi *op1, long op2);
/* returns 1 if op1 > 0, 0 if op1 = 0, and -1 if op1 < 0 */
int mpi_sgn(mpi *op);
/* fastly swaps contents of op1 and op2 */
void mpi_swap(mpi *op1, mpi *op2);
/* returns 1 if op fits in a signed long int, 0 otherwise */
int mpi_fiti(mpi *op);
/* converts mp integer to long int
* to know if the value will fit, call mpi_fiti */
long mpi_geti(mpi *op);
/* convert mp integer to double */
double mpi_getd(mpi *op);
/* returns exact number of characters to represent mp integer
* in given base, excluding sign and ending null character.
* base must be in the range 2 to 36 */
unsigned long mpi_getsize(mpi *op, int base);
/* returns pointer to string with representation of mp integer
* if str is not NULL, it must have enough space to store integer
* representation, if NULL a newly allocated string is returned.
* base must be in the range 2 to 36 */
char *mpi_getstr(char *str, mpi *op, int base);
/* RATIO FUNCTIONS */
#define mpr_num(op) (&((op)->num))
#define mpr_den(op) (&((op)->den))
/* initialize op and set it to 0/1 */
void mpr_init(mpr *op);
/* clear memory associated to op */
void mpr_clear(mpr *op);
/* set rop to the value of op */
void mpr_set(mpr *rop, mpr *op);
/* set rop to num/den */
void mpr_seti(mpr *rop, long num, long den);
/* set rop to the value of d */
void mpr_setd(mpr *rop, double d);
/* initialize rop to number representation in str in the given base.
* leading zeros are skipped.
* if sign present, it is processed.
* base must be in the range 2 to 36. */
void mpr_setstr(mpr *rop, char *str, int base);
/* remove common factors of op */
void mpr_canonicalize(mpr *op);
/* adds two mp rationals */
void mpr_add(mpr *rop, mpr *op1, mpr *op2);
/* adds op1 and op2 */
void mpr_addi(mpr *rop, mpr *op1, long op2);
/* subtracts two mp rationals */
void mpr_sub(mpr *rop, mpr *op1, mpr *op2);
/* subtracts op2 from op1 */
void mpr_subi(mpr *rop, mpr *op1, long op2);
/* multiply two mp rationals */
void mpr_mul(mpr *rop, mpr *op1, mpr *op2);
/* multiply op1 by op2 */
void mpr_muli(mpr *rop, mpr *op1, long op2);
/* divide two mp rationals */
void mpr_div(mpr *rop, mpr *op1, mpr *op2);
/* divides op1 by op2 */
void mpr_divi(mpr *rop, mpr *op1, long op2);
/* sets rop to 1/op */
void mpr_inv(mpr *rop, mpr *op);
/* sets rop to -op */
void mpr_neg(mpr *rop, mpr *op);
/* sets rop to the absolute value of op */
void mpr_abs(mpr *rop, mpr *op);
/* compares op1 and op2
* returns >0 if op1 > op2, 0 if op1 = op2, and <0 if op1 < op2 */
int mpr_cmp(mpr *op1, mpr *op2);
/* mpr_cmp with a long integer operand */
int mpr_cmpi(mpr *op1, long op2);
/* compares absolute value of op1 and op2
* returns >0 if abs(op1) > abs(op2), 0 if abs(op1) = abs(op2),
* and <0 if abs(op1) < abs(op2) */
int mpr_cmpabs(mpr *op1, mpr *op2);
/* mpr_cmpabs with a long integer operand */
int mpr_cmpabsi(mpr *op1, long op2);
/* fastly swaps contents of op1 and op2 */
void mpr_swap(mpr *op1, mpr *op2);
/* returns 1 if op fits in a signed long int, 0 otherwise */
int mpr_fiti(mpr *op);
/* convert mp rational to double */
double mpr_getd(mpr *op);
/* returns pointer to string with representation of mp rational
* if str is not NULL, it must have enough space to store rational
* representation, if NULL a newly allocated string is returned.
* base must be in the range 2 to 36 */
char *mpr_getstr(char *str, mpr *op, int base);
#endif /* __mp_h_ */

1660
app/xedit/lisp/mp/mpi.c Normal file
View file

File diff suppressed because it is too large Load diff

436
app/xedit/lisp/mp/mpr.c Normal file
View file

@ -0,0 +1,436 @@
/*
* Copyright (c) 2002 by The XFree86 Project, Inc.
*
* Permission is hereby granted, free of charge, to any person obtaining a
* copy of this software and associated documentation files (the "Software"),
* to deal in the Software without restriction, including without limitation
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
* and/or sell copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included in
* all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
* THE XFREE86 PROJECT BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF
* OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
* SOFTWARE.
*
* Except as contained in this notice, the name of the XFree86 Project shall
* not be used in advertising or otherwise to promote the sale, use or other
* dealings in this Software without prior written authorization from the
* XFree86 Project.
*
* Author: Paulo César Pereira de Andrade
*/
/* $XFree86$ */
#include "mp.h"
/*
* TODO:
* Implement a fast gcd and divexact for integers, so that this code
* could be changed to do intermediary calculations faster using smaller
* numbers.
*/
/*
* Prototypes
*/
/* do the hard work of mpr_add and mpr_sub */
static void mpr_addsub(mpr *rop, mpr *op1, mpr *op2, int sub);
/* do the hard work of mpr_cmp and mpr_cmpabs */
static int mpr_docmp(mpr *op1, mpr *op2, int sign);
/*
* Implementation
*/
void
mpr_init(mpr *op)
{
op->num.digs = mp_malloc(sizeof(BNS));
op->num.sign = 0;
op->num.size = op->num.alloc = 1;
op->num.digs[0] = 0;
op->den.digs = mp_malloc(sizeof(BNS));
op->den.sign = 0;
op->den.size = op->den.alloc = 1;
op->den.digs[0] = 1;
}
void
mpr_clear(mpr *op)
{
op->num.sign = 0;
op->num.size = op->num.alloc = 0;
mp_free(op->num.digs);
op->den.sign = 0;
op->den.size = op->den.alloc = 0;
mp_free(op->den.digs);
}
void
mpr_set(mpr *rop, mpr *op)
{
if (rop != op) {
mpi_set(mpr_num(rop), mpr_num(op));
mpi_set(mpr_den(rop), mpr_den(op));
}
}
void
mpr_seti(mpr *rop, long num, long den)
{
mpi_seti(mpr_num(rop), num);
mpi_seti(mpr_den(rop), den);
}
void
mpr_setd(mpr *rop, double d)
{
double val, num;
int e, sign;
/* initialize */
if (d < 0) {
sign = 1;
val = -d;
}
else {
sign = 0;
val = d;
}
val = frexp(val, &e);
while (modf(val, &num) != 0.0 && val <= DBL_MAX / 2.0) {
--e;
val *= 2.0;
}
if (e >= 0) {
mpi_setd(mpr_num(rop), d);
mpi_seti(mpr_den(rop), 1);
}
else {
mpi_setd(mpr_num(rop), sign ? -num : num);
mpi_setd(mpr_den(rop), ldexp(1.0, -e));
}
}
void
mpr_setstr(mpr *rop, char *str, int base)
{
char *slash = strchr(str, '/');
mpi_setstr(mpr_num(rop), str, base);
if (slash != NULL)
mpi_setstr(mpr_den(rop), slash + 1, base);
else
mpi_seti(mpr_den(rop), 1);
}
void
mpr_canonicalize(mpr *op)
{
mpi gcd;
memset(&gcd, '\0', sizeof(mpi));
mpi_gcd(&gcd, mpr_num(op), mpr_den(op));
if (mpi_cmpabsi(&gcd, 1)) {
mpi_div(mpr_num(op), mpr_num(op), &gcd);
mpi_div(mpr_den(op), mpr_den(op), &gcd);
}
if (op->den.sign) {
op->num.sign = !op->num.sign;
op->den.sign = 0;
}
mpi_clear(&gcd);
}
void
mpr_add(mpr *rop, mpr *op1, mpr *op2)
{
mpr_addsub(rop, op1, op2, 0);
}
void
mpr_addi(mpr *rop, mpr *op1, long op2)
{
mpi prod;
memset(&prod, '\0', sizeof(mpi));
mpi_muli(&prod, mpr_den(op1), op2);
mpi_add(mpr_num(rop), mpr_num(op1), &prod);
mpi_clear(&prod);
}
void
mpr_sub(mpr *rop, mpr *op1, mpr *op2)
{
mpr_addsub(rop, op1, op2, 1);
}
void
mpr_subi(mpr *rop, mpr *op1, long op2)
{
mpi prod;
memset(&prod, '\0', sizeof(mpi));
mpi_muli(&prod, mpr_den(op1), op2);
mpi_sub(mpr_num(rop), mpr_num(op1), &prod);
mpi_clear(&prod);
}
static void
mpr_addsub(mpr *rop, mpr *op1, mpr *op2, int sub)
{
mpi prod1, prod2;
memset(&prod1, '\0', sizeof(mpi));
memset(&prod2, '\0', sizeof(mpi));
mpi_mul(&prod1, mpr_num(op1), mpr_den(op2));
mpi_mul(&prod2, mpr_num(op2), mpr_den(op1));
if (sub)
mpi_sub(mpr_num(rop), &prod1, &prod2);
else
mpi_add(mpr_num(rop), &prod1, &prod2);
mpi_clear(&prod1);
mpi_clear(&prod2);
mpi_mul(mpr_den(rop), mpr_den(op1), mpr_den(op2));
}
void
mpr_mul(mpr *rop, mpr *op1, mpr *op2)
{
/* check if temporary storage is required */
if (op1 == op2 && rop == op1) {
mpi prod;
memset(&prod, '\0', sizeof(mpi));
mpi_mul(&prod, mpr_num(op1), mpr_num(op2));
mpi_mul(mpr_den(rop), mpr_den(op1), mpr_den(op2));
mpi_set(mpr_num(rop), &prod);
mpi_clear(&prod);
}
else {
mpi_mul(mpr_num(rop), mpr_num(op1), mpr_num(op2));
mpi_mul(mpr_den(rop), mpr_den(op1), mpr_den(op2));
}
}
void
mpr_muli(mpr *rop, mpr *op1, long op2)
{
mpi_muli(mpr_num(rop), mpr_num(op1), op2);
}
void
mpr_div(mpr *rop, mpr *op1, mpr *op2)
{
/* check if temporary storage is required */
if (op1 == op2 && rop == op1) {
mpi prod;
memset(&prod, '\0', sizeof(mpi));
mpi_mul(&prod, mpr_num(op1), mpr_den(op2));
mpi_mul(mpr_den(rop), mpr_num(op2), mpr_den(op1));
mpi_set(mpr_num(rop), &prod);
mpi_clear(&prod);
}
else {
mpi_mul(mpr_num(rop), mpr_num(op1), mpr_den(op2));
mpi_mul(mpr_den(rop), mpr_num(op2), mpr_den(op1));
}
}
void
mpr_divi(mpr *rop, mpr *op1, long op2)
{
mpi_muli(mpr_den(rop), mpr_den(op1), op2);
}
void
mpr_inv(mpr *rop, mpr *op)
{
if (rop == op)
mpi_swap(mpr_num(op), mpr_den(op));
else {
mpi_set(mpr_num(rop), mpr_den(op));
mpi_set(mpr_den(rop), mpr_num(op));
}
}
void
mpr_neg(mpr *rop, mpr *op)
{
mpi_neg(mpr_num(rop), mpr_num(op));
mpi_set(mpr_den(rop), mpr_den(op));
}
void
mpr_abs(mpr *rop, mpr *op)
{
if (mpr_num(op)->sign)
mpi_neg(mpr_num(rop), mpr_num(op));
else
mpi_set(mpr_num(rop), mpr_num(op));
/* op may not be canonicalized */
if (mpr_den(op)->sign)
mpi_neg(mpr_den(rop), mpr_den(op));
else
mpi_set(mpr_den(rop), mpr_den(op));
}
int
mpr_cmp(mpr *op1, mpr *op2)
{
return (mpr_docmp(op1, op2, 1));
}
int
mpr_cmpi(mpr *op1, long op2)
{
int cmp;
mpr rat;
mpr_init(&rat);
mpi_seti(mpr_num(&rat), op2);
cmp = mpr_docmp(op1, &rat, 1);
mpr_clear(&rat);
return (cmp);
}
int
mpr_cmpabs(mpr *op1, mpr *op2)
{
return (mpr_docmp(op1, op2, 0));
}
int
mpr_cmpabsi(mpr *op1, long op2)
{
int cmp;
mpr rat;
mpr_init(&rat);
mpi_seti(mpr_num(&rat), op2);
cmp = mpr_docmp(op1, &rat, 0);
mpr_clear(&rat);
return (cmp);
}
static int
mpr_docmp(mpr *op1, mpr *op2, int sign)
{
int cmp, neg;
mpi prod1, prod2;
neg = 0;
if (sign) {
/* if op1 is negative */
if (mpr_num(op1)->sign ^ mpr_den(op1)->sign) {
/* if op2 is positive */
if (!(mpr_num(op2)->sign ^ mpr_den(op2)->sign))
return (-1);
else
neg = 1;
}
/* if op2 is negative */
else if (mpr_num(op2)->sign ^ mpr_den(op2)->sign)
return (1);
/* else same sign */
}
/* if denominators are equal, compare numerators */
if (mpi_cmpabs(mpr_den(op1), mpr_den(op2)) == 0) {
cmp = mpi_cmpabs(mpr_num(op1), mpr_num(op2));
if (cmp == 0)
return (0);
if (sign && neg)
return (cmp < 0 ? 1 : -1);
return (cmp);
}
memset(&prod1, '\0', sizeof(mpi));
memset(&prod2, '\0', sizeof(mpi));
/* "divide" op1 by op2
* if result is smaller than 1, op1 is smaller than op2 */
mpi_mul(&prod1, mpr_num(op1), mpr_den(op2));
mpi_mul(&prod2, mpr_num(op2), mpr_den(op1));
cmp = mpi_cmpabs(&prod1, &prod2);
mpi_clear(&prod1);
mpi_clear(&prod2);
if (sign && neg)
return (cmp < 0 ? 1 : -1);
return (cmp);
}
void
mpr_swap(mpr *op1, mpr *op2)
{
if (op1 != op2) {
mpr swap;
memcpy(&swap, op1, sizeof(mpr));
memcpy(op1, op2, sizeof(mpr));
memcpy(op2, &swap, sizeof(mpr));
}
}
int
mpr_fiti(mpr *op)
{
return (mpi_fiti(mpr_num(op)) && mpi_fiti(mpr_den(op)));
}
double
mpr_getd(mpr *op)
{
return (mpi_getd(mpr_num(op)) / mpi_getd(mpr_den(op)));
}
char *
mpr_getstr(char *str, mpr *op, int base)
{
int len;
if (str == NULL) {
len = mpi_getsize(mpr_num(op), base) + mpr_num(op)->sign + 1 +
mpi_getsize(mpr_den(op), base) + mpr_den(op)->sign + 1;
str = mp_malloc(len);
}
(void)mpi_getstr(str, mpr_num(op), base);
len = strlen(str);
str[len] = '/';
(void)mpi_getstr(str + len + 1, mpr_den(op), base);
return (str);
}