[PATCH] fix '^' in Interval

[Qian Yun]

(1) -> x := interval(-2.0,-1.0)$Interval Float

(1) [- 2.0, - 1.0]
Type: Interval(Float)
(2) -> x^(2/1)

(2) [0.0, 4.0000000000_000000001]
Type: Interval(Float)
(3) -> x^2

(3) [0.9999999999_9999999999, 4.0000000000_000000001]
Type: Interval(Float)


Notice that for the fractional power (2/1), the answer is wrong.

- Qian

diff --git a/src/algebra/interval.spad b/src/algebra/interval.spad
index 180b1c7b..8bd84ec1 100644
--- a/src/algebra/interval.spad
+++ b/src/algebra/interval.spad
@@ -429,7 +429,7 @@
error "fractional power only defined for x > 0"
even?(numer(n)) =>
hi < 0 =>
- interval(0, lo^n)
+ interval(hi^n, lo^n)
interval(0, max(lo^n, hi^n))
interval(lo^n, hi^n)
interval(lo^n, hi^n)

[Ralf Hemmecke]

Well, it depends what "interval" wants to express.

[0.0, 4.0] is correct in the sense that all values that lie in x their
square lies in [0.0, 4.0]. Maybe it is not the closest interval that is
possible. Given that the implementation of interval(a,b) rounds up and
down accordingly, your patch should be safe.

Ralf

[Waldek Hebisch]

On Thu, May 07, 2026 at 05:26:02PM +0800, Qian Yun wrote:
> (1) -> x := interval(-2.0,-1.0)$Interval Float
>
> (1) [- 2.0, - 1.0]
> Type: Interval(Float)
> (2) -> x^(2/1)
>
> (2) [0.0, 4.0000000000_000000001]
> Type: Interval(Float)
> (3) -> x^2
>
> (3) [0.9999999999_9999999999, 4.0000000000_000000001]
> Type: Interval(Float)
>
>
> Notice that for the fractional power (2/1), the answer is wrong.

As Ralf noted the anwer is not strictly speaking wrong, just
suboptimal. AFAICS one could get better intervals in other
case, by adding something like:

0 < lo => interval(lo^n, hi^n)

0 is only needed for case when lo <= 0 <= hi

> - Qian
>
> diff --git a/src/algebra/interval.spad b/src/algebra/interval.spad
> index 180b1c7b..8bd84ec1 100644
> --- a/src/algebra/interval.spad
> +++ b/src/algebra/interval.spad
> @@ -429,7 +429,7 @@
> error "fractional power only defined for x > 0"
> even?(numer(n)) =>
> hi < 0 =>
> - interval(0, lo^n)
> + interval(hi^n, lo^n)
> interval(0, max(lo^n, hi^n))
> interval(lo^n, hi^n)
> interval(lo^n, hi^n)

--
Waldek Hebisch

[Qian Yun]

lo <= 0 =>
even?(denom(n)) =>

error "fractional power only defined for x > 0"
even?(numer(n)) =>
hi < 0 =>

interval(hi^n, lo^n) -- lo < hi <=0
interval(0, max(lo^n, hi^n)) -- lo <=0 < hi
interval(lo^n, hi^n) -- lo <=0, odd power
interval(lo^n, hi^n) -- lo > 0


I think there is no problem in other branch?

Also, I just found out that Interval supports negative fractional
power, but not negative integer power.

AIntervalCategory does not support both.

I guess that Interval did not want to deal with inverse of
[a,b] where a<0<b?

- Qian

[Waldek Hebisch]

On Mon, May 11, 2026 at 11:41:39AM +0800, Qian Yun wrote:
> lo <= 0 =>
> even?(denom(n)) =>
> error "fractional power only defined for x > 0"
> even?(numer(n)) =>
> hi < 0 =>
> interval(hi^n, lo^n) -- lo < hi <=0
> interval(0, max(lo^n, hi^n)) -- lo <=0 < hi
> interval(lo^n, hi^n) -- lo <=0, odd power
> interval(lo^n, hi^n) -- lo > 0
>
>
> I think there is no problem in other branch?

Right, I overlooked to 'lo <= 0' test above.

