!standard A.20(0) 16-12-19 AI12-0208-1/00 !class Amendment 16-12-19 !status work item 16-12-19 !status received 16-09-27 !priority Low !difficulty Medium !subject Predefined bignum support !summary Define a "bignum" package. [Editor's note: I surely hope that we have a better name for it than "bignum", which is as non-Ada a name as possible (two words with no underscore and an unnecessary abbreviation.] !problem Some applications need larger numbers than Standard.Integer. All Ada compilers have this capability in order to implement static expressions; shouldn't some such package be available to Ada users as well? (Yes.) !proposal (See Summary.) !wording ** TBD **. !discussion !ASIS No ASIS effect (assuming this is ONLY a library). !ACATS test An ACATS C-Test is needed to check that the new capabilities are supported. !appendix From: Steve Baird Sent: Tuesday, September 27, 2016 4:09 PM professor at U. of Utah: blog.regehr.org/archives/1401 Regehr says: In most programming languages, the default integer type should be a bignum: an arbitrary-precision integer that allocates more space when needed. Efficient bignum libraries exist and most integers never end up needing more than one machine word anyway, except in domains like crypto. Nobody is suggesting changing how Standard.Integer works for Ada, but a language-defined Bignum package (presumably supporting Rationals as well as Integers) would be a step in the right direction. It seems like the same arguments which were used (correctly, IMO) to justify adding predefined container packages to the language also apply here. As Tuck phrased it in a private message: portability and more capability "out of the box." Does some de facto standard already exist? **************************************************************** From: Bob Duff Sent: Tuesday, September 27, 2016 4:32 PM > Nobody is suggesting changing how Standard.Integer works But somebody might suggest that things like "type T is range 1..10**100;" should be supported by all Ada compilers. > It seems like the same arguments which were used (correctly, IMO) to > justify adding predefined container packages to the language also > apply here. As Tuck phrased it in a private message: > portability and more capability "out of the box." Plus the fact that all Ada compilers have to support that functionality at compile time, but can't provide it to their users in a portable way at run time. > Does some de facto standard already exist? For C and C++, yes. For Ada, no. For Common Lisp, Java, C#, and many others, a de jure standard exists. **************************************************************** From: Randy Brukardt Sent: Wednesday, September 28, 2016 12:49 PM > Does some de facto standard already exist? No. I could be convinced to contribute RR's Univmath package as a starting point for discussion. **************************************************************** From: Jean-Pierre Rosen Sent: Thursday, September 29, 2016 12:22 AM There are several packages available, see http://bignumber.chez.com/index.html **************************************************************** From: Randy Brukardt Sent: Thursday, September 29, 2016 12:28 PM Surely, like containers there are as many Bignum packages as there are Ada programmers (much like containers - everybody has one). But is someone putting them into RM format?? That's what it means to "contribute" a package here. **************************************************************** From: John Barnes Sent: Thursday, September 29, 2016 2:05 PM I see there has been chatter on big number packages. I wrote such a package many years ago. I was intending to write a book called Fun with Ada using big examples of Ada 83 programs. But it got overtaken by events such as having to write the book on Ada 95. But I kept the package, used some child stuff from Ada 95 but otherwise left it alone, I still use it for dabbling with large prime numbers and so on. I think it is based on base 10,000 which will run on a 16 bit machine and is easy for conversion for printing. But I fear that agreeing on something might be tricky. **************************************************************** From: Florian Schanda Sent: Friday, September 30, 2016 2:35 AM > But I kept the package, used some child stuff from Ada 95 but > otherwise left it alone, I still use it for dabbling with large prime > numbers and so on. I think it is based on base 10,000 which will run > on a 16 bit machine and is easy for conversion for printing. Generally, these days, you would probably want to stick to largest power-of- two as printing these is not a massive concern but performance is. :) Anyway, I think whatever we come up with, it should be possible to implement it via a binding to GMP [https://gmplib.org] which is more or less the gold standard for arbitrary precision arithmetic. Of course, some runtime may wish to have a more verifiable implementation... So, I think there are two requirements we should make sure to fulfil: 1. the api should be amenable to static analysis and formal verification 2. the api should make it easy to bind to gmp (Not saying this list is exhaustive.) I just want to avoid people starting from various in-house and private projects; its probably a good idea instead to start from established libraries. **************************************************************** From: Steve Baird Sent: Friday, September 30, 2016 12:27 PM > So, I think there are two > requirements we should make sure to fulfil: > > 1. the api should be amenable to static analysis and formal verification > 2. the api should make it easy to bind to gm It is also at least possible that we'll want something similar to what we have with the containers, where we have one version for use in situations where controlled types and dynamic storage allocation are ok and another for use in other situations. **************************************************************** From: Jean-Pierre Rosen Sent: Friday, September 30, 2016 2:44 PM Hmmm... bounded and unbounded bignums? **************************************************************** From: Tucker Taft Sent: Friday, September 30, 2016 3:52 PM Perhaps: "big enough nums, already..." **************************************************************** From: Steve Baird Sent: Tuesday, December 12, 2017 7:20 PM I thought I'd take a look at how Java and C++ do bignums to see if there are any ideas there worth incorporating. My going-in idea is to have two packages with similar specs; one has "Capacity" discriminants and the other is implemented using dynamic storage allocation of some sort (e.g., controlled types and allocators). Like the bounded/unbounded versions of the containers. C++ doesn't really have a standard for bignums, but the GCC/GMP stuff looks pretty similar to what I expected. Java, however, surprised me (note that I am far from a Java expert so it could be that I am just confused here). The Java big-real spec doesn't have Numerator and Denominator functions which yield big-ints. The Java type seems to be named BigDecimal. BigDecimal is implemented as a single big-int value accompanied by two ints (Scale and Precision), at least according to stackoverflow.com/questions/10655254/how-bigdecimal-is-implemented Which leads to my question: If Ada defined a spec where the intended implementation for bigreal is clearly two bigints (one for numerator, one for denominator), would this result in lots of "I coded up the same algorithm in Ada and Java and performance was a lot worse in Ada" horror stories? Apparently BigDecimal lets you have, in effect, a lot of decimal digits but the value "one third" still cannot be represented exactly. Why did the Java folks do it that way? It seems like you lose a lot of value if you can't exactly represent, for example, one third. But perhaps most folks don't care about that functionality and the performance/functionality tradeoff chosen by Java is closer to what most folks want. Opinions? Opinions about Java are of some interest, but what I really want is opinions about what we should do in Ada. p.s. Note that the current plan for this AI is to add one or more new predefined packages but no changes to language rules. In particular, numeric literals for a non-numeric type is the topic of another AI. **************************************************************** From: Tucker Taft Sent: Wednesday, December 13, 2017 9:15 AM We want rational, not decimal, computations, I believe. So I would ignore Java's BigDecimal. A different and interesting capability is true "real" arithmetic, which works for transcendentals, etc. It is intriguing, but probably not what people really want. I'll send the PDF for an article by Hans Boehm about "real" arithmetic separately, since it will probably not make it through the SPAM filter! **************************************************************** From: Randy Brukardt Sent: Wednesday, December 13, 2017 10:57 AM > We want rational, not decimal, computations, I believe. So I would > ignore Java's BigDecimal. Isn't that Steve's question? Ada compiler vendors use rational computations since that is required by the ACATS (it's necessary that 1/3 /= 0.33333333333333333333333333333). But is that the best choice for the Ada community? I don't know. > A different and interesting capability is true "real" > arithmetic, which works for transcendentals, etc. It is intriguing, > but probably not what people really want. > > I'll send the PDF for an article by Hans Boehm about "real" > arithmetic separately, since it will probably not make it through the > SPAM filter! It might be too large for the list as well. If so, I can post it in the Grab Bag if you send it directly to me. **************************************************************** From: Edmond Schonberg Sent: Wednesday, December 13, 2017 1:36 PM > I'll send the PDF for an article by Hans Boehm about "real" arithmetic > separately, since it will probably not make it through the SPAM filter! Both the rational representation and Boehm’s approach require arbitrary precision integer arithmetic, so the spec of that new package is straightforward. The papers describing the implementation of Boehm’s approach claim that it is much more efficient than working on rationals, where numerator and denominator grow very rapidly, while the other method only computes required bits. I have no idea whether numerical analysts use this method, From the literature it seems to be of interest to number theorists. **************************************************************** From: Randy Brukardt Sent: Wednesday, December 13, 2017 4:49 PM > > It might be too large for the list as well. If so, I can post it in > > the Grab Bag if you send it directly to me. > > It was too big. By about 4 Megabytes. :-) Since the article is copyrighted, I put it in the private part of the website. Find it at: http://www.ada-auth.org/standards/private/real_arithmetic-boehm.pdf Ed suggested privately: > Additional details on the underlying model and its implementation in: > http://keithbriggs.info/documents/xr-paper2.pdf **************************************************************** From: Randy Brukardt Sent: Wednesday, December 13, 2017 5:10 PM ... > Both the rational representation and Boehm's approach require > arbitrary precision integer arithmetic, so the spec of that new > package is straightforward. > The papers describing the implementation of Boehm's approach claim > that it is much more efficient than working on rationals, where > numerator and denominator grow very rapidly, while the other method > only computes required bits. I have no idea whether numerical analysts > use this method, From the literature it seems to be of interest to > number theorists. The problem with the Boehm method is that it requires specifying those "required bits", which seems problematic in a programming environment. One could do it with a form of type declaration (Ada's "digits" seems to be the right general idea), but that doesn't make much sense in a library form. Boehm's actual use gets those on the fly (by interacting with the user), and they also use a rational representation as a backup. So it seems that a rational representation is going to show up somewhere. I've repeatedly had the fever dream of a class-wide bignum base type, something like: package Root_Bignum is type Root_Bignum_Type is abstract tagged null record; function Bignum (Val : in Long_Float) return Root_Bignum_Type is abstract; -- To get the effect of literals and conversions. function "+" (Left, Right : in Root_Bignum_Type) return Root_Bignum_Type is abstract; -- And all of the rest. function Expected_Digits return Natural is abstract; -- The number of digits supported by this type; 0 is returned -- if the number is essentially infinite. -- And probably some other queries. end Root_Bignum; And then there could be multiple specific implementations with different performance characteristics, everything from Long_Float itself thru infinite rational representations. This would allow one to create most of the interesting algorithms as class-wide operations (with implementations that could adjust to the characteristics of the underlying type), for instance: function Sqrt (Val : in Root_Bignum_Type'Class; Required_Digits : Natural := 0) return Root_Bignum_Type'Class; Here with "Required_Digits" specifies how many digits of result are needed. If 0, the value would be retrieved from the underlying representation. (Probably have to raise an exception if that gives "infinite".) Such a layout would also allow easy changing of representations, which probably would be needed for tuning purposes (most of these maths being slow, at least by conventional standards). This would have the clear advantage of avoiding being locked into a single form of Bignum math, when clearly there are other choices out there useful for particular purposes. **************************************************************** From: Randy Brukardt Sent: Wednesday, December 13, 2017 5:25 PM I said: ... > The problem with the Boehm method is that it requires specifying those > "required bits", which seems problematic in a programming environment. but then also noted: > function Sqrt (Val : in Root_Bignum_Type'Class; > Required_Digits : Natural := 0) return > Root_Bignum_Type'Class; essentially, recognizing that many non-terminating algorithms have to have some sort of termination criteria. For ease of use purposes, one would prefer to only specify numbers of digits if they're really needed (as in Sqrt or PI, etc.). But if there are going to be a lot of such operations, one would want to be able to specify that once. Hope that explains my thinking here. Also, a Bignum library needs a corresponding Text_IO library. And probably a custom version of GEF. (The Janus/Ada compiler library has most of these features, and they are used extensively.) **************************************************************** From: Steve Baird Sent: Friday, January 19, 2018 2:36 PM We have agreed that we want bignum support in the form of one or more predefined packages with no other language extensions (e.g., no new rules for numeric literals) as part of this AI. The general approach seems fairly clear, although there are a lot of details to decide (not the least of which are the choices for names). I think we want two forms, "vanilla" and "bounded" (analogous to, for example, Ada.Containers.Vectors and Ada.Containers.Bounded_Vectors). In one form, the two "big" numeric types (tentatively named Big_Integer and Big_Rational) are defined as undiscriminated types. In the second form, these types are discriminated with some sort of a capacity discriminant. The idea is that the first form is allowed to use dynamic storage allocation and controlled types in its implementation while the second form is not; the discriminant somehow indicates the set of representable values via some mapping (should this mapping be implementation dependent?). At a high level, we might have something like package Ada.Big_Numbers is -- empty spec like Ada.Containers package end; package Ada.Big_Numbers.Big_Integers is type Big_Integer is private; function GCD (Left, Right : Big_Integer) return Integer; function "+" (Arg : Some_Concrete_Integer_Type_TBD) return Big_Integer; ... ops for Big_Integer ... end Ada.Big_Numbers.Big_Integers. with Ada.Big_Numbers.Big_Integers; package Ada.Big_Numbers.Big_Rationals is use type Big_Integers.Big_Integer; type Big_Rational is private with Type_Invariant => Big_Rational = +0 or else Big_Integers.GCD (Big_Integers.Numerator (Big_Rational), Big_Integers.Denominator (Big_Rational)) = +1; function Numerator (Arg : Big_Rational) return Big_Integer; function Denominator (Arg : Big_Rational) return Big_Integer; function "/" (Num, Den : Big_Integer) return Big_Rational with Pre => Den /= +0; ... other ops for Big_Rational ... end Ada.Big_Numbers.Big_Rationals; package Ada.Big_Numbers.Bounded_Big_Integers is ... end; package Ada.Big_Numbers.Bounded_Big_Rationals is ... end; Questions/observations include: 1) Do we declare deferred constants, parameterless functions, or neither for things like Zero, One, and Two? 2) Which ops do we include? It seems obvious that we define at least the arithmetic and relational ops that are defined for any predefined integer (respectively float) type for Big_Integer (respectively, Big_Rational). What Pre/Postconditions are specified for these ops? These might involve subtype predicates. For example (suggested by Bob), do we want subtype Nonzero_Integer is Big_Integer with Predicate => Nonzero_Integer /= Zero; function "/" (X: Big_Integer; Y: Nonzero_Integer) return Big_Integer; -- similar for "mod", "rem". ? What other operations should be provided? - Conversion between Big_Int and what concrete integer types? I'd say define a type with range Min_Int .. Max_Int and provide conversion functions for that type. Also provide two generic conversion functions that take a generic formal signed/modular type. - Conversion between Big_Rational and what concrete integer or float types? Same idea. Conversion between a maximal floating point type and then a pair of conversion generics with formal float/fixed parameters. - What shortcuts do we provide (i.e., ops that can easily be built out of other ops)? Assignment procedures like Add (X, Y); -- X := X + Y or mixed-type operators whose only purpose is to spare users from having to write explicit conversion? 3) It seems clear that we don't want the bounded form of either package to "with" the unbounded form but we do want conversion functions for going between corresponding bounded and unbounded types. Perhaps these go in child units of the two bounded packages (those child units could then "with" the corresponding unbounded packages). Should streaming of the two forms be compatible as with vectors and bounded vectors? 4) We need an Assign procedure. In the unbounded case it can be just a wrapper for predefined assignment, but in the bounded case it has to deal with the case where the two arguments have different capacities. It's fairly obvious what to do in most cases, but what about assigning a Big_Rational value which cannot be represented exactly given the capacity of the target. Raise an exception or round? In either case, we probably want to provide a Round function that deterministically finds an approximation to a given value which can be represented as a value having a given capacity. This can be useful in the unbounded case just to save storage. Should this Round function be implementation-dependent? If not, then we might end up talking about convergents and semi-convergents in the Ada RM (or at least in the AARM), which would be somewhat odd (see shreevatsa.wordpress.com/2011/01/10/not-all-best-rational-approximations-are-the-convergents-of-the-continued-fraction ). I do not think we want to define Succ/Pred functions which take a Big_Rational and a capacity value. 5) We want to be sure that a binding to GNU/GMP is straightforward in the unbounded case. [Fortunately, that does not require using the same identifiers used in GNU/GMP (mpz_t and mpq_t).] See gmplib.org/manual for the GNU/GMP interfaces. 6) Do we want functions to describe the mapping between Capacity discriminant values and the associated set of representable values? For example, a function from a value (Big_Integer or Big_Rational) to the smallest capacity value that could be used to represent it. For Big_Integer there could presumably be Min and Max functions that take a capacity argument. For Big_Rational, it's not so clear. We could require, for example, that a given capacity value allows representing a given Big_Rational value if it is >= the sum of the capacity requirements of the Numerator and the Denominator. 7) Bob feels (and I agree) that the ARG should not formally approve any changes until we have experience with an implementation. At this point the ARG should be focused on providing informal guidance on this topic. Opinions? **************************************************************** From: Randy Brukardt Sent: Friday, January 19, 2018 10:18 PM ... > Questions/observations include: 0) Should Big_Integer and (especially) Big_Rational be visibly tagged? If so, then we can use prefix notation on functions like Numerator and Denominator. We could also consider deriving both versions (usual and bounded) from an abstract ancestor. > 1) Do we declare deferred constants, parameterless functions, > or neither for things like Zero, One, and Two? If tagged, I'll finally get an excuse to show why what I called "tag propagation" is necessary to implement the dispatching rules in 3.9.2. :-) (One has to consider a set of calls, not a single call, for determining the static or dynamic tag for dispatching. That's demonstratably necessary to process tagged expressions with constants or literals.) Anyway, the answer to this depends on whether there is a sufficiently short constructor -- and that really depends on whether Tucker invents a useful "literals for private type" AI. So I don't think this can be answered until we find out about that. > 2) Which ops do we include? It seems obvious that we define at least > the arithmetic and relational ops that are defined for any > predefined integer (respectively float) type for Big_Integer > (respectively, Big_Rational). > > What Pre/Postconditions are specified for these ops? > These might involve subtype predicates. > For example (suggested by Bob), do we want > > subtype Nonzero_Integer is Big_Integer with > Predicate => Nonzero_Integer /= Zero; > function "/" > (X: Big_Integer; Y: Nonzero_Integer) return Big_Integer; > -- similar for "mod", "rem". > > ? Shouldn't this predicate raise Constraint_Error rather than defaulting to Assertion_Error, to be more like the other numeric operations? Otherwise, I'm all in favor of this formulation. Note, however, that since the underlying type is likely to be controlled and thus tagged, this would require some changes to other rules; there is already an AI about that (AI12-0243-1). > What other operations should be provided? > - Conversion between Big_Int and what concrete integer types? > I'd say define a type with range Min_Int .. Max_Int > and provide conversion functions for that type. Also provide > two generic conversion functions that take a generic formal > signed/modular type. Sounds OK. > - Conversion between Big_Rational and what concrete integer or > float types? Same idea. Conversion between a maximal > floating point type and then a pair of conversion generics > with formal float/fixed parameters. Sounds OK again. > - What shortcuts do we provide (i.e., ops that can easily be > built out of other ops)? Assignment procedures like > Add (X, Y); -- X := X + Y > or mixed-type operators whose only purpose is to spare users > from having to write explicit conversion? The only reason for mixed type operators is to make literals available. But if one does those, then we can't add literals properly in the future (Ada.Strings.Unbounded is damaged by this). So I say no. I wouldn't bother with any other routines until at least such time as Bob :-) has built some ACATS tests. > 3) It seems clear that we don't want the bounded form of either > package to "with" the unbounded form but we do want conversion > functions for going between corresponding bounded and unbounded > types. Perhaps these go in child units of the two bounded packages > (those child units could then "with" the corresponding unbounded > packages). Alternatively, both could be derived from an abstract type, and a class-wide conversion provided. That would get rid of the empty package in your proposal. :-) > Should streaming of the two forms be compatible as with > vectors and bounded vectors? Yes. > 4) We need an Assign procedure. In the unbounded case it can be just > a wrapper for predefined assignment, but in the bounded case it > has to deal with the case where the two arguments have different > capacities. It's fairly obvious what to do in most cases, but what > about assigning a Big_Rational value which cannot be represented > exactly given the capacity of the target. Raise an exception or > round? I think I'd raise Capacity_Error. (Isn't that what the containers do?) Having exact math be silently non-exact seems like exactly (pun) the wrong thing to do. > In either case, we probably want to provide a Round function > that deterministically finds an approximation to a given > value which can be represented as a value having a given > capacity. This can be useful in the unbounded case just to save > storage. Should this Round function be implementation-dependent? > If not, then we might end up talking about convergents and > semi-convergents in the Ada RM (or at least in the AARM), > which would be somewhat odd (see > shreevatsa.wordpress.com/2011/01/10/not-all-best-rational-appr > oximations-are-the-convergents-of-the-continued-fraction > ). I do not think we want to define Succ/Pred functions which take > a Big_Rational and a capacity value. ??? I don't think Round (or any other operation) ought to be implementation-dependent, so I think it would need a real definition. Hopefully with "semi-convergents" or other terms that no one has heard of. ;-) > 5) We want to be sure that a binding to GNU/GMP is straightforward in > the unbounded case. [Fortunately, that does not require using the > same identifiers used in GNU/GMP (mpz_t and mpq_t).] > See gmplib.org/manual for the GNU/GMP interfaces. Makes sense. > 6) Do we want functions to describe the mapping between Capacity > discriminant values and the associated set of representable values? > For example, a function from a value (Big_Integer or Big_Rational) > to the smallest capacity value that could be used to represent it. > For Big_Integer there could presumably be Min and Max functions > that take a capacity argument. For Big_Rational, it's not so clear. > We could require, for example, that a given capacity value allows > representing a given Big_Rational value if it is >= the sum of > the capacity requirements of the Numerator and the Denominator. It seems that the Capacity needs to mean something to the end user, not just the compiler. So such functions seem necessary, but KISS for those!! > 7) Bob feels (and I agree) that the ARG should not formally approve any > changes until we have experience with an implementation. At this > point the ARG should be focused on providing informal guidance on > this topic. I agree that Bob should prototype these packages, including writing ACATS-style tests for them, so that we can put them into the Ada 2020 Standard. I'll put it on his action item list. ;-) Seriously, we already have an ARG rule that all Amendment AIs are supposed to include (some) ACATS tests, and we really should have a similar rule that proposed packages are prototyped as well. This is the assumed responsibility of an AI author, so if you can't get Bob to help, you're pretty much stuck, and need to do that before the AI could be assumed complete. OTOH, we haven't required that from any other AI author, so why start now?? (We really ought to, I don't have a very big budget to write Ada 2020 ACATS tests. Topic to discuss during the call?) **************************************************************** From: Jean-Pierre Rosen Sent: Saturday, January 20, 2018 12:22 AM > Questions/observations include: > [...] > I'd add: 8) IOs Should an IO package be associated to each of these bignums? Note that the issue of IO may influence the representation of of bignums: I once knew an implementation where each super-digit was limited to 1_000_000_000 (instead of the natural 2_147_483_647), just to avoid terribly inefficient IOs. **************************************************************** From: Tucker Taft Sent: Saturday, January 20, 2018 11:08 AM > ... > >> 1) Do we declare deferred constants, parameterless functions, >> or neither for things like Zero, One, and Two? > > If tagged, I'll finally get an excuse to show why what I called "tag > propagation" is necessary to implement the dispatching rules in 3.9.2. > :-) (One has to consider a set of calls, not a single call, for > determining the static or dynamic tag for dispatching. That's > demonstratably necessary to process tagged expressions with constants > or literals.) I agree that you have to do "tag propagation" to properly handle tag indeterminate calls. Has anyone claimed otherwise? > > Anyway, the answer to this depends on whether there is a sufficiently > short constructor -- and that really depends on whether Tucker invents > a useful "literals for private type" AI. So I don't think this can be > answered until we find out about that. I'm on it. ;-) **************************************************************** From: Randy Brukardt Sent: Saturday, January 20, 2018 7:29 PM > I agree that you have to do "tag propagation" to properly handle tag > indeterminate calls. Has anyone claimed otherwise? Not that I know of, but based on my compiler surveys, no one implements it other than Janus/Ada. Admittedly, I haven't checked this recently. I've long had a tagged Bignum-like package on my ACATS test to-construct list (because one needs usage-orientation for such tests) in order to test this rule. So far as I can tell, the ACATS doesn't currrently test cases like those that arise in Bignum: procedure Something (Val : in out Num'Class) is begin Val := + Zero; -- Zero gets the tag of Val, propagated through "+". declare Org : Num'Class := Val + (- One); -- Org and One get the tag of Val. begin ... end; end Something; I'll probably come up with more realistic-looking expressions for this test, but the idea should be obvious. (I'll have to test both static and dynamic binding, as well as tag indeterminate cases.) **************************************************************** From: John Barnes Sent: Monday, January 22, 2018 5:49 AM I wrote a bignum package in Ada 83 some 30 years ago. I did make some updates to use Ada 95, mainly child packages. I still use it for numerical stuff for courses at Oxford. Notable points perhaps. I did use a power of 10 for the base to ease IO. It was originally on a 16 bit machine. (386 perhaps). It still works on this horrid Windows 10. Not much faster than on my old XP laptop. I don't know what Windows 10 is doing. Obviously playing with itself - ridiculous. I provided constants Zero and One. I didn't think any others were necessary. Others were provided by eg Two: Number := Make-Number(2); I provided a package for subprograms Add, Sub, Mul, Div, Neg, Compare, Length, To_Number, To_Text, To_Integer. And a package for functions +. -, abs, *, / rem, mod, <, <=, >, >=, = And other packages for I/O. Long time ago. Certainly very useful. **************************************************************** From: Steve Baird Sent: Monday, January 22, 2018 12:33 PM > I'd add: > 8) IOs > Should an IO package be associated to each of these bignums? Good question. If we provide conversion functions to and from String then would any further IO support be needed? **************************************************************** From: Steve Baird Sent: Monday, January 22, 2018 1:24 PM > ... >> Questions/observations include: > > 0) Should Big_Integer and (especially) Big_Rational be visibly tagged? > > If so, then we can use prefix notation on functions like Numerator and > Denominator. We could also consider deriving both versions (usual and > bounded) from an abstract ancestor. If we go this way, then should this common ancestor be an interface type? I'd say yes. Does it then get all the same ops, so that the non-abstract ops declared for the Bounded and Unbounded types would all be overriding? Would this make the AI12-0243-ish issues any worse (consider the proposed Nonzero_Integer parameter subtype mentioned earlier)? I know these problems are bad enough already, but my question is whether this would make matters any worse. >> 2) Which ops do we include? It seems obvious that we define at least >> the arithmetic and relational ops that are defined for any >> predefined integer (respectively float) type for Big_Integer >> (respectively, Big_Rational). >> >> What Pre/Postconditions are specified for these ops? >> These might involve subtype predicates. >> For example (suggested by Bob), do we want >> >> subtype Nonzero_Integer is Big_Integer with >> Predicate => Nonzero_Integer /= Zero; >> function "/" >> (X: Big_Integer; Y: Nonzero_Integer) return Big_Integer; >> -- similar for "mod", "rem". >> >> ? > > Shouldn't this predicate raise Constraint_Error rather than defaulting > to Assertion_Error, to be more like the other numeric operations? Good point; I agree. >> 3) It seems clear that we don't want the bounded form of either >> package to "with" the unbounded form but we do want conversion >> functions for going between corresponding bounded and unbounded >> types. Perhaps these go in child units of the two bounded packages >> (those child units could then "with" the corresponding unbounded >> packages). > > Alternatively, both could be derived from an abstract type, and a > class-wide conversion provided. That would get rid of the empty > package in your proposal. :-) Could you provide a more detailed spec? I don't see how this would work, but I suspect that I'm misunderstanding your proposal. >> 4) We need an Assign procedure. In the unbounded case it can be just >> a wrapper for predefined assignment, but in the bounded case it >> has to deal with the case where the two arguments have different >> capacities. It's fairly obvious what to do in most cases, but what >> about assigning a Big_Rational value which cannot be represented >> exactly given the capacity of the target. Raise an exception or >> round? > > I think I'd raise Capacity_Error. (Isn't that what the containers do?) > Having exact math be silently non-exact seems like exactly (pun) the > wrong thing to do. Is it that simple? Suppose somebody wants large rationals (e.g., 2048-bit numerators and denominators) with rounding. It's not that they require exact arithmetic - they just want a lot more range/precision than what you get from Ada's numeric types. It may be that this is an unimportant corner case and you are right to dismiss it; I don't know. >> 6) Do we want functions to describe the mapping between Capacity >> discriminant values and the associated set of representable values? >> For example, a function from a value (Big_Integer or Big_Rational) >> to the smallest capacity value that could be used to represent it. >> For Big_Integer there could presumably be Min and Max functions >> that take a capacity argument. For Big_Rational, it's not so clear. >> We could require, for example, that a given capacity value allows >> representing a given Big_Rational value if it is >= the sum of >> the capacity requirements of the Numerator and the Denominator. > > It seems that the Capacity needs to mean something to the end user, > not just the compiler. So such functions seem necessary, but KISS for those!! Am I right in guessing that you'd like these functions to be portable (as opposed to being implementation-defined)? **************************************************************** From: Randy Brukardt Sent: Monday, January 22, 2018 3:41 PM > > I'd add: > > 8) IOs > > Should an IO package be associated to each of these bignums? > > Good question. > > If we provide conversion functions to and from String then would any > further IO support be needed? We currently have Text_IO nested packages or children for pretty much any type for which it makes sense to have text input-output, despite the fact that every such type has an Image function or the equivalent (To_String for unbounded strings). So I'd rather expect a Ada.Text_IO.BigNum_IO package. If we don't define it now, we will the next time around. (The Janus/Ada UnivMath package has a complete set of Text_IO packages, and they are heavily used. I believe they can output both rational and decimal representation for the universal_real type.) **************************************************************** From: Randy Brukardt Sent: Monday, January 22, 2018 3:36 PM > > Steve Baird writes: > > ... > >> Questions/observations include: > > > > 0) Should Big_Integer and (especially) Big_Rational be visibly tagged? > > > > If so, then we can use prefix notation on functions like Numerator > > and Denominator. We could also consider deriving both versions > > (usual and > > bounded) from an abstract ancestor. > > If we go this way, then should this common ancestor be an interface > type? I'd say yes. I suggested making it abstract so it could have some concrete operations if those made sense. But perhaps they don't make sense. > Does it then get all the same ops, so that the non-abstract ops > declared for the Bounded and Unbounded types would all be overriding? I would expect that the vast majority of operations are in the interface, so dispatching can be used, and one can write class-wide algorithms that work with any Bignum representation. Probably the capacity-specific operations would be left out. > Would this make the AI12-0243-ish issues any worse (consider the > proposed Nonzero_Integer parameter subtype mentioned earlier)? I know > these problems are bad enough already, but my question is whether this > would make matters any worse. It just makes a solution more urgent, but it doesn't change the issues any. ... > >> 3) It seems clear that we don't want the bounded form of either > >> package to "with" the unbounded form but we do want conversion > >> functions for going between corresponding bounded and unbounded > >> types. Perhaps these go in child units of the two bounded packages > >> (those child units could then "with" the corresponding unbounded > >> packages). > > > > Alternatively, both could be derived from an abstract type, and a > > class-wide conversion provided. That would get rid of the empty > > package in your proposal. :-) > > Could you provide a more detailed spec? I don't see how this would > work, but I suspect that I'm misunderstanding your proposal. I was thinking about including cross-cut operations in the spec, something like: type BigNum is abstract tagged with private; function Convert (Val : in Bignum'Class) return Bignum; but thinking about it now, I can't figure out how one would implement one of those. You'd probably have to have a concrete universal representation to make that work: function Convert (Val : in Bignum) return Universal_Big; function Convert (Val : in Universal_Big) return BigNum; but of course that would bring in the memory allocation/finalization issues that you are trying to avoid. So at this moment I'm thinking that direct conversions would have to be left out; you could generally do it through intermediary types like Max_Integer using Numerator/Demomonator. > >> 4) We need an Assign procedure. In the unbounded case it can be just > >> a wrapper for predefined assignment, but in the bounded case it > >> has to deal with the case where the two arguments have different > >> capacities. It's fairly obvious what to do in most cases, but what > >> about assigning a Big_Rational value which cannot be represented > >> exactly given the capacity of the target. Raise an exception or > >> round? > > > > I think I'd raise Capacity_Error. (Isn't that what the containers > > do?) Having exact math be silently non-exact seems like exactly > > (pun) the wrong thing to do. > > Is it that simple? Suppose somebody wants large rationals (e.g., > 2048-bit numerators and denominators) with rounding. > It's not that they require exact arithmetic - they just want a lot > more range/precision than what you get from Ada's numeric types. > It may be that this is an unimportant corner case and you are right to > dismiss it; I don't know. We're not trying to be all things to all people. I'd consider these "exact" math packages and treat them accordingly. If there is an abstract root, one can "easily" make a clone version that uses rounding if someone needs that. (Defining the rounding is hard, as you noted elsewhere.) > >> 6) Do we want functions to describe the mapping between Capacity > >> discriminant values and the associated set of representable values? > >> For example, a function from a value (Big_Integer or Big_Rational) > >> to the smallest capacity value that could be used to represent it. > >> For Big_Integer there could presumably be Min and Max functions > >> that take a capacity argument. For Big_Rational, it's not so clear. > >> We could require, for example, that a given capacity value allows > >> representing a given Big_Rational value if it is >= the sum of > >> the capacity requirements of the Numerator and the Denominator. > > > > It seems that the Capacity needs to mean something to the end user, > > not just the compiler. So such functions seem necessary, but KISS > > for those!! > > Am I right in guessing that you'd like these functions to be portable > (as opposed to being implementation-defined)? I think so; otherwise it rather defeats the purpose of language-defined packages (to provide the ultimate in portability). **************************************************************** From: Bob Duff Sent: Sunday, January 28, 2018 11:29 AM > Steve Baird writes: > ... > > Questions/observations include: > > 0) Should Big_Integer and (especially) Big_Rational be visibly tagged? Surely not. I think we want to be competetive (efficiency-wise) with all sorts of other languages, and taggedness will destroy that. Let's not have another "tampering" fiasco. > If so, then we can use prefix notation on functions like Numerator and > Denominator. I'm not a big fan of that feature, but if we want it, we should figure out how to do it for untagged types. >... We could also consider deriving both versions (usual and > bounded) from an abstract ancestor. Consider, ..., and reject. ;-) **************************************************************** From: Jeff Cousins Sent: Sunday, January 28, 2018 12:21 PM John Barnes wrote: I wrote a bignum package in Ada 83 some 30 years ago Would you be able to let us see the spec for this? **************************************************************** From: Randy Brukardt Sent: Sunday, January 28, 2018 9:13 PM > > 0) Should Big_Integer and (especially) Big_Rational be > visibly tagged? > > Surely not. I think we want to be competetive > (efficiency-wise) with all sorts of other languages, and taggedness > will destroy that. ??? Tags (as opposed to controlled types) add almost no overhead, especially in a case like unbounded Bignum which probably will have to be controlled anyway. (The only overhead of a tagged type is initializing the tag in the object.) So long as one uses a single specific type, everything is statically bound and the cost is essentially the same as an untagged type (again, especially as the underlying specific type most likely will be tagged and certainly will be large). I wasn't suggesting that we define any class-wide operations other than representation conversion (which should be rarely used in any case). Class-wide operations are the only operations that add overhead. > Let's not have another "tampering" fiasco. I'm still waiting for an example program showing this supposed "fiasco". No one has ever submitted one to the ARG. We've essentially been asked to believe this issue by repeated assertion. (And most tampering checks can be done at compile-time, with sufficient will.) If there was a fiasco here, it was that the goals of the containers did not include making them particularly fast. If they are then misused for high-performance code, one is going to get the expected disappointment. Perhaps we started with the wrong set of goals. > > If so, then we can use prefix notation on functions like Numerator > > and Denominator. > > I'm not a big fan of that feature, but if we want it, we should figure > out how to do it for untagged types. We've already discussed that in a different e-mail thread. It seems dangerous. > >... We could also consider deriving both versions (usual and > > bounded) from an abstract ancestor. > > Consider, ..., and reject. ;-) Again, why? We have a request for a "universal" numeric type, and the only sane way to provide that is with dispatching. Probably, we'll just forget that request, but it seems worth spending a bit of time to see if it makes sense. **************************************************************** From: John Barnes Sent: Tuesday, January 30, 2018 4:22 AM I am feverishly giving lectures on numbers at Oxford at the moment but I am trying to keep an eye on what the ARG is up to. Did you know that a new Mersenne prime was discovered on Boxing Day (26 December) 2017. It is 2**77232917 - 1 and has only 23,249,425 digits. Will the Bignum package cope with it? **************************************************************** From: Tucker Taft Sent: Tuesday, January 30, 2018 3:02 PM Yes, I noticed that new Mersenne prime as well. And 23 mega-digit is nothing for a modern iPhone. ;-) Just be sure to set aside a bit of extra time to print it out using the Image function. Except in base 2, of course, which I could do right now. Ready: 1111111111 [... 77,232,900 1's] 1111111! **************************************************************** From: Jeff Cousins Sent: Wednesday, January 31, 2018 6:31 AM [This is John Barnes' Bignum package and some test programs - Editor.] -- file books\fun\progs\numbers.ada -- Restructured using children -- Types and No_Of_Places in parent package -- 20-10-06 package Numbers is Max_Index: constant := 1000; subtype Index is Integer range 0 .. Max_Index; type Number(Max_Digits: Index := 1) is private; Zero, One: constant Number; Number_Error : exception; private Base_Exp: constant := 4; Base: constant := 10 ** Base_Exp; type Base_Digit is range -Base .. 2 * Base - 1; type Base_Digit_Array is array(Index range <>) of Base_Digit; type Number(Max_Digits: Index := 1) is record Sign: Integer := +1; Length: Index := 0; D: Base_Digit_Array(1..Max_Digits); end record; Zero: constant Number := (0, +1, 0, (others => 0)); One: constant Number := (1, +1, 1, (1 => 1)); function No_Of_Places(N: Number) return Integer; end Numbers; package body Numbers is function No_Of_Places(N: Number) return Integer is begin if N.Length = 0 then return 1; else return N.Length * Base_Exp; end if; end No_Of_Places; end Numbers; -------------------------------------------------------------------- package Numbers.Proc is subtype Index is Numbers.Index; subtype Number is Numbers.Number; Zero: Number renames Numbers.Zero; One: Number renames Numbers.One; Number_Error: exception renames Numbers.Number_Error; procedure Add(X, Y: Number; Z: out Number); procedure Sub(X, Y: Number; Z: out Number); procedure Mul(X, Y: Number; Z: out Number); procedure Div(X, Y: Number; Quotient, Remainder: out Number); procedure Neg(X: in out Number); function Compare(X, Y: Number) return Integer; function Length(N: Number) return Index; procedure To_Number(S: String; N: out Number); procedure To_Number(I: Integer; N: out Number); procedure To_Text(N: Number; S: out String); procedure To_Integer(N: Number; I: out Integer); end Numbers.Proc; -------------------------------------------------------------------- package Numbers.IO is Default_Width: Natural := 0; procedure Put(Item: Number; Width: Natural := Default_Width); procedure Get(Item: out Number); end Numbers.IO; -------------------------------------------------------------------- package Numbers.Func is subtype Number is Numbers.Number; Zero: Number renames Numbers.Zero; One: Number renames Numbers.One; Number_Error: exception renames Numbers.Number_Error; function "+" (X: Number) return Number; function "-" (X: Number) return Number; function "abs" (X: Number) return Number; function "+" (X, Y: Number) return Number; function "-" (X, Y: Number) return Number; function "*" (X, Y: Number) return Number; function "/" (X, Y: Number) return Number; function "rem" (X, Y: Number) return Number; function "mod" (X, Y: Number) return Number; function "**" (X: Number; N: Natural) return Number; function "<" (X, Y: Number) return Boolean; function "<=" (X, Y: Number) return Boolean; function ">" (X, Y: Number) return Boolean; function ">=" (X, Y: Number) return Boolean; function "=" (X, Y: Number) return Boolean; function Make_Number(S: String) return Number; function Make_Number(I: Integer) return Number; function String_Of(N: Number) return String; function Integer_Of(N: Number) return Integer; end Numbers.Func; -------------------------------------------------------------------- package body Numbers.Proc is Base_Squared: constant := Base**2; subtype Single is Base_Digit; type Double is range -Base_Squared .. 2*Base_Squared - 1; function Unsigned_Compare(X, Y: Number) return Integer is -- ignoring signs -- returns +1, 0 or -1 according as X >, = or < Y begin if X.Length > Y.Length then return +1; end if; if X.Length < Y.Length then return -1; end if; for I in reverse 1 .. X.Length loop if X.D(I) > Y.D(I) then return +1; end if; if X.D(I) < Y.D(I) then return -1; end if; end loop; return 0; -- the numbers are equal end Unsigned_Compare; function Compare(X, Y: Number) return Integer is -- returns +1, 0 or -1 according as X >, = or < Y begin if X.Sign /= Y.Sign then return X.Sign; end if; return Unsigned_Compare(X, Y) * X.Sign; end Compare; procedure Raw_Add(X, Y: Number; Z: out Number) is -- assumes X not smaller than Y Carry: Single := 0; Digit: Single; ZL: Index := X.Length; -- length of answer begin if Z.Max_Digits < ZL then raise Number_Error; -- Z not big enough to hold X end if; for I in 1 .. ZL loop Digit := X.D(I) + Carry; if I <= Y.Length then Digit := Digit + Y.D(I); end if; if Digit >= Base then Carry := 1; Digit := Digit - Base; else Carry := 0; end if; Z.D(I) := Digit; end loop; if Carry /= 0 then if ZL = Z.Max_Digits then raise Number_Error; -- too big to fit in Z end if; ZL := ZL + 1; Z.D(ZL) := Carry; end if; Z.Length := ZL; end Raw_Add; procedure Raw_Sub(X, Y: Number; Z: out Number) is -- assumes X not smaller than Y Carry: Single := 0; Digit: Single; ZL: Index := X.Length; -- length of answer begin if Z.Max_Digits < ZL then raise Number_Error; -- Z not big enough to hold X end if; for I in 1 .. ZL loop Digit := X.D(I) - Carry; if I <= Y.Length then Digit := Digit - Y.D(I); end if; if Digit < 0 then Carry := 1; Digit := Digit + Base; else Carry := 0; end if; Z.D(I) := Digit; end loop; while Z.D(ZL) = 0 loop -- SHOULD THIS NOT FAIL??? ZL := ZL - 1; -- remove leading zeroes end loop; Z.Length := ZL; end Raw_Sub; procedure Add(X, Y: Number; Z: out Number) is UCMPXY: Integer := Unsigned_Compare(X, Y); begin if X.Sign = Y.Sign then Z.Sign := X.Sign; if UCMPXY >= 0 then Raw_Add(X, Y, Z); else Raw_Add(Y, X, Z); -- reverse if Y larger end if; else if UCMPXY > 0 then Raw_Sub(X, Y, Z); Z.Sign := X.Sign; elsif UCMPXY < 0 then Raw_Sub(Y, X, Z); Z.Sign := -X.Sign; else -- answer is zero Z.Sign := +1; Z.Length := 0; end if; end if; end Add; procedure Sub(X, Y: Number; Z: out Number) is UCMPXY: Integer := Unsigned_Compare(X, Y); begin if X.Sign /= Y.Sign then Z.Sign := X.Sign; if UCMPXY >= 0 then Raw_Add(X, Y, Z); else Raw_Add(Y, X, Z); -- reverse if Y larger end if; else if UCMPXY > 0 then Raw_Sub(X, Y, Z); Z.Sign := X.Sign; elsif UCMPXY < 0 then Raw_Sub(Y, X, Z); Z.Sign := -X.Sign; else -- answer is zero Z.Sign := +1; Z.Length := 0; end if; end if; end Sub; procedure Neg(X: in out Number) is begin -- do nothing in zero case if X.Length = 0 then return; end if; X.Sign := -X.Sign; end Neg; function Length(N: Number) return Index is begin return N.Length; end Length; procedure Mul(X, Y: Number; Z: out Number) is Carry: Double; Digit: Double; ZL: Index; begin if Z.Max_Digits < X.Length + Y.Length then raise Number_Error; end if; if X.Length = 0 or Y.Length = 0 then -- zero case Z.Sign := +1; Z.Length := 0; return; end if; ZL := X.Length + Y.Length - 1; -- lower possible length of answer -- copy X to top of Z; so X and Z can be same array for I in reverse 1 .. X.Length loop Z.D(I + Y.Length) := X.D(I); end loop; declare -- initialise limits and length of cycle Z_Index: Index; Y_Index: Index; Initial_Z_Index: Index := Y.Length + 1; Initial_Y_Index: Index := 1; Cycle_Length: Index := 1; begin Carry := 0; for I in 1 .. ZL loop Digit := Carry; Carry := 0; Z_Index := Initial_Z_Index; Y_Index := Initial_Y_Index; for J in 1 .. Cycle_Length loop if Digit > Base_Squared then Digit := Digit - Base_Squared; Carry := Carry + Base; end if; Digit := Digit + Double(Z.D(Z_Index)) * Double(Y.D(Y_Index)); Z_Index := Z_Index + 1; Y_Index := Y_Index - 1; end loop; -- now adjust limits and length of cycle if I < Y.Length then Cycle_Length := Cycle_Length + 1; Initial_Y_Index := Initial_Y_Index + 1; else Initial_Z_Index := Initial_Z_Index + 1; end if; if I < X.Length then Cycle_Length := Cycle_Length + 1; end if; Cycle_Length := Cycle_length - 1; Carry := Carry + Digit / Base; Z.D(I) := Single(Digit mod Base); end loop; end; if Carry /= 0 then -- one more digit in answer ZL := ZL + 1; Z.D(ZL) := Single(Carry); end if; Z.Length := ZL; Z.Sign := X.Sign * Y.Sign; end Mul; procedure Div(X, Y: Number; Quotient, Remainder: out Number) is U: Number renames Quotient; V: Number renames Remainder; Digit, Scale, Carry: Double; U0, U1, U2: Double; V1, V2: Double; QD: Double; LOQ: constant Index := Y.Length; HIQ: constant Index := X.Length; QL: Index; RL: Index; QStart: Index; J : Index; begin if Y.Length = 0 then raise Number_Error; end if; if Quotient.Max_Digits < X.Length or Remainder.Max_Digits < Y.Length then raise Number_Error; end if; if X.Length < Y.Length then -- Quotient is definitely zero Quotient.Sign := +1; Quotient.Length := 0; Remainder.Sign := X.Sign; Remainder.Length := X.Length; for I in 1 .. X.Length loop Remainder.D(I) := X.D(I); end loop; return; end if; QL := X.Length - Y.Length + 1; RL := Y.Length; QStart := QL; -- compute normalizing factor Scale := Base/Double(Y.D(Y.Length)+1); -- scale X and copy to U Carry := 0; for I in 1 .. X.Length loop Digit := Double(X.D(I)) * Scale + Carry; Carry := Digit / Base; U.D(I) := Single(Digit mod Base); end loop; U0 := Carry; -- leading digit of dividend -- scale Y and copy to V Carry := 0; for I in 1 .. Y.Length loop Digit := Double(Y.D(I)) * Scale + Carry; Carry := Digit / Base; V.D(I) := Single(Digit mod Base); end loop; -- no further carry -- set V1 and V2 to first two digits of divisor V1 := Double(V.D(Y.Length)); if Y.Length > 1 then V2 := Double(V.D(Y.Length-1)); else V2 := 0; end if; -- now iterate over digits in answer -- with U0, U1 and U2 being first three digits of dividend for I in reverse LOQ .. HIQ loop U1 := Double(U.D(I)); if Y.Length > 1 then U2 := Double(U.D(I-1)); else U2 := 0; end if; -- now set initial estimate of digit in quotient if U0 = V1 then QD := Base - 1; else QD := (U0 * Base + U1) / V1; end if; -- now refine estimate by considering U2 also while V2*QD > (U0*Base+U1-QD*V1)*Base + U2 loop QD := QD - 1; end loop; -- QD is now correct digit or possibly one too big -- subtract QD times V from U Carry := 0; J := QStart; for I in 1 .. Y.Length loop Digit := Double(U.D(J)) - Carry - QD * Double(V.D(I)); if Digit < 0 then Carry := (-1-Digit) / Base + 1; Digit := Digit + Carry * Base; else Carry := 0; end if; U.D(J) := Single(Digit); J := J + 1; end loop; if Carry > U0 then -- estimate was too large declare Carry, Digit: Single; begin QD := QD - 1; Carry := 0; J := QStart; for I in 1 .. Y.Length loop Digit := U.D(J) + Carry + V.D(I); if Digit >= Base then Carry := 1; Digit := Digit - Base; else Carry := 0; end if; U.D(J) := Digit; J := J + 1; end loop; end; end if; -- QD is now the required digit U0 := Double(U.D(I)); U.D(I) := Single(QD); QStart := QStart - 1; end loop; -- delete possible leading zero in quotient if U.D(HIQ) = 0 then QL := QL - 1; end if; -- copy remainder into place and scale -- top digit is in U0 still Digit := U0; for I in reverse 2 .. RL loop Remainder.D(I) := Single(Digit/Scale); Carry := Digit mod Scale; Digit := Double(U.D(I-1)) + Carry * Base; end loop; Remainder.D(1) := Single(Digit/Scale); -- delete leading zeroes in remainder while RL > 0 and then Remainder.D(RL) = 0 loop RL := RL - 1; end loop; Remainder.Length := RL; if Remainder.Length = 0 then Remainder.Sign := +1; else Remainder.Sign := X.Sign; end if; -- slide quotient into place -- Quotient.D(1 .. QL) := U.D(LOQ .. HIQ); for I in 1 .. QL loop Quotient.D(I) := U.D(I + LOQ - 1); end loop; Quotient.Length := QL; if Quotient.Length = 0 then Quotient.Sign := +1; else Quotient.Sign := X.Sign * Y.Sign; end if; end Div; procedure To_Number(S: String; N: out Number) is NL: Index := 0; Place: Integer := 0; Is_A_Number: Boolean := False; Digit: Single := 0; Ch: Character; Last_I: Positive; Dig_Of: constant array (Character range '0' .. '9') of Single := (0, 1, 2, 3, 4, 5, 6, 7, 8, 9); begin N.Sign := +1; -- set default sign -- scan string from end for I in reverse S'Range loop Last_I := I; -- note how far we have got Ch := S(I); case Ch is when '0' .. '9' => -- add digit to number so far if Place = 0 then NL := NL + 1; if NL > N.Max_Digits then raise Number_Error; end if; end if; Digit := Digit + Dig_Of(Ch) * 10**Place; Place := Place + 1; if Place = Base_Exp then N.D(NL) := Digit; Digit := 0; Place := 0; end if; Is_A_Number := True; when '_' => -- underscore must be embedded in digits if not Is_A_Number then raise Number_Error; end if; Is_A_Number := False; when '+' | '-' | ' ' => -- lump so far must be a valid number if not Is_A_Number then raise Number_Error; end if; if Ch ='-' then N.Sign := -1; end if; exit; -- leave loop when others => raise Number_Error; end case; end loop; -- check we had a number if not Is_A_Number then raise Number_Error; end if; -- add the last digit if necessary if Place /= 0 then N.D(NL) := Digit; end if; -- check that any other characters are leading spaces for I in S'First .. Last_I - 1 loop if S(I) /= ' ' then raise Number_Error; end if; end loop; -- remove leading zeroes if any, beware zero case while NL > 0 and then N.D(NL) = 0 loop NL := NL - 1; end loop; N.Length := NL; end To_Number; procedure To_Number(I: Integer; N: out Number) is NL: Index := 0; II: Integer; begin if I = 0 then N.Sign := +1; N.Length := 0; return; end if; if I > 0 then II := I; N.Sign := +1; else II := -I; N.Sign := -1; end if; while II /= 0 loop NL := NL + 1; if NL > N.Max_Digits then raise Number_Error; end if; N.D(NL) := Single(II mod Base); II := II / Base; end loop; N.Length := NL; end To_Number; procedure To_Text(N: Number; S: out String) is SI: Natural := S'Last; Digit: Single; Char_Of: constant array (Single range 0 .. 9) of Character := "0123456789"; begin if N.Length = 0 then -- zero case if SI < 2 then raise Number_Error; end if; S(SI) := Char_Of(0); S(SI-1) := '+'; for I in 1 .. SI-2 loop S(I) := ' '; end loop; return; end if; if SI < Base_Exp * N.Length + 1 then raise Number_Error; end if; for I in 1 .. N.Length loop Digit := N.D(I); for J in 1 .. Base_Exp loop S(SI) := Char_Of(Digit mod 10); Digit := Digit / 10; SI := SI - 1; end loop; end loop; while S(SI + 1) = '0' loop SI := SI + 1; -- delete leading zeroes end loop; if N.Sign = +1 then S(SI) := '+'; else S(SI) := '-'; end if; for I in 1 .. SI - 1 loop S(I) := ' '; end loop; end To_Text; procedure To_Integer(N: Number; I: out Integer) is II: Integer := 0; begin for I in reverse 1 .. N.Length loop II := II * Base + Integer(N.D(I)); end loop; if N.Sign = -1 then II := -II; end if; I := II; end To_Integer; end Numbers.Proc; -------------------------------------------------------------------- with Ada.Text_IO; use Ada; with Numbers.Proc; package body Numbers.IO is use Proc; procedure Put(Item: Number; Width: Natural := Default_Width) is Block_Size: constant := 3; Places: Integer := No_Of_Places(Item); S: String(1 .. Places + 1); SP: Positive := 1; Before_Break: Integer; begin To_Text(Item, S); -- allow for leading spaces in S while S(SP) = ' ' loop SP := SP + 1; Places := Places - 1; end loop; -- now output leading spaces for padding if any for I in 1 .. Width - (Places + 1 + (Places - 1) / Block_Size) loop Text_IO.Put(' '); end loop; if S(SP) = '+' then S(SP) := ' '; end if; Text_IO.Put(S(SP)); -- output minus or space -- output digits with underscores every "Blocksize" Before_Break := (Places - 1) rem Block_Size + 1; for I in SP + 1 .. S'Last loop if Before_Break = 0 then Text_IO.Put('_'); Before_Break := Block_Size; end if; Text_IO.Put(S(I)); Before_Break := Before_Break - 1; end loop; end Put; procedure Get(Item: out Number) is -- declare string large enough to hold maximum value -- allows every other character to be an underscore! S: String(1 .. Base_Exp * Max_Index * 2); SP: Positive := 1; Places: Integer := 0; Ch: Character; EOL: Boolean; -- end of line begin -- loop for first digit or sign, skipping spaces loop Text_IO.Get(Ch); case Ch is when ' ' => null; when '+'| '-' => S(SP) := Ch; SP := SP + 1; exit; when '0' .. '9' => S(SP) := Ch; SP := SP + 1; Places := 1; exit; when others => raise Number_Error; end case; end loop; -- now accept only digits and underscores -- count the digits in Places -- stop on end of line or other character loop Text_IO.Look_Ahead(Ch, EOL); exit when EOL; Text_IO.Get(Ch); case Ch is when '0' .. '9' => S(SP) := Ch; SP := SP + 1; Places := Places + 1; when '_' => S(SP) := Ch; SP := SP + 1; when others => exit; end case; end loop; -- now declare a Number big enough -- note Item assumed unconstrained declare Result: Number((Places - 1)/Base_Exp + 1); begin To_Number(S(1 .. SP - 1), Result); Item := Result; end; end Get; end Numbers.IO; -------------------------------------------------------------------- with Numbers.Proc; package body Numbers.Func is use Proc; function "+" (X: Number) return Number is begin return X; end "+"; function "-" (X: Number) return Number is N: Number(X.Max_Digits); begin N := X; Neg(N); return N; end "-"; function "abs" (X: Number) return Number is begin if X < Zero then return -X; else return X; end if; end "abs"; function "+" (X, Y: Number) return Number is Z: Number(Index'Max(Length(X), Length(Y)) + 1); begin Add(X, Y, Z); return Z; end "+"; function "-" (X, Y: Number) return Number is Z: Number(Index'Max(Length(X), Length(Y)) + 1); begin Sub(X, Y, Z); return Z; end "-"; function "*" (X, Y: Number) return Number is Z: Number(Length(X) + Length(Y)); begin Mul(X, Y, Z); return Z; end "*"; function "/" (X, Y: Number) return Number is Q: Number(Length(X)); R: Number(Length(Y)); begin Div(X, Y, Q, R); return Q; end "/"; function "rem" (X, Y: Number) return Number is Q: Number(Length(X)); R: Number(Length(Y)); begin Div(X, Y, Q, R); return R; end "rem"; function "mod" (X, Y: Number) return Number is Q: Number(Length(X)); R: Number(Length(Y)); begin Div(X, Y, Q, R); if (X < Zero and Y > Zero) or (X > Zero and Y < Zero) then R := R + Y; end if; return R; end "mod"; function "**" (X: Number; N: Natural) return Number is Result: Number := One; Term: Number := X; M: Natural := N; begin loop if M rem 2 /= 0 then Result := Term * Result; end if; M := M / 2; exit when M = 0; Term := Term * Term; end loop; return Result; end "**"; function "<" (X, Y: Number) return Boolean is begin return Compare(X, Y) < 0; end "<"; function "<=" (X, Y: Number) return Boolean is begin return Compare(X, Y) <= 0; end "<="; function ">" (X, Y: Number) return Boolean is begin return Compare(X, Y) > 0; end ">"; function ">=" (X, Y: Number) return Boolean is begin return Compare(X, Y) >= 0; end ">="; function "=" (X, Y: Number) return Boolean is begin return Compare(X, Y) = 0; end "="; function Make_Number(S: String) return Number is Result: Number((S'Length - 1) / Base_Exp + 1); begin To_Number(S, Result); return Result; end Make_Number; function Make_Number(I: Integer) return Number is Base_Digits: Index := 0; II: Integer := abs I; Base: constant := 10 ** Base_Exp; begin -- loop to determine discriminant for result while II /= 0 loop Base_Digits := Base_Digits + 1; II := II / Base; end loop; declare Result: Number(Base_Digits); begin To_Number(I, Result); return Result; end; end Make_Number; function String_Of(N: Number) return String is Places: Integer := No_Of_Places(N); S: String(1 .. Places + 1); SP: Positive := 1; begin To_Text(N, S); -- allow for leading spaces in S while S(SP) = ' ' loop SP := SP + 1; end loop; return S(SP .. S'Last); end String_Of; function Integer_Of(N: Number) return Integer is Result: Integer; begin To_Integer(N, Result); return Result; end Integer_Of; end Numbers.Func; --- Numbers calculator --------------------------------------------- --with ID, Numbers.Func; use Numbers.Func; --package Numbers_ID is new ID(Number, Zero, One); --with ID.IO, Numbers.IO; use Numbers.IO; --package Numbers_ID.The_IO is new Numbers_ID.IO; --with Numbers_ID.The_IO; --with Numbers.Func; --with Calculator; --procedure Numcalc is -- new Calculator.Run(Numbers_ID, -- Numbers_ID.The_IO, -- Numbers.Func."+"); --- test 1 - powers of 11 and 99 ----------------------------------- with Numbers.Proc; use Numbers.Proc; with Ada.Text_IO; use Ada; procedure Test_11_99 is U: Number(2); procedure P(X: Number) is S: String(1 .. 150); begin To_Text(X, S); Text_IO.Put(S); Text_IO.New_Line; end P; procedure Power(U: Number; V: Integer) is W: Number(50); begin To_Number(1, W); for I in 1 .. V loop Mul(W, U, W); P(W); end loop; end Power; begin To_Number(11, U); Power(U, 50); To_Number(99, U); Power(U, 50); end; -- test 2 - Mersenne using procedural forms ------------------------ with Ada.Calendar; use Ada.Calendar; with Numbers.Proc; use Numbers.Proc; with Ada.Text_IO, Ada.Integer_Text_IO; use Ada.Text_IO, Ada.Integer_Text_IO; procedure Test2 is Nop: constant := 30; Loop_Start, Loop_End: Integer; T_Start, T_End: Time; Is_Prime: Boolean; MM: Number(50); Primes: array(1 .. Nop) of Integer := (3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59, 61,67,71,73,79,83,89,97,101,103,107,109,113,127); package Duration_IO is new Fixed_IO(Duration); use Duration_IO; procedure Lucas_Lehmer (P: Integer; Mersenne: out Number; Is_Prime: out Boolean) is Two: Number(1); M: Number(Mersenne.Max_Digits); L, W, Quotient: Number(M.Max_Digits*2); begin To_Number(2, Two); To_Number(4, L); To_Number(1, M); for I in 1 .. P loop Mul(M, Two, M); end loop; Sub(M, One, M); for I in 1 .. P-2 loop Mul(L, L, W); Sub(W, Two, W); Div(W, M, Quotient, L); -- L := (L**2 - Two) mod M; end loop; Is_Prime := Compare(L, Zero) = 0; Mersenne := M; end Lucas_Lehmer; procedure Put(X: Number) is S: String(1 .. 45); begin To_Text(X, S); Put(S); end Put; begin Put_Line("Start loop? "); Get(Loop_Start); Put_Line("End loop? "); Get(Loop_End); Put_Line("Mersenne Primes"); Put_Line(" Time P 2**P-1"); for I in Loop_Start .. Loop_End loop New_Line; T_Start := Clock; Lucas_Lehmer(Primes(I), MM, Is_Prime); T_End := Clock; New_Line; Put(T_End - T_Start, 2, 1); Put(Primes(I), 4); Put(" : "); Put(MM); if Is_Prime then Put(" is prime"); else Put(" is not prime"); end if; end loop; end Test2; --- test 5 - Mersenne using functional forms ----------------------- with Ada.Calendar; use Ada.Calendar; with Numbers.Func; use Numbers.Func; with Numbers.IO; use Numbers.IO; with Ada.Text_IO, Ada.Integer_Text_IO; use Ada.Text_IO, Ada.Integer_Text_IO; procedure Test5 is Nop: constant := 30; Loop_Start, Loop_End: Integer; T_Start, T_End: Time; Is_Prime: Boolean; LL, MM: Number; Primes: array(1 .. Nop) of Integer := (3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59, 61,67,71,73,79,83,89,97,101,103,107,109,113,127); package Duration_IO is new Fixed_IO(Duration); use Duration_IO; procedure Lucas_Lehmer (Q: Integer; Mersenne: out Number; Lout: out Number; Is_Prime: out Boolean) is Two: constant Number := Make_Number(2); M: constant Number := Two**Q - One; L: Number := Make_Number(4); begin for I in 1 .. Q-2 loop L := (L**2 - Two); -- mod M; -- mod M here is optional; Put(L); New_line(2); end loop; Is_Prime := L mod M = Zero; Lout := l; Mersenne := M; end Lucas_Lehmer; begin Put_Line("Start loop? "); Get(Loop_Start); Put_Line("End loop? "); Get(Loop_End); Put_Line("Mersenne Primes"); Put_Line(" Time P 2**P-1"); for I in Loop_Start .. Loop_End loop New_Line; T_Start := Clock; Lucas_Lehmer(Primes(I), MM, LL, Is_Prime); T_End := Clock; New_Line; Put(T_End - T_Start, 2, 1); Put(Primes(I), 4); Put(" : "); Put(MM, 20); if Is_Prime then Put(" is prime"); New_line; put((MM+One)/Make_number(2)*MM, 30); Put(" is perfect"); New_Line; else Put(" is not prime"); end if; New_line; -- comment out next four lines to avoid detail Put("L is "); Put(LL); Put(" equals "); Put(MM); Put(" times "); Put(LL/MM); if not Is_prime then New_Line; Put(" remainder = "); Put(LL mod MM); end if; -- end of comment end loop; skip_Line(2); end Test5; --- test 6 - GCD and Mersenne -------------------------------------- with Ada.Text_IO, Ada.Integer_Text_IO; use Ada.Text_IO, Ada.Integer_Text_IO; with Numbers.Proc; use Numbers.Proc; procedure Test6 is subtype CNumber is Number(50); M1, M2, M3: CNumber; P1, P2, P3: Integer; G, H: CNumber; XX, YY, QQ, ZZ: CNumber; Two: CNumber; N: Cnumber; Start, Stop: Integer; procedure GCD(X, Y: CNumber; Z: out CNumber) is begin XX := X; YY := Y; while Compare(YY, Zero) /= 0 loop Div(XX, YY, QQ, ZZ); XX := YY; YY := ZZ; end loop; Z := XX; end GCD; procedure Mersenne(P: Integer; M: in out CNumber) is begin To_Number(2, Two); To_Number(1, N); for I in 1 .. P loop Mul(N, Two, N); end loop; Sub(N, One, N); M := N; end Mersenne; begin To_Number(2, Two); Put_Line("Start? "); Get(Start); Put_Line("Stop? "); Get(Stop); for I in Start .. Stop loop P1 := 2*I + 1; Mersenne(P1, M1); for J in 1 .. I -1 loop P2 := 2*J + 1; Mersenne(P2, M2); GCD(M1, M2, G); for K in 1 .. I loop P3 := 2*K + 1; Mersenne(P3, M3); if Compare(M3, G) > 0 then exit; end if; if Compare(M3, G) =0 then New_Line; Put(P1); Put(P2); Put(P3); exit; end if; end loop; end loop; end loop; end Test6; ------- test 7 multipication with Ada.Text_IO, Ada.Integer_Text_IO; use Ada.Text_IO, Ada.Integer_Text_IO; with Numbers.Proc; use Numbers.Proc; with Numbers.IO; use Numbers.IO; procedure Test7 is X: Number; Y: Number; Z: Number(50); begin Put_line("Multiplier test"); Put("X = "); Get(X); Put("Y = "); Get(Y); Mul(X, Y, Z); Put_Line("product is "); Put(Z); New_Line(2); end Test7; **************************************************************** From: Jeff Cousins Sent: Wednesday, January 31, 2018 6:17 AM [This some additional test programs for John Barnes' Bignum package; the Numbers package was duplicated at the front, which I removed. - Editor.] --package Monitor is -- C1, C2, C3, C4: Integer; --end; package Primes is Max: constant := 20000; Prime: array (1..Max) of Integer; pragma Elaborate_Body; end; package body Primes is N: Integer:= 2; Index: Integer := Prime'First; Found: Boolean := False; begin -- initialization part to build prime table Prime(Prime'First) := 2; loop N := N+1; Found := True; for I in Prime'First .. Index loop if n rem Prime(I) = 0 then -- divides by existing prime Found := False; exit; end if; end loop; if Found then -- found a new prime Index := Index+1; Prime(Index) := n; exit when Index = Max; end if; end loop; end Primes; Package Squares is Square_Digits: Integer := 4; Square_Ten_Power: Integer := 10**Square_Digits; Poss_Last_Digits: array(0..Square_Ten_Power-1) of Boolean := (others => False); pragma Elaborate_Body; end; package body Squares is begin -- make Poss_Last_Digits array for I in 1 .. Square_Ten_Power loop Poss_Last_Digits(I*I mod Square_Ten_Power) := True; end loop; end Squares; with Numbers.Func; use Numbers.Func; -- with Monitor; procedure Square_root(XX: in Number; Try: Number; X: out Number; R: out Number) is -- X is the largest X such that X*X + R equals XX with R >= zero -- Try is initial guess K, KK: Number; DK: Number; Three: Number := Make_Number(3); begin k := Try; loop KK := K*K; DK := (XX - KK)/(K+K); -- iterate until nearly there if abs DK < Three then -- nearly there if XX < KK then K := XX/K; -- ensure K*K is less than XX end if; exit; end if; -- do another iterate K := K+DK; -- Monitor.c2 := Monitor.C2 + 1; end loop; -- now loop from below loop KK := K*K; -- Monitor.C2 := Monitor.C2 + 1; if KK >= XX then if KK = XX then X := K; R := Zero; return; end if; X := K - One; R := XX - X*X; return; end if; K := K+One; end loop; end Square_Root; with Square_Root; with Numbers.Func; use Numbers.Func; -- with Monitor; with Squares; use Squares; procedure Fermat_Long(N: in Number; Min_Prime: in Number; P, Q: out Number) is -- we know that factors up to Min_Prime have been removed, so max square to try is -- roughly (N/Min_Prime + Min_Prime)/2; we add 2 for luck Two: Number := Make_Number(2); Number_Square_Digits: Number := Make_Number(Square_Digits); Number_Square_Ten_Power: Number := Make_Number(Square_Ten_Power); X: Number; Y: Number; R: Number; K: Number; DK: Number; N2, X2, K2: Integer; Last_Digits: Integer; Try: Number := One; Max_square : Number := (N/Min_Prime + Min_Prime)/ Two + Two; begin Square_Root(N, One, X, R); if R = Zero then -- N was a perfect square P := X; Q := X; return; end if; -- Monitor.C1 := Monitor.C2; -- Monitor.C2 := 0; K := X*X-n; DK := X+X+One; N2 := Integer_Of(N rem Number_Square_Ten_Power); X2 := Integer_Of(X rem Number_Square_Ten_Power); loop -- Monitor.C3 := Monitor.C3 + 1; X := X + One; if X > Max_Square then -- must be prime P := N; Q := One; return; end if; X2 := (X2 + 1) rem Square_Ten_Power; K := K + DK; -- omit if DK not used DK := DK + Two; -- K := X*X-N; -- omit if DK used K2 := (X2*X2-N2) mod Square_Ten_Power; -- Last_Digits := Integer_Of(K rem Number_Square_Ten_Power); Last_Digits := K2; if Poss_Last_Digits(Last_Digits) then -- Monitor.C4 := Monitor.C4 + 1; Square_Root(K, Try, Y, R); if R = Zero then -- X*X-N was a perfect square P := X+Y; Q := X-Y; return; end if; Try := Y; end if; end loop; exception when others => p := Zero; Q := Zero; end Fermat_long; with Numbers.IO; use Numbers.IO; with Numbers.Func; use Numbers.Func; with Ada.Text_IO; use Ada.Text_IO; with Ada.Integer_Text_IO; use Ada.Integer_Text_IO; -- with Monitor; with Primes; use Primes; with Fermat_Long; with Ada.Calendar; use Ada.Calendar; procedure Do_Fermat_Long is NN, PP, QQ: Number; T_Start, T_End: Time; package Duration_IO is new Fixed_IO(Duration); use Duration_IO; begin Put_Line("Welcome to Fermat's method (multilength)"); loop <> Begin New_Line; Put("Insert number N = "); Get(NN); exception when others => Put_Line("Not a number or too big"); Skip_Line; goto Again; end; exit when NN = Zero; -- check to see if N divides by a known prime for i in Prime'Range loop loop if NN rem make_Number(Prime(I)) = Zero then Put("N divides by "); Put(Prime(I), 0); Put_Line(" so removing factor"); NN := NN / Make_Number(Prime(I)); else exit; end if; end loop; end loop; if nN rem Make_number(4) = Make_Number(2) then Put_Line("Algorithm fails on odd multiples of 2, so halving N"); NN := NN/Make_Number(2); Put("New N = "); Put(NN, 0); New_Line; end if; if NN = One then Put_Line("all factors removed"); else -- Monitor.C1 := 0; -- Monitor.C2 := 0; -- Monitor.C3 := 0; -- Monitor.C4 := 0; T_Start := Clock; Fermat_Long(NN, Make_Number(Prime(Max)), PP, QQ); T_End := Clock; if PP = Zero and QQ = zero then Put_Line("Failed internally"); goto Again; end if; Put("Two factors are "); Put(PP, 0); Put(" "); Put(QQ, 0); New_Line; Put(T_End-T_Start, 2, 1); New_Line; -- Put(Monitor.C1, 9); Put(Monitor.C2, 9); Put(Monitor.C3, 9); Put(Monitor.C4, 9); New_Line; end if; end loop; Put_line("Goodbye"); Skip_Line(2); end Do_Fermat_Long; **************************************************************** From: Randy Brukardt Sent: Wednesday, January 31, 2018 5:03 PM >Sending in three parts as the original message was too big for the mail >server. Just so the rest of you know, "Part 3" was an executable Windows program which we've decided not to distribute at all (it's too big for the list, and executables for specific computers seem out-of-bounds for this list anyway, given all of us know how to operate our favorite Ada compiler and probably several non-favorite compilers as well). So don't look for the non-existent part 3. ****************************************************************