Re: [boost] Is there interest in Computable Calculus?

Subject: Re: [boost] Is there interest in Computable Calculus?
From: Damian Vicino (damian.vicino_at_[hidden])
Date: 2016-02-05 12:48:52

Thanks for all the comments Alex, we are still in an early stage and every input is useful for us. Sorry for the long delay in the reply.

I'm adding some context of our current work before going to the replies.

My motivation for the first implementation of this library was in the domain of Discrete-Event Simulation (DES). Looking for a representation for Time variables.
The main problem of using Floating point or other approximated data types in (some) DES is that errors of the approximations could not be bounded.

In the case of Time, only 3 operations were interesting for me, (+, <, =).

Each number was represented by a vector <c0, c1, c2, ….>. The vector could be interpreted as the sum(<c0, c1, c2, c3….> x <1, i1, i2, i3…>) where cx are rational coefficients and ix are irrational constants.
The c values could be declared on the construction of a new number, or obtained as result of an operation. The i1 constants are declared as Template parameters, each of them is a generation algorithm defined in operator().

In this context, ensuring termination was enforced in two ways. First, the set of ix constants used in each simulation is “linearly independent over Q”, and second by a MAX_EXPANSION_ALLOWED constant which forces the throw of an exception according to what user considers “too much” execution, or works as safe guard for errors in choosing the set of constants.

We started our work with using only pi, e^pi, and some particular roots of rational numbers (e.g. sq_rt(2)). Afterward,  we expanded to parametric sets of algorithms, e.g. roots of rational numbers, for them, we have to deal with conversions of operations where operating with a particular root result is related to a new irrational constant or may become rational.

> What do you mean by restricting the set of numbers? Restricting the
> expressibility for the digit generating functions; e.g., providing some
> basic reals (pi, e, ...) and the only other way to generate new ones is
> by using arithmetic operations?

Yes. 
But, leaving the set open to be introduced by the user as template parameters. For example, e and pi is not proven to be linearly independent set of constants over Q, so it should not be used. However, pi and sq_rt(2), or  e and sq_rt(2) are both valid sets of irrational constants.

> 
> In that setting, equality might indeed be decidable, but not by the
> algorithm you proposed. But I think you have to restrain to symbolic
> manipulations, most likely by finding some normal form for arithmetic
> expressions and comparing them. (Btw, that's what I meant by Computer
> Algebra System in my earlier message.)

Yes, It is a CAS then :)
> 
> For clarity, is the goal infinite precision equality?

Actually, my motivation was infinite precision of < comparison, but equality was required for getting that (at least using this solution).

> 
>>> Boost.Multiprecision also has infinite precision (in theory). The difference of your suggestion would be that as you do lazy-evaluation, you don't have to specify the precision a priori.
>> 
>> Here, multi precision needs the whole representation to be available for every operation. In the case of comparing using algorithms, you only have representation for the current digit dependencies. 
>> For example, if you have compared the first few digits in an operation and they are equal, then for comparing the following thousands you don’t need them (unless you require them to compare the following digits, but thats a problem with the algorithm used for the number, not with the approach).
>> For example the binary number composed by all zeros, but positions p where p is power of 2. Only requires a counter for deciding each number. And the limitation of how many digits can be expanded is imposed by time, almost no space usage at all. It is a tradeoff between time/space.
> 
> I don't really understand. I see that for certain numbers the approach
> might work -- if, together with the number, you also provide
> implementations for equality and all the other relations and operations.

I don’t see how you distinguish (in multi precision) a very large rational from an irrational. 
Then, you cannot tell when digits are exhausted was because of an approximation of the true rational value. Also, you will operate with every digit in every operation. 
In our described approach, equality and addition only work with the coefficients with no requirement of a large representation (0.1pi == 0.1 pi only checks 0.1==0.1, and not every digit used if represented as MP, which we expect to be a lot of digits).

In our case, expansion is only used for < comparison (only after checking != is true).
> 
> Let's make it more specific: You can generate pi, for example, by
> different series, some converging faster than others. Do you mean that
> by generator? Then, of course, the speed of the generators is the basis
> for the performance, but that's not my point.
> 
Yes, I meant that. 

> Maybe I don't understand, so let's consider an example: Take any natural
> i and put y = x_i, where x_0=pi and x_i=((x_{i-1}+1)-1 for natural i>0.
> How do you implement deciding y == pi?

We take a set of constants including pi, for example, <1, pi, sqrt2> (the one is always implicitly included for representing rationals).

Something like this:

class pi : public irrational_constant{
…
};
class sqrt2 : public irrational_constant{
…
};

real<pi, sqrt2> r_x = {0, 1, 0}; // zero rational + one pi + zero square root of 2
real<pi, sqrt2> r_1 = {1, 0, 0}; //one rational + zero pi + zero square root of 2

int i = 20; //or any other number
for (int j=0;j < i; j++){
     r_x += r_1; //This adds the coefficients, takes three rational additions
     r_x -= r_1; //This adds the coefficients, takes three rational additions
     assert(r_x == real<pi, sqrt2>(0, 1, 0)); //this is always true, and only requires 3 equality operations, comparing the coefficients.
     assert(! (r_x != real<pi, sqrt2>(0, 1, 0))); //this is also always true, and only requires 3 equality operations, comparing the coefficients.
}

The <  operations is way more complex to describe. For that one, we expand each component multiplied by the constant using the algorithms provided in the books I referenced in earlier emails.
Also, operations as product and division are more fragile to enforce uniqueness of representation. For those operations, having linear independence over Q is not enough.