> Also, I just found out that Interval supports negative fractional
> power, but not negative integer power.
>
> AIntervalCategory does not support both.
>
> I guess that Interval did not want to deal with inverse of
> [a,b] where a<0<b?

Yes, we do not want this: this is moral equivalent of division
by 0.

--
Waldek Hebisch

[Qian Yun]

On 5/11/26 11:49 AM, Waldek Hebisch wrote:
>>
>> I guess that Interval did not want to deal with inverse of
>> [a,b] where a<0<b?
>
> Yes, we do not want this: this is moral equivalent of division
> by 0.
>

So for the "^" with negative Fraction Integer argument,
shall we give an error for all intervals, or just for
intervals that contains 0?

- Qian

[Ralf Hemmecke]

> Also, I just found out that Interval supports negative fractional
> power, but not negative integer power.

> AIntervalCategory does not support both.

In AIntervalCategory, negative powers could be added, but only
conditionnally in the same place where

inv(x : %) : % ==
contains?(x, 0$R) => error "inv: interval contains 0"
interval(inv sup x, inv inf x)

is implemented.

For fractional powers, I am a bit hesitant. I introduced
AIntervalCategory to be able to compute with intervals also for more
general rings than just those that inherit from FloatingPointSystem.
However, I don't think it would be a problem to add the default
implementation

if R has Join(OrderedRing, Field, RadicalCategory) then
((x : %) ^ (n : Fraction Integer)) : % == ...

to AIntervalCategory.
Luckily, FloatingPointSystem fulfills this condition.

Ralf

[Waldek Hebisch]

Negative fractional power should be mostly eqivalent to computing
inverse and then positive power, but should avoid double
rounding impled by composition of operations. Negative
interval can not be raised to power 1/2, and it makes
sense to completely forbid raising negative numbers to
fractional powers (but current version allows this). But
if both lower bound and upper bound are positive, then negative
fractional power is well defined.

--
Waldek Hebisch

[Qian Yun]

The implementation can be done as such.

About Exports: "^ : (%, Integer) -> %" is naturally expressed in
DivisionRing. I think it's ok for Interval to have this category.
The original authors intention is to stop at Field?

- Qian

[Waldek Hebisch]

Note that Interval requires FloatingPointSystem, that is Field.
And that was clearly the intent. Also, sound use of Interval
requires directed rounding. If you want more general domains,
then you get different thing, like AInterval.

--
Waldek Hebisch

[Qian Yun]

On 5/13/26 10:07 AM, Waldek Hebisch wrote:
>>
>> About Exports: "^ : (%, Integer) -> %" is naturally expressed in
>> DivisionRing. I think it's ok for Interval to have this category.
>> The original authors intention is to stop at Field?
>
> Note that Interval requires FloatingPointSystem, that is Field.
> And that was clearly the intent. Also, sound use of Interval
> requires directed rounding. If you want more general domains,
> then you get different thing, like AInterval.
>

Just to be clear, I am worried about consistency here:

If Interval supports negative fraction power, then it should
support negative integer power.

If Interval supports negative integer power (inverse),
then it's natural to support division.

Where do we draw the line to stop?

- Qian

[Ralf Hemmecke]

On 5/13/26 04:19, Qian Yun wrote:
> If Interval supports negative fraction power, then it should
> support negative integer power.
>
> If Interval supports negative integer power (inverse),
> then it's natural to support division.
>
> Where do we draw the line to stop?

Usefulness.

I implemented AInteval, because we didn't have it before and I needed
interval computation for rational numbers.
I have no problem with adding ^:(%,Integer)->% (and default
implementation in AIntervalCategory) (division is already there) and
generalizing from Field to DivisionRing.

However, do you have a concrete usecase for intervals in a DivisionRing?

Ralf

[Qian Yun]

Oh! AInterval supports Fraction. I never realized that.

There's no concrete usecase for ^ and / in Interval,
as I said, just concerns about consistency.

- Qian

On 5/13/26 2:38 PM, 'Ralf Hemmecke' via FriCAS - computer algebra system
wrote:

[Ralf Hemmecke]

On 5/13/26 09:10, Qian Yun wrote:
> Oh! AInterval supports Fraction. I never realized that.

