mp.h revision 5dfecf96
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: xc/programs/xedit/lisp/mp/mp.h,v 1.5tsi Exp $ */ 31 32#include <stdio.h> 33#include <math.h> 34#ifdef sun 35#include <ieeefp.h> 36#endif 37#include <float.h> 38#include <stdlib.h> 39#include <limits.h> 40#include <ctype.h> 41#include <string.h> 42 43#ifndef __mp_h_ 44#define __mp_h_ 45 46#ifdef __GNUC__ 47#define INLINE __inline__ 48#else 49#define INLINE /**/ 50#endif 51 52/* this normally is better for multiplication and also 53 * simplify addition loops putting the larger value first */ 54#define MP_SWAP(op1, op2, len1, len2) { \ 55 BNS *top = op1; \ 56 BNI tlen = len1; \ 57 \ 58 op1 = op2; \ 59 len1 = len2; \ 60 op2 = top; \ 61 len2 = tlen; \ 62} 63 64/* 65 * At least this length to use Karatsuba multiplication method 66 */ 67#define KARATSUBA 32 68 69/* 70 * At least this length to use Toom multiplication method 71 */ 72#define TOOM 128 73 74#if ULONG_MAX > 4294967295UL 75 /* sizeof(long) == 8 and sizeof(int) == 4 */ 76# define BNI unsigned long 77# define BNS unsigned int 78# define MINSLONG 0x8000000000000000UL 79# define CARRY 0x100000000 80# define LMASK 0xffffffff00000000UL 81# define SMASK 0x00000000ffffffffUL 82# define BNIBITS 64 83# define BNSBITS 32 84# ifndef LONG64 85# define LONG64 86# endif 87#else 88 /* sizeof(long) == 4 and sizeof(short) == 2 */ 89# define BNI unsigned long 90# define BNS unsigned short 91# define MINSLONG 0x80000000UL 92# define CARRY 0x10000 93# define LMASK 0xffff0000UL 94# define SMASK 0x0000ffffUL 95# define BNIBITS 32 96# define BNSBITS 16 97#endif 98 99#ifdef MAX 100#undef MAX 101#endif 102#define MAX(a, b) ((a) > (b) ? (a) : (b)) 103 104#ifdef MIN 105#undef MIN 106#endif 107#define MIN(a, b) ((a) < (b) ? (a) : (b)) 108 109/* 110 * Types 111 */ 112typedef struct _mpi { 113 unsigned int size : 31; 114 unsigned int sign : 1; 115 BNI alloc; 116 BNS *digs; /* LSF format */ 117} mpi; 118 119typedef struct _mpr { 120 mpi num; 121 mpi den; 122} mpr; 123 124typedef void *(*mp_malloc_fun)(size_t); 125typedef void *(*mp_calloc_fun)(size_t, size_t); 126typedef void *(*mp_realloc_fun)(void*, size_t); 127typedef void (*mp_free_fun)(void*); 128 129/* 130 * Prototypes 131 */ 132/* GENERIC FUNCTIONS */ 133 /* memory allocation wrappers */ 134void *mp_malloc(size_t size); 135void *mp_calloc(size_t nmemb, size_t size); 136void *mp_realloc(void *pointer, size_t size); 137void mp_free(void *pointer); 138mp_malloc_fun mp_set_malloc(mp_malloc_fun); 139mp_calloc_fun mp_set_calloc(mp_calloc_fun); 140mp_realloc_fun mp_set_realloc(mp_realloc_fun); 141mp_free_fun mp_set_free(mp_free_fun); 142 143 /* adds op1 and op2, stores result in rop 144 * rop must pointer to at least len1 + len2 + 1 elements 145 * rop can be either op1 or op2 */ 146long mp_add(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2); 147 148 /* subtracts op2 from op1, stores result in rop 149 * rop must pointer to at least len1 + len2 elements 150 * op1 must be >= op2 151 * rop can be either op1 or op2 */ 152long mp_sub(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2); 153 154 /* shift op to the left shift bits 155 * rop must have enough storage for result 156 * rop can be op */ 157long mp_lshift(BNS *rop, BNS *op, BNI len, long shift); 158 159 /* shift op to the right shift bits 160 * shift must be positive 161 * rop can be op */ 162long mp_rshift(BNS *rop, BNS *op, BNI len, long shift); 163 164 /* use simple generic multiplication method 165 * rop cannot be the same as op1 or op2 166 * rop must be zeroed 167 * op1 can be op2 */ 168long mp_base_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2); 169 170 /* use Karatsuba method 171 * MIN(len1, len2) must be larger than (MAX(len1, len2) + 1) >> 1 172 * MAX(len1, len2) should be at least 2 173 * rop cannot be the same as op1 or op2 174 * rop must be zeroed 175 * op1 can be op2 */ 176long mp_karatsuba_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2); 177 178 /* use Toom method 179 * len1 / 3 should be equal to len2 / 3 180 * len1 / 3 should be at least 1 181 * rop cannot be the same as op1 or op2 182 * rop must be zeroed 183 * op1 can be op2 */ 184long mp_toom_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2); 185 186 /* chooses the available multiplication methods based on it's input 187 * rop must be a pointer to len1 + len2 elements 188 * rop cannot be the same as op1 or op2 189 * rop must be zeroed 190 * op1 can be op2 */ 191long mp_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2); 192 193/* INTEGER FUNCTIONS */ 194 /* initialize op and set it to 0 */ 195void mpi_init(mpi *op); 196 197 /* clear memory associated to op */ 198void mpi_clear(mpi *op); 199 200 /* set rop to the value of op */ 201void mpi_set(mpi *rop, mpi *op); 202 203 /* set rop to the value of si */ 204void mpi_seti(mpi *rop, long si); 205 206 /* set rop to the floor(fabs(d)) */ 207void mpi_setd(mpi *rop, double d); 208 209 /* initialize rop to number representation in str in the given base. 210 * leading zeros are skipped. 211 * if sign present, it is processed. 212 * base must be in the range 2 to 36. */ 213void mpi_setstr(mpi *rop, char *str, int base); 214 215 /* adds two mp integers */ 216void mpi_add(mpi *rop, mpi *op1, mpi *op2); 217 218 /* adds op1 and op2 */ 219void mpi_addi(mpi *rop, mpi *op1, long op2); 220 221 /* subtracts two mp integers */ 222void mpi_sub(mpi *rop, mpi *op1, mpi *op2); 223 224 /* subtracts op2 from op1 */ 225void mpi_subi(mpi *rop, mpi *op1, long op2); 226 227 /* multiply two mp integers */ 228void mpi_mul(mpi *rop, mpi *op1, mpi *op2); 229 230 /* multiply op1 by op2 */ 231void mpi_muli(mpi *rop, mpi *op1, long op2); 232 233 /* divides num by den and sets rop to result */ 234void mpi_div(mpi *rop, mpi *num, mpi *den); 235 236 /* divides num by den and sets rop to the remainder */ 237void mpi_rem(mpi *rop, mpi *num, mpi *den); 238 239 /* divides num by den, sets quotient to qrop and remainder to rrop 240 * qrop is truncated towards zero. 241 * qrop and rrop are optional 242 * qrop and rrop cannot be the same variable */ 243void mpi_divqr(mpi *qrop, mpi *rrop, mpi *num, mpi *den); 244 245 /* divides num by then and stores result in rop */ 246void mpi_divi(mpi *rop, mpi *num, long den); 247 248 /* divides num by den and returns remainder */ 249long mpi_remi(mpi *num, long den); 250 251 /* divides num by den 252 * stores quotient in qrop and returns remainder */ 253long mpi_divqri(mpi *qrop, mpi *num, long den); 254 255 /* sets rop to num modulo den */ 256void mpi_mod(mpi *rop, mpi *num, mpi *den); 257 258 /* returns num modulo den */ 259long mpi_modi(mpi *num, long den); 260 261 /* sets rop to the greatest common divisor of num and den 262 * result is always positive */ 263void mpi_gcd(mpi *rop, mpi *num, mpi *den); 264 265 /* sets rop to the least common multiple of num and den 266 * result is always positive */ 267void mpi_lcm(mpi *rop, mpi *num, mpi *den); 268 269 /* sets rop to op raised to exp */ 270void mpi_pow(mpi *rop, mpi *op, unsigned long exp); 271 272 /* sets rop to the integer part of the nth root of op. 273 * returns 1 if result is exact, 0 otherwise */ 274int mpi_root(mpi *rop, mpi *op, unsigned long nth); 275 276 /* sets rop to the integer part of the square root of op. 277 * returns 1 if result is exact, 0 otherwise */ 278int mpi_sqrt(mpi *rop, mpi *op); 279 280 /* bit shift, left if shift positive, right if negative 281 * a fast way to multiply and divide by powers of two */ 282void mpi_ash(mpi *rop, mpi *op, long shift); 283 284 /* sets rop to op1 logand op2 */ 285void mpi_and(mpi *rop, mpi *op1, mpi *op2); 286 287 /* sets rop to op1 logior op2 */ 288void mpi_ior(mpi *rop, mpi *op1, mpi *op2); 289 290 /* sets rop to op1 logxor op2 */ 291void mpi_xor(mpi *rop, mpi *op1, mpi *op2); 292 293 /* sets rop to one's complement of op */ 294void mpi_com(mpi *rop, mpi *op); 295 296 /* sets rop to -op */ 297void mpi_neg(mpi *rop, mpi *op); 298 299 /* sets rop to the absolute value of op */ 300void mpi_abs(mpi *rop, mpi *op); 301 302 /* compares op1 and op2 303 * returns >0 if op1 > op2, 0 if op1 = op2, and <0 if op1 < op2 */ 304int mpi_cmp(mpi *op1, mpi *op2); 305 306 /* mpi_cmp with a long integer operand */ 307int mpi_cmpi(mpi *op1, long op2); 308 309 /* compares absolute value of op1 and op2 310 * returns >0 if abs(op1) > abs(op2), 0 if abs(op1) = abs(op2), 311 * and <0 if abs(op1) < abs(op2) */ 312int mpi_cmpabs(mpi *op1, mpi *op2); 313 314 /* mpi_cmpabs with a long integer operand */ 315int mpi_cmpabsi(mpi *op1, long op2); 316 317 /* returns 1 if op1 > 0, 0 if op1 = 0, and -1 if op1 < 0 */ 318int mpi_sgn(mpi *op); 319 320 /* fastly swaps contents of op1 and op2 */ 321void mpi_swap(mpi *op1, mpi *op2); 322 323 /* returns 1 if op fits in a signed long int, 0 otherwise */ 324int mpi_fiti(mpi *op); 325 326 /* converts mp integer to long int 327 * to know if the value will fit, call mpi_fiti */ 328long mpi_geti(mpi *op); 329 330 /* convert mp integer to double */ 331double mpi_getd(mpi *op); 332 333 /* returns exact number of characters to represent mp integer 334 * in given base, excluding sign and ending null character. 335 * base must be in the range 2 to 36 */ 336unsigned long mpi_getsize(mpi *op, int base); 337 338 /* returns pointer to string with representation of mp integer 339 * if str is not NULL, it must have enough space to store integer 340 * representation, if NULL a newly allocated string is returned. 341 * base must be in the range 2 to 36 */ 342char *mpi_getstr(char *str, mpi *op, int base); 343 344/* RATIO FUNCTIONS */ 345#define mpr_num(op) (&((op)->num)) 346#define mpr_den(op) (&((op)->den)) 347 348 /* initialize op and set it to 0/1 */ 349void mpr_init(mpr *op); 350 351 /* clear memory associated to op */ 352void mpr_clear(mpr *op); 353 354 /* set rop to the value of op */ 355void mpr_set(mpr *rop, mpr *op); 356 357 /* set rop to num/den */ 358void mpr_seti(mpr *rop, long num, long den); 359 360 /* set rop to the value of d */ 361void mpr_setd(mpr *rop, double d); 362 363 /* initialize rop to number representation in str in the given base. 364 * leading zeros are skipped. 365 * if sign present, it is processed. 366 * base must be in the range 2 to 36. */ 367void mpr_setstr(mpr *rop, char *str, int base); 368 369 /* remove common factors of op */ 370void mpr_canonicalize(mpr *op); 371 372 /* adds two mp rationals */ 373void mpr_add(mpr *rop, mpr *op1, mpr *op2); 374 375 /* adds op1 and op2 */ 376void mpr_addi(mpr *rop, mpr *op1, long op2); 377 378 /* subtracts two mp rationals */ 379void mpr_sub(mpr *rop, mpr *op1, mpr *op2); 380 381 /* subtracts op2 from op1 */ 382void mpr_subi(mpr *rop, mpr *op1, long op2); 383 384 /* multiply two mp rationals */ 385void mpr_mul(mpr *rop, mpr *op1, mpr *op2); 386 387 /* multiply op1 by op2 */ 388void mpr_muli(mpr *rop, mpr *op1, long op2); 389 390 /* divide two mp rationals */ 391void mpr_div(mpr *rop, mpr *op1, mpr *op2); 392 393 /* divides op1 by op2 */ 394void mpr_divi(mpr *rop, mpr *op1, long op2); 395 396 /* sets rop to 1/op */ 397void mpr_inv(mpr *rop, mpr *op); 398 399 /* sets rop to -op */ 400void mpr_neg(mpr *rop, mpr *op); 401 402 /* sets rop to the absolute value of op */ 403void mpr_abs(mpr *rop, mpr *op); 404 405 /* compares op1 and op2 406 * returns >0 if op1 > op2, 0 if op1 = op2, and <0 if op1 < op2 */ 407int mpr_cmp(mpr *op1, mpr *op2); 408 409 /* mpr_cmp with a long integer operand */ 410int mpr_cmpi(mpr *op1, long op2); 411 412 /* compares absolute value of op1 and op2 413 * returns >0 if abs(op1) > abs(op2), 0 if abs(op1) = abs(op2), 414 * and <0 if abs(op1) < abs(op2) */ 415int mpr_cmpabs(mpr *op1, mpr *op2); 416 417 /* mpr_cmpabs with a long integer operand */ 418int mpr_cmpabsi(mpr *op1, long op2); 419 420 /* fastly swaps contents of op1 and op2 */ 421void mpr_swap(mpr *op1, mpr *op2); 422 423 /* returns 1 if op fits in a signed long int, 0 otherwise */ 424int mpr_fiti(mpr *op); 425 426 /* convert mp rational to double */ 427double mpr_getd(mpr *op); 428 429 /* returns pointer to string with representation of mp rational 430 * if str is not NULL, it must have enough space to store rational 431 * representation, if NULL a newly allocated string is returned. 432 * base must be in the range 2 to 36 */ 433char *mpr_getstr(char *str, mpr *op, int base); 434 435#endif /* __mp_h_ */ 436