2 Copyright (C) 2001-2006, The Perl Foundation.
7 src/bignum.c - A decimal arithmetic library Parrot
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.
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
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?">
26 I<On this occasion the kitten only purred: and it was impossible to
27 guess whether it meant `yes' or `no'.>
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.
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
41 Access digits, macros assume length given.
51 #include <string.h> /* XXX:ajg fixme later*/
53 /* * This lot wants to go in a (bignum specific seperate header file * */
57 =item C<BN_setd(BIGNUM*, pos, value)>
59 Set digit at C<pos> (zero is lsd) to C<value>.
61 =item C<int BN_getd(BIGNUM*, pos)>
63 Get value of digit at C<pos>.
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)
78 =item C<CHECK_OVERFLOW(bn, incr, context)>
80 If increasing the exponent of C<bn> by C<incr> will cause overflow (as
81 decided by C<elimit>), returns true.
87 #define CHECK_OVERFLOW(bn, incr, context) \
88 ((context)->elimit - ((bn)->expn + (bn)->digits -1) < (incr) ? 1 : 0)
92 =item C<CHECK_UNDERFLOW(bn, decrement, context)>
94 If subtracting C<decrement> (a positive number) from the exponent
95 of C<bn> would cause underflow, returns true.
101 #define CHECK_UNDERFLOW(bn, decr, context) \
102 ((context)->elimit + ((bn)->expn + (bn)->digits -1) > (decr) ? 0 : 1)
114 True if C<bn> is +Infinity or -Infinity.
120 #define am_INF(bn) ((bn)->flags & BN_INF_FLAG)
126 True if C<bn> is either a quiet or signalling NaN.
132 #define am_NAN(bn) ((bn)->flags & (BN_sNAN_FLAG | BN_qNAN_FLAG))
138 True if C<bn> is a signalling NaN.
144 #define am_sNAN(bn) ((bn)->flags & BN_sNAN_FLAG)
150 True if C<bn> is a quiet NaN.
156 #define am_qNAN(bn) ((bn)->flags & BN_qNAN_FLAG)
158 /* For the sake of debugging */
159 char* BN_lazydbprint(BIGNUM
* foo
) {
161 BN_to_scientific_string(foo
, &s
);
165 /* Internal functions + types */
166 typedef enum { /* Indicate to idivide when to stop */
172 /* Used to restore INTENT(IN) arguments to functions */
173 typedef struct BN_SAVE_PREC
{
179 BN_imultiply(PINTD_ BIGNUM
* result
, BIGNUM
*one
, BIGNUM
*two
,
180 BN_CONTEXT
*context
);
182 BN_idivide(PINTD_ BIGNUM
* result
, BIGNUM
*one
, BIGNUM
*two
,
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
,
203 int BN_set_verybig(PINTD_ BIGNUM
* bn
, BN_CONTEXT
*context
);
209 =head2 Creation and Memory Management Functions
214 BN_new(PINTD_ INTVAL length)>
216 Create a new C<BIGNUM>. C<length> is number of I<decimal> digits
217 required. The bignumber will be equal to zero.
224 BN_new(PINTD_ INTVAL length
) {
226 bn
= (BIGNUM
*)BN_alloc(PINT_
sizeof (BIGNUM
));
228 BN_EXCEPT(PINT_ BN_INSUFFICIENT_STORAGE
, "Cannot allocate new BigNum");
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");
246 BN_grow(PINTD_ BIGNUM *in, INTVAL length)>
248 Grows bn so that it can contain C<length> I<decimal> digits, does not
249 modify the value of the bignumber.
256 BN_grow(PINTD_
NOTNULL(BIGNUM
*in
), INTVAL length
) {
258 if (length
<= in
->nibs
* BN_D_PER_NIB
) {
261 if (length
> BN_MAX_DIGITS
) {
262 BN_EXCEPT(PINT_ BN_OVERFLOW
, "Attempt to grow BIGNUM beyond limits");
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");
276 BN_destroy(PINTD_ BIGNUM *bn)>
278 Frees all the memory used by the BIGNUM.
285 BN_destroy(PINTD_ BIGNUM
*bn
) {
286 PARROT_ASSERT(bn
!=NULL
);
288 BN_free(PINT_ bn
->buffer
);
296 BN_create_context(PINTD_ INTVAL precision)>
298 Creates a new context object, with specified I<precision>, other fields
299 are initialised as follows:
301 elimit = BN_HARD_EXPN_LIMIT (defined during configure)
302 rounding = ROUND_HALF_UP
305 traps = Division by zero, invalid operation, overflow, underflow
306 and rounded are enabled.
307 Lost digits and inexact are disabled.
309 The context object can be destroyed with C<free()>.
316 BN_create_context(PINTD_ INTVAL precision
) {
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");
325 result
->precision
= precision
;
327 result
->elimit
= BN_HARD_EXPN_LIMIT
;
328 result
->rounding
= ROUND_HALF_UP
;
329 result
->extended
= 1;
331 result
->traps
= (BN_F_DIVISION_BY_ZERO
|
332 BN_F_INVALID_OPERATION
|
342 BN_set_digit(PINT_ BIGNUM* bn, INTVAL pos, INTVAL value)>
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
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
) {
359 PARROT_ASSERT(value
< 10);
360 PARROT_ASSERT(value
> -1);
361 BN_setd(bn
, pos
, value
);
362 if (bn
->digits
< pos
+1) {
371 BN_get_digit(PINTD_ BIGNUM* bn, INTVAL pos)>
373 Get the value of the decimal digit at C<pos>, returns -1 if C<pos> is
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
);
389 =item C<int BN_set_inf(PINTD_ BIGNUM* bn)>
391 =item C<int BN_set_qNAN(PINTD_ BIGNUM* bn)>
393 =item C<int BN_set_sNAN(PINTD_ BIGNUM* bn)>
395 Sets its argument to appropriate value.
397 Infinity is represented as having zero digits, an undefined exponent
398 and private C<flags> set to C<BN_inf_FLAGS>.
400 sNAN is represented as having zero digits, an undefined exponent, an
401 undefined sign and both qNAN and sNAN bits set.
403 qNAN is represented as having zero digits, an undefined exponent
404 and only the qNAN bit set.
410 int BN_set_inf(PINTD_ BIGNUM
* bn
) {
411 PARROT_ASSERT(bn
!= NULL
);
413 bn
->flags
= (bn
->flags
& (~(UINTVAL
)255)) | BN_INF_FLAG
;
417 int BN_set_qNAN(PINTD_ BIGNUM
* bn
) {
418 PARROT_ASSERT(bn
!= NULL
);
420 bn
->flags
= (bn
->flags
& (~(UINTVAL
)255)) | BN_qNAN_FLAG
;
424 int BN_set_sNAN(PINTD_ BIGNUM
* bn
) {
425 PARROT_ASSERT(bn
!= NULL
);
427 bn
->flags
= (bn
->flags
& (~(UINTVAL
)255)) | BN_qNAN_FLAG
| BN_sNAN_FLAG
;
434 BN_set_verybig(PINTD_ BIGNUM* bn, BN_CONTEXT *context)>
436 Used when an operation has overflowed, sets C<bn> according to
437 C<< context->rounding >> and the sign of C<bn>:
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
450 BN_set_verybig(PINTD_ BIGNUM
* bn
, BN_CONTEXT
*context
) {
451 int massive
= 0; /* 0 => inf, 1=> 99999999999 etc...*/
452 switch (context
->rounding
) {
454 case ROUND_HALF_EVEN
:
457 if (!bn
->sign
) massive
= 1;
460 if (bn
->sign
) massive
= 1;
466 BN_EXCEPT(PINT_ BN_INVALID_OPERATION
,
467 "Unknown rounding during overflow");
469 if (context
->precision
> 0 && massive
) {
471 BN_grow(PINT_ bn
, context
->precision
);
472 for (i
= 0; i
< context
->precision
; i
++) {
475 bn
->digits
= context
->precision
;
476 bn
->expn
= context
->elimit
- context
->precision
+ 1;
479 BN_set_inf(PINT_ bn
);
486 BN_copy(PINTD_ BIGNUM* one, BIGNUM* two)>
488 Copies two into one, returning one for convenience.
495 BN_copy(PINTD_ BIGNUM
* one
, BIGNUM
* two
) {
496 PARROT_ASSERT(one
!= NULL
); PARROT_ASSERT(two
!= NULL
);
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
;
512 BN_new_from_int(PINTD_ INTVAL value)>
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.
523 BN_new_from_int(PINTD_ INTVAL value
) {
526 new = BN_new(PINT_ BN_D_PER_INT
);
533 current
= value
% 10;
534 BN_setd(new, i
, current
);
546 BN_PRINT_DEBUG(BIGNUM *bn, char* mesg)>
548 Dump the bignum for testing, along with a little message.
555 BN_PRINT_DEBUG(BIGNUM
*bn
, char* mesg
) {
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);
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");
575 BN_nonfatal(PINTD_ BN_CONTEXT *context, BN_EXCEPTIONS except, char *msg)>
577 When an exceptional condition occurs after which execution could
578 continue. If context specifies that death occurs, then so be it.
585 BN_nonfatal(PINTD_ BN_CONTEXT
*context
, BN_EXCEPTIONS except
, const char *msg
) {
586 /* See extended standard for details */
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
);
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
);
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
);
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
);
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
);
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
);
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
);
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
);
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
);
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
);
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
);
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
);
674 /* something was rounded */
675 context
->flags
|= BN_F_ROUNDED
;
676 if (context
->traps
& (BN_F_ROUNDED
)) {
677 BN_EXCEPT(PINT_ except
, msg
);
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
);
688 BN_EXCEPT(PINT_ BN_INVALID_OPERATION
, "An unknown error occurred");
696 BN_exception(PINTD_ BN_EXCEPTIONS exception, const char* message)>
698 Throw `exception'. Should be accessed via a C<BN_EXCEPT> macro, this
699 version is provided until Parrot exceptions are sorted out properly.
706 BN_exception(PINTD_ BN_EXCEPTIONS exception
, const char* message
) {
707 printf("Exception %d %s\n", exception
, message
);
714 BN_to_scientific_string(PINTD_ BIGNUM*bn, char **dest)>
716 Converts bn into a scientific representation, stored in dest.
719 BN_to_engineering_string(PINTD_ BIGNUM*bn, char **dest)>
721 Converts C<*bn> into a engineering representation, stored in C<**dest>.
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.
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.
729 Memory pointed to by C<dest> is not freed by this function.
737 BN_to_scientific_string(PINTD_ BIGNUM
*bn
, char **dest
) {
738 BN_to_scieng_string(PINT_ bn
, dest
, 0);
741 BN_to_engineering_string(PINTD_ BIGNUM
*bn
, char **dest
) {
742 BN_to_scieng_string(PINT_ bn
, dest
, 1);
748 BN_to_scieng_string(PINTD_ BIGNUM* bn, char **dest, int eng)>
750 Does the heavy string handling work, C<eng> defines the conversion to
758 BN_to_scieng_string(PINTD_ BIGNUM
* bn
, char **dest
, int eng
) {
760 INTVAL adj_exp
= 0; /* as bn->expn is relative to 0th digit */
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
);
766 BN_EXCEPT(PINT_ BN_INSUFFICIENT_STORAGE
,
767 "Cannot create buffer to hold output string");
772 /* Do we have a special value? */
777 else { /* must be signalling */
784 if (bn
->sign
) *cur
++ = '-';
785 strcpy(cur
, "Infinity");
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) {
796 /*pad with zeros if appropriate, plonk a point where we want it */
797 if (bn
->digits
+ bn
->expn
<= 0) {
800 if (bn
->digits
+bn
->expn
<0) *cur
++ = '.';
801 for (i
=1; i
<= -(bn
->digits
+ bn
->expn
); i
++) *cur
++ = '0';
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
);
809 else { /* Use exponential notation, different for sci and eng */
810 if (bn
->sign
) *cur
++ = '-'; /* We don't prefix '+' */
815 deficit
= (-adj_exp
) % 3;
820 else if (deficit
== 2) {
826 deficit
= adj_exp
% 3;
829 /* so, d = 0. x.yyyy, d=1 xx.y d=2 xxx.yyy special case if xxx*/
831 *cur
++ = '0' + BN_getd(bn
,bn
->digits
-1);
832 if (deficit
== 0) *cur
++ = '.';
834 if (bn
->digits
== 1) {
836 if (deficit
== 1) *cur
++ = '.';
839 else if (bn
->digits
== 2) {
840 *cur
++ = '0' + BN_getd(bn
,bn
->digits
-2);
841 if (deficit
== 1) *cur
++ = '.';
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
++ = '.';
850 for (cur_dig
=bn
->digits
-4; cur_dig
>-1; cur_dig
--) {
851 *cur
++ = '0' + BN_getd(bn
, cur_dig
);
856 sprintf(cur
, "%+i", adj_exp
);
859 else { /* Scientific */
860 if (bn
->digits
== 1) { /* We don't want 1.E+7 */
861 *cur
++ = '0'+ BN_getd(bn
, 0);
863 else { /* We have x.xE */
864 *cur
++ = '0' + BN_getd(bn
, bn
->digits
-1);
866 for (cur_dig
= bn
->digits
-2; cur_dig
> -1; cur_dig
--) {
867 *cur
++ = '0' + BN_getd(bn
, cur_dig
);
870 *cur
++ = 'E'; /* Eza Good, Eza Good */
871 sprintf(cur
, "%+i", adj_exp
);
881 BN_from_string(PINTD_ char* s2, BN_CONTEXT *context)>
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.
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.
893 Does not yet check for exponent overflow.
900 BN_from_string(PINTD_
char* s2
, BN_CONTEXT
*context
) {
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 '.' */
914 int non_zero_digits
= 0; /* have we seen *any* digits */
915 int seen_plus
= 0; /* was number prefixed with '+' */
920 temp
= BN_new(PINT_
1); /* We store coeff reversed in temp */
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 */
929 if (seen_dot
) fake_exponent
--;
936 BN_grow(PINT_ temp
, pos
+10);
937 BN_setd(temp
, pos
, digit
);
943 else if (isdigit((unsigned char)*s2
) && in_exp
) {
944 exponent
= 10 * exponent
+ (*s2
- '0'); /*XXX: overflow check */
946 else if (!in_number
) {
947 /* we've not yet seen any digits */
949 if (seen_plus
|| negative
|| seen_dot
) {
950 if (!context
->extended
) {
951 BN_EXCEPT(PINT_ BN_CONVERSION_SYNTAX
,
952 "Incorrect number format");
960 else if (*s2
== '.') {
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");
973 seen_plus
= 1; /* be very quiet */
975 else if (context
->extended
) {
976 if (*s2
== 'i' || *s2
== 'I') {
977 if (!strncasecmp("inf", s2
, 4)) { /* We include the \0 */
979 /* For certain, restricted values of infinity */
982 else if (!strncasecmp("infinity", s2
, 9)) {
991 else if (*s2
== 'n' || *s2
== 'N') {
992 qNAN
= 1; /* Don't need to check, as default.. */
995 else if (*s2
== 's' || *s2
== 'S') {
996 if (!strncasecmp("snan", s2
, 5)) {
1007 } /* don't know, not in extended mode... */
1009 BN_EXCEPT(PINT_ BN_CONVERSION_SYNTAX
,
1010 "Incorrect number format");
1014 /* we've seen some digits, are we done yet? */
1015 if (!seen_dot
&& *s2
== '.' && !in_exp
) {
1018 else if (!seen_e
&& (*s2
== 'e' || *s2
== 'E')) {
1022 else if (seen_e
&& !exp_sign
) {
1026 else if (*s2
== '-') {
1030 if (!context
->extended
) {
1031 BN_EXCEPT(PINT_ BN_CONVERSION_SYNTAX
,
1032 "Incorrect number format");
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");
1049 s2
++; /* rinse, lather... */
1052 if (!(qNAN
|| sNAN
|| infinity
)) {
1053 if (in_number
&& !pos
) { /* Only got zeros */
1055 BN_setd(temp
, 0, 0);
1058 if (pos
==0) { /* This includes ".e+20" */
1059 if (!context
->extended
) {
1060 BN_EXCEPT(PINT_ BN_CONVERSION_SYNTAX
, "no digits in string");
1068 result
= BN_new(pos
+1);
1070 /* copy reversed string of digits backwards into result */
1071 if (!(qNAN
|| sNAN
|| infinity
)) { /* Normal */
1074 for (i
=0; i
< temp
->digits
; i
++) {
1075 BN_setd(result
, i
, BN_getd(temp
, temp
->digits
-i
-1));
1078 result
->sign
= negative
;
1079 result
->digits
= pos
;
1080 if (exp_sign
== -1) {
1081 result
->expn
= fake_exponent
- exponent
;
1084 result
->expn
= fake_exponent
+ exponent
;
1087 else { /* Special */
1089 BN_set_inf(PINT_ result
);
1090 result
->sign
= negative
;
1093 BN_set_sNAN(PINT_ result
);
1096 BN_set_qNAN(PINT_ result
);
1101 BN_destroy(PINT_ temp
);
1102 BN_really_zero(PINT_ result
,context
->extended
);
1110 BN_strip_lead_zeros(PINTD_ BIGNUM* bn, BN_CONTEXT *context)>
1112 Removes any zeros before the msd and after the lsd.
1119 BN_strip_lead_zeros(PINTD_ BIGNUM
* bn
, BN_CONTEXT
*context
) {
1122 if (bn
->digits
== 0) return 0; /* Cannot "fail" with special nums */
1126 while (0==BN_getd(bn
, msd
) && msd
> 0) {
1130 bn
->digits
-= bn
->digits
-1 - msd
;
1136 BN_strip_tail_zeros(PINTD_ BIGNUM *bn, BN_CONTEXT *context)>
1138 Removes trailing zeros and increases the exponent appropriately.
1139 Does not remove zeros before the decimal point.
1146 BN_strip_tail_zeros(PINTD_ BIGNUM
*bn
, BN_CONTEXT
*context
) {
1151 while (0==BN_getd(bn
, lsd
)) {
1154 if (bn
->expn
>= 0) {
1155 lsd
= 0; /* units column */
1157 else if (bn
->expn
+ lsd
> 0) {
1160 for (i
=0; i
< bn
->digits
-lsd
; i
++) {
1161 BN_setd(bn
, i
, BN_getd(bn
, i
+lsd
));
1164 if (CHECK_OVERFLOW(bn
, lsd
, context
)) {
1165 BN_EXCEPT(PINT_ BN_OVERFLOW
, "overflow when striping zeros");
1174 BN_make_integer(PINTD_ BIGNUM* bn, BN_CONTEXT* context)>
1176 Convert the number to a plain integer I<if> precision such that this
1184 BN_make_integer(PINTD_ BIGNUM
* bn
, BN_CONTEXT
* context
) {
1186 if (bn
->expn
> 0 && bn
->digits
+ bn
->expn
<= context
->precision
) {
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
));
1192 for (i
=0; i
<bn
->expn
; i
++) {
1195 bn
->digits
+= bn
->expn
;
1199 return; /* XXX: fixed precision, bigints */
1206 BN_really_zero(PINTD_ BIGNUM* bn, int allow_neg_zero)>
1208 Sets any number which should be zero to a canonical zero.
1210 To check if a number is equal to zero, use C<BN_is_zero()>.
1217 BN_really_zero(PINTD_ BIGNUM
* bn
, int allow_neg_zero
) {
1219 if (bn
->digits
== 0) return;
1220 for (i
=0; i
< bn
->digits
; i
++)
1221 if (BN_getd(bn
, i
) != 0) return;
1225 if (!allow_neg_zero
) bn
->sign
= 0;
1232 BN_round(PINTD_ BIGNUM *bn, BN_CONTEXT* context)>
1234 Rounds C<*bn> according to C<*context>.
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
,
1246 BN_set_sNAN(PINT_ bn
);
1250 BN_iround(PINT_ bn
, context
);
1258 BN_iround(PINTD_ BIGNUM *bn, BN_CONTEXT* context)>
1260 Rounds victim according to context.
1262 Round assumes that any leading zeros are significant (after an
1263 addition operation, for instance).
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
1272 eg. for 1.234567E+3 with rounding of C<ROUND_DOWN>
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
1286 BN_iround(PINTD_ BIGNUM
*bn
, BN_CONTEXT
* context
) {
1287 PARROT_ASSERT(bn
!= NULL
);
1288 PARROT_ASSERT(context
!= NULL
);
1290 if (bn
->digits
== 0) {
1291 return 0; /* rounding special values always works */
1294 if (context
->precision
< 1) { /* Rounding a BigInt or fixed */
1295 BN_round_as_integer(PINT_ bn
, context
);
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 */
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");
1315 if (context
->rounding
== ROUND_DOWN
) {
1316 return BN_round_down(PINT_ bn
, context
);
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
);
1323 return BN_round_down(PINT_ bn
, context
);
1326 else if (context
->rounding
== ROUND_HALF_EVEN
) {
1328 if (BN_getd(bn
, bn
->digits
- context
->precision
-1) > 5) {
1329 return BN_round_up(PINT_ bn
, context
);
1331 else if (BN_getd(bn
, bn
->digits
- context
->precision
-1) < 5) {
1332 return BN_round_down(PINT_ bn
, context
);
1335 INTVAL i
= bn
->digits
- context
->precision
-2;
1338 if (BN_getd(bn
, i
) != 0) {
1339 return BN_round_up(PINT_ bn
, context
);
1344 switch (BN_getd(bn
, bn
->digits
-context
->precision
)) {
1350 return BN_round_down(PINT_ bn
, context
);
1352 return BN_round_up(PINT_ bn
, context
);
1357 else if (context
->rounding
== ROUND_CEILING
) {
1360 return BN_round_down(PINT_ bn
, context
);
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
);
1367 return BN_round_down(PINT_ bn
, context
);
1369 else if (context
->rounding
== ROUND_FLOOR
) {
1372 return BN_round_down(PINT_ bn
, context
);
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
);
1379 return BN_round_down(PINT_ bn
, context
);
1381 BN_EXCEPT(PINT_ BN_INVALID_OPERATION
, "Unknown rounding attempted");
1389 BN_round_up(PINTD_ BIGNUM *bn, BN_CONTEXT* context)>
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>.
1400 BN_round_up(PINTD_ BIGNUM
*bn
, BN_CONTEXT
* context
) {
1403 /* Do a cheap num += 1E+(num->expn) */
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);
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");
1417 bn
->digits
= context
->precision
+1;
1418 return BN_iround(PINT_ bn
, context
);
1421 INTVAL extra
= bn
->digits
- context
->precision
;
1422 if (CHECK_OVERFLOW(bn
, extra
, context
)) {
1423 BN_EXCEPT(PINT_ BN_OVERFLOW
, "overflow while rounding");
1426 bn
->digits
= context
->precision
;
1434 BN_round_down(PINT_ BIGNUM *bn, BN_CONTEXT* context)>
1436 Truncates the coefficient of C<bn> to have C<precision> digits. Do
1437 not call this function with non-positive precision.
1444 BN_round_down(PINT_ BIGNUM
*bn
, BN_CONTEXT
* context
) {
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
);
1452 extra
= bn
->digits
- context
->precision
;
1453 if (CHECK_OVERFLOW(bn
, extra
, context
)) {
1454 BN_EXCEPT(PINT_ BN_OVERFLOW
, "overflow while rounding");
1457 bn
->digits
= context
->precision
;
1465 BN_round_as_integer(PINTD_ BIGNUM *bn, BN_CONTEXT *context)>
1467 C<precision> must be less than one. This rounds so that C<expn> is at
1468 least C<precision>. Name is slightly misleading.
1475 BN_round_as_integer(PINTD_ BIGNUM
*bn
, BN_CONTEXT
*context
) {
1477 BN_CONTEXT temp_context
;
1479 if (bn
->expn
< context
->precision
) {
1481 /* Are we losing information? */
1483 i
< (context
->precision
- bn
->expn
) && i
<bn
->digits
;
1485 if (BN_getd(bn
, i
) != 0) {
1486 BN_nonfatal(PINT_ context
, BN_LOST_DIGITS
,
1487 "digits lost while rounding");
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);
1502 BN_iround(PINT_ bn
, &temp_context
);
1505 BN_iround(PINT_ bn
, &temp_context
);
1507 BN_really_zero(PINT_ bn
, context
->extended
);
1509 /* XXX: if using warning flags on context, | with temp context here */
1519 =head2 Arithmetic operations
1521 Operations are performed like this:
1527 Both operands are rounded to have no more than C<< context->precision >>
1532 The operation is computed.
1534 =item Rounding of result
1536 The result is then rounded to context->precision digits.
1538 =item Conversion to zero and integerisation
1540 If the result is equal to zero, it is made exactly zero.
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.
1547 The general form for all arithmetic operations is:
1549 void BN_operation(result, one, two, context)
1561 BN_arith_setup(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
1562 BN_CONTEXT *context, BN_SAVE_PREC* restore)>
1564 Rounds one and two ready for arithmetic operation.
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
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.
1575 If overflow or underflow occurs during rounding, the numbers will be
1576 modified to the appropriate representation and will not be restorable.
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
);
1590 restore
->one
= *one
;
1591 restore
->two
= *two
;
1598 BN_arith_cleanup(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
1599 BN_CONTEXT *context, BN_SAVE_PREC* restore)>
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
1610 BN_arith_cleanup(PINTD_ BIGNUM
* result
, BIGNUM
*one
, BIGNUM
*two
,
1611 BN_CONTEXT
*context
, BN_SAVE_PREC
* restore
) {
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
;
1621 one
->digits
= restore
->one
.digits
;
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
;
1629 two
->digits
= restore
->two
.digits
;
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
);
1642 BN_strip_lead_zeros(PINT_ result
, context
);
1643 BN_really_zero(PINT_ result
, context
->extended
);
1644 BN_make_integer(PINT_ result
, context
);
1650 BN_align(PINTD_ BIGNUM* one, BIGNUM* two)>
1652 Adds zero digits so that decimal points of each number are at the same
1660 BN_align(PINTD_ BIGNUM
* one
, BIGNUM
* two
) {
1664 diff
= one
->expn
- two
->expn
;
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
) {
1675 BN_grow(PINT_ two
, one
->digits
);
1676 for (i
=two
->digits
; i
<one
->digits
; i
++) {
1679 two
->digits
= one
->digits
;
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. */
1688 /* force smallest exponent in two */
1696 if (one
->digits
+ diff
< two
->digits
) {
1697 final
= two
->digits
;
1700 final
= one
->digits
+ diff
;
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
++)
1709 /* Add zeros to start of one */
1710 for (i
=one
->digits
+ diff
; i
< final
; i
++)
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
));
1717 /* Set end of one to zeros */
1718 for (i
=0; i
< diff
; i
++)
1721 one
->digits
= two
->digits
= final
;
1722 one
->expn
= two
->expn
;
1729 BN_add(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two, BN_CONTEXT *context)>
1731 Adds one to two, returning value in result.
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
,
1747 BN_set_qNAN(PINT_ result
);
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
);
1759 BN_set_inf(PINT_ result
);
1760 result
->sign
= one
->sign
;
1764 /* So we've only got one infinity... */
1765 BN_set_inf(PINT_ result
);
1766 result
->sign
= am_INF(one
) ? one
->sign
: two
->sign
;
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
)) {
1774 BN_setd(result
, 0, 0);
1775 if (one
->sign
& two
->sign
) {
1778 else if (context
->rounding
== ROUND_FLOOR
&& (one
->sign
^ two
->sign
)) {
1788 BN_arith_setup(PINT_ result
, one
, two
, context
, &restore
);
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
);
1794 else if (one
->sign
&& two
->sign
) { /* -a + -b = -(a+b) */
1795 BN_iadd(PINT_ result
, one
, two
, context
);
1798 else if (two
->sign
) { /* a + -b = (a-b) */
1799 BN_isubtract(PINT_ result
, one
, two
, context
);
1801 else { /* a + b = (a+b) */
1802 BN_iadd(PINT_ result
, one
, two
, context
);
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
)) {
1816 BN_iadd(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
1817 BN_CONTEXT *context)>
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.
1828 BN_iadd(PINTD_ BIGNUM
* result
, BIGNUM
*one
, BIGNUM
*two
,
1829 BN_CONTEXT
*context
) {
1833 /* Make sure we don't do work we don't need, or add precision where
1835 if (BN_is_zero(PINT_ one
, context
)) {
1836 BN_copy(PINT_ result
, two
);
1840 else if (BN_is_zero(PINT_ two
, context
)) {
1841 BN_copy(PINT_ result
, one
);
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 */
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
);
1858 for (i
= i
; i
< context
->precision
; i
++) {
1859 BN_setd(result
, context
->precision
- i
, 0);
1861 BN_setd(result
, 0, 1);
1862 result
->digits
= context
->precision
+ 1;
1864 result
->expn
= one
->expn
+ one
->digits
- context
->precision
- 1;
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
);
1874 for (i
= i
; i
< context
->precision
; i
++) {
1875 BN_setd(result
, context
->precision
- i
, 0);
1877 BN_setd(result
, 0, 1);
1878 result
->digits
= context
->precision
+ 1;
1880 result
->expn
= two
->expn
+ two
->digits
- context
->precision
- 1;
1884 /* Ok, we can't be lazy, we'll have to do it all ourselves */
1885 BN_align(PINT_ one
, two
);
1887 BN_grow(PINT_ result
, one
->digits
+ 1);
1890 for (i
=0; i
< one
->digits
; i
++) {
1891 carry
+= BN_getd(one
, i
) + BN_getd(two
, i
);
1893 BN_setd(result
, i
, dig
);
1897 BN_setd(result
, i
, carry
);
1898 result
->digits
= i
+1;
1904 result
->expn
= one
->expn
;
1911 BN_subtract(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
1912 BN_CONTEXT *context)>
1914 Subtracts C<*two> from C<*one>, returning value in C<*result>.
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");
1931 BN_set_qNAN(PINT_ result
);
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
);
1943 BN_set_inf(PINT_ result
);
1944 result
->sign
= one
->sign
;
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
));
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
)) {
1958 BN_setd(result
, 0, 0);
1959 if (one
->sign
&& !two
->sign
&& context
->extended
) {
1962 else if (context
->rounding
== ROUND_FLOOR
&&
1963 (one
->sign
== two
->sign
) &&
1964 context
->extended
) {
1973 BN_arith_setup(PINT_ result
, one
, two
, context
, &restore
);
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
);
1980 else if (one
->sign
&& two
->sign
) { /* -a - -b = (b-a) */
1981 BN_isubtract(PINT_ result
, two
, one
, context
);
1983 else if (two
->sign
) { /* a - -b = (a+b) */
1984 BN_iadd(PINT_ result
, one
, two
, context
);
1986 else { /* a - b = (a-b) */
1987 BN_isubtract(PINT_ result
, one
, two
, context
);
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
)) {
2001 BN_isubtract(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
2002 BN_CONTEXT *context)>
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.
2013 BN_isubtract(PINTD_ BIGNUM
* result
, BIGNUM
*one
, BIGNUM
*two
,
2014 BN_CONTEXT
*context
) {
2017 /* Make sure we don't do work we don't need, or add precision where
2019 if (BN_is_zero(PINT_ one
, context
)) {
2020 BN_copy(PINT_ result
, two
);
2024 else if (BN_is_zero(PINT_ two
, context
)) {
2025 BN_copy(PINT_ result
, one
);
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
2034 if (context
->precision
> -1 &&
2035 one
->expn
> two
->expn
+ context
->precision
+ 1) {
2036 BN_grow(PINT_ result
, context
->precision
+ 1);
2038 for (i
=0; i
<one
->digits
; i
++) {
2039 carry
= carry
+ BN_getd(one
, i
);
2041 BN_setd(result
, context
->precision
- one
->digits
+ i
+ 1,
2046 BN_setd(result
, context
->precision
- one
->digits
+ i
+ 1,
2051 for (i
=0; i
<context
->precision
+ 1 - one
->digits
; i
++) {
2052 BN_setd(result
, i
, 9);
2054 result
->expn
= one
->expn
+ one
->digits
- context
->precision
-1;
2056 result
->digits
= context
->precision
+ 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);
2063 for (i
=0; i
<two
->digits
; i
++) {
2064 carry
= carry
+ BN_getd(two
, i
);
2066 BN_setd(result
, context
->precision
- two
->digits
+ i
+ 1,
2071 BN_setd(result
, context
->precision
- two
->digits
+ i
+ 1,
2076 for (i
=0; i
<context
->precision
+ 1 - two
->digits
; i
++) {
2077 BN_setd(result
, i
, 9);
2079 result
->expn
= two
->expn
+ two
->digits
- context
->precision
-1;
2081 result
->digits
= context
->precision
+ 1;
2086 BN_align(PINT_ one
, two
);
2088 /* as a-b == -(b-a), we find larger of
2089 a and b and make sure it is in one */
2091 for (i
=one
->digits
-1; i
>-1; i
--) {
2092 carry
= BN_getd(one
, i
) - BN_getd(two
, i
);
2096 if (!carry
) { /* a==b*/
2099 BN_setd(result
, 0,0);
2102 else if (carry
< 0) { /* b > a */
2103 BN_isubtract(PINT_ result
, two
, one
, context
);
2107 BN_grow(PINT_ result
, one
->digits
+ 1);
2110 for (i
=0; i
<one
->digits
; i
++) {
2111 carry
= carry
+ BN_getd(one
, i
) - BN_getd(two
,i
);
2113 BN_setd(result
, i
, 10+carry
);
2117 BN_setd(result
, i
, carry
);
2122 PARROT_ASSERT(!carry
); /* as to get here a > b*/
2124 result
->digits
= one
->digits
;
2125 result
->expn
= one
->expn
;
2134 BN_plus(PINTD_ BIGNUM* result, BIGNUM *one, BN_CONTEXT *context)>
2136 Perform unary C<+> on C<*one>. Does all the rounding and what have you.
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
,
2149 BN_set_qNAN(PINT_ result
);
2152 else { /* Infinity */
2153 BN_set_inf(PINT_ result
);
2154 result
->sign
= one
->sign
;
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
);
2168 BN_minus(PINTD_ BIGNUM* result, BIGNUM *one, BN_CONTEXT *context)>
2170 Perform unary C<-> (minus) on C<*one>. Does all the rounding and what
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
,
2184 BN_set_qNAN(PINT_ result
);
2187 else { /* Infinity */
2188 BN_set_inf(PINT_ result
);
2189 result
->sign
= 1 & (1 ^ one
->sign
);
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
);
2204 BN_compare(PINT_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
2205 BN_CONTEXT *context)>
2207 Numerically compares C<*one> and C<*two>, storing the result (as a
2208 BIGNUM) in C<*result>.
2210 result = 1 => one > two
2211 result = -1 => two > one
2212 result = 0 => one == two
2219 BN_compare(PINT_ BIGNUM
* result
, BIGNUM
*one
, BIGNUM
*two
,
2220 BN_CONTEXT
*context
) {
2223 /* Special values */
2224 if (one
->digits
== 0 || two
->digits
==0) {
2226 if (am_NAN(one
) || am_NAN(two
)) {
2227 if (am_sNAN(one
) || am_sNAN(two
)) {
2228 BN_nonfatal(PINT_ context
, BN_INVALID_OPERATION
,
2231 BN_set_qNAN(PINT_ result
);
2234 /* Otherwise at least one of the operands is an infinity */
2235 if (one
->sign
!= two
->sign
) {
2236 cmp
= one
->sign
? -1 : 1;
2238 else if (am_INF(one
) && am_INF(two
)) {
2241 else if (am_INF(one
)) {
2242 cmp
= one
->sign
? -1 : 1;
2245 cmp
= one
->sign
? 1 : -1;
2249 BN_arith_setup(PINT_ result
, one
, two
, context
, NULL
);
2250 cmp
= BN_comp(PINT_ one
, two
, context
);
2256 BN_setd(result
, 0, 0);
2260 BN_setd(result
, 0, 1);
2264 BN_setd(result
, 0, 1);
2267 BN_arith_cleanup(PINT_ result
, one
, two
, context
, NULL
);
2273 BN_multiply(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
2274 BN_CONTEXT *context)>
2276 Multiplies C<*one> and C<*two>, storing the result in C<*result>.
2283 BN_multiply(PINTD_ BIGNUM
* result
, BIGNUM
*one
, BIGNUM
*two
,
2284 BN_CONTEXT
*context
)
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");
2292 BN_set_qNAN(PINT_ result
);
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
);
2304 /* (anything but 0) * Inf => Inf */
2305 BN_set_inf(PINT_ result
);
2306 result
->sign
= 1 & (one
->sign
^ two
->sign
);
2310 BN_arith_setup(PINT_ result
, one
, two
, context
, NULL
);
2312 BN_imultiply(PINT_ result
, one
, two
, context
);
2314 BN_strip_lead_zeros(PINT_ result
, context
);
2315 BN_arith_cleanup(PINT_ result
, one
, two
, context
, NULL
);
2322 BN_imultiply(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
2323 BN_CONTEXT *context)>
2325 Multiplication without the rounding and other set up.
2332 BN_imultiply(PINTD_ BIGNUM
* result
, BIGNUM
*one
, BIGNUM
*two
,
2333 BN_CONTEXT
*context
)
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);
2343 /* make sure largest coeff is in one */
2344 if (one
->digits
< two
->digits
) {
2350 /* multiply element by element */
2351 for (i
=0; i
<two
->digits
; i
++) {
2352 dig
= BN_getd(two
, i
);
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);
2360 BN_setd(result
, i
+j
, carry
);
2364 /* extend if there's still stuff to take care of */
2366 result
->digits
= one
->digits
+ two
->digits
+ 1;
2368 else if (BN_getd(result
, one
->digits
+ two
->digits
- 1)) {
2369 result
->digits
= one
->digits
+ two
->digits
;
2372 result
->digits
= one
->digits
+ two
->digits
- 1;
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");
2380 if (i
< -context
->elimit
) {
2381 BN_EXCEPT(PINT_ BN_UNDERFLOW
, "underflow in multiplication");
2385 result
->sign
= 1 & (one
->sign
^ two
->sign
);
2392 BN_divide(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
2393 BN_CONTEXT *context)>
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()>.
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
2410 BN_divide(PINTD_ BIGNUM
* result
, BIGNUM
*one
, BIGNUM
*two
,
2411 BN_CONTEXT
*context
) {
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
,
2420 BN_set_qNAN(PINT_ result
);
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
);
2430 if (BN_is_zero(PINT_ two
, context
)) {
2431 BN_nonfatal(PINT_ context
, BN_DIVISION_BY_ZERO
,
2432 "Inf / 0 in divide");
2434 BN_set_inf(PINT_ result
);
2435 result
->sign
= 1 & (one
->sign
^ two
->sign
);
2438 /* so we're left with x / Inf */
2441 BN_setd(result
, 0, 0);
2442 result
->sign
= 1 & (one
->sign
^ two
->sign
);
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
);
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
);
2461 /* We're left with the case that only arg 1 is zero */
2462 if (BN_is_zero(PINT_ one
, context
)) {
2465 BN_setd(result
, 0 ,0);
2466 result
->sign
= 1 & (one
->sign
^ two
->sign
);
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
);
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 */
2477 BN_strip_lead_zeros(PINT_ result
, context
);
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
);
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)) {
2500 BN_round_down(PINT_ result
, context
);
2503 BN_round_up(PINT_ result
, context
);
2507 BN_nonfatal(PINT_ context
, BN_INEXACT
,
2508 "Loss of precision in divide");
2509 BN_round_up(PINT_ result
, context
);
2513 if (BN_getd(result
, 0) !=0) {
2514 BN_nonfatal(PINT_ context
, BN_INEXACT
,
2515 "Loss of precision in divide");
2517 else if (!BN_is_zero(PINT_ result
, context
)) {
2518 BN_nonfatal(PINT_ context
, BN_INEXACT
,
2519 "Loss of precision in divide");
2521 BN_round_down(PINT_ result
, context
);
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
);
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");
2537 BN_round_down(PINT_ result
, context
);
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
);
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
);
2550 BN_round_down(PINT_ result
, context
);
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
);
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");
2566 BN_round_down(PINT_ result
, context
);
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
);
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
);
2579 BN_round_down(PINT_ result
, context
);
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
);
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");
2600 context
->traps
= save_lost
;
2601 context
->flags
= (context
->flags
& ~(unsigned char)BN_F_LOST_DIGITS
)
2602 | (flags_save
& BN_F_LOST_DIGITS
);
2605 BN_really_zero(PINT_ result
, context
->extended
);
2607 BN_strip_tail_zeros(PINT_ result
, context
);
2609 BN_make_integer(PINT_ result
, context
);
2611 /* Remove trailing zeros if positive exponent */
2612 if (result
->expn
> 0) {
2614 for (i
=0; i
<result
->digits
; i
++) {
2615 if (BN_getd(result
, i
) != 0) break;
2618 for (j
=i
; j
<result
->digits
; j
++) {
2619 BN_setd(result
, j
-i
, BN_getd(result
, j
));
2622 if (CHECK_OVERFLOW(result
, i
, context
)) {
2623 BN_EXCEPT(PINT_ BN_OVERFLOW
, "overflow in divide");
2628 BN_destroy(PINT_ rem
);
2634 BN_divide_integer(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
2635 BN_CONTEXT *context)>
2637 Places the integer part of C<*one> divided by C<*two> into C<*result>.
2644 BN_divide_integer(PINTD_ BIGNUM
* result
, BIGNUM
*one
, BIGNUM
*two
,
2645 BN_CONTEXT
*context
) {
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");
2654 BN_set_qNAN(PINT_ result
);
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
);
2664 if (BN_is_zero(PINT_ two
, context
)) {
2665 BN_nonfatal(PINT_ context
, BN_DIVISION_BY_ZERO
,
2666 "Inf / 0 in divide-integer");
2668 BN_set_inf(PINT_ result
);
2669 result
->sign
= 1 & (one
->sign
^ two
->sign
);
2672 /* so we're left with x / Inf */
2675 BN_setd(result
, 0, 0);
2676 result
->sign
= 1 & (one
->sign
^ two
->sign
);
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
);
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
);
2695 /* We're left with the case that only arg 1 is zero */
2696 if (BN_is_zero(PINT_ one
, context
)) {
2699 BN_setd(result
, 0 ,0);
2700 result
->sign
= 1 & (one
->sign
^ two
->sign
);
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
);
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
);
2718 BN_destroy(PINT_ rem
);
2719 if (result
->expn
!= 0) {
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
));
2726 for (i
=0; i
<result
->expn
; i
++) {
2727 BN_setd(result
, i
,0);
2729 result
->digits
+= result
->expn
;
2733 BN_really_zero(PINT_ result
, context
->extended
);
2734 BN_make_integer(PINT_ result
, context
);
2740 BN_remainder(PINTD_ BIGNUM* result, BIGNUM *one, BIGNUM *two,
2741 BN_CONTEXT *context)>
2743 Places the remainder from divide-integer (above) into C<*result>.
2750 BN_remainder(PINTD_ BIGNUM
* result
, BIGNUM
*one
, BIGNUM
*two
,
2751 BN_CONTEXT
*context
) {
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");
2761 BN_set_qNAN(PINT_ result
);
2764 /* Infinities, first cover Inf rem x and Inf rem Inf */
2766 BN_nonfatal(PINT_ context
, BN_INVALID_OPERATION
,
2767 "x rem Inf in remainder");
2768 BN_set_qNAN(PINT_ result
);
2771 /* now cover x rem Inf => 0 */
2774 BN_setd(result
, 0,0);
2775 result
->sign
= one
->sign
;
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
);
2786 BN_nonfatal(PINT_ context
, BN_INVALID_OPERATION
,
2787 "x rem 0 in remainder");
2788 BN_set_qNAN(PINT_ result
);
2792 if (BN_is_zero(PINT_ one
, context
)) {
2795 BN_setd(result
, 0, 0);
2799 BN_arith_setup(PINT_ result
, one
, two
, context
, NULL
);
2801 BN_idivide(PINT_ fake
, one
, two
, context
, BN_DIV_REMAIN
, result
);
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
);
2815 BN_destroy(PINT_ fake
);
2817 result
->sign
= one
->sign
;
2819 BN_arith_cleanup(PINT_ result
, one
, two
, context
, NULL
);
2824 =item C<BN_idivide(PINT_ BIGNUM *result, BIGNUM *one, BIGNUM *two,
2825 BN_CONTEXT *context, BN_DIV_ENUM operation, BIGNUM *rem)>
2827 Does the heavy work for the various division wossnames.
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
;
2838 /* We assume we've been given something to divide */
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*/
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 */
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)");
2859 if (newexpn
< -context
->elimit
) {
2860 BN_EXCEPT(PINT_ BN_UNDERFLOW
, "underflow in divide (1)");
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;
2870 for (j
=1; j
<=10;j
++) {
2873 BN_imultiply(PINT_ result
, t1
, two
, context
);
2874 cmp
= BN_comp(PINT_ result
, div
, context
);
2876 BN_setd(t2
, value
, j
);
2879 j
= j
+1; /* for multiply below */
2884 if (j
==1 && value
== 0) break; /* don't collect leading 0s */
2885 BN_setd(t2
, value
, j
-1);
2894 BN_imultiply(PINT_ result
, t1
, two
, context
);
2895 BN_isubtract(PINT_ rem
, div
, result
, context
);
2898 /* Are we done yet? */
2899 if (value
&& rem
->digits
==1 && BN_getd(rem
, 0)==0) {
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) {
2912 if (t1
->expn
== context
->precision
-1) break;
2914 if (operation
== BN_DIV_DIVINT
|| operation
== BN_DIV_REMAIN
) {
2915 if (t1
->expn
==0) break;
2917 if (CHECK_UNDERFLOW(t1
, 1, context
)) {
2918 BN_EXCEPT(PINT_ BN_UNDERFLOW
, "underflow in divide (2)");
2921 if (divided
) BN_copy(PINT_ div
, rem
);
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
));
2928 if (t2
->digits
== 0||(!divided
&&!value
)) {
2930 BN_setd(result
, 0, 0);
2934 result
->digits
= t2
->digits
;
2935 result
->sign
= 1&(one
->sign
^ s2
);
2936 result
->expn
= t1
->expn
; /* We know this is fine, from above */
2939 rem
->sign
= 1&(one
->sign
^ s2
);
2941 BN_destroy(PINT_ t1
);
2942 BN_destroy(PINT_ t2
);
2943 BN_destroy(PINT_ div
);
2951 BN_comp(PINTD_ BIGNUM *one, BIGNUM *two, BN_CONTEXT* context)>
2953 Comparison with no rounding etc.
2960 BN_comp(PINTD_ BIGNUM
*one
, BIGNUM
*two
, BN_CONTEXT
* context
) {
2964 BN_strip_lead_zeros(PINT_ one
, context
);
2965 BN_strip_lead_zeros(PINT_ two
, context
);
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 */
2971 return one
->sign
? -1 : 1;
2973 else if (one
->expn
+ one
->digits
> two
->expn
+ two
->digits
) {
2974 return one
->sign
? -1 : 1;
2976 else if (one
->expn
+ one
->digits
< two
->expn
+ two
->digits
) {
2977 return one
->sign
? 1 : -1;
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
;
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
;
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
;
2998 return one
->sign
? -cmp
: cmp
;
3006 BN_power(PINTD_ BIGNUM* result, BIGNUM* bignum,
3007 BIGNUM* expn, BN_CONTEXT* context)>
3009 Calculate C<result> = C<bignum> to the power of C<*expn>;
3016 BN_power(PINTD_ BIGNUM
* result
, BIGNUM
* bignum
,
3017 BIGNUM
* expn
, BN_CONTEXT
* context
)
3019 BN_arith_setup(PINT_ result
, bignum
, expn
, context
, NULL
);
3020 BN_arith_cleanup(PINT_ result
, bignum
, expn
, context
, NULL
);
3026 BN_rescale(PINTD_ BIGNUM* result, BIGNUM* one, BIGNUM* two,
3027 BN_CONTEXT* context)>
3029 Rescales C<*one> to have an exponent of C<*two>.
3036 BN_rescale(PINTD_ BIGNUM
* result
, BIGNUM
* one
, BIGNUM
* two
,
3037 BN_CONTEXT
* context
) {
3039 unsigned char lost
= context
->traps
;
3040 context
->traps
&= ~(unsigned char)BN_F_LOST_DIGITS
;
3042 BN_arith_setup(PINT_ result
, one
, two
, context
, NULL
);
3044 expn
= BN_to_int(PINT_ two
, context
);
3046 context
->traps
= lost
;
3048 BN_arith_cleanup(PINT_ result
, one
, two
, context
, NULL
);
3054 BN_to_int(PINT_ BIGNUM* bn, BN_CONTEXT* context)>
3056 Converts the bignum into an integer, raises overflow if an exact
3057 representation cannot be created.
3064 BN_to_int(PINT_ BIGNUM
* bn
, BN_CONTEXT
* context
) {
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
;
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");
3077 if (bn
->expn
< 0 && bn
->expn
+ bn
->digits
> BN_D_PER_INT
) {
3079 BN_EXCEPT(PINT_ BN_OVERFLOW
, "bignum too large to fit in an int");
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");
3089 if (bn
->digits
+ bn
->expn
> context
->precision
&& context
->precision
> 0) {
3090 BN_EXCEPT(PINT_ BN_LOST_DIGITS
, "digits lost in conv -> int");
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
);
3098 for (i
=0; i
<bn
->expn
; i
++) result
= result
* 10;
3101 for (i
=bn
->digits
-1; i
>-1-bn
->expn
; i
--) {
3102 result
= result
* 10 + BN_getd(bn
, i
);
3106 return bn
->sign
? -result
: result
;
3112 BN_is_zero(BIGNUM* foo, BN_CONTEXT* context)>
3114 Returns a boolean value indicating whether C<*foo> is zero.
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) {
3137 This is currently not used yet. Parrot has no BigNum support yet.
3139 Parrot string playing, exception raising
3143 F<docs/docs/pdds/draft/pdd14_bignum.pod>,
3144 L<https://rt.perl.org/rt3/Ticket/Display.html?id=36330>
3152 * c-file-style: "parrot"
3154 * vim: expandtab shiftwidth=4: