Mercurial > hg > xemacs-beta

comparison src/number-mp.c @ 1983:9c872f33ecbe

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[xemacs-hg @ 2004-04-05 22:49:31 by james] Add bignum, ratio, and bigfloat support.
author james
date Mon, 05 Apr 2004 22:50:11 +0000
parents
children 4f5e0297b73f
comparison
equal deleted inserted replaced
1982:a748951fd4fb 1983:9c872f33ecbe
1 /* Numeric types for XEmacs using the MP library.
2 Copyright (C) 2004 Jerry James.
3
4 This file is part of XEmacs.
5
6 XEmacs is free software; you can redistribute it and/or modify it
7 under the terms of the GNU General Public License as published by the
8 Free Software Foundation; either version 2, or (at your option) any
9 later version.
10
11 XEmacs is distributed in the hope that it will be useful, but WITHOUT
12 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14 for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with XEmacs; see the file COPYING. If not, write to
18 the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
19 Boston, MA 02111-1307, USA. */
20
21 /* Synched up with: Not in FSF. */
22
23 #include <config.h>
24 #include <limits.h>
25 #include <math.h>
26 #include "lisp.h"
27
28 static MINT *bignum_bytesize, *bignum_long_sign_bit, *bignum_one, *bignum_two;
29 MINT *bignum_zero, *intern_bignum;
30 MINT *bignum_min_int, *bignum_max_int, *bignum_max_uint;
31 MINT *bignum_min_long, *bignum_max_long, *bignum_max_ulong;
32 short div_rem;
33
34 char *
35 bignum_to_string (bignum b, int base)
36 {
37 REGISTER unsigned int i;
38 unsigned int bufsize = 128U, index = 0U;
39 int sign;
40 char *buffer = xnew_array (char, 128), *retval;
41 MINT *quo = MP_ITOM (0);
42 short rem;
43
44 /* FIXME: signal something if base is < 2 or doesn't fit into a short. */
45
46 /* Save the sign for later */
47 sign = MP_MCMP (b, bignum_zero);
48
49 if (sign == 0)
50 {
51 XREALLOC_ARRAY (buffer, char, 2);
52 buffer[0] = '0';
53 buffer[1] = '\0';
54 return buffer;
55 }
56 /* Copy abs(b) into quo for destructive modification */
57 else if (sign < 0)
58 MP_MSUB (bignum_zero, b, quo);
59 else
60 MP_MOVE (b, quo);
61
62 quo = MP_ITOM (0);
63
64 /* Loop over the digits of b (in BASE) and place each one into buffer */
65 for (i = 0U; MP_MCMP(quo, bignum_zero) > 0; i++)
66 {
67 MP_SDIV (quo, base, quo, &rem);
68 if (index == bufsize)
69 {
70 bufsize <<= 1;
71 XREALLOC_ARRAY (buffer, char, bufsize);
72 }
73 buffer[index++] = rem < 10 ? rem + '0' : rem - 10 + 'a';
74 }
75 MP_MFREE (quo);
76
77 /* Reverse the digits, maybe add a minus sign, and add a null terminator */
78 bufsize = index + (sign < 0 ? 1 : 0) + 1;
79 retval = xnew_array (char, bufsize);
80 if (sign < 0)
81 {
82 retval[0] = '-';
83 i = 1;
84 }
85 else
86 i = 0;
87 for (; i < bufsize - 1; i++)
88 retval[i] = buffer[--index];
89 retval[bufsize - 1] = '\0';
90 xfree (buffer, char *);
91 return retval;
92 }
93
94 #define BIGNUM_TO_TYPE(type,accumtype) do { \
95 MP_MULT (b, quo, quo); \
96 for (i = 0U; i < sizeof(type); i++) \
97 { \
98 MP_SDIV (quo, 256, quo, &rem); \
99 retval |= ((accumtype) rem) << (8 * i); \
100 } \
101 MP_MFREE (quo); \
102 } while (0)
103
104 int
105 bignum_to_int (bignum b)
106 {
107 short rem, sign;
108 unsigned int retval = 0;
109 REGISTER unsigned int i;
110 MINT *quo;
111
112 sign = MP_MCMP (b, bignum_zero) < 0 ? -1 : 1;
113 quo = MP_ITOM (sign);
114 BIGNUM_TO_TYPE (int, unsigned int);
115 return ((int) retval) * sign;
116 }
117
118 unsigned int
119 bignum_to_uint (bignum b)
120 {
121 short rem;
122 unsigned int retval = 0U;
123 REGISTER unsigned int i;
124 MINT *quo;
125
126 quo = MP_ITOM (MP_MCMP (b, bignum_zero) < 0 ? -1 : 1);
127 BIGNUM_TO_TYPE (unsigned int, unsigned int);
128 return retval;
129 }
130
131 long
132 bignum_to_long (bignum b)
133 {
134 short rem, sign;
135 unsigned long retval = 0L;
136 REGISTER unsigned int i;
137 MINT *quo;
138
139 sign = MP_MCMP (b, bignum_zero) < 0 ? -1 : 1;
140 quo = MP_ITOM (sign);
141 BIGNUM_TO_TYPE (long, unsigned long);
142 return ((long) retval) * sign;
143 }
144
145 unsigned long
146 bignum_to_ulong (bignum b)
147 {
148 short rem;
149 unsigned long retval = 0UL;
150 REGISTER unsigned int i;
151 MINT *quo;
152
153 quo = MP_ITOM (MP_MCMP (b, bignum_zero) < 0 ? -1 : 1);
154 BIGNUM_TO_TYPE (unsigned long, unsigned long);
155 return retval;
156 }
157
158 double
159 bignum_to_double (bignum b)
160 {
161 short rem, sign;
162 double retval = 0.0;
163 REGISTER unsigned int i;
164 MINT *quo;
165
166 sign = MP_MCMP (b, bignum_zero) < 0 ? -1 : 1;
167 quo = MP_ITOM (sign);
168 MP_MULT (b, quo, quo);
169 for (i = 0U; MP_MCMP(quo, bignum_zero) > 0; i++)
170 {
171 MP_SDIV (quo, 256, quo, &rem);
172 retval += rem * i * 256;
173 }
174 MP_MFREE (quo);
175 return retval * sign;
176 }
177
178 static short
179 char_to_number (char c)
180 {
181 if (c >= '0' && c <= '9')
182 return c - '0';
183 if (c >= 'a' && c <= 'z')
184 return c - 'a' + 10;
185 if (c >= 'A' && c <= 'Z')
186 return c - 'A' + 10;
187 return -1;
188 }
189
190 int
191 bignum_set_string (bignum b, const char *s, int base)
192 {
193 MINT *mbase;
194 short digit;
195
196 if (base == 0)
197 {
198 if (s[0] == '0' && (s[1] == 'x' || s[1] == 'X'))
199 {
200 base = 16;
201 s += 2;
202 }
203 else if (*s == '0')
204 {
205 base = 8;
206 s++;
207 }
208 else
209 base = 10;
210 }
211
212 /* FIXME: signal something if base is < 2 or doesn't fit into a short. */
213
214 mbase = MP_ITOM ((short) base);
215 MP_MOVE (bignum_zero, b);
216
217 for (digit = char_to_number (*s); digit >= 0 && digit < base;
218 digit = char_to_number (*++s))
219 {
220 MINT *temp;
221
222 MP_MULT (b, mbase, b);
223 temp = MP_ITOM (digit);
224 MP_MADD (b, temp, b);
225 MP_MFREE (temp);
226 }
227
228 return (digit >= 0) ? -1 : 0;
229 }
230
231 void
232 bignum_set_long (MINT *b, long l)
233 {
234 /* Negative l is hard, not least because -LONG_MIN == LONG_MIN. We pretend
235 that l is unsigned, then subtract off the amount equal to the sign bit. */
236 bignum_set_ulong (b, (unsigned long) l);
237 if (l < 0L)
238 MP_MSUB (b, bignum_long_sign_bit, b);
239 }
240
241 void
242 bignum_set_ulong (bignum b, unsigned long l)
243 {
244 REGISTER unsigned int i;
245 MINT *multiplier = MP_ITOM (1);
246
247 MP_MOVE (bignum_zero, b);
248 for (i = 0UL; l > 0UL; l >>= 8, i++)
249 {
250 MINT *temp = MP_ITOM ((short) (l & 255));
251 MP_MULT (multiplier, temp, temp);
252 MP_MADD (b, temp, b);
253 MP_MULT (multiplier, bignum_bytesize, multiplier);
254 MP_MFREE (temp);
255 }
256 MP_MFREE (multiplier);
257 }
258
259 void
260 bignum_set_double (bignum b, double d)
261 {
262 REGISTER unsigned int i;
263 int negative = (d < 0) ? 1 : 0;
264 MINT *multiplier = MP_ITOM (1);
265
266 MP_MOVE (bignum_zero, b);
267 if (negative)
268 d = -d;
269 for (i = 0UL; d > 0.0; d /= 256, i++)
270 {
271 MINT *temp = MP_ITOM ((short) fmod (d, 256.0));
272 MP_MULT (multiplier, temp, temp);
273 MP_MADD (b, temp, b);
274 MP_MULT (multiplier, bignum_bytesize, multiplier);
275 MP_MFREE (temp);
276 }
277 MP_MFREE (multiplier);
278 if (negative)
279 MP_MSUB (bignum_zero, b, b);
280 }
281
282 /* Return nonzero if b1 is exactly divisible by b2 */
283 int
284 bignum_divisible_p (bignum b1, bignum b2)
285 {
286 int retval;
287 MINT *rem = MP_ITOM (0);
288 MP_MDIV (b1, b2, intern_bignum, rem);
289 retval = (MP_MCMP (rem, bignum_zero) == 0);
290 MP_MFREE (rem);
291 return retval;
292 }
293
294 void bignum_ceil (bignum quotient, bignum N, bignum D)
295 {
296 MP_MDIV (N, D, quotient, intern_bignum);
297 if (MP_MCMP (intern_bignum, bignum_zero) > 0 &&
298 MP_MCMP (quotient, bignum_zero) > 0)
299 MP_MADD (quotient, bignum_one, quotient);
300 }
301
302 void bignum_floor (bignum quotient, bignum N, bignum D)
303 {
304 MP_MDIV (N, D, quotient, intern_bignum);
305 if (MP_MCMP (intern_bignum, bignum_zero) > 0 &&
306 MP_MCMP (quotient, bignum_zero) < 0)
307 MP_MSUB (quotient, bignum_one, quotient);
308 }
309
310 /* RESULT = N to the POWth power */
311 void
312 bignum_pow (bignum result, bignum n, unsigned long pow)
313 {
314 MP_MOVE (bignum_one, result);
315 for ( ; pow > 0UL; pow--)
316 MP_MULT (result, n, result);
317 }
318
319 /* lcm(b1,b2) = b1 * b2 / gcd(b1, b2) */
320 void
321 bignum_lcm (bignum result, bignum b1, bignum b2)
322 {
323 MP_MULT (b1, b2, result);
324 MP_GCD (b1, b2, intern_bignum);
325 MP_MDIV (result, intern_bignum, result, intern_bignum);
326 }
327
328 /* FIXME: We can't handle negative args, so right now we just make them
329 positive before doing anything else. How should we really handle negative
330 args? */
331 #define bignum_bit_op(result, b1, b2, op) \
332 REGISTER unsigned int i; \
333 MINT *multiplier = MP_ITOM (1), *n1 = MP_ITOM (0), *n2 = MP_ITOM (0); \
334 \
335 if (MP_MCMP (bignum_zero, b1) > 0) \
336 MP_MSUB (bignum_zero, b1, n1); \
337 else \
338 MP_MOVE (b1, n1); \
339 if (MP_MCMP (bignum_zero, b2) > 0) \
340 MP_MSUB (bignum_zero, b2, n2); \
341 else \
342 MP_MOVE (b2, n2); \
343 \
344 MP_MOVE (bignum_zero, result); \
345 \
346 for (i = 0UL; MP_MCMP (bignum_zero, n1) < 0 && \
347 MP_MCMP (bignum_zero, n2) < 0; i++) \
348 { \
349 short byte1, byte2; \
350 MINT *temp; \
351 \
352 MP_SDIV (n1, 256, n1, &byte1); \
353 MP_SDIV (n2, 256, n2, &byte2); \
354 temp = MP_ITOM (byte1 op byte2); \
355 MP_MULT (multiplier, temp, temp); \
356 MP_MADD (result, temp, result); \
357 MP_MULT (multiplier, bignum_bytesize, multiplier); \
358 MP_MFREE (temp); \
359 } \
360 MP_MFREE (n2); \
361 MP_MFREE (n1); \
362 MP_MFREE (multiplier)
363
364 void
365 bignum_and (bignum result, bignum b1, bignum b2)
366 {
367 bignum_bit_op (result, b1, b2, &);
368 }
369
370 void
371 bignum_ior (bignum result, bignum b1, bignum b2)
372 {
373 bignum_bit_op (result, b1, b2, |);
374 }
375
376 void
377 bignum_xor (bignum result, bignum b1, bignum b2)
378 {
379 bignum_bit_op (result, b1, b2, ^);
380 }
381
382 /* NOT is not well-defined for bignums ... where do you stop flipping bits?
383 We just flip until we see the last one. This is probably a bad idea. */
384 void
385 bignum_not (bignum result, bignum b)
386 {
387 REGISTER unsigned int i;
388 MINT *multiplier = MP_ITOM (1), *n = MP_ITOM (0);
389
390 if (MP_MCMP (bignum_zero, b) > 0)
391 MP_MSUB (bignum_zero, b, n);
392 else
393 MP_MOVE (b, n);
394
395 MP_MOVE (bignum_zero, result);
396
397 for (i = 0UL; MP_MCMP (bignum_zero, n) < 0; i++)
398 {
399 short byte;
400 MINT *temp;
401
402 MP_SDIV (n, 256, n, &byte);
403 temp = MP_ITOM (~byte);
404 MP_MULT (multiplier, temp, temp);
405 MP_MADD (result, temp, result);
406 MP_MULT (multiplier, bignum_bytesize, multiplier);
407 MP_MFREE (temp);
408 }
409 MP_MFREE (n);
410 MP_MFREE (multiplier);
411 }
412
413 void
414 bignum_setbit (bignum b, unsigned long bit)
415 {
416 bignum_pow (intern_bignum, bignum_two, bit);
417 bignum_ior (b, b, intern_bignum);
418 }
419
420 /* This is so evil, even I feel queasy. */
421 void
422 bignum_clrbit (bignum b, unsigned long bit)
423 {
424 MINT *num = MP_ITOM (0);
425
426 /* See if the bit is already set, and subtract it off if not */
427 MP_MOVE (b, intern_bignum);
428 bignum_pow (num, bignum_two, bit);
429 bignum_ior (intern_bignum, intern_bignum, num);
430 if (MP_MCMP (b, intern_bignum) == 0)
431 MP_MSUB (b, num, b);
432 MP_MFREE (num);
433 }
434
435 int
436 bignum_testbit (bignum b, unsigned long bit)
437 {
438 bignum_pow (intern_bignum, bignum_two, bit);
439 bignum_and (intern_bignum, b, intern_bignum);
440 return MP_MCMP (intern_bignum, bignum_zero);
441 }
442
443 void
444 bignum_lshift (bignum result, bignum b, unsigned long bits)
445 {
446 bignum_pow (intern_bignum, bignum_two, bits);
447 MP_MULT (b, intern_bignum, result);
448 }
449
450 void
451 bignum_rshift (bignum result, bignum b, unsigned long bits)
452 {
453 bignum_pow (intern_bignum, bignum_two, bits);
454 MP_MDIV (b, intern_bignum, result, intern_bignum);
455 }
456
457 void bignum_random_seed(unsigned long seed)
458 {
459 /* FIXME: Implement me */
460 }
461
462 void bignum_random(bignum result, bignum limit)
463 {
464 /* FIXME: Implement me */
465 MP_MOVE (bignum_zero, result);
466 }
467
468 void
469 init_number_mp ()
470 {
471 REGISTER unsigned int i;
472
473 bignum_zero = MP_ITOM (0);
474 bignum_one = MP_ITOM (1);
475 bignum_two = MP_ITOM (2);
476
477 /* intern_bignum holds throwaway values from macro expansions in
478 number-mp.h. Its value is immaterial. */
479 intern_bignum = MP_ITOM (0);
480
481 /* bignum_bytesize holds the number of bits in a byte. */
482 bignum_bytesize = MP_ITOM (256);
483
484 /* bignum_long_sign_bit holds an adjustment for negative longs. */
485 bignum_long_sign_bit = MP_ITOM (256);
486 for (i = 1UL; i < sizeof (long); i++)
487 MP_MULT (bignum_bytesize, bignum_long_sign_bit, bignum_long_sign_bit);
488
489 /* The MP interface only supports turning short ints into MINTs, so we have
490 to set these the hard way. */
491
492 bignum_min_int = MP_ITOM (0);
493 bignum_set_long (bignum_min_int, INT_MIN);
494
495 bignum_max_int = MP_ITOM (0);
496 bignum_set_long (bignum_max_int, INT_MAX);
497
498 bignum_max_uint = MP_ITOM (0);
499 bignum_set_ulong (bignum_max_uint, UINT_MAX);
500
501 bignum_min_long = MP_ITOM (0);
502 bignum_set_long (bignum_min_long, LONG_MIN);
503
504 bignum_max_long = MP_ITOM (0);
505 bignum_set_long (bignum_max_long, LONG_MAX);
506
507 bignum_max_ulong = MP_ITOM (0);
508 bignum_set_ulong (bignum_max_ulong, ULONG_MAX);
509 }