However, as Waldek mentioned, the default implementation of
AIntervalCategory should be overridden in Interval, because applying the
up/down rounding for both inv and * might lead to a bigger interval than
necessary.

Additionally, Interval(R) does deliberately not inherit the
implementation of AInterval. For me the reason was the qinterval
definition in AInterval. As far as I understand (Waldek please correct
me, if I am wrong), that qinterval would be used in all the definitions
coming from the defaults in AIntervalCategory. If we then had Interval
inheriting the implementation of AInterval, then, (even though Interval
defines its own qinterval), the implementation of abs would be the one
from AInterval (and thus using the definition of qinterval from
AInterval, not Interval).

I simply did not think of ^: (%, Integer) -> % for AIntervalCategory.
I agree, adding it.

Ralf

[Qian Yun]

Regarding qinterval, documentation in AIntervalCategory says
++ qinterval(inf,sup) creates a new interval without checking the
++ ordering on the elements.

Implementation in Interval says
qinterval(i : R, s : R) : % ==
[roundDown i, roundUp s]


I think in Interval we should have "qinterval(i, s) = [i, s]",
and have a local function that uses "[roundDown i, roundUp s]".
And choose to use one or the other wisely depending on the situation.

- Qian

On 5/13/26 3:44 PM, 'Ralf Hemmecke' via FriCAS - computer algebra system

[Ralf Hemmecke]

> Implementation in Interval says
> qinterval(i : R, s : R) : % ==
> [roundDown i, roundUp s]

> I think in Interval we should have "qinterval(i, s) = [i, s]",
> and have a local function that uses "[roundDown i, roundUp s]".
> And choose to use one or the other wisely depending on the situation.

Honestly, when I saw the qinterval implementation in Interval(R) for the
first time, I was a bit puzzled. Why should one ever have to round up or
down if you i and s are given on input and that is the stuff one should
use? I didn't want to touch that stuff, so I left Interval unchanged.

Now looking more closely at what roundUp/roundDown does, I believe the
problem the original author wanted to solve was that we might be in
Float (software floats where the number of bits may change in the course
of the computation). I don't want to invest time to understand the
details now, but as you see in float.spad, the precision is global for
the whole domain.

My guess is, qinterval(i,s) was supposed to return an interval
iv:=qinterval(i,s) such that if i<=x<=s with the precision at creation
time, it should also hold that contains?(iv,x) at any later point in
time when the precision might have changed. I do not currently see where
this can be violated, but I don't even know whether the above condition
is the condition the original author cared about. when roundDown/roundUp
was chosen for qinterval.

Ralf

[Waldek Hebisch]

On Wed, May 13, 2026 at 05:15:49PM +0800, Qian Yun wrote:
> Regarding qinterval, documentation in AIntervalCategory says
> ++ qinterval(inf,sup) creates a new interval without checking the
> ++ ordering on the elements.
>
> Implementation in Interval says
> qinterval(i : R, s : R) : % ==
> [roundDown i, roundUp s]
>
>
> I think in Interval we should have "qinterval(i, s) = [i, s]",
> and have a local function that uses "[roundDown i, roundUp s]".
> And choose to use one or the other wisely depending on the situation.

Let me say that relation between AInterval and Interval is not
a easy one, they have quite different motivation. Motivation
for Interval is to have exact computations, that is true result
of computation should be contained in the resulting interval.
By neccesity size of such interval depends on effort that we
ready to spent on it. Large intervals say very litte about
result, so we do not care so much if they are somewhat bigger
than necessary. OTOH we would like to limit growth of small
intervals.

As I wrote we want that result of operations contain true
result. I must say that it is not clear to me if this property
holds for Interval. Namely, "^" for Interval looks fishy,
in general floating point power can have higher error than
simple operations and then correction done by 'roundUp' and
'roundDown' is too small to ensure correctness. Due to this
I consider all "complex" operations in Interval to be highly
suspect, and hence of limited use.

--
Waldek Hebisch

[Qian Yun]

Consider this proposal:

In Interval, rename "qinterval" to "rinterval" (for round
interval), and rename "exactInterval" to "qinterval".

Yes, there are other inaccuracies in Interval, that
can be fixed later.

I propose this change to make the interface clear and
consistent.

- Qian