cgit.sxemacs.org Git - sxemacs/blob - src/ent/ent-gmp.c

   1 /*
   2   ent-gmp.c -- Numeric types for SXEmacs
   3   Copyright (C) 2004 Jerry James
   4   Copyright (C) 2004, 2005, 2006 Sebastian Freundt
   5
   6   Author:  Jerry James
   7   Backport:  Sebastian Freundt
   8
   9 This file is part of SXEmacs
  10
  11 SXEmacs is free software: you can redistribute it and/or modify
  12 it under the terms of the GNU General Public License as published by
  13 the Free Software Foundation, either version 3 of the License, or
  14 (at your option) any later version.
  15
  16 SXEmacs is distributed in the hope that it will be useful,
  17 but WITHOUT ANY WARRANTY; without even the implied warranty of
  18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  19 GNU General Public License for more details.
  20
  21 You should have received a copy of the GNU General Public License
  22 along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  23
  24
  25 #include <config.h>
  26 #include <limits.h>
  27 #include <math.h>
  28 #include "lisp.h"
  29 #include "sysproc.h"    /* For qxe_getpid */
  30
  31 #include "ent-gmp.h"
  32
  33 static mpf_t float_print_min, float_print_max;
  34 gmp_randstate_t random_state;
  35
  36 bigz ent_scratch_bigz;
  37 bigq ent_scratch_bigq;
  38 bigf ent_scratch_bigf;
  39
  40 static ase_nullary_operation_f Qent_mpz_zero, Qent_mpz_one;
  41 static ase_nullary_operation_f Qent_mpq_zero, Qent_mpq_one;
  42 static ase_nullary_operation_f Qent_mpf_zero, Qent_mpf_one;
  43
  44 \f
  45 #define yrealloc_array(ptr, type, len)          \
  46         ((void)(ptr = (type*)yrealloc(ptr, (len) * sizeof(type))))
  47
  48
  49 \f
  50 /************************* Big Rational Integers ****************************/
  51 static void
  52 bigz_print (Lisp_Object obj, Lisp_Object printcharfun, int SXE_UNUSED(escapeflag))
  53 {
  54         Bufbyte *bstr;
  55
  56         bstr = (Bufbyte*)bigz_to_string(XBIGZ_DATA(obj), 10);
  57         write_c_string((char*)bstr, printcharfun);
  58         xfree(bstr);
  59         bstr = (Bufbyte *)NULL;
  60 }
  61
  62 static int
  63 bigz_equal (Lisp_Object obj1, Lisp_Object obj2, int SXE_UNUSED(depth))
  64 {
  65         return bigz_eql(XBIGZ_DATA(obj1), XBIGZ_DATA(obj2));
  66 }
  67
  68 static unsigned long
  69 bigz_hash (Lisp_Object obj, int SXE_UNUSED(depth))
  70 {
  71         return (unsigned long)bigz_hashcode(XBIGZ_DATA(obj));
  72 }
  73
  74 static const struct lrecord_description bigz_description[] = {
  75         { XD_OPAQUE_DATA_PTR, offsetof(Lisp_Bigz, data) },
  76         { XD_END }
  77 };
  78
  79 DEFINE_BASIC_LRECORD_IMPLEMENTATION("bigz", bigz,
  80                                     NULL, bigz_print, NULL,
  81                                     bigz_equal, bigz_hash,
  82                                     bigz_description, Lisp_Bigz);
  83
  84 \f
  85 /************************** Rational Integer Fractions **********************/
  86 static void
  87 bigq_print (Lisp_Object obj, Lisp_Object printcharfun, int SXE_UNUSED(escapeflag))
  88 {
  89         Bufbyte *rstr;
  90
  91         rstr = (Bufbyte*)bigq_to_string(XBIGQ_DATA(obj), 10);
  92         write_c_string((char *)rstr, printcharfun);
  93         xfree(rstr);
  94         rstr = (Bufbyte *)NULL;
  95         return;
  96 }
  97
  98 static int
  99 bigq_equal (Lisp_Object obj1, Lisp_Object obj2, int SXE_UNUSED(depth))
 100 {
 101         return bigq_eql(XBIGQ_DATA(obj1), XBIGQ_DATA(obj2));
 102 }
 103
 104 static unsigned long
 105 bigq_hash (Lisp_Object obj, int SXE_UNUSED(depth))
 106 {
 107         return bigq_hashcode(XBIGQ_DATA(obj));
 108 }
 109
 110 static const struct lrecord_description bigq_description[] = {
 111         { XD_OPAQUE_DATA_PTR, offsetof (Lisp_Bigq, data) },
 112         { XD_END }
 113 };
 114
 115 DEFINE_BASIC_LRECORD_IMPLEMENTATION("bigq", bigq,
 116                                     NULL, bigq_print, NULL,
 117                                     bigq_equal, bigq_hash,
 118                                     bigq_description, Lisp_Bigq);
 119
 120 \f
 121 /********************************** Bigfs ***********************************/
 122 static void
 123 bigf_print(Lisp_Object obj, Lisp_Object printcharfun, int SXE_UNUSED(escapeflag))
 124 {
 125         Bufbyte *fstr = bigf_to_string(XBIGF_DATA(obj), 10);
 126         write_c_string((char*)fstr, printcharfun);
 127         xfree(fstr);
 128         fstr = (Bufbyte *)NULL;
 129         return;
 130 }
 131
 132 static int
 133 bigf_equal(Lisp_Object obj1, Lisp_Object obj2, int SXE_UNUSED(depth))
 134 {
 135         return bigf_eq(XBIGF_DATA(obj1), XBIGF_DATA(obj2));
 136 }
 137
 138 static unsigned long
 139 bigf_hash(Lisp_Object obj, int SXE_UNUSED(depth))
 140 {
 141         return bigf_hashcode(XBIGF_DATA(obj));
 142 }
 143
 144 static const struct lrecord_description bigf_description[] = {
 145         { XD_OPAQUE_DATA_PTR, offsetof(Lisp_Bigf, data) },
 146         { XD_END }
 147 };
 148
 149 DEFINE_BASIC_LRECORD_IMPLEMENTATION("bigf", bigf,
 150                                     NULL, bigf_print, NULL,
 151                                     bigf_equal, bigf_hash,
 152                                     bigf_description, Lisp_Bigf);
 153
 154 DEFUN("bigf-get-precision", Fbigf_get_precision, 1, 1, 0, /*
 155 Return the precision of bigf F as an integer.
 156 */
 157        (f))
 158 {
 159         CHECK_BIGF(f);
 160         return make_integer((signed long)XBIGF_GET_PREC(f));
 161 }
 162
 163 DEFUN("bigf-set-precision", Fbigf_set_precision, 2, 2, 0, /*
 164 Set the precision of F, a bigf, to PRECISION, a nonnegative integer.
 165 The new precision of F is returned.  Note that the return value may differ
 166 from PRECISION if the underlying library is unable to support exactly
 167 PRECISION bits of precision.
 168 */
 169        (f, precision))
 170 {
 171         unsigned long prec;
 172
 173         CHECK_BIGF(f);
 174         if (INTP(precision)) {
 175                 prec = (XINT(precision) <= 0)
 176                         ? 1UL : (unsigned long)XINT(precision);
 177         }
 178 #ifdef HAVE_MPZ
 179         else if (BIGZP(precision)) {
 180                 prec = bigz_fits_ulong_p(XBIGZ_DATA(precision))
 181                         ? bigz_to_ulong(XBIGZ_DATA(precision))
 182                         : UINT_MAX;
 183         }
 184 #endif  /* HAVE_MPZ */
 185         else {
 186                 dead_wrong_type_argument(Qintegerp, f);
 187                 return Qnil;
 188         }
 189
 190         XBIGF_SET_PREC(f, prec);
 191         return Fbigf_get_precision(f);
 192 }
 193
 194
 195 \f
 196 Bufbyte *
 197 bigf_to_string(mpf_t f, int base)
 198 {
 199         mp_exp_t expt;
 200         Bufbyte *str = (Bufbyte*)mpf_get_str(NULL, &expt, base, 0, f);
 201         const int sign = mpf_sgn(f);
 202         const int neg = (sign < 0) ? 1 : 0;
 203         int len = strlen((char *)str) + 1;  /* Count the null terminator */
 204
 205         /* Move digits down to insert a radix point */
 206         if (expt <= 0) {
 207                 /* We need room for a radix point and leading zeroes */
 208                 const int space = -expt + 2;
 209                 xrealloc_array(str, Bufbyte, len + space);
 210                 memmove(&str[space + neg], &str[neg], len - neg);
 211                 memset(&str[neg], '0', space);
 212                 str[neg + 1] = '.';
 213                 len += space;
 214         } else if (expt < len) {
 215                 /* We just need room for a radix point */
 216                 xrealloc_array(str, Bufbyte, len + 1);
 217                 memmove(&str[expt + neg + 1],
 218                         &str[expt + neg],
 219                         len - (expt + neg));
 220                 str[expt + neg] = '.';
 221                 len++;
 222         } else {
 223                 /* We need room for trailing zeroes */
 224                 xrealloc_array(str, Bufbyte, expt + 1);
 225                 memset(&str[len-1], '0', expt+2-len);
 226                 str[expt] = '\0';
 227                 len = expt + 1;
 228         }
 229 #if 0
 230         /* never want this here */
 231         /* Computerized scientific notation */
 232         /* We need room for a radix point, format identifier, and exponent */
 233         const int space = (expt < 0)
 234                 ? (int)(log(-expt) / log(base)) + 3
 235                 : (int)(log(expt) / log(base)) + 2;
 236         xrealloc_array(str, Bufbyte, len + space);
 237         memmove(&str[neg + 2], &str[neg + 1], len - neg);
 238         str[len + 1] = 'l';
 239         sprintf ((char *)&str[len + 2], "%ld", expt);
 240 }
 241 #endif
 242         return str;
 243 }
 244
 245 \f
 246 /* reader funs */
 247
 248 Lisp_Object
 249 read_bigz_string(const char *cp, int base)
 250 {
 251         Lisp_Object result;
 252         bigz bz;
 253
 254         bigz_init(bz);
 255
 256         /* MPZ bigz_set_string has no effect
 257          * with initial + sign */
 258         if (*cp == '+') {
 259                 cp++;
 260         }
 261
 262         bigz_set_string(bz, (const char*)cp, base);
 263         result = make_bigz_bz(bz);
 264
 265         bigz_fini(bz);
 266         return result;
 267 }
 268
 269
 270 Lisp_Object read_bigq_string(char *cp)
 271 {
 272         bigq bq;
 273         Lisp_Object result;
 274
 275         bigq_init(bq);
 276
 277         /* The GMP version of bigq_set_string (mpq_set_str) has the following
 278            limitations:
 279            - If p starts with a '+' sign, it does nothing; i.e., it leaves its
 280            ratio argument untouched.
 281            - If p has a '+' sign after the '/' (e.g., 300/+400), it sets the
 282            numerator from the string, but *leaves the denominator unchanged*.
 283            - If p has trailing nonnumeric characters, it sets the numerator from
 284            the string, but leaves the denominator unchanged.
 285            - If p has more than one '/', (e.g., 1/2/3), then it sets the
 286            numerator from the string, but leaves the denominator unchanged.
 287
 288            Therefore, move p past any leading '+' signs, temporarily drop a null
 289            after the numeric characters we are trying to convert, and then put
 290            the nulled character back afterward.  I am not going to fix problem
 291            #2; just don't write ratios that look like that. */
 292         if (*cp == '+') {
 293                 cp++;
 294         }
 295
 296         bigq_set_string(bq, cp, 0);
 297         bigq_canonicalize(bq);
 298
 299         result = ent_mpq_downgrade_maybe(bq);
 300
 301         bigq_fini(bq);
 302         return result;
 303 }
 304
 305 Lisp_Object read_bigf_string(char *cp)
 306 {
 307         bigf bf;
 308         Lisp_Object result;
 309
 310         bigf_init_prec(bf, internal_get_precision(Qnil));
 311
 312         /* The GMP version of bigf_set_string (mpf_set_str)
 313            has the following limitation: if p starts with a '+'
 314            sign, it does nothing; i.e., it leaves its bigfloat
 315            argument untouched.
 316            Therefore, move p past any leading '+' signs. */
 317         if (*cp == '+') {
 318                 cp++;
 319         }
 320
 321         bigf_set_string(bf, cp, 0);
 322         result = make_bigf_bf(bf);
 323
 324         bigf_fini(bf);
 325         return result;
 326 }
 327
 328 /* bigz ops */
 329 static inline int
 330 _ent_mpz_zerop(bigz n)
 331 {
 332         return (bigz_sign(n) == 0);
 333 }
 334
 335 static int
 336 ent_mpz_zerop(Lisp_Object l)
 337 {
 338         return _ent_mpz_zerop(XBIGZ_DATA(l));
 339 }
 340
 341 static inline int
 342 _ent_mpz_onep(bigz n)
 343 {
 344         return (bigz_fits_long_p(n) && bigz_to_long(n) == 1L);
 345 }
 346
 347 static int
 348 ent_mpz_onep(Lisp_Object l)
 349 {
 350         return _ent_mpz_onep(XBIGZ_DATA(l));
 351 }
 352
 353 static inline int
 354 _ent_mpz_unitp(bigz n)
 355 {
 356         return (bigz_fits_long_p(n) &&
 357                 (bigz_to_long(n) == 1L || bigz_to_long(n) == -1L));
 358 }
 359
 360 static int
 361 ent_mpz_unitp(Lisp_Object l)
 362 {
 363         return _ent_mpz_unitp(XBIGZ_DATA(l));
 364 }
 365
 366 static Lisp_Object
 367 ent_sum_BIGZ_T(Lisp_Object l, Lisp_Object r)
 368 {
 369         bigz_add(ent_scratch_bigz, XBIGZ_DATA(l), XBIGZ_DATA(r));
 370         return ent_mpz_downgrade_maybe(ent_scratch_bigz);
 371 }
 372 static Lisp_Object
 373 ent_sum_BIGZ_T_INT_T(Lisp_Object l, Lisp_Object r)
 374 {
 375         bigz_set_long(ent_scratch_bigz, ent_int(r));
 376         bigz_add(ent_scratch_bigz, XBIGZ_DATA(l), ent_scratch_bigz);
 377         return ent_mpz_downgrade_maybe(ent_scratch_bigz);
 378 }
 379 static Lisp_Object
 380 ent_sum_INT_T_BIGZ_T(Lisp_Object l, Lisp_Object r)
 381 {
 382         return ent_sum_BIGZ_T_INT_T(r, l);
 383 }
 384 #ifdef HAVE_FPFLOAT
 385 static Lisp_Object
 386 ent_sum_BIGZ_T_FLOAT_T(Lisp_Object l, Lisp_Object r)
 387 {
 388         return __ent_binop_lift_1(
 389                 ASE_BINARY_OP_SUM, BIGZ_T, l, FLOAT_T, r, NULL);
 390 }
 391 static Lisp_Object
 392 ent_sum_FLOAT_T_BIGZ_T(Lisp_Object l, Lisp_Object r)
 393 {
 394         return __ent_binop_lift_2(
 395                 ASE_BINARY_OP_SUM, FLOAT_T, l, BIGZ_T, r, NULL);
 396 }
 397 #endif
 398
 399 static Lisp_Object
 400 ent_diff_BIGZ_T(Lisp_Object l, Lisp_Object r)
 401 {
 402         bigz_sub(ent_scratch_bigz, XBIGZ_DATA(l), XBIGZ_DATA(r));
 403         return ent_mpz_downgrade_maybe(ent_scratch_bigz);
 404 }
 405 static Lisp_Object
 406 ent_diff_BIGZ_T_INT_T(Lisp_Object l, Lisp_Object r)
 407 {
 408         bigz_set_long(ent_scratch_bigz, ent_int(r));
 409         bigz_sub(ent_scratch_bigz, XBIGZ_DATA(l), ent_scratch_bigz);
 410         return ent_mpz_downgrade_maybe(ent_scratch_bigz);
 411 }
 412 static Lisp_Object
 413 ent_diff_INT_T_BIGZ_T(Lisp_Object l, Lisp_Object r)
 414 {
 415         bigz_set_long(ent_scratch_bigz, ent_int(l));
 416         bigz_sub(ent_scratch_bigz, ent_scratch_bigz, XBIGZ_DATA(r));
 417         return ent_mpz_downgrade_maybe(ent_scratch_bigz);
 418 }
 419 #ifdef HAVE_FPFLOAT
 420 static Lisp_Object
 421 ent_diff_BIGZ_T_FLOAT_T(Lisp_Object l, Lisp_Object r)
 422 {
 423         return __ent_binop_lift_1(
 424                 ASE_BINARY_OP_DIFF, BIGZ_T, l, FLOAT_T, r, NULL);
 425 }
 426 static Lisp_Object
 427 ent_diff_FLOAT_T_BIGZ_T(Lisp_Object l, Lisp_Object r)
 428 {
 429         return __ent_binop_lift_2(
 430                 ASE_BINARY_OP_DIFF, FLOAT_T, l, BIGZ_T, r, NULL);
 431 }
 432 #endif
 433
 434 static Lisp_Object
 435 ent_neg_BIGZ_T(Lisp_Object l)
 436 {
 437         bigz_neg(ent_scratch_bigz, XBIGZ_DATA(l));
 438         return make_bigz_bz(ent_scratch_bigz);
 439 }
 440
 441 static Lisp_Object
 442 ent_prod_BIGZ_T(Lisp_Object l, Lisp_Object r)
 443 {
 444         bigz_mul(ent_scratch_bigz, XBIGZ_DATA(l), XBIGZ_DATA(r));
 445         return ent_mpz_downgrade_maybe(ent_scratch_bigz);
 446 }
 447 static Lisp_Object
 448 ent_prod_BIGZ_T_INT_T(Lisp_Object l, Lisp_Object r)
 449 {
 450         bigz_set_long(ent_scratch_bigz, ent_int(r));
 451         bigz_mul(ent_scratch_bigz, XBIGZ_DATA(l), ent_scratch_bigz);
 452         return ent_mpz_downgrade_maybe(ent_scratch_bigz);
 453 }
 454 static Lisp_Object
 455 ent_prod_INT_T_BIGZ_T(Lisp_Object l, Lisp_Object r)
 456 {
 457         return ent_prod_BIGZ_T_INT_T(r, l);
 458 }
 459 #ifdef HAVE_FPFLOAT
 460 static Lisp_Object
 461 ent_prod_BIGZ_T_FLOAT_T(Lisp_Object l, Lisp_Object r)
 462 {
 463         return __ent_binop_lift_1(
 464                 ASE_BINARY_OP_PROD, BIGZ_T, l, FLOAT_T, r, NULL);
 465 }
 466 static Lisp_Object
 467 ent_prod_FLOAT_T_BIGZ_T(Lisp_Object l, Lisp_Object r)
 468 {
 469         return __ent_binop_lift_2(
 470                 ASE_BINARY_OP_PROD, FLOAT_T, l, BIGZ_T, r, NULL);
 471 }
 472 #endif
 473
 474 static Lisp_Object
 475 ent_div_BIGZ_T(Lisp_Object l, Lisp_Object r)
 476 {
 477         if (bigz_sign(XBIGZ_DATA(r)) == 0) {
 478                 int lsgn = bigz_sign(XBIGZ_DATA(l));
 479                 if (lsgn > 0)
 480                         return make_indef(POS_INFINITY);
 481                 else if (lsgn < 0)
 482                         return make_indef(NEG_INFINITY);
 483                 else
 484                         return make_indef(NOT_A_NUMBER);
 485         }
 486         bigz_div(ent_scratch_bigz, XBIGZ_DATA(l), XBIGZ_DATA(r));
 487         return ent_mpz_downgrade_maybe(ent_scratch_bigz);
 488 }
 489 static Lisp_Object
 490 ent_div_BIGZ_T_INT_T(Lisp_Object l, Lisp_Object r)
 491 {
 492         if (ent_int(r) == 0) {
 493                 int lsgn = bigz_sign(XBIGZ_DATA(l));
 494                 if (lsgn > 0)
 495                         return make_indef(POS_INFINITY);
 496                 else if (lsgn < 0)
 497                         return make_indef(NEG_INFINITY);
 498                 else
 499                         return make_indef(NOT_A_NUMBER);
 500         }
 501
 502         bigz_set_long(ent_scratch_bigz, ent_int(r));
 503         bigz_div(ent_scratch_bigz, XBIGZ_DATA(l), ent_scratch_bigz);
 504         return ent_mpz_downgrade_maybe(ent_scratch_bigz);
 505 }
 506 static Lisp_Object
 507 ent_div_INT_T_BIGZ_T(Lisp_Object l, Lisp_Object r)
 508 {
 509         if (bigz_sign(XBIGZ_DATA(r)) == 0) {
 510                 EMACS_INT rl = ent_int(l);
 511                 if (rl > 0)
 512                         return make_indef(POS_INFINITY);
 513                 else if (rl < 0)
 514                         return make_indef(NEG_INFINITY);
 515                 else
 516                         return make_indef(NOT_A_NUMBER);
 517         }
 518
 519         bigz_set_long(ent_scratch_bigz, ent_int(l));
 520         bigz_div(ent_scratch_bigz, ent_scratch_bigz, XBIGZ_DATA(r));
 521         return ent_mpz_downgrade_maybe(ent_scratch_bigz);
 522 }
 523 #ifdef HAVE_FPFLOAT
 524 static Lisp_Object
 525 ent_div_BIGZ_T_FLOAT_T(Lisp_Object l, Lisp_Object r)
 526 {
 527         return __ent_binop_lift_1(
 528                 ASE_BINARY_OP_DIV, BIGZ_T, l, FLOAT_T, r, NULL);
 529 }
 530 static Lisp_Object
 531 ent_div_FLOAT_T_BIGZ_T(Lisp_Object l, Lisp_Object r)
 532 {
 533         return __ent_binop_lift_2(
 534                 ASE_BINARY_OP_DIV, FLOAT_T, l, BIGZ_T, r, NULL);
 535 }
 536 #endif
 537
 538 static Lisp_Object
 539 ent_quo_BIGZ_T(Lisp_Object l, Lisp_Object r)
 540 {
 541         if (bigz_sign(XBIGZ_DATA(r)) == 0) {
 542                 int lsgn = bigz_sign(XBIGZ_DATA(l));
 543                 if (lsgn > 0)
 544                         return make_indef(POS_INFINITY);
 545                 else if (lsgn < 0)
 546                         return make_indef(NEG_INFINITY);
 547                 else
 548                         return make_indef(NOT_A_NUMBER);
 549         }
 550
 551         bigq_set_bigz_bigz(ent_scratch