summaryrefslogtreecommitdiff
path: root/src/mp.h
blob: 634bb0147757d3f48208b3731edc2b3d1170d000 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
/*
 * Copyright (C) 1987-2008 Sun Microsystems, Inc. All Rights Reserved.
 * Copyright (C) 2008-2011 Robert Ancell.
 * 
 * This program is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License as published by the Free Software
 * Foundation, either version 2 of the License, or (at your option) any later
 * version. See http://www.gnu.org/copyleft/gpl.html the full text of the
 * license.
 */

/*  This maths library is based on the MP multi-precision floating-point
 *  arithmetic package originally written in FORTRAN by Richard Brent,
 *  Computer Centre, Australian National University in the 1970's.
 *
 *  It has been converted from FORTRAN into C using the freely available
 *  f2c translator, available via netlib on research.att.com.
 *
 *  The subsequently converted C code has then been tidied up, mainly to
 *  remove any dependencies on the libI77 and libF77 support libraries.
 *
 *  FOR A GENERAL DESCRIPTION OF THE PHILOSOPHY AND DESIGN OF MP,
 *  SEE - R. P. BRENT, A FORTRAN MULTIPLE-PRECISION ARITHMETIC
 *  PACKAGE, ACM TRANS. MATH. SOFTWARE 4 (MARCH 1978), 57-70.
 *  SOME ADDITIONAL DETAILS ARE GIVEN IN THE SAME ISSUE, 71-81.
 *  FOR DETAILS OF THE IMPLEMENTATION, CALLING SEQUENCES ETC. SEE
 *  THE MP USERS GUIDE.
 */

#ifndef MP_H
#define MP_H

#include <stdbool.h>
#include <stdint.h>
#include <glib.h>

/* Size of the multiple precision values */
#define MP_SIZE 1000

/* Base for numbers */
#define MP_BASE 10000

/* Object for a high precision floating point number representation
 *
 * x = sign * (MP_BASE^(exponent-1) + MP_BASE^(exponent-2) + ...)
 */
typedef struct
{
   /* Sign (+1, -1) or 0 for the value zero */
   int sign, im_sign;

   /* Exponent (to base MP_BASE) */
   int exponent, im_exponent;

   /* Normalized fraction */
   int fraction[MP_SIZE], im_fraction[MP_SIZE];
} MPNumber;

typedef enum
{
    MP_RADIANS,
    MP_DEGREES,
    MP_GRADIANS
} MPAngleUnit;

/* Returns error string or NULL if no error */
// FIXME: Global variable
const char *mp_get_error(void);

/* Clear any current error */
void mp_clear_error(void);

/* Returns:
 *  0 if x == y
 * <0 if x < y
 * >0 if x > y
 */
int    mp_compare_mp_to_mp(const MPNumber *x, const MPNumber *y);

/* Return true if the value is x == 0 */
bool   mp_is_zero(const MPNumber *x);

/* Return true if x < 0 */
bool   mp_is_negative(const MPNumber *x);

/* Return true if x is integer */
bool   mp_is_integer(const MPNumber *x);

/* Return true if x is a positive integer */
bool   mp_is_positive_integer(const MPNumber *x);

/* Return true if x is a natural number (an integer ≥ 0) */
bool   mp_is_natural(const MPNumber *x);

/* Return true if x has an imaginary component */
bool   mp_is_complex(const MPNumber *x);

/* Return true if x == y */
bool   mp_is_equal(const MPNumber *x, const MPNumber *y);

/* Return true if x ≥ y */
bool   mp_is_greater_equal(const MPNumber *x, const MPNumber *y);

/* Return true if x > y */
bool   mp_is_greater_than(const MPNumber *x, const MPNumber *y);

/* Return true if x ≤ y */
bool   mp_is_less_equal(const MPNumber *x, const MPNumber *y);

/* Return true if x < y */
bool   mp_is_less_than(const MPNumber *x, const MPNumber *y);

/* Sets z = |x| */
void   mp_abs(const MPNumber *x, MPNumber *z);

/* Sets z = Arg(x) */
void   mp_arg(const MPNumber *x, MPAngleUnit unit, MPNumber *z);

