src/mp.h

   1
   2 /*  Copyright (c) 1987-2008 Sun Microsystems, Inc. All Rights Reserved.
   3  *  Copyright (c) 2008-2009 Robert Ancell
   4  *
   5  *  This program is free software; you can redistribute it and/or modify
   6  *  it under the terms of the GNU General Public License as published by
   7  *  the Free Software Foundation; either version 2, or (at your option)
   8  *  any later version.
   9  *
  10  *  This program is distributed in the hope that it will be useful, but
  11  *  WITHOUT ANY WARRANTY; without even the implied warranty of
  12  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  13  *  General Public License for more details.
  14  *
  15  *  You should have received a copy of the GNU General Public License
  16  *  along with this program; if not, write to the Free Software
  17  *  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  18  *  02111-1307, USA.
  19  */
  20
  21 /*  This maths library is based on the MP multi-precision floating-point
  22  *  arithmetic package originally written in FORTRAN by Richard Brent,
  23  *  Computer Centre, Australian National University in the 1970's.
  24  *
  25  *  It has been converted from FORTRAN into C using the freely available
  26  *  f2c translator, available via netlib on research.att.com.
  27  *
  28  *  The subsequently converted C code has then been tidied up, mainly to
  29  *  remove any dependencies on the libI77 and libF77 support libraries.
  30  *
  31  *  FOR A GENERAL DESCRIPTION OF THE PHILOSOPHY AND DESIGN OF MP,
  32  *  SEE - R. P. BRENT, A FORTRAN MULTIPLE-PRECISION ARITHMETIC
  33  *  PACKAGE, ACM TRANS. MATH. SOFTWARE 4 (MARCH 1978), 57-70.
  34  *  SOME ADDITIONAL DETAILS ARE GIVEN IN THE SAME ISSUE, 71-81.
  35  *  FOR DETAILS OF THE IMPLEMENTATION, CALLING SEQUENCES ETC. SEE
  36  *  THE MP USERS GUIDE.
  37  */
  38
  39 #ifndef MP_H
  40 #define MP_H
  41
  42 #include <stdbool.h>
  43 #include <stdint.h>
  44 #include <glib.h>
  45
  46 /* Size of the multiple precision values */
  47 #define MP_SIZE 1000
  48
  49 /* Base for numbers */
  50 #define MP_BASE 10000
  51
  52 /* Object for a high precision floating point number representation
  53  *
  54  * x = sign * (MP_BASE^(exponent-1) + MP_BASE^(exponent-2) + ...)
  55  */
  56 typedef struct
  57 {
  58    /* Sign (+1, -1) or 0 for the value zero */
  59    int sign, im_sign;
  60
  61    /* Exponent (to base MP_BASE) */
  62    int exponent, im_exponent;
  63
  64    /* Normalized fraction */
  65    int fraction[MP_SIZE], im_fraction[MP_SIZE];
  66 } MPNumber;
  67
  68 typedef enum
  69 {
  70     MP_RADIANS,
  71     MP_DEGREES,
  72     MP_GRADIANS
  73 } MPAngleUnit;
  74
  75 /* Returns error string or NULL if no error */
  76 // FIXME: Global variable
  77 const char *mp_get_error(void);
  78
  79 /* Clear any current error */
  80 void mp_clear_error(void);
  81
  82 /* Returns:
  83  *  0 if x == y
  84  * <0 if x < y
  85  * >0 if x > y
  86  */
  87 int    mp_compare_mp_to_mp(const MPNumber *x, const MPNumber *y);
  88
  89 /* Return true if the value is x == 0 */
  90 bool   mp_is_zero(const MPNumber *x);
  91
  92 /* Return true if x < 0 */
  93 bool   mp_is_negative(const MPNumber *x);
  94
  95 /* Return true if x is integer */
  96 bool   mp_is_integer(const MPNumber *x);
  97
  98 /* Return true if x is a positive integer */
  99 bool   mp_is_positive_integer(const MPNumber *x);
 100
 101 /* Return true if x is a natural number (an integer ≥ 0) */
 102 bool   mp_is_natural(const MPNumber *x);
 103
 104 /* Return true if x has an imaginary component */
 105 bool   mp_is_complex(const MPNumber *x);
 106
 107 /* Return true if x == y */
 108 bool   mp_is_equal(const MPNumber *x, const MPNumber *y);
 109
 110 /* Return true if x ≥ y */
 111 bool   mp_is_greater_equal(const MPNumber *x, const MPNumber *y);
 112
 113 /* Return true if x > y */
 114 bool   mp_is_greater_than(const MPNumber *x, const MPNumber *y);
 115
 116 /* Return true if x ≤ y */
 117 bool   mp_is_less_equal(const MPNumber *x, const MPNumber *y);
 118
 119 /* Return true if x < y */
 120 bool   mp_is_less_than(const MPNumber *x, const MPNumber *y);
 121
 122 /* Sets z = |x| */
 123 void   mp_abs(const MPNumber *x, MPNumber *z);
 124
 125 /* Sets z = Arg(x) */
 126 void   mp_arg(const MPNumber *x, MPAngleUnit unit, MPNumber *z);
 127
 128 /* Sets z = ‾̅x */
 129 void   mp_conjugate(const MPNumber *x, MPNumber *z);
 130
 131 /* Sets z = Re(x) */
 132 void   mp_real_component(const MPNumber *x, MPNumber *z);
 133
 134 /* Sets z = Im(x) */
 135 void   mp_imaginary_component(const MPNumber *x, MPNumber *z);
 136
 137 /* Sets z = −x */
 138 void   mp_invert_sign(const MPNumber *x, MPNumber *z);
 139
 140 /* Sets z = x + y */
 141 void   mp_add(const MPNumber *x, const MPNumber *y, MPNumber *z);
 142
 143 /* Sets z = x + y */
 144 void   mp_add_integer(const MPNumber *x, int64_t y, MPNumber *z);
 145
 146 /* Sets z = x + numerator ÷ denominator */
 147 void   mp_add_fraction(const MPNumber *x, int64_t numerator, int64_t denominator, MPNumber *z);
 148
 149 /* Sets z = x − y */
 150 void   mp_subtract(const MPNumber *x, const MPNumber *y, MPNumber *z);
 151
 152 /* Sets z = x × y */
 153 void   mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z);
 154
 155 /* Sets z = x × y */
 156 void   mp_multiply_integer(const MPNumber *x, int64_t y, MPNumber *z);
 157
 158 /* Sets z = x × numerator ÷ denominator */
 159 void   mp_multiply_fraction(const MPNumber *x, int64_t numerator, int64_t denominator, MPNumber *z);
 160
 161 /* Sets z = x ÷ y */
 162 void   mp_divide(const MPNumber *x, const MPNumber *y, MPNumber *z);
 163
 164 /* Sets z = x ÷ y */
 165 void   mp_divide_integer(const MPNumber *x, int64_t y, MPNumber *z);
 166
 167 /* Sets z = 1 ÷ x */
 168 void   mp_reciprocal(const MPNumber *, MPNumber *);
 169
 170 /* Sets z = sgn(x) */
 171 void   mp_sgn(const MPNumber *x, MPNumber *z);
 172
 173 void   mp_integer_component(const MPNumber *x, MPNumber *z);
 174
 175 /* Sets z = x mod 1 */
 176 void   mp_fractional_component(const MPNumber *x, MPNumber *z);
 177
 178 /* Sets z = {x} */
 179 void   mp_fractional_part(const MPNumber *x, MPNumber *z);
 180
 181 /* Sets z = ⌊x⌋ */
 182 void   mp_floor(const MPNumber *x, MPNumber *z);
 183
 184 /* Sets z = ⌈x⌉ */
 185 void   mp_ceiling(const MPNumber *x, MPNumber *z);
 186
 187 /* Sets z = [x] */
 188 void   mp_round(const MPNumber *x, MPNumber *z);
 189
 190 /* Sets z = ln x */
 191 void   mp_ln(const MPNumber *x, MPNumber *z);
 192
 193 /* Sets z = log_n x */
 194 void   mp_logarithm(int64_t n, const MPNumber *x, MPNumber *z);
 195
 196 /* Sets z = π */
 197 void   mp_get_pi(MPNumber *z);
 198
 199 /* Sets z = e */
 200 void   mp_get_eulers(MPNumber *z);
 201
 202 /* Sets z = i (√−1) */
 203 void   mp_get_i(MPNumber *z);
 204
 205 /* Sets z = n√x */
 206 void   mp_root(const MPNumber *x, int64_t n, MPNumber *z);
 207
 208 /* Sets z = √x */
 209 void   mp_sqrt(const MPNumber *x, MPNumber *z);
 210
 211 /* Sets z = x! */
 212 void   mp_factorial(const MPNumber *x, MPNumber *z);
 213
 214 /* Sets z = x mod y */
 215 void   mp_modulus_divide(const MPNumber *x, const MPNumber *y, MPNumber *z);
 216
 217 /* Sets z = x^y */
 218 void   mp_xpowy(const MPNumber *x, const MPNumber *y, MPNumber *z);
 219
 220 /* Sets z = x^y */
 221 void   mp_xpowy_integer(const MPNumber *x, int64_t y, MPNumber *z);
 222
 223 /* Sets z = e^x */
 224 void   mp_epowy(const MPNumber *x, MPNumber *z);
 225
 226 /* Returns a list of all prime factors in x as MPNumbers */
 227 GList* mp_factorize(const MPNumber *x);
 228
 229 /* Sets z = x */
 230 void   mp_set_from_mp(const MPNumber *x, MPNumber *z);
 231
 232 /* Sets z = x */
 233 void   mp_set_from_float(float x, MPNumber *z);
 234
 235 /* Sets z = x */
 236 void   mp_set_from_double(double x, MPNumber *z);
 237
 238 /* Sets z = x */
 239 void   mp_set_from_integer(int64_t x, MPNumber *z);
 240
 241 /* Sets z = x */
 242 void   mp_set_from_unsigned_integer(uint64_t x, MPNumber *z);
 243
 244 /* Sets z = numerator ÷ denominator */
 245 void   mp_set_from_fraction(int64_t numerator, int64_t denominator, MPNumber *z);
 246
 247 /* Sets z = r(cos theta + i sin theta) */
 248 void   mp_set_from_polar(const MPNumber *r, MPAngleUnit unit, const MPNumber *theta, MPNumber *z);
 249
 250 /* Sets z = x + iy */
 251 void   mp_set_from_complex(const MPNumber *x, const MPNumber *y, MPNumber *z);
 252
 253 /* Sets z to be a uniform random number in the range [0, 1] */
 254 void   mp_set_from_random(MPNumber *z);
 255
 256 /* Sets z from a string representation in 'text'.
 257  * Returns true on success.
 258  */
 259 bool   mp_set_from_string(const char *text, int default_base, MPNumber *z);
 260
 261 /* Returns x as a native single-precision floating point number */
 262 float  mp_cast_to_float(const MPNumber *x);
 263
 264 /* Returns x as a native double-precision floating point number */
 265 double mp_cast_to_double(const MPNumber *x);
 266
 267 /* Returns x as a native integer */
 268 int64_t mp_cast_to_int(const MPNumber *x);
 269
 270 /* Returns x as a native unsigned integer */
 271 uint64_t mp_cast_to_unsigned_int(const MPNumber *x);
 272
 273 /* Sets z = sin x */
 274 void   mp_sin(const MPNumber *x, MPAngleUnit unit, MPNumber *z);
 275
 276 /* Sets z = cos x */
 277 void   mp_cos(const MPNumber *x, MPAngleUnit unit, MPNumber *z);
 278
 279 /* Sets z = tan x */
 280 void   mp_tan(const MPNumber *x, MPAngleUnit unit, MPNumber *z);
 281
 282 /* Sets z = sin⁻¹ x */
 283 void   mp_asin(const MPNumber *x, MPAngleUnit unit, MPNumber *z);
 284
 285 /* Sets z = cos⁻¹ x */
 286 void   mp_acos(const MPNumber *x, MPAngleUnit unit, MPNumber *z);
 287
 288 /* Sets z = tan⁻¹ x */
 289 void   mp_atan(const MPNumber *x, MPAngleUnit unit, MPNumber *z);
 290
 291 /* Sets z = sinh x */
 292 void   mp_sinh(const MPNumber *x, MPNumber *z);
 293
 294 /* Sets z = cosh x */
 295 void   mp_cosh(const MPNumber *x, MPNumber *z);
 296
 297 /* Sets z = tanh x */
 298 void   mp_tanh(const MPNumber *x, MPNumber *z);
 299
 300 /* Sets z = sinh⁻¹ x */
 301 void   mp_asinh(const MPNumber *x, MPNumber *z);
 302
 303 /* Sets z = cosh⁻¹ x */
 304 void   mp_acosh(const MPNumber *x, MPNumber *z);
 305
 306 /* Sets z = tanh⁻¹ x */
 307 void   mp_atanh(const MPNumber *x, MPNumber *z);
 308
 309 /* Returns true if x is cannot be represented in a binary word of length 'wordlen' */
 310 bool   mp_is_overflow(const MPNumber *x, int wordlen);
 311
 312 /* Sets z = boolean AND for each bit in x and z */
 313 void   mp_and(const MPNumber *x, const MPNumber *y, MPNumber *z);
 314
 315 /* Sets z = boolean OR for each bit in x and z */
 316 void   mp_or(const MPNumber *x, const MPNumber *y, MPNumber *z);
 317
 318 /* Sets z = boolean XOR for each bit in x and z */
 319 void   mp_xor(const MPNumber *x, const MPNumber *y, MPNumber *z);
 320
 321 /* Sets z = boolean XNOR for each bit in x and z for word of length 'wordlen' */
 322 void   mp_xnor(const MPNumber *x, const MPNumber *y, int wordlen, MPNumber *z);
 323
 324 /* Sets z = boolean NOT for each bit in x and z for word of length 'wordlen' */
 325 void   mp_not(const MPNumber *x, int wordlen, MPNumber *z);
 326
 327 /* Sets z = x masked to 'wordlen' bits */
 328 void   mp_mask(const MPNumber *x, int wordlen, MPNumber *z);
 329
 330 /* Sets z = x shifted by 'count' bits.  Positive shift increases the value, negative decreases */
 331 void   mp_shift(const MPNumber *x, int count, MPNumber *z);
 332
 333 /* Sets z to be the ones complement of x for word of length 'wordlen' */
 334 void   mp_ones_complement(const MPNumber *x, int wordlen, MPNumber *z);
 335
 336 /* Sets z to be the twos complement of x for word of length 'wordlen' */
 337 void   mp_twos_complement(const MPNumber *x, int wordlen, MPNumber *z);
 338
 339 #endif /* MP_H */