mCRL2
Searching...
No Matches
big_numbers.h
Go to the documentation of this file.
1// Author(s): Jan Friso Groote
2// Copyright: see the accompanying file COPYING or copy at
3// https://github.com/mCRL2org/mCRL2/blob/master/COPYING
4//
6// (See accompanying file LICENSE_1_0.txt or copy at
8//
9
18#ifndef MCRL2_UTILITIES_BIG_NUMBERS_H
19#define MCRL2_UTILITIES_BIG_NUMBERS_H
20
23#include <algorithm>
24#include <limits>
25#include <string>
26
27// Prototype.
28namespace mcrl2
29{
30namespace utilities
31{
32class big_natural_number;
33
34inline std::string pp(const big_natural_number& l);
35
36} // namespace utilities
37} // namespace mcrl2
38
39
40namespace mcrl2
41{
42namespace utilities
43{
44namespace detail
45{
46
47 // Calculate <carry,result>:=n1+n2+carry. The carry can be either 0 or 1, both
48 // at the input and the output.
49 inline std::size_t add_single_number(const std::size_t n1, const std::size_t n2, std::size_t& carry)
50 {
51 assert(carry<=1);
52 std::size_t result=n1+n2+carry;
53 if (carry==0)
54 {
55 if (result<n1)
56 {
57 carry=1;
58 }
59 }
60 else // carry==1
61 {
62 if (result>n1)
63 {
64 carry=0;
65 }
66 }
67 return result;
68 }
69
70 // Calculate <carry,result>:=n1-n2-carry. The carry can be either 0 or 1, both
71 // at the input and the output. If the carry is 1, this indicated that 1 must be subtracted.
72 inline std::size_t subtract_single_number(const std::size_t n1, const std::size_t n2, std::size_t& carry)
73 {
74 assert(carry<=1);
75 std::size_t result=n1-n2-carry;
76 if (carry==0)
77 {
78 if (result>n1)
79 {
80 carry=1;
81 }
82 }
83 else // carry==1
84 {
85 if (result<n1)
86 {
87 carry=0;
88 }
89 }
90 return result;
91 }
92
93 // Calculate <carry,result>:=n1*n2+carry, where the lower bits of the calculation
94 // are stored in the result, and the higher bits are stored in carry.
95 inline std::size_t multiply_single_number(const std::size_t n1, const std::size_t n2, std::size_t& multiplication_carry)
96 {
97 // TODO: It is more concise and efficient to use 128bit machine calculation when available.
98 const int no_of_bits_per_digit=std::numeric_limits<std::size_t>::digits;
99
100 // split input numbers into no_of_bits_per_digit/2 digits
101 std::size_t n1ls = n1 & ((1LL<<(no_of_bits_per_digit/2))-1);
102 std::size_t n1ms = n1 >> (no_of_bits_per_digit/2);
103 std::size_t n2ls = n2 & ((1LL<<(no_of_bits_per_digit/2))-1);
104 std::size_t n2ms = n2 >> (no_of_bits_per_digit/2);
105
106 // First calculate the result of the least significant no_of_bits_per_digit.
107 std::size_t local_carry=0;
109 std::size_t cumulative_carry=local_carry;
110 local_carry=0;
112 cumulative_carry=cumulative_carry+local_carry;
113 local_carry=0;
115 cumulative_carry=cumulative_carry+local_carry;
116
117 // Now calculate the result of the most significant no_of_bits_per_digit.
118 multiplication_carry=cumulative_carry;
119 multiplication_carry=multiplication_carry+((n1ms*n2ls)>>(no_of_bits_per_digit/2));
120 multiplication_carry=multiplication_carry+((n1ls*n2ms)>>(no_of_bits_per_digit/2));
121 multiplication_carry=multiplication_carry+n1ms*n2ms;
122
123 return result;
124 }
125
126 // Calculate <result,remainder>:=(remainder * 2^64 + p) / q assuming the result
127 // fits in 64 bits. More concretely, q>remainder.
128 inline std::size_t divide_single_number(const std::size_t p, const std::size_t q, std::size_t& remainder)
129 {
130 assert(q>remainder);
131 const int no_of_bits_per_digit=std::numeric_limits<std::size_t>::digits;
132
133 // Split input numbers into no_of_bits_per_digit/2 digits.
134 // First get the least significant part.
135 std::size_t pms = (p >> (no_of_bits_per_digit/2)) + (remainder << (no_of_bits_per_digit/2));
136
137 // First divide the most significant part by q.
138 std::size_t resultms = pms/q;
139 assert((resultms >> (no_of_bits_per_digit/2)) == 0);
140 std::size_t remainderms = pms%q;
141
142 // Now obtain the least significant part.
143 std::size_t pls = (p & ((1LL<<(no_of_bits_per_digit/2))-1)) + (remainderms << (no_of_bits_per_digit/2));
144
145 // Second divide the least significant part by q.
146 std::size_t resultls = pls/q;
147 remainder = pls%q;
148 assert((resultls >> (no_of_bits_per_digit/2)) == 0);
149
150 return resultls + (resultms << (no_of_bits_per_digit/2));
151 }
152} // namespace detail
153
154class big_natural_number;
155inline std::string pp(const big_natural_number& l);
156
158{
159 friend std::hash<big_natural_number>;
160 friend inline void swap(big_natural_number& x, big_natural_number& y);
161
162 protected:
163 // Numbers are stored as std::size_t words, with the most significant number last.
164 // Note that the number representation is not unique. Numbers have no trailing
165 // zero's, i.e., this->back()!=0 (if this->size()>0). Therefore their representation is unique.
166 std::vector<std::size_t> m_number;
167
168 /* Multiply the current number by n and add the carry */
169 void multiply_by(std::size_t n, std::size_t carry)
170 {
171 for(std::size_t& i: m_number)
172 {
174 }
175 if (carry)
176 {
177 /* Add an extra digit with the carry */
178 m_number.push_back(carry);
179 }
182 }
183
184 // Remove zero's at the end of m_number vector that are
185 // leading and irrelevant zero's in the number representation.
187 {
188 for( ; m_number.size()>0 && m_number.back()==0 ; )
189 {
190 m_number.pop_back();
191 }
192 }
193
194 // Check that the number is well defined, in the sense
195 // that there are no most significant digits that are zero.
196 void is_well_defined() const
197 {
198 assert(m_number.size()==0 || m_number.back()!=0);
199 }
200
201 // \brief This functions prints a number in internal represenation. This function is useful and only meant for debugging.
202 void print_number(const std::string& s) const
203 {
204 std::cerr << s << " " << m_number.size() << "\n";
205 for(std::size_t i: m_number)
206 {
207 std::cerr << i << " ";
208 }
209 std::cerr << "\n---------------------\n";
210 }
211
212 public:
213 // Default constructor. The value is 0 by default.
215 {
216 }
217
218 // Constructor based on a std::size_t. The value of the number will be n.
219 explicit big_natural_number(const std::size_t n)
220 {
221 if (n>0)
222 {
223 m_number.push_back(n);
224 }
226 }
227
228 // Constructor based on a string. The string is a digital number.
229 big_natural_number(const std::string& s)
230 {
231 for(char c: s)
232 {
233 if (!isdigit(c))
234 {
235 throw mcrl2::runtime_error("Number " + s + " contains symbol '" + c + "' which is not a digit.");
236 }
237 if (c>='0')
238 {
239 multiply_by(10,std::size_t(c-'0'));
240 }
241 }
243 }
244
248 bool is_zero() const
249 {
251 return m_number.size()==0;
252 }
253
257 bool is_number(std::size_t n) const
258 {
260 if (n==0)
261 {
262 return m_number.size()==0;
263 }
264 return m_number.size()==1 && m_number.front()==n;
265 }
266
270 void clear()
271 {
273 m_number.clear();
274 }
275
279 explicit operator std::size_t() const
280 {
282 if (m_number.size()>1)
283 {
284 throw mcrl2::runtime_error("It is not possible to transform a big natural number into a machine size number if it does not fit.");
285 }
286 if (m_number.size()==0)
287 {
288 return 0;
289 }
290 return m_number.front();
291 }
292
293 /* \brief Standard comparison operator.
294 */
295 bool operator==(const big_natural_number& other) const
296 {
298 other.is_well_defined();
299 return m_number==other.m_number;
300 }
301
302 /* \brief Standard comparison operator.
303 */
304 bool operator!=(const big_natural_number& other) const
305 {
306 return !this->operator==(other);
307 }
308
309 /* \brief Standard comparison operator.
310 */
311 bool operator<(const big_natural_number& other) const
312 {
314 other.is_well_defined();
315 if (m_number.size()<other.m_number.size())
316 {
317 return true;
318 }
319 if (m_number.size()>other.m_number.size())
320 {
321 return false;
322 }
323 assert(m_number.size()==other.m_number.size());
324 std::vector<std::size_t>::const_reverse_iterator j=other.m_number.rbegin();
325 for(std::vector<std::size_t>::const_reverse_iterator i=m_number.rbegin(); i!=m_number.rend(); ++i, ++j)
326 {
327 if (*i < *j)
328 {
329 return true;
330 }
331 if (*i > *j)
332 {
333 return false;
334 }
335 }
336 // The numbers are equal.
337 assert(m_number==other.m_number);
338 return false;
339 }
340
341 /* \brief Standard comparison operator.
342 */
343 bool operator<=(const big_natural_number& other) const
344 {
345 return !this->operator>(other);
346 }
347
348 /* \brief Standard comparison operator.
349 */
350 bool operator>(const big_natural_number& other) const
351 {
352 return other.operator<(*this);
353 }
354
355 /* \brief Standard comparison operator.
356 */
357 bool operator>=(const big_natural_number& other) const
358 {
359 return other.operator<=(*this);
360 }
361
362 /* Divide the current number by n. If there is a remainder return it. */
363 std::size_t divide_by(std::size_t n)
364 {
365 std::size_t remainder=0;
366 for(std::vector<std::size_t>::reverse_iterator i=m_number.rbegin(); i!=m_number.rend(); ++i)
367 {
368 *i=detail::divide_single_number(*i,n,remainder);
369 }
372 return remainder;
373 }
374
375 // Add the argument to this big natural number
377 {
379 other.is_well_defined();
380
381 // big_natural_number result;
382 // result.m_number.reserve((std::max)(m_number.size(),other.m_number.size()));
383 std::size_t carry=0;
384 for(std::size_t i=0; i < (std::max)(m_number.size(),other.m_number.size()); ++i)
385 {
386 if (i>=m_number.size())
387 {
389 }
390 else if (i>=other.m_number.size())
391 {
393 }
394 else
395 {
397 }
398 }
399 if (carry>0)
400 {
401 m_number.push_back(carry);
402 }
404 }
405
406
407 /* \brief Standard addition operator.
408 */
410 {
412 other.is_well_defined();
413
414 big_natural_number result=*this;
416 return result;
417 }
418
419 /* \brief Standard subtraction.
420 \detail Subtract other from this number. Throws an exception if
421 the result is negative and cannot be represented.
422 */
423 void subtract(const big_natural_number& other)
424 {
426 other.is_well_defined();
427
428 assert(m_number.size()==0 || m_number.back()!=0);
429 assert(other.m_number.size()==0 || other.m_number.back()!=0);
430 if (m_number.size()<other.m_number.size())
431 {
432 throw mcrl2::runtime_error("Subtracting numbers " + pp(*this) + " and " + pp(other) + " yields a non representable negative result (1).");
433 }
434 std::size_t carry=0;
435 for(std::size_t i=0; i<m_number.size(); ++i)
436 {
437 if (i>=other.m_number.size())
438 {
440 }
441 else
442 {
444 }
445 }
446 if (carry>0)
447 {
448 throw mcrl2::runtime_error("Subtracting numbers " + pp(*this) + " and " + pp(other) + " yields a non representable negative result (2).");
449 }
452 }
453
454 /* \brief Standard subtraction operator. Throws an exception if the result
455 * is negative and cannot be represented.
456 */
458 {
460 other.is_well_defined();
461 big_natural_number result=*this;
462 result.subtract(other);
463 return result;
464 }
465
466
467
468 /* \brief Efficient multiplication operator that does not declare auxiliary vectors.
469 \detail Initially result must be zero. At the end: result equals (*this)*other+result.
470 The calculation_buffer does not need to be initialised.
471 */
472 void multiply(const big_natural_number& other,
473 big_natural_number& result,
474 big_natural_number& calculation_buffer_for_multiplicand) const
475 {
477 other.is_well_defined();
478 std::size_t offset=0;
479 for(std::size_t digit: other.m_number)
480 {
481 // Move n2 to the multiplicand and multiply it with base^offset.
482 calculation_buffer_for_multiplicand.clear();
483 for(std::size_t i=0; i< offset; ++i) { calculation_buffer_for_multiplicand.m_number.push_back(0); }
484 for(std::size_t n: m_number)
485 {
486 calculation_buffer_for_multiplicand.m_number.push_back(n);
487 }
488
489 // Multiply the multiplicand with digit.
490 calculation_buffer_for_multiplicand.multiply_by(digit,0);
491 // Add the multiplicand to the result.
493 offset++;
494 }
495 result.is_well_defined();
496 }
497
498 /* \brief Standard multiplication operator. This is not very efficient as it declares two temporary vectors.
499 */
501 {
502 big_natural_number result, buffer;
503 multiply(other,result,buffer);
504 return result;
505 }
506
507 // This is an auxiliary function getting the n-th digit,
508 // where digit 0 is the least significant one.
509 // It is not necessary that n is in range.
510 /* std::size_t getdigit(const std::size_t n) const
511 {
512 if (n>=m_number.size())
513 {
514 return 0;
515 }
516 return m_number[n];
517 } */
518
519 /* \brief Efficient divide operator that does not declare auxiliary vectors.
520 \detail Initially result must be zero. At the end: result equals (*this)*other+result.
521 The calculation_buffer does not need to be initialised.
522 The algorithm uses standard "primary school" division, except that
523 the digits in this case are 64 bits numbers. The calculation is tricky
524 as the most significant digit of this may be a remaining digit from
525 the previous calculation.
526 */
527 void div_mod(const big_natural_number& other,
528 big_natural_number& result,
529 big_natural_number& remainder,
530 big_natural_number& calculation_buffer_divisor) const
531 {
533 other.is_well_defined();
534 assert(!other.is_zero());
535
536 if (m_number.size()==1 && other.m_number.size()==1)
537 {
538 std::size_t n=m_number.front()/other.m_number.front(); // Calculate div.
539 if (n==0)
540 {
541 result.clear();
542 }
543 else
544 {
545 result.m_number.resize(1);
546 result.m_number[0]=n;
547 }
548 n=m_number.front() % other.m_number.front(); // Calculate mod.
549 if (n==0)
550 {
551 remainder.clear();
552 }
553 else
554 {
555 remainder.m_number.resize(1);
556 remainder.m_number[0]=n;
557 }
558 result.is_well_defined();
559 remainder.is_well_defined();
560 return;
561 }
562
563 // TODO: The procedure below works bitwise, as no efficient division algorithm has yet
564 // been implemented. A natural candidate is the algorithm in "Per Brinch Hansen, Multiple
565 // length division revisited: A tour of the minefield. Software practice and experience 24,
566 // 579-601, 1994.
567 result.clear();
568 remainder=*this;
569
570 if (m_number.size()<other.m_number.size())
571 {
572 result.is_well_defined();
573 remainder.is_well_defined();
574 return;
575 }
576 const int no_of_bits_per_digit=std::numeric_limits<std::size_t>::digits;
577
578 big_natural_number divisor;
579 calculation_buffer_divisor.clear();
580 for(std::size_t i=0; i< (1+m_number.size())-other.m_number.size(); ++i) { calculation_buffer_divisor.m_number.push_back(0); }
581
582 // Place 0 digits at least significant position of the calculation_buffer_divisor to make it of comparable length as the remainder.
583 for(std::size_t i: other.m_number)
584 {
585 calculation_buffer_divisor.m_number.push_back(i);
586 }
587 calculation_buffer_divisor.remove_significant_digits_that_are_zero();
588
589 for(std::size_t i=0; i<=no_of_bits_per_digit*(m_number.size()-other.m_number.size()+1); ++i)
590 {
591 if (remainder<calculation_buffer_divisor)
592 {
593 // We cannot subtract the calculation_buffer_divisor from the remainder.
594 result.multiply_by(2,0); // result=result*2
595 }
596 else
597 {
598 // We subtract the calculation_buffer_divisor from the remainder.
599 result.multiply_by(2,1); // result=result*2 + 1
600 remainder.subtract(calculation_buffer_divisor);
601 }
602 calculation_buffer_divisor.divide_by(2); // Shift the calculation_buffer_divisor one bit to the left.
603 }
604
606 result.is_well_defined();
607 remainder.is_well_defined();
608 }
609
610 /* void div_mod(const big_natural_number& other,
611 big_natural_number& result,
612 big_natural_number& remainder,
613 big_natural_number& calculation_buffer_subtractor) const
614 {
615 is_well_defined();
616 other.is_well_defined();
617 assert(!other.is_zero());
618
619 result.clear();
620 remainder=*this;
621
622 if (m_number.size()<other.m_number.size())
623 {
624 return;
625 }
626this->print_number("div_mod: this: ");
627other.print_number("div_mod: other: ");
628 std::size_t remaining_digit=0;
629 for(std::size_t n=remainder.m_number.size()-1; n>=other.m_number.size(); n--)
630 {
631std::cerr << "Wat is n " << n << "\n";
632 std::size_t r=remainder.getdigit(n+1);
633 std::size_t divisor=detail::divide_single_number(
634 remainder.getdigit(n),
635 other.m_number.back(),
636 r);
637
638result.print_number("div_mod: result: ");
639remainder.print_number("div_mod: remainder: ");
640std::cerr << "dmwhile: divisor: " << divisor <<"\n";
641calculation_buffer_subtractor.is_well_defined();
642 calculation_buffer_subtractor.m_number=std::vector<std::size_t>(n-other.m_number.size(),0); Inefficient, want constructie
643 for(std::size_t i: other.m_number)
644 {
645 calculation_buffer_subtractor.m_number.push_back(i);
646 }
647 calculation_buffer_subtractor.multiply_by(divisor,0);
648calculation_buffer_subtractor.print_number("div_mod: subtractor: ");
649 if (remainder<calculation_buffer_subtractor)
650 {
651 divisor=divisor-1;
652std::cerr << "dmwhile: divisor adapted: " << divisor <<"\n";
653 calculation_buffer_subtractor.m_number=std::vector<std::size_t>(n-other.m_number.size(),0);
654 for(std::size_t i: other.m_number)
655 {
656 calculation_buffer_subtractor.m_number.push_back(i);
657 }
658 calculation_buffer_subtractor.multiply_by(divisor,0);
659 }
660 result.m_number.push_back(divisor);
661remainder.print_number("remainder before assert");
662calculation_buffer_subtractor.print_number("subtractor before assert");
663 assert(remainder>=calculation_buffer_subtractor);
664 std::size_t old_remainder_size=remainder.m_number.size();
665 remainder.subtract(calculation_buffer_subtractor);
666 remaining_digit=(remaining_digit?
667 remainder.m_number.size()+1==old_remainder_size:
668 remainder.m_number.size()==old_remainder_size);
669 }
670result.print_number("div_mod: result voor swap: ");
671 if (result.m_number.size()==0)
672 {
673 return;
674 }
675 // Result must be reverted.
676 std::size_t begin = 0;
677 std::size_t end = result.m_number.size()-1;
678 while (begin<end)
679 {
680 std::swap(result.m_number[begin],result.m_number[end]);
681 begin++;
682 end--;
683 }
684 result.remove_significant_digits_that_are_zero();
685result.print_number("div_mod: result na swap: ");
686 result.is_well_defined();
687 remainder.is_well_defined();
688 } */
689
690 /* \brief Standard division operator. This is currently implemented by a bit wise subtraction. Can be optimized by a 64 bit calculation.
691 \detail. This routine is not particularly efficient as it declares three temporary vectors.
692 */
694 {
695 // Division by zero is not allowed.
696 if (other.is_zero())
697 {
698 throw mcrl2::runtime_error("Division by zero.");
699 }
700 // Zero divided by something is zero.
701 if (is_zero())
702 {
703 return *this;
704 }
705
706 // Often numbers only consist of one digit. Deal with this using machine division.
707 if (m_number.size()==1 && other.m_number.size()==1)
708 {
709 return big_natural_number(m_number.front()/other.m_number.front());
710 }
711
712 // Otherwise do a multiple digit division.
713 big_natural_number result, remainder, buffer;
714 div_mod(other,result,remainder,buffer);
715 return result;
716 }
717
718 /* \brief Standard modulo operator. This is currently implemented by a bit wise subtraction. Can be optimized by a 64 bit calculation.
719 \detail. This routine is not particularly efficient as it declares three temporary vectors.
720 */
722 {
723 // Modulo zero is yields the value itself.
724 // Zero modulo something is zero.
725 if (other.is_zero() || is_zero())
726 {
727 return *this;
728 }
729
730 // Often numbers only consist of one digit. Deal with this using machine division.
731 if (m_number.size()==1 && other.m_number.size()==1)
732 {
733 return big_natural_number(m_number.front()%other.m_number.front());
734 }
735
736 big_natural_number result, remainder, buffer;
737 div_mod(other,result,remainder,buffer);
738 return remainder;
739
740 }
741};
742
743inline std::ostream& operator<<(std::ostream& ss, const big_natural_number& l)
744{
745 static big_natural_number n; // This code is not re-entrant, but it avoids declaring a vector continuously.
746 n=l;
747 std::string s; // This string contains the number in reverse ordering.
748 for( ; !n.is_zero() ; )
749 {
750 std::size_t remainder = n.divide_by(10); /* This divides n by 10. */
751 s.push_back(static_cast<char>('0'+remainder));
752 }
753 if (s.empty())
754 {
755 ss << "0";
756 }
757 else
758 {
759 for(std::string::const_reverse_iterator i=s.rbegin(); i!=s.rend(); ++i)
760 {
761 ss << *i;
762 }
763 }
764 return ss;
765}
766
767
771{
772 x.m_number.swap(y.m_number);
773}
774
775} // namespace utilities
776} // namespace mcrl2
777
778namespace std
779{
780
781template <>
782struct hash< mcrl2::utilities::big_natural_number >
783{
785 {
786 hash<std::vector<std::size_t> > hasher;
787 return hasher(n.m_number);
788 }
789};
790
791
792} // namespace std
793
794namespace mcrl2
795{
796namespace utilities
797{
798/* \brief A pretty print operator on action labels, returning it as a string.
799*/
800inline std::string pp(const big_natural_number& l)
801{
802 std::stringstream s;
803 s << l;
804 return s.str();
805}
806} // namespace utilities
807} // namespace mcrl2
808
809
810#endif // MCRL2_UTILITIES_BIG_NUMBERS_H
811
812
Standard exception class for reporting runtime errors.
Definition exception.h:27
void div_mod(const big_natural_number &other, big_natural_number &result, big_natural_number &remainder, big_natural_number &calculation_buffer_divisor) const
bool is_number(std::size_t n) const
Returns whether this number equals a number of std::size_t.
big_natural_number(const std::size_t n)
big_natural_number operator*(const big_natural_number &other) const
big_natural_number(const std::string &s)
bool operator>=(const big_natural_number &other) const
void print_number(const std::string &s) const
std::vector< std::size_t > m_number
big_natural_number operator+(const big_natural_number &other) const
bool operator==(const big_natural_number &other) const
std::size_t divide_by(std::size_t n)
big_natural_number operator%(const big_natural_number &other) const
big_natural_number operator-(const big_natural_number &other) const
friend void swap(big_natural_number &x, big_natural_number &y)
bool operator<=(const big_natural_number &other) const
big_natural_number operator/(const big_natural_number &other) const
void clear()
Sets the number to zero.
bool operator!=(const big_natural_number &other) const
bool is_zero() const
Returns whether this number equals zero.
void multiply(const big_natural_number &other, big_natural_number &result, big_natural_number &calculation_buffer_for_multiplicand) const
bool operator<(const big_natural_number &other) const
void subtract(const big_natural_number &other)
bool operator>(const big_natural_number &other) const
void multiply_by(std::size_t n, std::size_t carry)
Exception classes for use in libraries and tools.
This file contains a specialisation for hashes on pairs. This is not a part of the standard,...
std::size_t multiply_single_number(const std::size_t n1, const std::size_t n2, std::size_t &multiplication_carry)
Definition big_numbers.h:95
std::size_t divide_single_number(const std::size_t p, const std::size_t q, std::size_t &remainder)
std::size_t add_single_number(const std::size_t n1, const std::size_t n2, std::size_t &carry)
Definition big_numbers.h:49
std::size_t subtract_single_number(const std::size_t n1, const std::size_t n2, std::size_t &carry)
Definition big_numbers.h:72
void swap(big_natural_number &x, big_natural_number &y)