/* Sets z = ‾̅x */
void   mp_conjugate(const MPNumber *x, MPNumber *z);

/* Sets z = Re(x) */
void   mp_real_component(const MPNumber *x, MPNumber *z);

/* Sets z = Im(x) */
void   mp_imaginary_component(const MPNumber *x, MPNumber *z);

/* Sets z = −x */
void   mp_invert_sign(const MPNumber *x, MPNumber *z);

/* Sets z = x + y */
void   mp_add(const MPNumber *x, const MPNumber *y, MPNumber *z);

/* Sets z = x + y */
void   mp_add_integer(const MPNumber *x, int64_t y, MPNumber *z);

/* Sets z = x + numerator ÷ denominator */
void   mp_add_fraction(const MPNumber *x, int64_t numerator, int64_t denominator, MPNumber *z);

/* Sets z = x − y */
void   mp_subtract(const MPNumber *x, const MPNumber *y, MPNumber *z);

/* Sets z = x × y */
void   mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z);

/* Sets z = x × y */
void   mp_multiply_integer(const MPNumber *x, int64_t y, MPNumber *z);

/* Sets z = x × numerator ÷ denominator */
void   mp_multiply_fraction(const MPNumber *x, int64_t numerator, int64_t denominator, MPNumber *z);

/* Sets z = x ÷ y */
void   mp_divide(const MPNumber *x, const MPNumber *y, MPNumber *z);

/* Sets z = x ÷ y */
void   mp_divide_integer(const MPNumber *x, int64_t y, MPNumber *z);

/* Sets z = 1 ÷ x */
void   mp_reciprocal(const MPNumber *, MPNumber *);

/* Sets z = sgn(x) */
void   mp_sgn(const MPNumber *x, MPNumber *z);

void   mp_integer_component(const MPNumber *x, MPNumber *z);

/* Sets z = x mod 1 */
void   mp_fractional_component(const MPNumber *x, MPNumber *z);

/* Sets z = {x} */
void   mp_fractional_part(const MPNumber *x, MPNumber *z);

/* Sets z = ⌊x⌋ */
void   mp_floor(const MPNumber *x, MPNumber *z);

/* Sets z = ⌈x⌉ */
void   mp_ceiling(const MPNumber *x, MPNumber *z);

/* Sets z = [x] */
void   mp_round(const MPNumber *x, MPNumber *z);

/* Sets z = ln x */
void   mp_ln(const MPNumber *x, MPNumber *z);

/* Sets z = log_n x */
void   mp_logarithm(int64_t n, const MPNumber *x, MPNumber *z);

/* Sets z = π */
void   mp_get_pi(MPNumber *z);

/* Sets z = e */
void   mp_get_eulers(MPNumber *z);

/* Sets z = i (√−1) */
void   mp_get_i(MPNumber *z);

/* Sets z = n√x */
void   mp_root(const MPNumber *x, int64_t n, MPNumber *z);

/* Sets z = √x */
void   mp_sqrt(const MPNumber *x, MPNumber *z);

/* Sets z = x! */
void   mp_factorial(const MPNumber *x, MPNumber *z);

/* Sets z = x mod y */
void   mp_modulus_divide(const MPNumber *x, const MPNumber *y, MPNumber *z);

/* Sets z = x^y */
void   mp_xpowy(const MPNumber *x, const MPNumber *y, MPNumber *z);

/* Sets z = x^y */
void   mp_xpowy_integer(const MPNumber *x, int64_t y, MPNumber *z);

/* Sets z = e^x */
void   mp_epowy(const MPNumber *x, MPNumber *z);

/* Returns a list of all prime factors in x as MPNumbers */
GList* mp_factorize(const MPNumber *x);

/* Sets z = x */
void   mp_set_from_mp(const MPNumber *x, MPNumber *z);

/* Sets z = x */
void   mp_set_from_float(float x, MPNumber *z);

/* Sets z = x */
void   mp_set_from_double(double x, MPNumber *z);

/* Sets z = x */
void   mp_set_from_integer(int64_t x, MPNumber *z);

/* Sets z = x */
void   mp_set_from_unsigned_integer(uint64_t x, MPNumber *z);

/* Sets z = numerator ÷ denominator */
void   mp_set_from_fraction(int64_t numerator, int64_t denominator, MPNumber *z);

/* Sets z = r(cos theta + i sin theta) */
void   mp_set_from_polar(const MPNumber *r, MPAngleUnit unit, const MPNumber *theta, MPNumber *z);

/* Sets z = x + iy */
void   mp_set_from_complex(const MPNumber *x, const MPNumber *y, MPNumber *z);

/* Sets z to be a uniform random number in the range [0, 1] */
void   mp_set_from_random(MPNumber *z);

/* Sets z from a string representation in 'text'.
 * Returns true on success.
 */
bool   mp_set_from_string(const char *text, int default_base, MPNumber *z);

/* Returns x as a native single-precision floating point number */
float  mp_cast_to_float(const MPNumber *x);

/* Returns x as a native double-precision floating point number */
double mp_cast_to_double(const MPNumber *x);

/* Returns x as a native integer */
int64_t mp_cast_to_int(const MPNumber *x);

/* Returns x as a native unsigned integer */
uint64_t mp_cast_to_unsigned_int(const MPNumber *x);

/* Sets z = sin x */
void   mp_sin(const MPNumber *x, MPAngleUnit unit, MPNumber *z);

/* Sets z = cos x */
void   mp_cos(const MPNumber *x, MPAngleUnit unit, MPNumber *z);

/* Sets z = tan x */
void   mp_tan(const MPNumber *x, MPAngleUnit unit, MPNumber *z);

/* Sets z = sin⁻¹ x */
void   mp_asin(const MPNumber *x, MPAngleUnit unit, MPNumber *z);

/* Sets z = cos⁻¹ x */
void   mp_acos(const MPNumber *x, MPAngleUnit unit, MPNumber *z);

/* Sets z = tan⁻¹ x */
void   mp_atan(const MPNumber *x, MPAngleUnit unit, MPNumber *z);

/* Sets z = sinh x */
void   mp_sinh(const MPNumber *x, MPNumber *z);

/* Sets z = cosh x */
void   mp_cosh(const MPNumber *x, MPNumber *z);

/* Sets z = tanh x */
void   mp_tanh(const MPNumber *x, MPNumber *z);

/* Sets z = sinh⁻¹ x */
void   mp_asinh(const MPNumber *x, MPNumber *z);

/* Sets z = cosh⁻¹ x */
void   mp_acosh(const MPNumber *x, MPNumber *z);

/* Sets z = tanh⁻¹ x */
void   mp_atanh(const MPNumber *x, MPNumber *z);

/* Returns true if x is cannot be represented in a binary word of length 'wordlen' */
bool   mp_is_overflow(const MPNumber *x, int wordlen);

/* Sets z = boolean AND for each bit in x and z */
void   mp_and(const MPNumber *x, const MPNumber *y, MPNumber *z);

/* Sets z = boolean OR for each bit in x and z */
void   mp_or(const MPNumber *x, const MPNumber *y, MPNumber *z);

/* Sets z = boolean XOR for each bit in x and z */
void   mp_xor(const MPNumber *x, const MPNumber *y, MPNumber *z);

/* Sets z = boolean XNOR for each bit in x and z for word of length 'wordlen' */
void   mp_xnor(const MPNumber *x, const MPNumber *y, int wordlen, MPNumber *z);

/* Sets z = boolean NOT for each bit in x and z for word of length 'wordlen' */
void   mp_not(const MPNumber *x, int wordlen, MPNumber *z);

/* Sets z = x masked to 'wordlen' bits */
void   mp_mask(const MPNumber *x, int wordlen, MPNumber *z);

/* Sets z = x shifted by 'count' bits.  Positive shift increases the value, negative decreases */
void   mp_shift(const MPNumber *x, int count, MPNumber *z);

/* Sets z to be the ones complement of x for word of length 'wordlen' */
void   mp_ones_complement(const MPNumber *x, int wordlen, MPNumber *z);

/* Sets z to be the twos complement of x for word of length 'wordlen' */
void   mp_twos_complement(const MPNumber *x, int wordlen, MPNumber *z);

#endif /* MP_H */