diff options
Diffstat (limited to 'lisp/mp/mp.h')
| -rw-r--r-- | lisp/mp/mp.h | 435 |
1 files changed, 435 insertions, 0 deletions
diff --git a/lisp/mp/mp.h b/lisp/mp/mp.h new file mode 100644 index 0000000..88f1b24 --- /dev/null +++ b/lisp/mp/mp.h @@ -0,0 +1,435 @@ +/* + * 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.6 2003/01/12 03:55:51 tsi 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 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 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_ */ |