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