src/bignum.c

   1 /*
   2 Copyright (C) 2001-2006, The Perl Foundation.
   3 $Id$
   4
   5 =head1 NAME
   6
   7 src/bignum.c - A decimal arithmetic library Parrot
   8
   9 =head1 DESCRIPTION
  10
  11 This code is intended for inclusion in the parrot project, and also
  12 for backporting into Perl5 (as a CPAN module).  Any patches to this
  13 code will likely find their way back to the Mother Ship, as it were.
  14
  15 There is a good deal of scope for improving the speed of this code,
  16 modifications are encouraged as long as the extended regression tests
  17 continue to pass.
  18 Alex Gough, 2002
  19
  20 I<It was a very inconvenient habit of kittens (Alice had once made the
  21 remark) that, whatever you say to them, they always purr.  "If they
  22 would only purr for `yes', and mew for `no', or any rule of that sort",
  23 she had said, "so that one could keep up a conversation!  But how can
  24 you talk with a person if they always say the same thing?">
  25
  26 I<On this occasion the kitten only purred: and it was impossible to
  27 guess whether it meant `yes' or `no'.>
  28
  29 =head2 When in parrot
  30
  31 When the library is used within parrot, all calls expect an additional
  32 first argument of an interpreter, for the purposes of memory allocation,
  33 some internal macros do not (getd/setd and CHECK(O|U)FLOW.
  34
  35 If you're being useful and inserting proper rapid fillins, start
  36 with the C<BN_i*> methods, but make sure any errors can still be
  37 thrown.
  38
  39 =head2 Macros
  40
  41 Access digits, macros assume length given.
  42
  43 =over 4
  44
  45 =cut
  46
  47 */
  48
  49 #include <stdio.h>
  50 #include "bignum.h"
  51 #include <string.h> /* XXX:ajg fixme later*/
  52
  53 /* * This lot wants to go in a (bignum specific seperate header file * */
  54
  55 /*
  56
  57 =item C<BN_setd(BIGNUM*, pos, value)>
  58
  59 Set digit at C<pos> (zero is lsd) to C<value>.
  60
  61 =item C<int BN_getd(BIGNUM*, pos)>
  62
  63 Get value of digit at C<pos>.
  64
  65 =cut
  66
  67 */
  68 #define BN_setd(x, y, z) \
  69     ((x)->buffer[(y) / BN_D_PER_NIB] = \
  70      ((z) << ((y) % BN_D_PER_NIB)*4) | \
  71      ((x)->buffer[(y) / BN_D_PER_NIB] & ~(15<< ((y) % BN_D_PER_NIB)*4)))
  72 #define BN_getd(x, y) \
  73     (((x)->buffer[(y) / BN_D_PER_NIB] >> \
  74      ((y) % BN_D_PER_NIB)*4) & 15)
  75
  76 /*
  77
  78 =item C<CHECK_OVERFLOW(bn, incr, context)>
  79
  80 If increasing the exponent of C<bn> by C<incr> will cause overflow (as
  81 decided by C<elimit>), returns true.
  82
  83 =cut
  84
  85 */
  86
  87 #define CHECK_OVERFLOW(bn, incr, context) \
  88     ((context)->elimit - ((bn)->expn + (bn)->digits -1) < (incr) ? 1 : 0)
  89
  90 /*
  91
  92 =item C<CHECK_UNDERFLOW(bn, decrement, context)>
  93
  94 If subtracting C<decrement> (a positive number) from the exponent
  95 of C<bn> would cause underflow, returns true.
  96
  97 =cut
  98
  99 */
 100
 101 #define CHECK_UNDERFLOW(bn, decr, context) \
 102     ((context)->elimit + ((bn)->expn + (bn)->digits -1) > (decr) ? 0 : 1)
 103
 104 /*
 105
 106 =back
 107
 108 Special Values
 109
 110 =over 4
 111
 112 =item C<am_INF(bn)>
 113
 114 True if C<bn> is +Infinity or -Infinity.
 115
 116 =cut
 117
 118 */
 119
 120 #define am_INF(bn)  ((bn)->flags & BN_INF_FLAG)
 121
 122 /*
 123
 124 =item C<am_NAN(bn)>
 125
 126 True if C<bn> is either a quiet or signalling NaN.
 127
 128 =cut
 129
 130 */
 131
 132 #define am_NAN(bn)  ((bn)->flags & (BN_sNAN_FLAG | BN_qNAN_FLAG))
 133
 134 /*
 135
 136 =item C<am_sNAN(bn)>
 137
 138 True if C<bn> is a signalling NaN.
 139
 140 =cut
 141
 142 */
 143
 144 #define am_sNAN(bn) ((bn)->flags & BN_sNAN_FLAG)
 145
 146 /*
 147
 148 =item C<am_qNAN(bn)>
 149
 150 True if C<bn> is a quiet NaN.
 151
 152 =cut
 153
 154 */
 155
 156 #define am_qNAN(bn) ((bn)->flags & BN_qNAN_FLAG)
 157
 158 /*
 159
 160 =item C<char* BN_lazydbprint(BIGNUM* foo)>
 161
 162 For the sake of debugging
 163
 164 =cut
 165
 166 */
 167
 168 char* BN_lazydbprint(BIGNUM* foo)
 169 {
 170     char*s;
 171     BN_to_scientific_string(foo, &s);
 172     return s;
 173 }
 174
 175 /* Internal functions + types */
 176 typedef enum {   /* Indicate to idivide when to stop */
 177     BN_DIV_DIVIDE,
 178     BN_DIV_DIVINT,
 179     BN_DIV_REMAIN
 180 } BN_DIV_ENUM;
 181
 182 /* Used to restore INTENT(IN) arguments to functions */
 183 typedef struct BN_SAVE_PREC {
 184     BIGNUM one;
 185     BIGNUM two;
 186 } BN_SAVE_PREC;
 187
 188 int
 189 BN_imultiply(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
 190               BN_CONTEXT *context);
 191 int
 192 BN_idivide(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
 193             BN_CONTEXT *context,
 194             BN_DIV_ENUM operation, BIGNUM* rem);
 195 int BN_iround(PINTD_ BIGNUM *bn, BN_CONTEXT* context);
 196 INTVAL BN_to_scieng_string(PINTD_ BIGNUM* bn, char **dest, int eng);
 197 int BN_strip_lead_zeros(PINTD_ BIGNUM* victim, BN_CONTEXT*);
 198 int BN_strip_tail_zeros(PINTD_ BIGNUM* victim, BN_CONTEXT*);
 199 int BN_round_up(PINTD_ BIGNUM *victim, BN_CONTEXT* context);
 200 int BN_round_down(PINTD_ BIGNUM *victim, BN_CONTEXT* context);
 201 int BN_make_integer(PINTD_ BIGNUM* bn, BN_CONTEXT* context);
 202 int BN_arith_setup(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
 203                    BN_CONTEXT *context, BN_SAVE_PREC* restore);
 204 int BN_arith_cleanup(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
 205                      BN_CONTEXT *context, BN_SAVE_PREC* restore);
 206 int BN_iadd(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
 207             BN_CONTEXT *context);
 208 int BN_isubtract(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
 209                  BN_CONTEXT *context);
 210 int BN_align(PINTD_ BIGNUM* one, BIGNUM* two);
 211 INTVAL BN_nonfatal(PINTD_ BN_CONTEXT *context, BN_EXCEPTIONS except,
 212                    const char *msg);
 213 int BN_set_verybig(PINTD_ BIGNUM* bn, BN_CONTEXT *context);
 214
 215 /*
 216
 217 =back
 218
 219 =head2 Creation and Memory Management Functions
 220
 221 =over 4
 222
 223 =item C<BIGNUM*
 224 BN_new(PINTD_ INTVAL length)>
 225
 226 Create a new C<BIGNUM>.  C<length> is number of I<decimal> digits
 227 required. The bignumber will be equal to zero.
 228
 229 =cut
 230
 231 */
 232
 233 BIGNUM*
 234 BN_new(PINTD_ INTVAL length)
 235 {
 236     BIGNUM* bn;
 237     bn = (BIGNUM*)BN_alloc(PINT_ sizeof (BIGNUM));
 238     if (NULL == bn) {
 239         BN_EXCEPT(PINT_ BN_INSUFFICIENT_STORAGE, "Cannot allocate new BigNum");
 240     }
 241     bn->nibs = 1 + length / BN_D_PER_NIB;
 242     bn->buffer = (BN_NIB*)BN_alloc(PINT_ sizeof (BN_NIB) * bn->nibs);
 243     if (NULL == bn->buffer) {
 244         BN_EXCEPT(PINT_ BN_INSUFFICIENT_STORAGE, "Cannot allocate new BigNum");
 245     }
 246     bn->sign = 0;
 247     bn->expn = 0;
 248     bn->digits = 1;
 249     bn->flags = 0;
 250
 251     return bn;
 252 }
 253
 254 /*
 255
 256 =item C<void
 257 BN_grow(PINTD_ NOTNULL(BIGNUM *in), INTVAL length)>
 258
 259 Grows bn so that it can contain C<length> I<decimal> digits, does not
 260 modify the value of the bignumber.
 261
 262 =cut
 263
 264 */
 265
 266 void
 267 BN_grow(PINTD_ NOTNULL(BIGNUM *in), INTVAL length)
 268 {
 269     PARROT_ASSERT(in);
 270     if (length <= in->nibs * BN_D_PER_NIB) {
 271         return;
 272     }
 273     if (length > BN_MAX_DIGITS) {
 274         BN_EXCEPT(PINT_ BN_OVERFLOW, "Attempt to grow BIGNUM beyond limits");
 275     }
 276     in->nibs = 1+ length / BN_D_PER_NIB;
 277     in->buffer = (BN_NIB*)BN_realloc(PINT_ in->buffer,
 278                                      sizeof (BN_NIB) *(in->nibs));
 279     if (NULL==in->buffer) {
 280         BN_EXCEPT(PINT_ BN_INSUFFICIENT_STORAGE, "Cannot grow BIGNUM");
 281     }
 282     return;
 283 }
 284
 285 /*
 286
 287 =item C<void
 288 BN_destroy(PINTD_ BIGNUM *bn)>
 289
 290 Frees all the memory used by the BIGNUM.
 291
 292 =cut
 293
 294 */
 295
 296 void
 297 BN_destroy(PINTD_ BIGNUM *bn)
 298 {
 299     PARROT_ASSERT(bn!=NULL);
 300
 301     BN_free(PINT_ bn->buffer);
 302     BN_free(PINT_ bn);
 303     return;
 304 }
 305
 306 /*
 307
 308 =item C<BN_CONTEXT*
 309 BN_create_context(PINTD_ INTVAL precision)>
 310
 311 Creates a new context object, with specified I<precision>, other fields
 312 are initialised as follows:
 313
 314  elimit   = BN_HARD_EXPN_LIMIT (defined during configure)
 315  rounding = ROUND_HALF_UP
 316  extended = 1
 317  flags    = 0
 318  traps    = Division by zero, invalid operation, overflow, underflow
 319             and rounded are enabled.
 320             Lost digits and inexact are disabled.
 321
 322 The context object can be destroyed with C<free()>.
 323
 324 =cut
 325
 326 */
 327
 328 BN_CONTEXT*
 329 BN_create_context(PINTD_ INTVAL precision)
 330 {
 331     BN_CONTEXT *result;
 332
 333     result = (BN_CONTEXT*)BN_alloc(PINT_ sizeof (BN_CONTEXT));
 334     if (result == NULL) {
 335         BN_EXCEPT(PINT_ BN_INSUFFICIENT_STORAGE,
 336                   "Cannot create a context");
 337     }
 338
 339     result->precision = precision;
 340
 341     result->elimit = BN_HARD_EXPN_LIMIT;
 342     result->rounding = ROUND_HALF_UP;
 343     result->extended = 1;
 344     result->flags = 0;
 345     result->traps = (BN_F_DIVISION_BY_ZERO  |
 346                       BN_F_EXCEPTION_INVALID_OPERATION |
 347                       BN_F_OVERFLOW          |
 348                       BN_F_ROUNDED           |
 349                       BN_F_UNDERFLOW);
 350     return result;
 351 }
 352
 353 /*
 354
 355 =item C<INTVAL
 356 BN_set_digit(PINT_ BIGNUM* bn, INTVAL pos, INTVAL value)>
 357
 358 Sets digit at C<pos> (zero based) to C<value>. Number is grown if digits
 359 > allocated space are accessed, but intermediate digits will have
 360 undefined values. If C<pos> is beyond C<digits> then C<digits> is also
 361 updated.
 362
 363 =cut
 364
 365 */
 366
 367 INTVAL
 368 BN_set_digit(PINT_ BIGNUM* bn, INTVAL pos, INTVAL value)
 369 {
 370     PARROT_ASSERT(bn != NULL);
 371     if (pos > bn->nibs * BN_D_PER_NIB) {
 372         BN_grow(bn, pos);
 373     }
 374     PARROT_ASSERT(value < 10);
 375     PARROT_ASSERT(value > -1);
 376     BN_setd(bn, pos, value);
 377     if (bn->digits < pos+1) {
 378         bn->digits = pos+1;
 379     }
 380     return value;
 381 }
 382
 383 /*
 384
 385 =item C<INTVAL
 386 BN_get_digit(PINTD_ BIGNUM* bn, INTVAL pos)>
 387
 388 Get the value of the decimal digit at C<pos>, returns -1 if C<pos> is
 389 out of bounds.
 390
 391 =cut
 392
 393 */
 394
 395 INTVAL
 396 BN_get_digit(PINTD_ BIGNUM* bn, INTVAL pos)
 397 {
 398     PARROT_ASSERT(bn != NULL);
 399     if (pos > bn->digits || pos < 0) return -1;
 400     return BN_getd(bn, pos);
 401 }
 402
 403 /*
 404
 405 =item C<int BN_set_inf(PINTD_ BIGNUM* bn)>
 406
 407 =item C<int BN_set_qNAN(PINTD_ BIGNUM* bn)>
 408
 409 =item C<int BN_set_sNAN(PINTD_ BIGNUM* bn)>
 410
 411 Sets its argument to appropriate value.
 412
 413 Infinity is represented as having zero digits, an undefined exponent
 414 and private C<flags> set to C<BN_inf_FLAGS>.
 415
 416 sNAN is represented as having zero digits, an undefined exponent, an
 417 undefined sign and both qNAN and sNAN bits set.
 418
 419 qNAN is represented as having zero digits, an undefined exponent
 420 and only the qNAN bit set.
 421
 422 =cut
 423
 424 */
 425
 426 int BN_set_inf(PINTD_ BIGNUM* bn)
 427 {
 428     PARROT_ASSERT(bn != NULL);
 429     bn->digits = 0;
 430     bn->flags = (bn->flags & (~(UINTVAL)255)) | BN_INF_FLAG;
 431     return;
 432 }
 433
 434 int BN_set_qNAN(PINTD_ BIGNUM* bn)
 435 {
 436     PARROT_ASSERT(bn != NULL);
 437     bn->digits = 0;
 438     bn->flags = (bn->flags & (~(UINTVAL)255)) | BN_qNAN_FLAG;
 439     return;
 440 }
 441
 442 int BN_set_sNAN(PINTD_ BIGNUM* bn)
 443 {
 444     PARROT_ASSERT(bn != NULL);
 445     bn->digits = 0;
 446     bn->flags = (bn->flags & (~(UINTVAL)255)) | BN_qNAN_FLAG | BN_sNAN_FLAG;
 447     return;
 448 }
 449
 450 /*
 451
 452 =item C<int
 453 BN_set_verybig(PINTD_ BIGNUM* bn, BN_CONTEXT *context)>
 454
 455 Used when an operation has overflowed, sets C<bn> according to
 456 C<< context->rounding >> and the sign of C<bn>:
 457
 458  ROUND_HALF_UP, ROUND_HALF_EVEN => sign Infinity
 459  ROUND_DOWN => sign, largest finite number in given precision (or Inf, if
 460                  infinite precision is specified)
 461  ROUND_CEILING => same as round down, if sign is 1, +Inf otherwise
 462  ROUND_FLOOR => same as round down, if sign is 0, -Inf otherwise
 463
 464 =cut
 465
 466 */
 467
 468 int
 469 BN_set_verybig(PINTD_ BIGNUM* bn, BN_CONTEXT *context)
 470 {
 471     int massive = 0; /* 0 => inf, 1=> 99999999999 etc...*/
 472     switch (context->rounding) {
 473       case ROUND_HALF_UP:
 474       case ROUND_HALF_EVEN:
 475         break;
 476       case ROUND_FLOOR:
 477         if (!bn->sign) massive = 1;
 478         break;
 479       case ROUND_CEILING:
 480         if (bn->sign) massive = 1;
 481         break;
 482       case ROUND_DOWN:
 483         massive = 1;
 484         break;
 485       default:
 486         BN_EXCEPT(PINT_ BN_EXCEPTION_INVALID_OPERATION,
 487                   "Unknown rounding during overflow");
 488     }
 489     if (context->precision > 0 && massive) {
 490         INTVAL i;
 491         BN_grow(PINT_ bn, context->precision);
 492         for (i = 0; i< context->precision; i++) {
 493             BN_setd(bn, i, 9);
 494         }
 495         bn->digits = context->precision;
 496         bn->expn = context->elimit - context->precision + 1;
 497     }
 498     else {
 499         BN_set_inf(PINT_ bn);
 500     }
 501 }
 502
 503 /*
 504
 505 =item C<BIGNUM*
 506 BN_copy(PINTD_ BIGNUM* one, BIGNUM* two)>
 507
 508 Copies two into one, returning one for convenience.
 509
 510 =cut
 511
 512 */
 513
 514 BIGNUM*
 515 BN_copy(PINTD_ BIGNUM* one, BIGNUM* two)
 516 {
 517     PARROT_ASSERT(one != NULL); PARROT_ASSERT(two != NULL);
 518
 519     BN_grow(PINT_ two, one->digits);
 520     memcpy((void*)one->buffer, (void*)two->buffer,
 521            (1+two->digits / BN_D_PER_NIB) * sizeof (BN_NIB));
 522     one->flags &= ~(UINTVAL)0xF;
 523     one->flags |= two->flags & 0xF;
 524     one->digits = two->digits;
 525     one->expn = two->expn;
 526     one->sign = two->sign;
 527     return one;
 528 }
 529
 530 /*
 531
 532 =item C<BIGNUM*
 533 BN_new_from_int(PINTD_ INTVAL value)>
 534
 535 Create a new bignum from a (signed) integer value (C<INTVAL>)
 536 We assume that the implementation limits are somewhat larger than
 537 those required to store a single integer into a bignum.
 538
 539 =cut
 540
 541 */
 542
 543 BIGNUM*
 544 BN_new_from_int(PINTD_ INTVAL value)
 545 {
 546     BIGNUM *new;
 547     int i, current;
 548     new = BN_new(PINT_ BN_D_PER_INT);
 549     if (value < 0) {
 550         new->sign = 1;
 551         value = -value;
 552     }
 553     i = 0;
 554     while (value) {
 555         current = value % 10;
 556         BN_setd(new, i, current);
 557         value = value / 10;
 558         i++;
 559     }
 560     new->digits = i;
 561     new->expn = 0;
 562     return new;
 563 }
 564
 565 /*
 566
 567 =item C<void
 568 BN_PRINT_DEBUG(BIGNUM *bn, char* mesg)>
 569
 570 Dump the bignum for testing, along with a little message.
 571
 572 =cut
 573
 574 */
 575
 576 void
 577 BN_PRINT_DEBUG(BIGNUM *bn, char* mesg)
 578 {
 579     INTVAL i;
 580     printf("%s: nibs %i digits %i sign %i expn %i \n", mesg,
 581            bn->nibs, bn->digits, bn->sign, bn->expn);
 582     if (bn->digits == 0) {
 583         printf("Special value, flags: %x", bn->flags & 127);
 584     }
 585     else {
 586         for (i=bn->digits-1; i>-1; i--) {
 587             printf("%d", BN_getd(bn, i));
 588             if (!(i%5)) printf(" ");
 589             if (!(i%70)) printf("\n");
 590         }
 591     }
 592     printf("\n");
 593 }
 594
 595 /*
 596
 597 =item C<INTVAL
 598 BN_nonfatal(PINTD_ BN_CONTEXT *context, BN_EXCEPTIONS except, const char *msg)>
 599
 600 When an exceptional condition occurs after which execution could
 601 continue.  If context specifies that death occurs, then so be it.
 602
 603 =cut
 604
 605 */
 606
 607 INTVAL
 608 BN_nonfatal(PINTD_ BN_CONTEXT *context, BN_EXCEPTIONS except, const char *msg)
 609 {
 610     /* See extended standard for details */
 611     switch (except) {
 612     case BN_CONVERSION_OVERFLOW :
 613         /* Asked to hold coeff|expn too large value */
 614         context->flags |= BN_F_OVERFLOW | BN_F_INEXACT | BN_F_ROUNDED;
 615         if (context->traps & (BN_F_OVERFLOW | BN_F_INEXACT | BN_F_ROUNDED)) {
 616             BN_EXCEPT(PINT_ except, msg);
 617         }
 618         break;
 619     case BN_CONVERSION_SYNTAX:
 620         /* string not conforming to numeric form */
 621         context->flags |= BN_F_EXCEPTION_INVALID_OPERATION;
 622         if (context->traps & (BN_F_EXCEPTION_INVALID_OPERATION)) {
 623             BN_EXCEPT(PINT_ except, msg);
 624         }
 625         break;
 626     case BN_CONVERSION_UNDERFLOW:
 627         /* expn of string too small to be held  */
 628         context->flags |= BN_F_UNDERFLOW;
 629         if (context->traps & (BN_F_UNDERFLOW)) {
 630             BN_EXCEPT(PINT_ except, msg);
 631         }
 632         break;
 633     case BN_DIVISION_BY_ZERO:
 634         /* dividend of div/div-int or pow zero  */
 635         context->flags |= BN_F_DIVISION_BY_ZERO;
 636         if (context->traps & (BN_F_DIVISION_BY_ZERO)) {
 637             BN_EXCEPT(PINT_ except, msg);
 638         }
 639         break;
 640     case BN_DIVISION_IMPOSSIBLE:
 641         /* integer result of div-int or rem > precision */
 642         context->flags |= BN_F_EXCEPTION_INVALID_OPERATION;
 643         if (context->traps & (BN_F_EXCEPTION_INVALID_OPERATION)) {
 644             BN_EXCEPT(PINT_ except, msg);
 645         }
 646         break;
 647     case BN_DIVISION_UNDEFINED:
 648         /* div by zero with zero on top also */
 649         context->flags |= BN_F_EXCEPTION_INVALID_OPERATION;
 650         if (context->traps & (BN_F_EXCEPTION_INVALID_OPERATION)) {
 651             BN_EXCEPT(PINT_ except, msg);
 652         }
 653         break;
 654     case BN_INEXACT:
 655         /* some sort of rounding: with loss of information */
 656         context->flags |= BN_F_INEXACT;
 657         if (context->traps & (BN_F_INEXACT)) {
 658             BN_EXCEPT(PINT_ except, msg);
 659         }
 660         break;
 661     case BN_INSUFFICIENT_STORAGE:
 662         /* not enough space to hold intermediate results */
 663         /* Often this will be raised directly (ie. fatally) */
 664         context->flags |= BN_F_EXCEPTION_INVALID_OPERATION;
 665         if (context->traps & (BN_F_EXCEPTION_INVALID_OPERATION)) {
 666             BN_EXCEPT(PINT_ except, msg);
 667         }
 668         break;
 669     case BN_INVALID_CONTEXT:
 670         /* context given was not valid (unknown round) */
 671         context->flags |= BN_F_EXCEPTION_INVALID_OPERATION;
 672         if (context->traps & (BN_F_EXCEPTION_INVALID_OPERATION)) {
 673             BN_EXCEPT(PINT_ except, msg);
 674         }
 675         break;
 676     case BN_EXCEPTION_INVALID_OPERATION:
 677         /* operation which is not valid */
 678         context->flags |= BN_F_EXCEPTION_INVALID_OPERATION;
 679         if (context->traps & (BN_F_EXCEPTION_INVALID_OPERATION)) {
 680             BN_EXCEPT(PINT_ except, msg);
 681         }
 682         break;
 683     case BN_LOST_DIGITS:
 684         /* digits lost in rounding  */
 685         context->flags |= BN_F_LOST_DIGITS;
 686         if (context->traps & (BN_F_LOST_DIGITS)) {
 687             BN_EXCEPT(PINT_ except, msg);
 688         }
 689         break;
 690     case BN_OVERFLOW:
 691         /* expn becomes larger than max allowed */
 692         context->flags |= BN_F_OVERFLOW | BN_F_ROUNDED | BN_F_INEXACT;
 693         if (context->traps & (BN_F_OVERFLOW | BN_F_ROUNDED | BN_F_INEXACT)) {
 694             BN_EXCEPT(PINT_ except, msg);
 695         }
 696         break;
 697     case BN_ROUNDED:
 698         /* something was rounded */
 699         context->flags |= BN_F_ROUNDED;
 700         if (context->traps & (BN_F_ROUNDED)) {
 701             BN_EXCEPT(PINT_ except, msg);
 702         }
 703         break;
 704     case BN_UNDERFLOW:
 705         /* expn becomes smaller than min allowed */
 706         context->flags |= BN_F_UNDERFLOW | BN_F_ROUNDED | BN_F_INEXACT;
 707         if (context->traps &(BN_F_UNDERFLOW | BN_F_ROUNDED | BN_F_INEXACT)) {
 708             BN_EXCEPT(PINT_ except, msg);
 709         }
 710         break;
 711     default:
 712         BN_EXCEPT(PINT_ BN_EXCEPTION_INVALID_OPERATION, "An unknown error occurred");
 713     }
 714
 715 }
 716
 717 /*
 718
 719 =item C<void
 720 BN_exception(PINTD_ BN_EXCEPTIONS exception, const char* message)>
 721
 722 Throw `exception'. Should be accessed via a C<BN_EXCEPT> macro, this
 723 version is provided until Parrot exceptions are sorted out properly.
 724
 725 =cut
 726
 727 */
 728
 729 void
 730 BN_exception(PINTD_ BN_EXCEPTIONS exception, const char* message)
 731 {
 732     printf("Exception %d %s\n", exception, message);
 733     exit(EXIT_SUCCESS);
 734 }
 735
 736 /*
 737
 738 =item C<INTVAL
 739 BN_to_scientific_string(PINTD_ BIGNUM*bn, char **dest)>
 740
 741 Converts bn into a scientific representation, stored in dest.
 742
 743 =item C<INTVAL
 744 BN_to_engineering_string(PINTD_ BIGNUM*bn, char **dest)>
 745
 746 Converts C<*bn> into a engineering representation, stored in C<**dest>.
 747
 748 These functions return C<char*> strings only, parrot may want to
 749 reimplement these so that locales and the like are nicely coped with.
 750
 751 Any reimplementation should be in a seperate file, this section of
 752 the main file can be C<#ifdef>ed out if this is done.
 753
 754 Memory pointed to by C<dest> is not freed by this function.
 755
 756 =cut
 757
 758 */
 759
 760
 761 INTVAL
 762 BN_to_scientific_string(PINTD_ BIGNUM*bn, char **dest)
 763 {
 764     BN_to_scieng_string(PINT_ bn, dest, 0);
 765 }
 766
 767 INTVAL
 768 BN_to_engineering_string(PINTD_ BIGNUM*bn, char **dest)
 769 {
 770     BN_to_scieng_string(PINT_ bn, dest, 1);
 771 }
 772
 773 /*
 774
 775 =item C<INTVAL
 776 BN_to_scieng_string(PINTD_ BIGNUM* bn, char **dest, int eng)>
 777
 778 Does the heavy string handling work, C<eng> defines the conversion to
 779 perform.
 780
 781 =cut
 782
 783 */
 784
 785 INTVAL
 786 BN_to_scieng_string(PINTD_ BIGNUM* bn, char **dest, int eng)
 787 {
 788     char* cur;
 789     INTVAL adj_exp = 0; /* as bn->expn is relative to 0th digit */
 790     INTVAL cur_dig = 0;
 791     PARROT_ASSERT(bn!=NULL);
 792     /* Special values have digits set to zero, so this should be enough */
 793     *dest = (char*)BN_alloc(PINT_ bn->digits + 5 + BN_D_PER_INT);
 794     if (dest == NULL) {
 795         BN_EXCEPT(PINT_ BN_INSUFFICIENT_STORAGE,
 796                   "Cannot create buffer to hold output string");
 797     }
 798
 799     cur = *dest;
 800
 801     /* Do we have a special value? */
 802     if (am_NAN(bn)) {
 803         if (am_qNAN(bn)) {
 804             strcpy(cur, "NaN");
 805         }
 806         else { /* must be signalling */
 807             strcpy(cur, "sNaN");
 808         }
 809         return 1;
 810     }
 811
 812     if (am_INF(bn)) {
 813         if (bn->sign) *cur++ = '-';
 814         strcpy(cur, "Infinity");
 815         return 1;
 816     }
 817
 818     /* Nope, nothing out of the ordinary, full steam ahead! */
 819     adj_exp = bn->digits + bn->expn -1;
 820     /* For values near to zero, we do not use exponential notation */
 821     if (bn->expn <= 0 && adj_exp >= -6) {
 822         if (bn->sign) {
 823             *cur++ = '-';
 824         }
 825         /*pad with zeros if appropriate, plonk a point where we want it */
 826         if (bn->digits + bn->expn <= 0) {
 827             int i;
 828             *cur++ = '0';
 829             if (bn->digits+bn->expn<0) *cur++ = '.';
 830             for (i=1; i <= -(bn->digits + bn->expn); i++) *cur++ = '0';
 831         }
 832         for (cur_dig = bn->digits-1; cur_dig >-1; cur_dig--) {
 833             if (1+cur_dig + bn->expn == 0) *cur++ ='.';
 834             *cur++ = '0' + BN_getd(bn, cur_dig);
 835         }
 836         *cur = 0;
 837     }
 838     else { /* Use exponential notation, different for sci and eng */
 839         if (bn->sign) *cur++ = '-'; /* We don't prefix '+' */
 840
 841         if (eng) {
 842             int deficit;
 843             if (adj_exp < 0) {
 844                 deficit = (-adj_exp) % 3;
 845                 if (deficit == 1) {
 846                     deficit = 2;
 847                     adj_exp -= 2;
 848                 }
 849                 else if (deficit == 2) {
 850                     deficit = 1;
 851                     adj_exp -=1;
 852                 }
 853             }
 854             else {
 855                 deficit = adj_exp % 3;
 856                 adj_exp -= deficit;
 857             }
 858             /* so, d = 0. x.yyyy, d=1 xx.y d=2 xxx.yyy special case if xxx*/
 859
 860             *cur++ = '0' + BN_getd(bn, bn->digits-1);
 861             if (deficit == 0) *cur++ = '.';
 862
 863             if (bn->digits == 1) {
 864                 *cur++ = '0';
 865                 if (deficit == 1) *cur++ = '.';
 866                 *cur++ = '0';
 867             }
 868             else if (bn->digits == 2) {
 869                 *cur++ = '0' + BN_getd(bn, bn->digits-2);
 870                 if (deficit == 1) *cur++ = '.';
 871                 *cur++ = '0';
 872             }
 873             else {
 874                 *cur++ = '0' + BN_getd(bn, bn->digits-2);
 875                 if (deficit == 1) *cur++ = '.';
 876                 *cur++ = '0' + BN_getd(bn, bn->digits-3);
 877                 if (bn->digits != 3 && deficit == 2) *cur++ = '.';
 878
 879                 for (cur_dig=bn->digits-4; cur_dig>-1; cur_dig--) {
 880                     *cur++ = '0' + BN_getd(bn, cur_dig);
 881                 }
 882             }
 883
 884             *cur++ = 'E';
 885             sprintf(cur, "%+i", adj_exp);
 886
 887         }
 888         else { /* Scientific */
 889             if (bn->digits == 1) { /* We don't want 1.E+7 */
 890                 *cur++ = '0'+ BN_getd(bn, 0);
 891             }
 892             else { /* We have x.xE */
 893                 *cur++ = '0' + BN_getd(bn, bn->digits-1);
 894                 *cur++ = '.';
 895                 for (cur_dig = bn->digits-2; cur_dig > -1; cur_dig--) {
 896                     *cur++ = '0' + BN_getd(bn, cur_dig);
 897                 }
 898             }
 899             *cur++ = 'E'; /* Eza Good, Eza Good */
 900             sprintf(cur, "%+i", adj_exp);
 901         }
 902     }
 903
 904     return 0;
 905 }
 906
 907 /*
 908
 909 =item C<BIGNUM*
 910 BN_from_string(PINTD_ char* s2, BN_CONTEXT *context)>
 911
 912 Convert a scientific string to a BIGNUM.  This function deals entirely
 913 with common-or-garden C byte strings, so the library can work
 914 anywhere.  Another version will be eventually required to cope with
 915 the parrot string fun.
 916
 917 This is the Highly Pedantic string conversion.  If C<context> has
 918 C<extended> as a true value, then the full range of extended number is
 919 made available, and any string which does not match the numeric syntax
 920 is converted to a quiet NaN.
 921
 922 Does not yet check for exponent overflow.
 923
 924 =cut
 925
 926 */
 927
 928 BIGNUM*
 929 BN_from_string(PINTD_ char* s2, BN_CONTEXT *context)
 930 {
 931     BIGNUM *result;
 932     BIGNUM *temp;
 933
 934     INTVAL pos = 0;             /* current digit in buffer */
 935     int negative = 0;        /* is it negative */
 936     int seen_dot = 0;        /* have we seen a '.' */
 937     int seen_e = 0;          /* have we seen an 'E' or 'e' */
 938     int exp_sign = 0;        /* is the exponent negative */
 939     int in_exp = 0;          /* are we reading exponent digits */
 940     int in_number = 0;       /* are we reading coeff digits */
 941     INTVAL exponent = 0;        /* the exponent */
 942     INTVAL fake_exponent = 0;   /* adjustment for digits after a '.' */
 943     INTVAL i = 0;
 944     int non_zero_digits = 0; /* have we seen *any* digits */
 945     int seen_plus = 0;       /* was number prefixed with '+' */
 946     int infinity =0;
 947     int qNAN = 0;
 948     int sNAN = 0;
 949
 950     temp = BN_new(PINT_ 1);     /* We store coeff reversed in temp */
 951
 952     while (*s2) { /* charge through the string */
 953         if (isdigit((unsigned char)*s2) && !in_exp) {
 954             /* We're somewhere in the main string of numbers */
 955             int digit = *s2 - '0'; /* byte me! */
 956             if (digit ==0 && !non_zero_digits) { /* ignore leading zeros */
 957                 in_number = 1;
 958                 s2++;
 959                 if (seen_dot) fake_exponent--;
 960                 continue;
 961             }
 962             else {
 963                 non_zero_digits = 1;
 964             }
 965             in_number = 1;
 966             BN_grow(PINT_ temp, pos+10);
 967             BN_setd(temp, pos, digit);
 968             pos++;
 969             if (seen_dot) {
 970                 fake_exponent--;
 971             }
 972         }
 973         else if (isdigit((unsigned char)*s2) && in_exp) {
 974             exponent = 10 * exponent + (*s2 - '0'); /*XXX: overflow check */
 975         }
 976         else if (!in_number) {
 977             /* we've not yet seen any digits */
 978             if (*s2 == '-') {
 979                 if (seen_plus || negative || seen_dot) {
 980                     if (!context->extended) {
 981                         BN_EXCEPT(PINT_ BN_CONVERSION_SYNTAX,
 982                                   "Incorrect number format");
 983                     }
 984                     else {
 985                         qNAN = 1; break;
 986                     }
 987                 }
 988                 negative = 1;
 989             }
 990             else if (*s2 == '.') {
 991                 seen_dot = 1;
 992             }
 993             else if (*s2 == '+') {
 994                 if (seen_plus || negative || seen_dot) {
 995                     if (!context->extended) {
 996                         BN_EXCEPT(PINT_ BN_CONVERSION_SYNTAX,
 997                                   "Incorrect number format");
 998                     }
 999                     else {
1000                         qNAN = 1; break;
1001                     }
1002                 }
1003                 seen_plus = 1; /* be very quiet */
1004             }
1005             else if (context->extended) {
1006                 if (*s2 == 'i' || *s2 == 'I') {
1007                     if (!strncasecmp("inf", s2, 4)) { /* We include the \0 */
1008                         infinity = 1;
1009                         /* For certain, restricted values of infinity */
1010                         break;
1011                     }
1012                     else if (!strncasecmp("infinity", s2, 9)) {
1013                         infinity = 1;
1014                         break;
1015                     }
1016                     else {
1017                         qNAN = 1;
1018                         break;
1019                     }
1020                 }
1021                 else if (*s2 == 'n' || *s2 == 'N') {
1022                     qNAN = 1; /* Don't need to check, as default.. */
1023                     break;
1024                 }
1025                 else if (*s2 == 's' || *s2 == 'S') {
1026                     if (!strncasecmp("snan", s2, 5)) {
1027                         sNAN = 1;
1028                         break;
1029                     }
1030                     else {
1031                         qNAN = 1;
1032                         break;
1033                     }
1034                 }
1035                 qNAN = 1;
1036                 break;
1037             } /* don't know, not in extended mode... */
1038             else {
1039                 BN_EXCEPT(PINT_ BN_CONVERSION_SYNTAX,
1040                           "Incorrect number format");
1041             }
1042         }
1043         else {
1044             /* we've seen some digits, are we done yet? */
1045             if (!seen_dot && *s2 == '.' && !in_exp) {
1046                 seen_dot = 1;
1047             }
1048             else if (!seen_e && (*s2 == 'e' || *s2 == 'E')) {
1049                 seen_e = 1;
1050                 in_exp = 1;
1051             }
1052             else if (seen_e && !exp_sign) {
1053                 if (*s2 == '+') {
1054                     exp_sign = 1;
1055                 }
1056                 else if (*s2 == '-') {
1057                     exp_sign = -1;
1058                 }
1059                 else {
1060                     if (!context->extended) {
1061                         BN_EXCEPT(PINT_ BN_CONVERSION_SYNTAX,
1062                                   "Incorrect number format");
1063                     }
1064                     else {
1065                         qNAN = 1; break;
1066                     }
1067                 }
1068             }
1069             else { /* We fall through here if we don't recognise something */
1070                 if (!context->extended) {
1071                     BN_EXCEPT(PINT_ BN_CONVERSION_SYNTAX,
1072                               "c Incorrect number format");
1073                 }
1074                 else {
1075                     qNAN = 1; break;
1076                 }
1077             }
1078         }
1079         s2++; /* rinse, lather... */
1080     }
1081
1082     if (!(qNAN || sNAN || infinity)) {
1083         if (in_number && !pos) { /* Only got zeros */
1084             pos = 1;
1085             BN_setd(temp, 0, 0);
1086         }
1087
1088         if (pos==0) { /* This includes ".e+20" */
1089             if (!context->extended) {
1090                 BN_EXCEPT(PINT_ BN_CONVERSION_SYNTAX, "no digits in string");
1091             }
1092             else {
1093                 qNAN = 1;
1094             }
1095         }
1096     }
1097
1098     result = BN_new(pos+1);
1099
1100     /* copy reversed string of digits backwards into result */
1101     if (!(qNAN || sNAN || infinity)) { /* Normal */
1102         temp->digits = pos;
1103
1104         for (i=0; i< temp->digits; i++) {
1105             BN_setd(result, i, BN_getd(temp, temp->digits-i-1));
1106         }
1107
1108         result->sign = negative;
1109         result->digits = pos;
1110         if (exp_sign == -1) {
1111             result->expn = fake_exponent - exponent;
1112         }
1113         else {
1114             result->expn = fake_exponent + exponent;
1115         }
1116     }
1117     else { /* Special */
1118         if (infinity) {
1119             BN_set_inf(PINT_ result);
1120             result->sign = negative;
1121         }
1122         else if (sNAN) {
1123             BN_set_sNAN(PINT_ result);
1124         }
1125         else {
1126             BN_set_qNAN(PINT_ result);
1127         }
1128     }
1129
1130
1131     BN_destroy(PINT_ temp);
1132     BN_really_zero(PINT_ result, context->extended);
1133     return result;
1134 }
1135
1136
1137 /*
1138
1139 =item C<int
1140 BN_strip_lead_zeros(PINTD_ BIGNUM* bn, BN_CONTEXT *context)>
1141
1142 Removes any zeros before the msd and after the lsd.
1143
1144 =cut
1145
1146 */
1147
1148 int
1149 BN_strip_lead_zeros(PINTD_ BIGNUM* bn, BN_CONTEXT *context)
1150 {
1151     INTVAL msd, i;
1152
1153     if (bn->digits == 0) return 0; /* Cannot "fail" with special nums */
1154
1155     msd = bn->digits-1;
1156
1157     while (0==BN_getd(bn, msd) && msd > 0) {
1158         msd--;
1159     }
1160
1161     bn->digits -= bn->digits-1 - msd;
1162 }
1163
1164 /*
1165
1166 =item C<int
1167 BN_strip_tail_zeros(PINTD_ BIGNUM *bn, BN_CONTEXT *context)>
1168
1169 Removes trailing zeros and increases the exponent appropriately.
1170 Does not remove zeros before the decimal point.
1171
1172 =cut
1173
1174 */
1175
1176 int
1177 BN_strip_tail_zeros(PINTD_ BIGNUM *bn, BN_CONTEXT *context)
1178 {
1179     INTVAL lsd, i;
1180
1181     lsd = 0;
1182
1183     while (0==BN_getd(bn, lsd)) {
1184         lsd++;
1185     }
1186     if (bn->expn >= 0) {
1187         lsd = 0; /* units column */
1188     }
1189     else if (bn->expn + lsd > 0) {
1190         lsd = -bn->expn;
1191     }
1192     for (i=0; i< bn->digits -lsd; i++) {
1193         BN_setd(bn, i, BN_getd(bn, i+lsd));
1194     }
1195
1196     if (CHECK_OVERFLOW(bn, lsd, context)) {
1197         BN_EXCEPT(PINT_ BN_OVERFLOW, "overflow when striping zeros");
1198     }
1199     bn->expn += lsd;
1200     bn->digits -= lsd;
1201 }
1202
1203 /*
1204
1205 =item C<int
1206 BN_make_integer(PINTD_ BIGNUM* bn, BN_CONTEXT* context)>
1207
1208 Convert the number to a plain integer I<if> precision such that this
1209 is possible.
1210
1211 =cut
1212
1213 */
1214
1215 int
1216 BN_make_integer(PINTD_ BIGNUM* bn, BN_CONTEXT* context)
1217 {
1218     /* Normal bignum */
1219     if (bn->expn > 0 && bn->digits + bn->expn <= context->precision) {
1220         INTVAL i;
1221         BN_grow(PINT_ bn, context->precision);
1222         for (i=bn->digits-1; i>= 0; i--) {
1223             BN_setd(bn, i+bn->expn, BN_getd(bn, i));
1224         }
1225         for (i=0; i<bn->expn; i++) {
1226             BN_setd(bn, i, 0);
1227         }
1228         bn->digits += bn->expn;
1229         bn->expn = 0;
1230     }
1231     else {
1232         return; /* XXX: fixed precision, bigints */
1233     }
1234 }
1235
1236 /*
1237
1238 =item C<int
1239 BN_really_zero(PINTD_ BIGNUM* bn, int allow_neg_zero)>
1240
1241 Sets any number which should be zero to a canonical zero.
1242
1243 To check if a number is equal to zero, use C<BN_is_zero()>.
1244
1245 =cut
1246
1247 */
1248
1249 int
1250 BN_really_zero(PINTD_ BIGNUM* bn, int allow_neg_zero)
1251 {
1252     INTVAL i;
1253     if (bn->digits == 0) return;
1254     for (i=0; i< bn->digits; i++)
1255         if (BN_getd(bn, i) != 0) return;
1256
1257     bn->digits = 1;
1258     bn->expn = 0;
1259     if (!allow_neg_zero) bn->sign = 0;
1260     return;
1261 }
1262
1263 /*
1264
1265 =item C<void
1266 BN_round(PINTD_ BIGNUM *bn, BN_CONTEXT* context)>
1267
1268 Rounds C<*bn> according to C<*context>.
1269
1270 =cut
1271
1272 */
1273
1274 void
1275 BN_round(PINTD_ BIGNUM *bn, BN_CONTEXT* context)
1276 {
1277     /* In exported version, must check for sNAN */
1278     if (bn->digits == 0 && am_sNAN(bn)) {
1279         BN_nonfatal(PINT_ context, BN_EXCEPTION_INVALID_OPERATION,
1280                     "sNaN in round");
1281         BN_set_sNAN(PINT_ bn);
1282         return;
1283     }
1284     else {
1285         BN_iround(PINT_ bn, context);
1286         return;
1287     }
1288 }
1289
1290 /*
1291
1292 =item C<int
1293 BN_iround(PINTD_ BIGNUM *bn, BN_CONTEXT* context)>
1294
1295 Rounds victim according to context.
1296
1297 Round assumes that any leading zeros are significant (after an
1298 addition operation, for instance).
1299
1300 If C<precision> is positive, the digit string is rounded to have no more
1301 than C<precision> digits.  If C<precision> is equal to zero, the number
1302 is treated as an integer, and any digits after the number's decimal
1303 point are removed.  If C<precision> is negative, the number is rounded
1304 so that there are no more than - C<precision> digits after the decimal
1305 point.
1306
1307 eg. for  1.234567E+3 with rounding of C<ROUND_DOWN>
1308
1309     precision:  4 =>  1.234E+3      1234
1310     precision:  6 =>  1.234567E+3   1234.56
1311     precision:  9 =>  1.234567E+3   1234.567
1312     precision:  0 =>  1234          1234
1313     precision: -1 =>  1.2345E+3     1234.5
1314     precision: -9 =>  1.234567E+3   1234.567
1315
1316 =cut
1317
1318 */
1319
1320 int
1321 BN_iround(PINTD_ BIGNUM *bn, BN_CONTEXT* context)
1322 {
1323     PARROT_ASSERT(bn!= NULL);
1324     PARROT_ASSERT(context != NULL);
1325
1326     if (bn->digits == 0) {
1327         return 0; /* rounding special values always works */
1328     }
1329
1330     if (context->precision < 1) { /* Rounding a BigInt or fixed */
1331         BN_round_as_integer(PINT_ bn, context);
1332         return;
1333     }
1334     /* Rounding a BigNum or sdaNum*/
1335     if (bn->digits > context->precision) {
1336         /* We're rounding, right... do we care? */
1337         BN_nonfatal(PINT_ context, BN_ROUNDED, "Argument rounded");
1338         { /* Check for lost digits */
1339             INTVAL digit;
1340             for (digit = 0; digit < bn->digits - context->precision; digit++) {
1341                 if (BN_getd(bn, digit) != 0) {
1342                     BN_nonfatal(PINT_ context, BN_LOST_DIGITS,
1343                                 "digits lost while rounding");
1344                     BN_nonfatal(PINT_ context, BN_INEXACT,
1345                                 "Loss of precision while rounding");
1346                     break;
1347                 }
1348             }
1349         }
1350
1351         if (context->rounding == ROUND_DOWN) {
1352             return BN_round_down(PINT_ bn, context);
1353         }
1354         else if (context->rounding == ROUND_HALF_UP) {
1355             if (BN_getd(bn, bn->digits - context->precision -1) > 4) {
1356                 return BN_round_up(PINT_ bn, context);
1357             }
1358             else {
1359                 return BN_round_down(PINT_ bn, context);
1360             }
1361         }
1362         else if (context->rounding == ROUND_HALF_EVEN) {
1363
1364             if (BN_getd(bn, bn->digits - context->precision -1) > 5) {
1365                 return BN_round_up(PINT_ bn, context);
1366             }
1367             else if (BN_getd(bn, bn->digits - context->precision -1) < 5) {
1368                 return BN_round_down(PINT_ bn, context);
1369             }
1370             else {
1371                 INTVAL i = bn->digits - context->precision -2;
1372                 if (i > -1) {
1373                     while (i>=0) {
1374                         if (BN_getd(bn, i) != 0) {
1375                             return BN_round_up(PINT_ bn, context);
1376                         }
1377                         i--;
1378                     }
1379                 }
1380                 switch (BN_getd(bn, bn->digits-context->precision)) {
1381                 case 0 :
1382                 case 2 :
1383                 case 4 :
1384                 case 6 :
1385                 case 8 :
1386                     return BN_round_down(PINT_ bn, context);
1387                 default:
1388                     return BN_round_up(PINT_ bn, context);
1389                 }
1390             }
1391
1392         }
1393         else if (context->rounding == ROUND_CEILING) {
1394             INTVAL i;
1395             if (bn->sign) {
1396                 return BN_round_down(PINT_ bn, context);
1397             }
1398             for (i = bn->digits - context->precision -1; i > -1; i--) {
1399                 if (BN_getd(bn, i) != 0) {
1400                     return BN_round_up(PINT_ bn, context);
1401                 }
1402             }
1403             return BN_round_down(PINT_ bn, context);
1404         }
1405         else if (context->rounding == ROUND_FLOOR) {
1406             INTVAL i;
1407             if (!bn->sign) {
1408                 return BN_round_down(PINT_ bn, context);
1409             }
1410             for (i = bn->digits - context->precision; i > -1; i--) {
1411                 if (BN_getd(bn, i) != 0) {
1412                     return BN_round_up(PINT_ bn, context);
1413                 }
1414             }
1415             return BN_round_down(PINT_ bn, context);
1416         }
1417         BN_EXCEPT(PINT_ BN_EXCEPTION_INVALID_OPERATION, "Unknown rounding attempted");
1418     }
1419     return;
1420 }
1421
1422 /*
1423
1424 =item C<int
1425 BN_round_up(PINTD_ BIGNUM *bn, BN_CONTEXT* context)>
1426
1427 Truncates coefficient of C<bn> to have C<precision> digits, then adds
1428 1 to the last digits and carries until done.  Do not call this
1429 function with non-positive values of C<precision>.
1430
1431 =cut
1432
1433 */
1434
1435 int
1436 BN_round_up(PINTD_ BIGNUM *bn, BN_CONTEXT* context)
1437 {
1438     INTVAL i, carry;
1439
1440     /* Do a cheap num += 1E+(num->expn) */
1441     carry = 1;
1442     for (i = bn->digits - context->precision; i< bn->digits; i++) {
1443         carry += BN_getd(bn, i);
1444         BN_setd(bn, i-bn->digits + context->precision, carry%10);
1445         carry = carry / 10;
1446     }
1447     if (carry) { /* We had 999999999 + 1, extend number */
1448         INTVAL extra = bn->digits - context->precision;
1449         BN_setd(bn, context->precision, carry);
1450         if (CHECK_OVERFLOW(bn, extra, context)) {
1451             BN_EXCEPT(PINT_ BN_OVERFLOW, "overflow while rounding");
1452         }
1453         bn->expn += extra;
1454         bn->digits = context->precision +1;
1455         return BN_iround(PINT_ bn, context);
1456     }
1457     else {
1458         INTVAL extra = bn->digits - context->precision;
1459         if (CHECK_OVERFLOW(bn, extra, context)) {
1460             BN_EXCEPT(PINT_ BN_OVERFLOW, "overflow while rounding");
1461         }
1462         bn->expn += extra;
1463         bn->digits = context->precision;
1464         return;
1465     }
1466 }
1467
1468 /*
1469
1470 =item C<int
1471 BN_round_down(PINT_ BIGNUM *bn, BN_CONTEXT* context)>
1472
1473 Truncates the coefficient of C<bn> to have C<precision> digits.  Do
1474 not call this function with non-positive precision.
1475
1476 =cut
1477
1478 */
1479
1480 int
1481 BN_round_down(PINT_ BIGNUM *bn, BN_CONTEXT* context)
1482 {
1483     INTVAL i =0;
1484     INTVAL extra;
1485
1486     for (i=0; i<context->precision; i++) {
1487         int temp = BN_getd(bn, i+bn->digits - context->precision);
1488         BN_setd(bn, i, temp);
1489     }
1490     extra = bn->digits - context->precision;
1491     if (CHECK_OVERFLOW(bn, extra, context)) {
1492         BN_EXCEPT(PINT_ BN_OVERFLOW, "overflow while rounding");
1493     }
1494     bn->expn += extra;
1495     bn->digits = context->precision;
1496
1497     return;
1498 }
1499
1500 /*
1501
1502 =item C<void
1503 BN_round_as_integer(PINTD_ BIGNUM *bn, BN_CONTEXT *context)>
1504
1505 C<precision> must be less than one.  This rounds so that C<expn> is at
1506 least C<precision>.  Name is slightly misleading.
1507
1508 =cut
1509
1510 */
1511
1512 void
1513 BN_round_as_integer(PINTD_ BIGNUM *bn, BN_CONTEXT *context)
1514 {
1515     INTVAL i;
1516     BN_CONTEXT temp_context;
1517
1518     if (bn->expn < context->precision) {
1519         {
1520             /* Are we losing information? */
1521             for (i=0;
1522                  i< (context->precision - bn->expn) && i<bn->digits;
1523                  i++) {
1524                 if (BN_getd(bn, i) != 0) {
1525                     BN_nonfatal(PINT_ context, BN_LOST_DIGITS,
1526                                 "digits lost while rounding");
1527                     break;
1528                 }
1529             }
1530         }
1531
1532         /* We'll cheat by passing a false context to the normal rounding.
1533            If "precision" < 1, we add a false zero to front and set p to 1 */
1534         temp_context = *context;
1535         temp_context.precision = bn->digits + bn->expn - context->precision;
1536         if (temp_context.precision < 1) {
1537             temp_context.precision = 1;
1538             BN_grow(bn, bn->digits + 1);
1539             BN_setd(bn, bn->digits, 0);
1540             bn->digits++;
1541             BN_iround(PINT_ bn, &temp_context);
1542         }
1543         else {
1544             BN_iround(PINT_ bn, &temp_context);
1545         }
1546         BN_really_zero(PINT_ bn, context->extended);
1547
1548         /* XXX: if using warning flags on context, | with temp context here */
1549     }
1550
1551     return;
1552 }
1553
1554 /*
1555
1556 =back
1557
1558 =head2  Arithmetic operations
1559
1560 Operations are performed like this:
1561
1562 =over 4
1563
1564 =item Rounding
1565
1566 Both operands are rounded to have no more than C<< context->precision >>
1567 digits.
1568
1569 =item Computation
1570
1571 The operation is computed.
1572
1573 =item Rounding of result
1574
1575 The result is then rounded to context->precision digits.
1576
1577 =item Conversion to zero and integerisation
1578
1579 If the result is equal to zero, it is made exactly zero.
1580
1581 Where the length of the coefficient + the exponent of the result is
1582 less than context->precision, the result is converted into an integer.
1583
1584 =back
1585
1586 The general form for all arithmetic operations is:
1587
1588     void BN_operation(result, one, two, context)
1589
1590 =cut
1591
1592 */
1593
1594
1595 /*
1596
1597 =over 4
1598
1599 =item C<int
1600 BN_arith_setup(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
1601                BN_CONTEXT *context, BN_SAVE_PREC* restore)>
1602
1603 Rounds one and two ready for arithmetic operation.
1604
1605 We assume that an operation might extend the digit buffer with zeros on
1606 either side, but not tamper with the actual digits of the number, we can
1607 then easily return the number to the correct (but still rounded)
1608 representation in _cleanup later
1609
1610 If you can promise that you will not modify the representation of one
1611 and two during your operation, then you may pass C<&restore> as a
1612 C<NULL> pointer to both setup and cleanup.
1613
1614 If overflow or underflow occurs during rounding, the numbers will be
1615 modified to the appropriate representation and will not be restorable.
1616
1617 =cut
1618
1619 */
1620
1621 int
1622 BN_arith_setup(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
1623                BN_CONTEXT *context, BN_SAVE_PREC* restore)
1624 {
1625     BN_strip_lead_zeros(PINT_ one, context);
1626     BN_strip_lead_zeros(PINT_ two, context);
1627     BN_iround(PINT_ one, context);
1628     BN_iround(PINT_ two, context);
1629     if (restore) {
1630         restore->one = *one;
1631         restore->two = *two;
1632     }
1633 }
1634
1635 /*
1636
1637 =item C<int
1638 BN_arith_cleanup(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
1639                  BN_CONTEXT *context, BN_SAVE_PREC* restore)>
1640
1641 Rounds C<result>, C<one>, C<two>, checks for zeroness and makes
1642 integers. Fixes C<one> and C<two> so they don't gain precision by
1643 mistake.
1644
1645 =cut
1646
1647 */
1648
1649 int
1650 BN_arith_cleanup(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
1651                  BN_CONTEXT *context, BN_SAVE_PREC* restore)
1652 {
1653     INTVAL i;
1654     unsigned char traps_save;
1655     unsigned char flags_save;
1656     if (restore && one->digits != restore->one.digits) {
1657         if (one->expn < restore->one.expn) {
1658             for (i=0; i<restore->one.digits; i++)
1659                 BN_setd(one, i, BN_getd(one, i+restore->one.expn - one->expn));
1660             one->expn = restore->one.expn;
1661         }
1662         one->digits = restore->one.digits;
1663     }
1664     if (restore && two->digits != restore->two.digits) {
1665         if (two->expn < restore->two.expn) {
1666             for (i=0; i<restore->two.digits; i++)
1667                 BN_setd(two, i, BN_getd(two, i+restore->two.expn - two->expn));
1668             two->expn = restore->two.expn;
1669         }
1670         two->digits = restore->two.digits;
1671     }
1672     /* We don't raise lost_digits after an operation, only beforehand,
1673        so we mask off the lost_digits handler to stop the error, and
1674        also clear any lost_digits flags*/
1675     traps_save = context->traps;
1676     flags_save = context->flags;
1677     context->traps &= ~(unsigned char)BN_F_LOST_DIGITS;
1678     BN_iround(PINT_ result, context);
1679     context->traps = traps_save;
1680     context->flags =  (context->flags & ~(unsigned char)BN_F_LOST_DIGITS)
1681                     | (flags_save & BN_F_LOST_DIGITS);
1682
1683     BN_strip_lead_zeros(PINT_ result, context);
1684     BN_really_zero(PINT_ result, context->extended);
1685     BN_make_integer(PINT_ result, context);
1686 }
1687
1688 /*
1689
1690 =item C<int
1691 BN_align(PINTD_ BIGNUM* one, BIGNUM* two)>
1692
1693 Adds zero digits so that decimal points of each number are at the same
1694 place.
1695
1696 =cut
1697
1698 */
1699
1700 int
1701 BN_align(PINTD_ BIGNUM* one, BIGNUM* two)
1702 {
1703     INTVAL i;
1704     INTVAL diff;
1705
1706     diff = one->expn - two->expn;
1707
1708     if (diff == 0) {
1709         /* The numbers have the same exponent, we merely need to extend
1710            the one with a shorter coeff length with zeros */
1711         if (one->digits < two->digits) {
1712             BIGNUM *temp = one;
1713             one = two;
1714             two = temp;
1715         }
1716
1717         BN_grow(PINT_ two, one->digits);
1718         for (i=two->digits; i<one->digits; i++) {
1719             BN_setd(two, i, 0);
1720         }
1721         two->digits = one->digits;
1722     }
1723     else {
1724         /* We need to pad both numbers to have the same number of digits
1725            the number with the most negative exponent only needs leading
1726            digits, while the number with the less negative expn may need
1727            both front and back padding, depending on its coeff length.
1728            Ideally we'll only move any digit once. */
1729         INTVAL final;
1730         /* force smallest exponent in two */
1731         if (diff < 0) {
1732             BIGNUM *temp = one;
1733             one = two;
1734             two = temp;
1735             diff = -diff;
1736         }
1737
1738         if (one->digits + diff < two->digits) {
1739             final = two->digits;
1740         }
1741         else {
1742             final = one->digits + diff;
1743         }
1744
1745         BN_grow(PINT_ one, final);
1746         BN_grow(PINT_ two, final);
1747         /* Add zeros to start of two */
1748         for (i=two->digits; i<final; i++)
1749             BN_setd(two, i, 0);
1750
1751         /* Add zeros to start of one */
1752         for (i=one->digits + diff; i< final; i++)
1753             BN_setd(one, i, 0);
1754
1755         /* Move one into new home */
1756         for (i=one->digits-1; i>-1; i--)
1757             BN_setd(one, i+diff, BN_getd(one, i));
1758
1759         /* Set end of one to zeros */
1760         for (i=0; i< diff; i++)
1761             BN_setd(one, i, 0);
1762
1763         one->digits = two->digits = final;
1764         one->expn = two->expn;
1765     }
1766
1767 }
1768 /*
1769
1770 =item C<void
1771 BN_add(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two, BN_CONTEXT *context)>
1772
1773 Adds one to two, returning value in result.
1774
1775 =cut
1776
1777 */
1778
1779 void
1780 BN_add(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two, BN_CONTEXT *context)
1781 {
1782     BN_SAVE_PREC restore;
1783     /* Special values */
1784     if (one->digits == 0 || two->digits == 0) {
1785         if (am_NAN(one) || am_NAN(two)) {
1786             if (am_sNAN(one) || am_sNAN(two)) {
1787                 BN_nonfatal(PINT_ context, BN_EXCEPTION_INVALID_OPERATION,
1788                             "sNAN in add");
1789             }
1790             BN_set_qNAN(PINT_ result);
1791             return;
1792         }
1793         /* Otherwise an infinity */
1794         if (am_INF(one) && am_INF(two)) {
1795             if (one->sign != two->sign) {
1796                 BN_nonfatal(PINT_ context, BN_EXCEPTION_INVALID_OPERATION,
1797                             "addition of +Inf and -Inf");
1798                 BN_set_qNAN(PINT_ result);
1799                 return;
1800             }
1801             else {
1802                 BN_set_inf(PINT_ result);
1803                 result->sign = one->sign;
1804                 return;
1805             }
1806         }
1807         /* So we've only got one infinity... */
1808         BN_set_inf(PINT_ result);
1809         result->sign = am_INF(one) ? one->sign : two->sign;
1810         return;
1811     }
1812
1813     /* Be careful to do 0 + -0 and -0 + 0 correctly */
1814     if (BN_is_zero(PINT_ one, context) && BN_is_zero(PINT_ two, context)) {
1815         result->digits = 1;
1816         result->expn = 0;
1817         BN_setd(result, 0, 0);
1818         if (one->sign & two->sign) {
1819             result->sign = 1;
1820         }
1821         else if (context->rounding == ROUND_FLOOR && (one->sign ^ two->sign)) {
1822             result->sign = 1;
1823         }
1824         else {
1825             result->sign = 0;
1826         }
1827         return;
1828     }
1829
1830
1831     BN_arith_setup(PINT_ result, one, two, context, &restore);
1832
1833     /* Do we mean add, or do we mean subtract? */
1834     if (one->sign && !two->sign) {             /* -a +  b =  (b-a) */
1835         BN_isubtract(PINT_ result, two, one, context);
1836     }
1837     else if (one->sign && two->sign) {         /* -a + -b = -(a+b) */
1838         BN_iadd(PINT_ result, one, two, context);
1839         result->sign = 1;
1840     }
1841     else if (two->sign) {                      /*  a + -b =  (a-b) */
1842         BN_isubtract(PINT_ result, one, two, context);
1843     }
1844     else {                                     /*  a +  b =  (a+b) */
1845         BN_iadd(PINT_ result, one, two, context);
1846     }
1847
1848     BN_arith_cleanup(PINT_ result, one, two, context, &restore);
1849     /* If using round_floor, need to make sure x + -x => -0 */
1850     if (context->rounding == ROUND_FLOOR && BN_is_zero(PINT_ result, context)
1851         && (one->sign ^ two->sign)) {
1852         result->sign = 1;
1853     }
1854 }
1855
1856 /*
1857
1858 =item C<int
1859 BN_iadd(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
1860               BN_CONTEXT *context)>
1861
1862 Adds together two aligned big numbers with coefficients of equal
1863 length.  Returns a result without reference to the signs of its
1864 arguments.  Cannot cope with special values.
1865
1866 =cut
1867
1868 */
1869
1870 int
1871 BN_iadd(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
1872               BN_CONTEXT *context)
1873 {
1874     INTVAL i;
1875     int carry, dig;
1876
1877     /* Make sure we don't do work we don't need, or add precision where
1878        it isn't wanted */
1879     if (BN_is_zero(PINT_ one, context)) {
1880         BN_copy(PINT_ result, two);
1881         result->sign = 0;
1882         return;
1883     }
1884     else if (BN_is_zero(PINT_ two, context)) {
1885         BN_copy(PINT_ result, one);
1886         result->sign = 0;
1887         return;
1888     }
1889
1890     /* Do the numbers overlap (within precision (and a bit) digits)?
1891        If not, we can simply use the first and round given the second
1892        by concatenating `01' to the result... Remember that we also know
1893        that neither is zero */
1894
1895     if (context->precision > -1 &&
1896         one->expn > two->expn + context->precision +1) {
1897         BN_grow(PINT_ result, context->precision + 1);
1898         for (i = 0; i < one->digits && i < context->precision; i++) {
1899             dig = BN_getd(one, one->digits - i -1);
1900             BN_setd(result, context->precision-i, dig);
1901         }
1902         for (i = i; i < context->precision; i++) {
1903             BN_setd(result, context->precision - i, 0);
1904         }
1905         BN_setd(result, 0, 1);
1906         result->digits = context->precision + 1;
1907         result->sign = 0;
1908         result->expn = one->expn + one->digits - context->precision - 1;
1909         return 0;
1910     } /* or two might be in the lead, but will he win by more than a length? */
1911     else if (context->precision > -1 &&
1912              two->expn > one->expn + context->precision + 1) {
1913         BN_grow(PINT_ result, context->precision + 1);
1914         for (i = 0; i < two->digits && i < context->precision; i++) {
1915             dig = BN_getd(two, two->digits -i-1);
1916             BN_setd(result, context->precision-i, dig);
1917         }
1918         for (i = i; i < context->precision; i++) {
1919             BN_setd(result, context->precision - i, 0);
1920         }
1921         BN_setd(result, 0, 1);
1922         result->digits = context->precision + 1;
1923         result->sign = 0;
1924         result->expn = two->expn + two->digits - context->precision - 1;
1925         return 0;
1926     }
1927
1928     /* Ok, we can't be lazy, we'll have to do it all ourselves */
1929     BN_align(PINT_ one, two);
1930
1931     BN_grow(PINT_ result, one->digits + 1);
1932
1933     carry = 0;
1934     for (i=0; i< one->digits; i++) {
1935         carry += BN_getd(one, i) + BN_getd(two, i);
1936         dig = carry % 10;
1937         BN_setd(result, i, dig);
1938         carry = carry / 10;
1939     }
1940     if (carry) {
1941         BN_setd(result, i, carry);
1942         result->digits = i+1;
1943     }
1944     else {
1945         result->digits = i;
1946     }
1947     result->sign = 0;
1948     result->expn = one->expn;
1949     return 0;
1950 }
1951
1952 /*
1953
1954 =item C<void
1955 BN_subtract(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
1956             BN_CONTEXT *context)>
1957
1958 Subtracts C<*two> from C<*one>, returning value in C<*result>.
1959
1960 =cut
1961
1962 */
1963
1964 void
1965 BN_subtract(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
1966             BN_CONTEXT *context)
1967 {
1968     BN_SAVE_PREC restore;
1969     /* Special values, like addition but careful with those signs eugene */
1970     if (one->digits == 0 || two->digits == 0) {
1971         if (am_NAN(one) || am_NAN(two)) {
1972             if (am_sNAN(one) || am_sNAN(two)) {
1973                 BN_nonfatal(PINT_ context, BN_EXCEPTION_INVALID_OPERATION,
1974                             "sNAN in subtract");
1975             }
1976             BN_set_qNAN(PINT_ result);
1977             return;
1978         }
1979         /* Otherwise an infinity */
1980         if (am_INF(one) && am_INF(two)) {
1981             if (one->sign == two->sign) {
1982                 BN_nonfatal(PINT_ context, BN_EXCEPTION_INVALID_OPERATION,
1983                             "subtraction of Inf and Inf");
1984                 BN_set_qNAN(PINT_ result);
1985                 return;
1986             }
1987             else {
1988                 BN_set_inf(PINT_ result);
1989                 result->sign = one->sign;
1990                 return;
1991             }
1992         }
1993         /* So we've only got one infinity... */
1994         BN_set_inf(PINT_ result);
1995         result->sign = am_INF(one) ? one->sign : (1 & (1 ^ two->sign));
1996         return;
1997     }
1998
1999     /* Be careful to do 0 + -0 and -0 + 0 correctly*/
2000     if (BN_is_zero(PINT_ one, context) && BN_is_zero(PINT_ two, context)) {
2001         result->digits = 1;
2002         result->expn = 0;
2003         BN_setd(result, 0, 0);
2004         if (one->sign && !two->sign && context->extended) {
2005             result->sign = 1;
2006         }
2007         else if (context->rounding == ROUND_FLOOR &&
2008                  (one->sign == two->sign) &&
2009                  context->extended) {
2010             result->sign = 1;
2011         }
2012         else {
2013             result->sign = 0;
2014         }
2015         return;
2016     }
2017
2018     BN_arith_setup(PINT_ result, one, two, context, &restore);
2019
2020     /* Do we mean subtract, or do we mean add? */
2021     if (one->sign && !two->sign) {         /* -a -  b = -(a+b) */
2022         BN_iadd(PINT_ result, one, two, context);
2023         result->sign = 1;
2024     }
2025     else if (one->sign && two->sign) {     /* -a - -b = (b-a)  */
2026         BN_isubtract(PINT_ result, two, one, context);
2027     }
2028     else if (two->sign) {                  /*  a - -b = (a+b)  */
2029         BN_iadd(PINT_ result, one, two, context);
2030     }
2031     else {                                 /*  a -  b = (a-b)  */
2032         BN_isubtract(PINT_ result, one, two, context);
2033     }
2034
2035     BN_arith_cleanup(PINT_ result, one, two, context, &restore);
2036     /* If using round_floor, need to make sure x + -x => -0 */
2037     if (context->rounding == ROUND_FLOOR && BN_is_zero(PINT_ result, context)
2038         && (one->sign == two->sign)) {
2039         result->sign = 1;
2040     }
2041 }
2042
2043 /*
2044
2045 =item C<int
2046 BN_isubtract(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
2047               BN_CONTEXT *context)>
2048
2049 Subtracts two from one, assumes both numbers have positive aligned
2050 coefficients of equal length.  Sets sign of result as appropriate.
2051 Cannot cope with special values.
2052
2053 =cut
2054
2055 */
2056
2057 int
2058 BN_isubtract(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
2059               BN_CONTEXT *context)
2060 {
2061     INTVAL i;
2062     int carry, dig;
2063     /* Make sure we don't do work we don't need, or add precision where
2064        it isn't wanted */
2065     if (BN_is_zero(PINT_ one, context)) {
2066         BN_copy(PINT_ result, two);
2067         result->sign = 1;
2068         return;
2069     }
2070     else if (BN_is_zero(PINT_ two, context)) {
2071         BN_copy(PINT_ result, one);
2072         result->sign = 0;
2073         return;
2074     }
2075
2076     /* Do the numbers fail to overlap?  If so, we can simplify the sum
2077        by taking a little bit, but essentially copying... We don't
2078        yet know which number is the bigger of the two, so need to do
2079        each by itself */
2080     if (context->precision > -1 &&
2081         one->expn > two->expn + context->precision + 1) {
2082         BN_grow(PINT_ result, context->precision + 1);
2083         carry = -1;
2084         for (i=0; i<one->digits; i++) {
2085             carry = carry + BN_getd(one, i);
2086             if (carry < 0) {
2087                 BN_setd(result, context->precision - one->digits + i + 1,
2088                         10+carry);
2089                 carry = -1;
2090             }
2091             else {
2092                 BN_setd(result, context->precision - one->digits + i + 1,
2093                         carry);
2094                 carry = 0;
2095             }
2096         }
2097         for (i=0; i<context->precision + 1 - one->digits; i++) {
2098             BN_setd(result, i, 9);
2099         }
2100         result->expn = one->expn + one->digits - context->precision -1;
2101         result->sign = 0;
2102         result->digits = context->precision + 1;
2103         return 1;
2104     } /* or, do we do [ickle] - b */
2105     else if (context->precision > -1 &&
2106              two->expn > one->expn + context->precision + 1) {
2107         BN_grow(PINT_ result, context->precision + 1);
2108         carry = -1;
2109         for (i=0; i<two->digits; i++) {
2110             carry = carry + BN_getd(two, i);
2111             if (carry < 0) {
2112                 BN_setd(result, context->precision - two->digits + i + 1,
2113                         10+carry);
2114                 carry = -1;
2115             }
2116             else {
2117                 BN_setd(result, context->precision - two->digits + i + 1,
2118                         carry);
2119                 carry = 0;
2120             }
2121         }
2122         for (i=0; i<context->precision + 1 - two->digits; i++) {
2123             BN_setd(result, i, 9);
2124         }
2125         result->expn = two->expn + two->digits - context->precision -1;
2126         result->sign = 1;
2127         result->digits = context->precision + 1;
2128         return 1;
2129     }
2130
2131
2132     BN_align(PINT_ one, two);
2133
2134     /* as a-b == -(b-a), we find larger of
2135        a and b and make sure it is in one */
2136     carry = 0;
2137     for (i=one->digits -1; i>-1; i--) {
2138         carry = BN_getd(one, i) - BN_getd(two, i);
2139         if (carry) break;
2140     }
2141
2142     if (!carry) { /* a==b*/
2143         result->digits = 1;
2144         result->sign = 0;
2145         BN_setd(result, 0, 0);
2146         return;
2147     }
2148     else if (carry < 0) { /* b > a */
2149         BN_isubtract(PINT_ result, two, one, context);
2150         result->sign = 1;
2151     }
2152     else {
2153         BN_grow(PINT_ result, one->digits + 1);
2154
2155         carry = 0;
2156         for (i=0; i<one->digits; i++) {
2157             carry = carry + BN_getd(one, i) - BN_getd(two, i);
2158             if (carry < 0) {
2159                 BN_setd(result, i, 10+carry);
2160                 carry = -1;
2161             }
2162             else {
2163                 BN_setd(result, i, carry);
2164                 carry = 0;
2165             }
2166         }
2167
2168         PARROT_ASSERT(!carry); /* as to get here a > b*/
2169
2170         result->digits = one->digits;
2171         result->expn = one->expn;
2172         result->sign = 0;
2173     }
2174
2175 }
2176
2177 /*
2178
2179 =item C<void
2180 BN_plus(PINTD_ BIGNUM* result, BIGNUM *one, BN_CONTEXT *context)>
2181
2182 Perform unary C<+> on C<*one>.  Does all the rounding and what have you.
2183
2184 =cut
2185
2186 */
2187
2188 void
2189 BN_plus(PINTD_ BIGNUM* result, BIGNUM *one, BN_CONTEXT *context)
2190 {
2191     /* Check for special values */
2192     if (one->digits ==0) {
2193         if (am_sNAN(one)) BN_nonfatal(PINT_ context, BN_EXCEPTION_INVALID_OPERATION,
2194                                       "sNAN in plus");
2195         if (am_NAN(one)) {
2196             BN_set_qNAN(PINT_ result);
2197             return;
2198         }
2199         else { /* Infinity */
2200             BN_set_inf(PINT_ result);
2201             result->sign = one->sign;
2202             return;
2203         }
2204     }
2205
2206     BN_arith_setup(PINT_ result, one, one, context, NULL);
2207     BN_copy(PINT_ result, one);
2208     BN_really_zero(PINT_ result, 0);
2209     BN_arith_cleanup(PINT_ result, one, one, context, NULL);
2210 }
2211
2212 /*
2213
2214 =item C<void
2215 BN_minus(PINTD_ BIGNUM* result, BIGNUM *one, BN_CONTEXT *context)>
2216
2217 Perform unary C<-> (minus) on C<*one>.  Does all the rounding and what
2218 have you.
2219
2220 =cut
2221
2222 */
2223
2224 void
2225 BN_minus(PINTD_ BIGNUM* result, BIGNUM *one, BN_CONTEXT *context)
2226 {
2227     /* Check for special values */
2228     if (one->digits ==0) {
2229         if (am_sNAN(one)) BN_nonfatal(PINT_ context, BN_EXCEPTION_INVALID_OPERATION,
2230                                       "sNAN in minus");
2231         if (am_NAN(one)) {
2232             BN_set_qNAN(PINT_ result);
2233             return;
2234         }
2235         else { /* Infinity */
2236             BN_set_inf(PINT_ result);
2237             result->sign = 1 & (1 ^ one->sign);
2238             return;
2239         }
2240     }
2241
2242     BN_arith_setup(PINT_ result, one, one, context, NULL);
2243     BN_copy(PINT_ result, one);
2244     result->sign = result->sign ? 0 : 1;
2245     BN_really_zero(PINT_ result, 0);
2246     BN_arith_cleanup(PINT_ result, one, one, context, NULL);
2247 }
2248
2249 /*
2250
2251 =item C<void
2252 BN_compare(PINT_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
2253             BN_CONTEXT *context)>
2254
2255 Numerically compares C<*one> and C<*two>, storing the result (as a
2256 BIGNUM) in C<*result>.
2257
2258     result = 1  => one > two
2259     result = -1 => two > one
2260     result = 0  => one == two
2261
2262 =cut
2263
2264 */
2265
2266 void
2267 BN_compare(PINT_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
2268             BN_CONTEXT *context)
2269 {
2270     INTVAL cmp;
2271
2272     /* Special values */
2273     if (one->digits == 0 || two->digits ==0) {
2274         /* NaN */
2275         if (am_NAN(one) || am_NAN(two)) {
2276             if (am_sNAN(one) || am_sNAN(two)) {
2277                 BN_nonfatal(PINT_ context, BN_EXCEPTION_INVALID_OPERATION,
2278                             "NaN in compare");
2279             }
2280             BN_set_qNAN(PINT_ result);
2281             return;
2282         }
2283         /* Otherwise at least one of the operands is an infinity */
2284         if (one->sign != two->sign) {
2285             cmp = one->sign ? -1 : 1;
2286         }
2287         else if (am_INF(one) && am_INF(two)) {
2288             cmp = 0;
2289         }
2290         else if (am_INF(one)) {
2291             cmp = one->sign ? -1 : 1;
2292         }
2293         else {
2294             cmp = one->sign ? 1 : -1;
2295         }
2296     }
2297     else {
2298         BN_arith_setup(PINT_ result, one, two, context, NULL);
2299         cmp = BN_comp(PINT_ one, two, context);
2300     }
2301     result->digits = 1;
2302     result->expn = 0;
2303
2304     if (cmp == 0) {
2305         BN_setd(result, 0, 0);
2306         result->sign = 0;
2307     }
2308     else if (cmp > 0) {
2309         BN_setd(result, 0, 1);
2310         result->sign = 0;
2311     }
2312     else {
2313         BN_setd(result, 0, 1);
2314         result->sign = 1;
2315     }
2316     BN_arith_cleanup(PINT_ result, one, two, context, NULL);
2317 }
2318
2319 /*
2320
2321 =item C<void
2322 BN_multiply(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
2323              BN_CONTEXT *context)>
2324
2325 Multiplies C<*one> and C<*two>, storing the result in C<*result>.
2326
2327 =cut
2328
2329 */
2330
2331 void
2332 BN_multiply(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
2333              BN_CONTEXT *context)
2334 {
2335     if (one->digits == 0 || two->digits == 0) {
2336         if (am_NAN(one) || am_NAN(two)) {
2337             if (am_sNAN(one) || am_sNAN(two)) {
2338                 BN_nonfatal(PINT_ context, BN_EXCEPTION_INVALID_OPERATION,
2339                             "sNAN in multiply");
2340             }
2341             BN_set_qNAN(PINT_ result);
2342             return;
2343         }
2344         /* We've got at least one infinity */
2345         /* 0 * Inf => NaN */
2346         if (BN_is_zero(PINT_ one, context) || BN_is_zero(PINT_ two, context)) {
2347             BN_nonfatal(PINT_ context, BN_EXCEPTION_INVALID_OPERATION,
2348                         "Attempt to multiply 0 and Infinity");
2349             BN_set_qNAN(PINT_ result);
2350             return;
2351         }
2352
2353         /* (anything but 0) * Inf => Inf */
2354         BN_set_inf(PINT_ result);
2355         result->sign = 1 & (one->sign ^ two->sign);
2356         return;
2357     }
2358
2359     BN_arith_setup(PINT_ result, one, two, context, NULL);
2360
2361     BN_imultiply(PINT_ result, one, two, context);
2362
2363     BN_strip_lead_zeros(PINT_ result, context);
2364     BN_arith_cleanup(PINT_ result, one, two, context, NULL);
2365
2366 }
2367
2368 /*
2369
2370 =item C<int
2371 BN_imultiply(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
2372               BN_CONTEXT *context)>
2373
2374 Multiplication without the rounding and other set up.
2375
2376 =cut
2377
2378 */
2379
2380 int
2381 BN_imultiply(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
2382               BN_CONTEXT *context)
2383 {
2384     INTVAL i, j;
2385     int carry, dig;
2386
2387     BN_grow(PINT_ result, one->digits + two->digits + 2);
2388     /* zero contents of result so that it can be used as intermediate */
2389     for (i=0; i<one->digits + two->digits +2; i++)
2390         BN_setd(result, i, 0);
2391
2392     /* make sure largest coeff is in one */
2393     if (one->digits < two->digits) {
2394         BIGNUM* temp = one;
2395         one = two;
2396         two = temp;
2397     }
2398
2399     /* multiply element by element */
2400     for (i=0; i<two->digits; i++) {
2401         dig = BN_getd(two, i);
2402         carry = 0;
2403         for (j=0; j<one->digits; j++) {
2404             carry += BN_getd(one, j) * dig + BN_getd(result, i+j);
2405             BN_setd(result, i+j, carry % 10);
2406             carry = carry / 10;
2407         }
2408         if (carry) {
2409             BN_setd(result, i+j, carry);
2410         }
2411     }
2412
2413     /* extend if there's still stuff to take care of */
2414     if (carry) {
2415         result->digits = one->digits + two->digits + 1;
2416     }
2417     else if (BN_getd(result, one->digits + two->digits - 1)) {
2418         result->digits = one->digits + two->digits;
2419     }
2420     else {
2421         result->digits = one->digits + two->digits - 1;
2422     }
2423
2424     i = one->expn + two->expn;
2425     /*XXX: use unsigned here to be sure? */
2426     if (i > context->elimit) {
2427         BN_EXCEPT(PINT_ BN_OVERFLOW, "overflow in multiplication");
2428     }
2429     if (i < -context->elimit) {
2430         BN_EXCEPT(PINT_ BN_UNDERFLOW, "underflow in multiplication");
2431     }
2432     result->expn = i;
2433
2434     result->sign = 1 & (one->sign ^ two->sign);
2435     return;
2436 }
2437
2438 /*
2439
2440 =item C<void
2441 BN_divide(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
2442           BN_CONTEXT *context)>
2443
2444 Divide two into one, storing up to C<precision> digits in result.
2445 Performs own rounding.  We also assume that this function B<will not
2446 be used> to produce a BigInt.  That is the job of C<divide_integer()>.
2447
2448 If you want to divide two integers to produce a float, you must do so
2449 with C<precision> greater than the number of significant digits in
2450 either operand.  If you want the result to be an integer or a numer
2451 with a fixed fractional part
2452
2453 =cut
2454
2455 */
2456
2457
2458 void
2459 BN_divide(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
2460           BN_CONTEXT *context)
2461 {
2462     BIGNUM* rem;
2463     /* Check for special values */
2464     if (one->digits == 0 || two->digits == 0) {
2465         if (am_NAN(one) || am_NAN(two)) {
2466             if (am_sNAN(one) || am_sNAN(two)) {
2467                 BN_nonfatal(PINT_ context, BN_EXCEPTION_INVALID_OPERATION,
2468                             "sNAN in divide");
2469             }
2470             BN_set_qNAN(PINT_ result);
2471             return;
2472         }
2473         if (am_INF(one) && am_INF(two)) {
2474             BN_nonfatal(PINT_ context, BN_EXCEPTION_INVALID_OPERATION,
2475                         "Inf / Inf in divide");
2476             BN_set_qNAN(PINT_ result);
2477             return;
2478         }
2479         if (am_INF(one)) {
2480             if (BN_is_zero(PINT_ two, context)) {
2481                 BN_nonfatal(PINT_ context, BN_DIVISION_BY_ZERO,
2482                             "Inf / 0 in divide");
2483             }
2484             BN_set_inf(PINT_ result);
2485             result->sign = 1 & (one->sign ^ two->sign);
2486             return;
2487         }
2488         /* so we're left with x / Inf */
2489         result->digits = 1;
2490         result->expn = 0;
2491         BN_setd(result, 0, 0);
2492         result->sign = 1 & (one->sign ^ two->sign);
2493         return;
2494     }
2495
2496     if (BN_is_zero(PINT_ two, context)) {
2497         if (BN_is_zero(PINT_ one, context)) {
2498             BN_nonfatal(PINT_ context, BN_DIVISION_UNDEFINED,
2499                         " 0 / 0 in divide");
2500             BN_set_qNAN(PINT_ result);
2501             return;
2502         }
2503
2504         BN_nonfatal(PINT_ context, BN_DIVISION_BY_ZERO,
2505                     "division by zero in divide");
2506         BN_set_inf(PINT_ result);
2507         result->sign = 1 & (one->sign ^ two->sign);
2508         return;
2509     }
2510
2511     /* We're left with the case that only arg 1 is zero */
2512     if (BN_is_zero(PINT_ one, context)) {
2513         result->expn = 0;
2514         result->digits = 1;
2515         BN_setd(result, 0, 0);
2516         result->sign = 1 & (one->sign ^ two->sign);
2517         return;
2518     }
2519
2520     rem = BN_new(PINT_ 1);
2521     BN_arith_setup(PINT_ result, one, two, context, NULL);
2522     BN_idivide(PINT_ result, one, two, context, BN_DIV_DIVIDE, rem);
2523
2524     /* Use rem to work out things like rounding here, we'll do our
2525        own clean up as it's all a little odd */
2526
2527     BN_strip_lead_zeros(PINT_ result, context);
2528
2529     /*XXX: write this rounding to cope with precision < 1 */
2530     if (context->rounding == ROUND_HALF_EVEN) {
2531         if (result->digits > context->precision) {
2532             BN_nonfatal(PINT_ context, BN_ROUNDED, "Rounded in divide");
2533             /* We collected precision + 1 digits... */
2534             BN_really_zero(PINT_ rem, context->extended);
2535             if (BN_getd(result, 0) > 5) {
2536                 BN_nonfatal(PINT_ context, BN_INEXACT,
2537                             "Loss of precision in divide");
2538                 BN_round_up(PINT_ result, context);
2539             }
2540             else if (BN_getd(result, 0) == 5) {
2541                 BN_nonfatal(PINT_ context, BN_INEXACT,
2542                             "Loss of precision in divide");
2543                 if (rem->digits == 1 && BN_getd(rem, 0)==0) {
2544                     switch (BN_getd(result, 1)) {
2545                     case 2:
2546                     case 4:
2547                     case 6:
2548                     case 8:
2549                     case 0:
2550                         BN_round_down(PINT_ result, context);
2551                         break;
2552                     default :
2553                         BN_round_up(PINT_ result, context);
2554                     }
2555                 }
2556                 else {
2557                     BN_nonfatal(PINT_ context, BN_INEXACT,
2558                                 "Loss of precision in divide");
2559                     BN_round_up(PINT_ result, context);
2560                 }
2561             }
2562             else {
2563                 if (BN_getd(result, 0) !=0) {
2564                     BN_nonfatal(PINT_ context, BN_INEXACT,
2565                                 "Loss of precision in divide");
2566                 }
2567                 else if (!BN_is_zero(PINT_ result, context)) {
2568                     BN_nonfatal(PINT_ context, BN_INEXACT,
2569                                 "Loss of precision in divide");
2570                 }
2571                 BN_round_down(PINT_ result, context);
2572             }
2573         }
2574     }
2575     else if (context->rounding == ROUND_CEILING) {
2576         if (result->digits > context->precision) {
2577             BN_nonfatal(PINT_ context, BN_ROUNDED,
2578                         "Rounded in divide");
2579             BN_really_zero(PINT_ rem, context->extended);
2580
2581             if (result->sign) {
2582                 if (BN_getd(result, 0) != 0 ||
2583                     !BN_is_zero(PINT_ result, context)) {
2584                     BN_nonfatal(PINT_ context, BN_INEXACT,
2585                                 "Loss of precision in divide");
2586                 }
2587                 BN_round_down(PINT_ result, context);
2588             }
2589             else if (BN_getd(result, 0) != 0) {
2590                 BN_nonfatal(PINT_ context, BN_INEXACT,
2591                             "Loss of precision in divide");
2592                 BN_round_up(PINT_ result, context);
2593             }
2594             else if (!BN_is_zero(PINT_ rem, context)) {
2595                 BN_nonfatal(PINT_ context, BN_INEXACT,
2596                             "Loss of precision in divide");
2597                 BN_round_up(PINT_ result, context);
2598             }
2599             else {
2600                 BN_round_down(PINT_ result, context);
2601             }
2602         }
2603     }
2604     else if (context->rounding == ROUND_FLOOR) {
2605         if (result->digits > context->precision) {
2606             BN_nonfatal(PINT_ context, BN_ROUNDED,
2607                         "Rounded in divide");
2608             BN_really_zero(PINT_ rem, context->extended);
2609
2610             if (!result->sign) {
2611                 if (BN_getd(result, 0) != 0 ||
2612                     !BN_is_zero(PINT_ result, context)) {
2613                     BN_nonfatal(PINT_ context, BN_INEXACT,
2614                                 "Loss of precision in divide");
2615                 }
2616                 BN_round_down(PINT_ result, context);
2617             }
2618             else if (BN_getd(result, 0) != 0) {
2619                 BN_nonfatal(PINT_ context, BN_INEXACT,
2620                             "Loss of precision in divide");
2621                 BN_round_up(PINT_ result, context);
2622             }
2623             else if (!BN_is_zero(PINT_ rem, context)) {
2624                 BN_nonfatal(PINT_ context, BN_INEXACT,
2625                             "Loss of precision in divide");
2626                 BN_round_up(PINT_ result, context);
2627             }
2628             else {
2629                 BN_round_down(PINT_ result, context);
2630             }
2631         }
2632     }
2633     else { /* Other roundings just need digits to play with */
2634         unsigned char save_lost;
2635         unsigned char flags_save;
2636         /* We don't warn on lost digits here, as is after an operation */
2637         save_lost = context->traps;
2638         context->traps &= ~(unsigned char)BN_F_LOST_DIGITS;
2639         flags_save = context->flags;
2640         BN_iround(PINT_ result, context);
2641
2642         /* We need to check the remainder here, as we might have
2643            passed "[digits we want]0[digits we've kept a secret]" into
2644            the rounding without knowing it*/
2645         if (!BN_is_zero(PINT_ rem, context)) {
2646             BN_nonfatal(PINT_ context, BN_INEXACT,
2647                         "Loss of precision in divide");
2648         }
2649
2650         context->traps = save_lost;
2651         context->flags =  (context->flags & ~(unsigned char)BN_F_LOST_DIGITS)
2652             | (flags_save & BN_F_LOST_DIGITS);
2653     }
2654
2655     BN_really_zero(PINT_ result, context->extended);
2656
2657     BN_strip_tail_zeros(PINT_ result, context);
2658
2659     BN_make_integer(PINT_ result, context);
2660
2661     /* Remove trailing zeros if positive exponent */
2662     if (result->expn > 0) {
2663         INTVAL i, j;
2664         for (i=0; i<result->digits; i++) {
2665             if (BN_getd(result, i) != 0) break;
2666         }
2667         if (i) {
2668             for (j=i; j<result->digits; j++) {
2669                 BN_setd(result, j-i, BN_getd(result, j));
2670             }
2671         }
2672         if (CHECK_OVERFLOW(result, i, context)) {
2673             BN_EXCEPT(PINT_ BN_OVERFLOW, "overflow in divide");
2674         }
2675         result->expn += i;
2676     }
2677
2678     BN_destroy(PINT_ rem);
2679 }
2680
2681 /*
2682
2683 =item C<void
2684 BN_divide_integer(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
2685                    BN_CONTEXT *context)>
2686
2687 Places the integer part of C<*one> divided by C<*two> into C<*result>.
2688
2689 =cut
2690
2691 */
2692
2693 void
2694 BN_divide_integer(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
2695                    BN_CONTEXT *context)
2696 {
2697     BIGNUM* rem;
2698     /* Check for special values (same as divide...) */
2699     if (one->digits == 0 || two->digits == 0) {
2700         if (am_NAN(one) || am_NAN(two)) {
2701             if (am_sNAN(one) || am_sNAN(two)) {
2702                 BN_nonfatal(PINT_ context, BN_EXCEPTION_INVALID_OPERATION,
2703                             "sNAN in divide-integer");
2704             }
2705             BN_set_qNAN(PINT_ result);
2706             return;
2707         }
2708         if (am_INF(one) && am_INF(two)) {
2709             BN_nonfatal(PINT_ context, BN_EXCEPTION_INVALID_OPERATION,
2710                         "Inf / Inf in divide-integer");
2711             BN_set_qNAN(PINT_ result);
2712             return;
2713         }
2714         if (am_INF(one)) {
2715             if (BN_is_zero(PINT_ two, context)) {
2716                 BN_nonfatal(PINT_ context, BN_DIVISION_BY_ZERO,
2717                             "Inf / 0 in divide-integer");
2718             }
2719             BN_set_inf(PINT_ result);
2720             result->sign = 1 & (one->sign ^ two->sign);
2721             return;
2722         }
2723         /* so we're left with x / Inf */
2724         result->digits = 1;
2725         result->expn = 0;
2726         BN_setd(result, 0, 0);
2727         result->sign = 1 & (one->sign ^ two->sign);
2728         return;
2729     }
2730
2731     if (BN_is_zero(PINT_ two, context)) {
2732         if (BN_is_zero(PINT_ one, context)) {
2733             BN_nonfatal(PINT_ context, BN_DIVISION_UNDEFINED,
2734                         " 0 / 0 in divide-integer");
2735             BN_set_qNAN(PINT_ result);
2736             return;
2737         }
2738
2739         BN_nonfatal(PINT_ context, BN_DIVISION_BY_ZERO,
2740                     "division by zero in divide-integer");
2741         BN_set_inf(PINT_ result);
2742         result->sign = 1 & (one->sign ^ two->sign);
2743         return;
2744     }
2745
2746     /* We're left with the case that only arg 1 is zero */
2747     if (BN_is_zero(PINT_ one, context)) {
2748         result->expn = 0;
2749         result->digits = 1;
2750         BN_setd(result, 0, 0);
2751         result->sign = 1 & (one->sign ^ two->sign);
2752         return;
2753     }
2754
2755     rem = BN_new(PINT_ 1);
2756     BN_arith_setup(PINT_ result, one, two, context, NULL);
2757     BN_idivide(PINT_ result, one, two, context, BN_DIV_DIVINT, rem);
2758
2759     BN_really_zero(PINT_ rem, context->extended);
2760     if (result->expn >0 && context->precision > 0 &&
2761         result->expn + result->digits > context->precision &&
2762         !(rem->digits == 0 && BN_getd(rem, 0) == 0)) {
2763         BN_nonfatal(PINT_ context, BN_DIVISION_IMPOSSIBLE,
2764                   "divide-integer requires more precision than available");
2765         BN_set_qNAN(PINT_ result);
2766         BN_destroy(PINT_ rem);
2767         return;
2768     }
2769     BN_destroy(PINT_ rem);
2770     if (result->expn != 0) {
2771         INTVAL i;
2772         BN_grow(PINT_ result, result->expn + result->digits);
2773         for (i=0; i<result->digits; i++) {
2774             BN_setd(result, result->expn + result->digits -1 -i,
2775                     BN_getd(result, result->digits - 1- i));
2776         }
2777         for (i=0; i<result->expn; i++) {
2778             BN_setd(result, i, 0);
2779         }
2780         result->digits += result->expn;
2781         result->expn = 0;
2782     }
2783
2784     BN_really_zero(PINT_ result, context->extended);
2785     BN_make_integer(PINT_ result, context);
2786 }
2787
2788 /*
2789
2790 =item C<void
2791 BN_remainder(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
2792               BN_CONTEXT *context)>
2793
2794 Places the remainder from divide-integer (above) into C<*result>.
2795
2796 =cut
2797
2798 */
2799
2800 void
2801 BN_remainder(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
2802               BN_CONTEXT *context)
2803 {
2804     BIGNUM* fake;
2805
2806     /* Check for special values */
2807     if (one->digits == 0 || two->digits == 0) {
2808         if (am_NAN(one) || am_NAN(two)) {
2809             if (am_sNAN(one) || am_sNAN(two)) {
2810                 BN_nonfatal(PINT_ context, BN_EXCEPTION_INVALID_OPERATION,
2811                             "sNAN in remainder");
2812             }
2813             BN_set_qNAN(PINT_ result);
2814             return;
2815         }
2816         /* Infinities, first cover Inf rem x and Inf rem Inf */
2817         if (am_INF(one)) {
2818             BN_nonfatal(PINT_ context, BN_EXCEPTION_INVALID_OPERATION,
2819                         "x rem Inf in remainder");
2820             BN_set_qNAN(PINT_ result);
2821             return;
2822         }
2823         /* now cover x rem Inf => 0 */
2824         result->expn = 0;
2825         result->digits = 1;
2826         BN_setd(result, 0, 0);
2827         result->sign = one->sign;
2828         return;
2829     }
2830
2831     if (BN_is_zero(PINT_ two, context)) {
2832         if (BN_is_zero(PINT_ one, context)) {
2833             BN_nonfatal(PINT_ context, BN_DIVISION_UNDEFINED,
2834                         "0 rem 0 in remainder");
2835             BN_set_qNAN(PINT_ result);
2836             return;
2837         }
2838         BN_nonfatal(PINT_ context, BN_EXCEPTION_INVALID_OPERATION,
2839                     "x rem 0 in remainder");
2840         BN_set_qNAN(PINT_ result);
2841         return;
2842     }
2843
2844     if (BN_is_zero(PINT_ one, context)) {
2845         result->digits = 1;
2846         result->sign = 0;
2847         BN_setd(result, 0, 0);
2848         return;
2849     }
2850
2851     BN_arith_setup(PINT_ result, one, two, context, NULL);
2852     fake = BN_new(1);
2853     BN_idivide(PINT_ fake, one, two, context, BN_DIV_REMAIN, result);
2854
2855     BN_really_zero(PINT_ result, context->extended);
2856     if (fake->expn >0 && context->precision > 0 &&
2857         fake->expn + result->digits > context->precision &&
2858         !(result->digits == 0 && BN_getd(result, 0) == 0)) {
2859         BN_nonfatal(PINT_ context, BN_DIVISION_IMPOSSIBLE,
2860                     "remainder requires more precision than available");
2861         BN_set_qNAN(PINT_ result);
2862         BN_destroy(PINT_ fake);
2863         return;
2864     }
2865
2866
2867     BN_destroy(PINT_ fake);
2868
2869     result->sign = one->sign;
2870
2871     BN_arith_cleanup(PINT_ result, one, two, context, NULL);
2872 }
2873
2874 /*
2875
2876 =item C<int BN_idivide(PINT_ BIGNUM *result, BIGNUM *one, BIGNUM *two,
2877             BN_CONTEXT *context, BN_DIV_ENUM operation, BIGNUM *rem)>
2878
2879 Does the heavy work for the various division wossnames.
2880
2881 =cut
2882
2883 */
2884 int
2885 BN_idivide(PINT_ BIGNUM *result, BIGNUM *one, BIGNUM *two,
2886             BN_CONTEXT *context, BN_DIV_ENUM operation, BIGNUM *rem)
2887 {
2888     INTVAL i, j, divided, newexpn;
2889     BIGNUM *div, *t1, *t2;
2890     int s2, value;
2891
2892     /* We assume we've been given something to divide */
2893
2894     /* Make some temporaries, set all signs to positive for simplicity */
2895     /* We use result as a temporary and store the reversed result in t2 */
2896     div = BN_new(PINT_ 1);
2897     BN_copy(PINT_ div, one);
2898     BN_copy(PINT_ rem, div); /* In case doing int-div and don't div*/
2899     div->sign = 0;
2900     t1 = BN_new(PINT_ 1);
2901     t2 = BN_new(PINT_ 1);
2902     t2->digits = 0; /* ok, as all internal */
2903     s2 = two->sign; /* store the sign of 2, as we set it +ve internally */
2904     two->sign = 0;
2905     result->digits = 1;
2906     rem->digits = 1;
2907
2908     /* First position to try to fill */
2909     newexpn = one->digits + one->expn - two->digits - two->expn;
2910     if (newexpn > context->elimit) {
2911         BN_EXCEPT(PINT_ BN_OVERFLOW, "overflow in divide (1)");
2912     }
2913     if (newexpn < -context->elimit) {
2914         BN_EXCEPT(PINT_ BN_UNDERFLOW, "underflow in divide (1)");
2915     }
2916     t1->expn = newexpn;
2917
2918     value = 0;
2919     for (;;) {
2920         if (!(t2->digits % 10)) BN_grow(PINT_ t2, t2->digits+11);
2921         if ((operation == BN_DIV_DIVINT || operation == BN_DIV_REMAIN) &&
2922             t1->expn < 0) break;
2923         divided = 0;
2924         for (j=1; j<=10;j++) {
2925             int cmp;
2926             BN_setd(t1, 0, j);
2927             BN_imultiply(PINT_ result, t1, two, context);
2928             cmp = BN_comp(PINT_ result, div, context);
2929             if (cmp ==0) {
2930                 BN_setd(t2, value, j);
2931                 t2->digits++;
2932                 value++;
2933                 j = j+1; /* for multiply below */
2934                 divided = 1;
2935                 break;
2936             }
2937             else if (cmp> 0) {
2938                 if (j==1 && value == 0) break; /* don't collect leading 0s */
2939                 BN_setd(t2, value, j-1);
2940                 t2->digits++;
2941                 value++;
2942                 divided = 1;
2943                 break;
2944             }
2945         }
2946         if (divided) {
2947             BN_setd(t1, 0, j-1);
2948             BN_imultiply(PINT_ result, t1, two, context);
2949             BN_isubtract(PINT_ rem, div, result, context);
2950         }
2951
2952         /* Are we done yet? */
2953         if (value && rem->digits ==1 && BN_getd(rem, 0)==0) {
2954             break;
2955         }
2956
2957         /* We collect one more digit than precision requires, then
2958            round in divide, if we're doing divint or rem then we terminate
2959            at the decimal point and return */
2960         if (context->precision > 0) {
2961             if (t2->digits == context->precision + 1) {
2962                 break;
2963             }
2964         }
2965         else {
2966             if (t1->expn == context->precision -1) break;
2967         }
2968         if (operation == BN_DIV_DIVINT|| operation == BN_DIV_REMAIN) {
2969             if (t1->expn ==0) break;
2970         }
2971         if (CHECK_UNDERFLOW(t1, 1, context)) {
2972             BN_EXCEPT(PINT_ BN_UNDERFLOW, "underflow in divide (2)");
2973         }
2974         t1->expn--;
2975         if (divided) BN_copy(PINT_ div, rem);
2976     }
2977
2978     /* Work out the sign and exponent of the result */
2979     for (i=0; i< t2->digits; i++) {
2980         BN_setd(result, i, BN_getd(t2, t2->digits - 1 -i));
2981     }
2982     if (t2->digits == 0||(!divided&&!value)) {
2983         result->digits = 1;
2984         BN_setd(result, 0, 0);
2985         result->sign = 0;
2986     }
2987     else {
2988         result->digits = t2->digits;
2989         result->sign = 1&(one->sign ^ s2);
2990         result->expn = t1->expn; /* We know this is fine, from above */
2991     }
2992     two->sign = s2;
2993     rem->sign = 1&(one->sign ^ s2);
2994
2995     BN_destroy(PINT_ t1);
2996     BN_destroy(PINT_ t2);
2997     BN_destroy(PINT_ div);
2998
2999     return; /* phew! */
3000 }
3001
3002 /*
3003
3004 =item C<INTVAL
3005 BN_comp(PINTD_ BIGNUM *one, BIGNUM *two, BN_CONTEXT* context)>
3006
3007 Comparison with no rounding etc.
3008
3009 =cut
3010
3011 */
3012
3013 INTVAL
3014 BN_comp(PINTD_ BIGNUM *one, BIGNUM *two, BN_CONTEXT* context)
3015 {
3016     INTVAL i, j;
3017     int cmp;
3018
3019     BN_strip_lead_zeros(PINT_ one, context);
3020     BN_strip_lead_zeros(PINT_ two, context);
3021
3022     if (one->sign != two->sign) {
3023         if (BN_is_zero(PINT_ one, context) && BN_is_zero(PINT_ two, context)) {
3024             return 0; /* as -0 == 0 */
3025         }
3026         return one->sign ? -1 : 1;
3027     }
3028     else if (one->expn + one->digits > two->expn + two->digits) {
3029         return one->sign ? -1 : 1;
3030     }
3031     else if (one->expn + one->digits < two->expn + two->digits) {
3032         return one->sign ? 1 : -1;
3033     }
3034     else { /* Same sign, same "size" */
3035         for (i=0; i<one->digits && i<two->digits; i++) {
3036             cmp = BN_getd(one, one->digits-1-i)
3037                       - BN_getd(two, two->digits-1-i);
3038             if (cmp) return one->sign ? -cmp : cmp;
3039         }
3040         if (!cmp) {
3041             if (i==one->digits) {
3042                 for (i=i; i<two->digits; i++) {
3043                     cmp = 0-BN_getd(two, two->digits-1-i);
3044                     if (cmp) return one->sign ? -cmp : cmp;
3045                 }
3046             }
3047             else if (i==two->digits) {
3048                 for (i=i; i<one->digits; i++) {
3049                     cmp = BN_getd(one, one->digits-1-i);
3050                     if (cmp) return one->sign ? -cmp : cmp;
3051                 }
3052             }
3053             return one->sign ? -cmp : cmp;
3054         }
3055     }
3056 }
3057
3058 /*
3059
3060 =item C<void
3061 BN_power(PINTD_ BIGNUM* result, BIGNUM* bignum,
3062               BIGNUM* expn, BN_CONTEXT* context)>
3063
3064 Calculate C<result> = C<bignum> to the power of C<*expn>;
3065
3066 =cut
3067
3068 */
3069
3070 void
3071 BN_power(PINTD_ BIGNUM* result, BIGNUM* bignum,
3072               BIGNUM* expn, BN_CONTEXT* context)
3073 {
3074     BN_arith_setup(PINT_ result, bignum, expn, context, NULL);
3075     BN_arith_cleanup(PINT_ result, bignum, expn, context, NULL);
3076 }
3077
3078 /*
3079
3080 =item C<void
3081 BN_rescale(PINTD_ BIGNUM* result, BIGNUM* one, BIGNUM* two,
3082                 BN_CONTEXT* context)>
3083
3084 Rescales C<*one> to have an exponent of C<*two>.
3085
3086 =cut
3087
3088 */
3089
3090 void
3091 BN_rescale(PINTD_ BIGNUM* result, BIGNUM* one, BIGNUM* two,
3092                 BN_CONTEXT* context)
3093 {
3094     INTVAL expn;
3095     unsigned char lost = context->traps;
3096     context->traps &= ~(unsigned char)BN_F_LOST_DIGITS;
3097
3098     BN_arith_setup(PINT_ result, one, two, context, NULL);
3099
3100     expn = BN_to_int(PINT_ two, context);
3101
3102     context->traps = lost;
3103
3104     BN_arith_cleanup(PINT_ result, one, two, context, NULL);
3105 }
3106
3107 /*
3108
3109 =item C<INTVAL
3110 BN_to_int(PINT_ BIGNUM* bn, BN_CONTEXT* context)>
3111
3112 Converts the bignum into an integer, raises overflow if an exact
3113 representation cannot be created.
3114
3115 =cut
3116
3117 */
3118
3119 INTVAL
3120 BN_to_int(PINT_ BIGNUM* bn, BN_CONTEXT* context)
3121 {
3122     INTVAL insig, i;
3123     INTVAL result = 0;
3124     INTVAL maxdigs = BN_D_PER_INT < context->precision ?
3125         BN_D_PER_INT : context->precision;
3126     if (context->precision < 0) maxdigs = BN_D_PER_INT;
3127
3128
3129     BN_strip_lead_zeros(PINT_ bn, context);
3130     /* Check for definite big as your head overflow */
3131     if (bn->expn >= 0 && bn->expn + bn->digits > BN_D_PER_INT) {
3132         BN_EXCEPT(PINT_ BN_OVERFLOW, "bignum too large to fit in an int");
3133     }
3134     if (bn->expn < 0 && bn->expn + bn->digits > BN_D_PER_INT) {
3135
3136         BN_EXCEPT(PINT_ BN_OVERFLOW, "bignum too large to fit in an int");
3137     }
3138
3139     /* if e>0, if we'll lose precision we'll also be too big, so lose
3140        above anyway.  On the other hand, with e<0, we can lose digits <
3141        . from this so need to check that we don't lose precision */
3142     if (bn->expn<0 && context->traps & BN_F_LOST_DIGITS) {
3143         BN_EXCEPT(PINT_ BN_LOST_DIGITS, "digits lost in conv -> int");
3144     }
3145
3146     if (bn->digits + bn->expn > context->precision && context->precision > 0) {
3147         BN_EXCEPT(PINT_ BN_LOST_DIGITS, "digits lost in conv -> int");
3148     }
3149
3150     /* luckily, we get to keep our digits, so let's get at 'em */
3151     if (bn->expn >= 0) {
3152         for (i = bn->digits-1; i>-1; i--) {
3153             result = result * 10 + BN_getd(bn, i);
3154         }
3155         for (i=0; i<bn->expn; i++) result = result * 10;
3156     }
3157     else {
3158         for (i=bn->digits-1; i>-1-bn->expn; i--) {
3159             result = result * 10 + BN_getd(bn, i);
3160         }
3161     }
3162
3163     return bn->sign ? -result : result;
3164 }
3165
3166 /*
3167
3168 =item C<INTVAL
3169 BN_is_zero(BIGNUM* foo, BN_CONTEXT* context)>
3170
3171 Returns a boolean value indicating whether C<*foo> is zero.
3172
3173 =cut
3174
3175 */
3176
3177 INTVAL
3178 BN_is_zero(BIGNUM* foo, BN_CONTEXT* context)
3179 {
3180     BN_really_zero(foo, context->extended);
3181     if (foo->digits == 1 && foo->expn == 0 && BN_getd(foo, 0) == 0) {
3182         return 1;
3183     }
3184     else {
3185         return 0;
3186     }
3187 }
3188
3189 /*
3190
3191 =back
3192
3193 =head1 TODO
3194
3195 This is currently not used yet. Parrot has no BigNum support yet.
3196
3197 Parrot string playing, exception raising
3198
3199 ==head1 SEE ALSO
3200
3201 F<docs/docs/pdds/draft/pdd14_bignum.pod>,
3202 L<https://rt.perl.org/rt3/Ticket/Display.html?id=36330>
3203
3204 =cut
3205
3206 */
3207
3208 /*
3209  * Local variables:
3210  *   c-file-style: "parrot"
3211  * End:
3212  * vim: expandtab shiftwidth=4:
3213  */