divexact/divides and sqrt/issquare

[wbhart]

This ticket is for discussion of the various issues related to divexact/divides and sqrt/issquare. @thofma @tthsqe12

We will not discuss the Euclidean division operator div and friends here. There is already a ticket for that issue, which is largely orthogonal.

Lets also not discuss whether 0 divides 0, which is also an orthogonal issue.

We will also not discuss square root mod composite moduli, which is beyond the scope of the current discussion.

Divexact

1. In Flint, for performance reasons, divexact does an exact division without checking or raising an exception if the division is not exact. Originally in AbstractAlgebra/Nemo divexact had the same meaning, however we eventually decided this is not good behaviour in the first instance for a CAS. This became most obvious when implementing interpolation where integer divexact must fail to indicate the interpolation has failed.

2. In AbstractAlgebra/Nemo we started implementing divides, which returns a tuple (flag, q) of a boolean to indicate whether the division was exact and a/the quotient in the case that it is and usually zero if not. We started implementing divides in Flint for the same purpose, though it isn't clear how far the rollout has gotten in C or Julia.

Square root

1. In Flint, sqrt computes the square root, usually returning 0/1 to indicate success/failure (an exception is integer sqrt and sqrt of integers mod small p, where for performance/interface reasons this is not done - Fredrik has done some work on an integer square root which tries to do a heuristic test to fail early before trying a sqrtrem). We also have the function is_square which is a simple predicate.

2. In AbstractAlgebra/Nemo we have sqrt which tries to return the square root, but raises an exception if the input is not square. We also have issquare, which is a simple predicate.

Issues

1. Integer sqrt and exact division and multivariate sqrt and exact division can be done faster/much faster if exactness is assumed but not tested.

2. Nemo uses Flint's scalar_div(exact) functions, which assume exactness, to implement ad hoc division of a polynomial/matrix by a coefficient. No test of exactness is performed.

3. There is no counterpart of divides for sqrt in AbstractAlgebra/Nemo though in some cases the underlying square root functions are implemented to support such an operation, in both Flint and AbstractAlgebra/Nemo.

4. divides does not start with the word "is", so one does not expect it to be a simple predicate returning true/false only. But it is hard to come up with a name for issquare that doesn't start with "is". Therefore one expects a simple predicate.

5. Some proposed names in AbstractAlgebra/Nemo include: sqrt_exact, issquare_with_sqrt, divexact -> div_exact, isdivisible (for the simple predicate).

6. There is a performance issue in returning a tuple (flag, q), from divides and the sqrt analogue, if one only needs the simple predicate, as an object q must always be created to hold the quotient/sqrt. This is particularly bad in cases where it is expected that exactness is unlikely and a fast heuristic test will reject negatives quickly.

7. Many of our divides functions do not do a fast heuristic test, but should.

Discussion

Mostly this ticket is about agreeing on an interface and rolling out the changes. I don't know of any cases where we don't know an algorithm, or how to implement whatever we decide on.

Note that in the worst case we end up with four functions in each case:

div_exact, divides, isdivisible, div_exact_unchecked

sqrt, ??, issquare, sqrt_unchecked

Of course div_exact could take an optional flag which allows checking to be turned off:

div_exact, divides, isdivisible

sqrt, ??, issquare

[wbhart]

One possible name for the analogue of divides could be "square" but I think I like "issquare_with_sqrt" better, even though it starts with "is".

[wbhart]

isdivisible_by may be better than isdivisible, as the latter has a specific interpretation in terms of rings.

[wbhart]

Also if we don't use sqrt_exact, I'm not too bothered about using div_exact instead of divexact.

We should also think about the interface for nth roots:

nth_root, isnth_root_with_root, isnth_root ??

[thofma]

We use the tuple return value quite extensively, also for other concepts, in Hecke. I have no objection changing it, but want to CC @fieker.

1. I am happy with the *_with_* approach, although I am sure this will not stop people from writing:

     if issquare(a); return sqrt(a); end;
   

(Would still be possible with issquare returning a tuple, but people would know that there is a function returning both, since they are using it.)

2. Here are my proposals for the names. Last two columns are the same functions, but with or without making it consistent with the rest.

Predicate only	Both	Return value or fail	Return value, assume exact
issquare(a)	issquare_with_sqrt(a)	sqrt(a)	sqrt(a, check = true)
ispower(a, n)	ispower_with_root(a, n)	root(a, n)	root(a, n, check = true)
isdivisible_by(a, b)	divides(a, b)	divexact(a, b)	divexact(a, b, check = true)

[wbhart]

The last two rows seem to be duplicates. Or are you simply giving two options?

[wbhart]

I guess is_power should be ispower under our convention of no underscore after "is".

[wbhart]

Naturally, another wrinkle is sqrt vs isqrt in Julia. However, we should assume we are taking about what the functions do inside AbstractAlgebra/Nemo. I'm not sure if it should be sqrt(:fmpz) or isqrt(::fmpz) outside Nemo.

[thofma]

Oh, right. I have updated the table. And yes, the last two are the same with different choices for the "Both" column. The divides does not really fit in, but I have a hard time coming up with a good consistent name.

[wbhart]

I strongly prefer divides out of the two options.

[wbhart]

It would have to be isdivisible_with_quotient, but this is just way too much typing. I think it is far more important for the simple predicates to be consistently named than the "both" functions.

But I'm ok with all the other suggestions.

[wbhart]

What would you call a function that checks if something is a power to any exponent (and returns the root and exponent). This is actually useful in factoring.

[wbhart]

Actually, we can still use ispower for that and just have a different signature which doesn't take the exponent.

[thofma]

Yes, I would call it ispower. We have this in Hecke for some types.

[fieker]

Just to clarify:

• I am "happy" with the names, the ".." indicates that I don't like having to change them, but a uniform scheme is important
• the assume exact, no fail: what is supposed to happen on illegal input? Any output? (sometimes VERY useful)?

[wbhart]

If you want a function which doesn't assume exact, use one of the other functions.

If you want the result checked, set checked=true.

I'm really not sure what you are asking. If the algorithm is given invalid input, anything could happen, including failure to terminate. That's what the check option is for, in case that is not an option for you.

[fieker]

On Thu, Feb 04, 2021 at 02:03:13AM -0800, wbhart wrote: If you want a function which doesn't assume exact, use one of the other functions. If you want the result checked, set checked=true. I'm really not sure what you are asking. If the algorithm is given invalid input, anything could happen, including failure to terminate. That's what the check option is for, in case that is not an option for you.

The best ispower(a::fmpz) uses ispower(a, p) to NOT check (way too slow) and return the correct result if a is a p-th power or unspecified garbage Hence the question

…

-- You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub: #862 (comment)

[wbhart]

The algorithm shouldn't be specified. I can't guarantee that some future algorithm is going to crash, not terminate, raise an exception or what upon invalid input. The point is, when the input is not checked, you cannot make any guarantees what will happen. If you need some other behaviour for a specific project, I guess you have to implement it.

[wbhart]

In the table above, is ispower(a) actually intended? Or should that be ispower(a, n)?

[thofma]

It is a typo. It should be ispower(a, n).

[fieker]

The existence of ispower(a) is independent?

[wbhart]

I just realised by "check=false" you probably meant a default argument, not a named argument, right?

[wbhart]

My branch for the divides related functions will be called divides_functions.

(Some docs will have to wait until the new documentation is merged, otherwise there will be massive merge conflicts.)

[wbhart]

I am not convinced check=false is the right default. We used to have the convention that because it is a computer algebra system the user needs to (by default) get feedback if their input was garbage.

Now we are effectively going through and systematically reversing that decision. And I think that has only happened because we wanted to rely on the fact that sqrt in GMP returns the truncated integer square root and setting check=true will break all those uses.

I'm literally changing hundreds of cases in the code to handle the default differently across all these functions. Admittedly we need to go through and check every individual usage of functions like divrem to see whether it should be checking or not, but the default check=false feels wrong to me.

Unfortunately I've already spent days changing stuff to check=false as the default. What worries me the most is the use of divexact in divides and then divides having fallbacks in terms of divexact. Both seem problematic to me.

Another thing that worries me is that I'm changing a lot of code like this:

flag, q = divides(a, b)
!flag && error("not an exact division"

to

if check
   flag, q = divides(a, b)
   !flag && error("not an exact division")
else
   q = divexact(a, b)
end

I would change it to

q = divexact(a, b; check=true)

except that our current fallbacks are back-to-front if I do this, and in some cases this occurs in code implementing divexact, etc. It's hard to then get the test code to pass.

Overall I think fixing divexact across the whole system is going to be an enormous amount of work. There are so many hacks to work around....

[tthsqe12]

default should be checked=true in my opinion

[wbhart]

Maybe I'm wrong. We have a lot of cases where we are taking the gcd of two things and then dividing by the gcd, in which case you maybe want the default to be false. I don't think we can win either way.

For the user, the default should be true though I reckon. Which probably means I have to redo all the work I just did over the last week.

[tthsqe12]

that is why I advocate for a gcd_with_cofactors function

[wbhart]

Yeah see in a lot of places we have a divides function implemented in full. Then we have something like the following:

function divexact(a::ResElem{T}, b::ResElem{T}; check::Bool=false) where {T <: RingElement}
   check_parent(a, b)
   fl, q = divides(a, b)
   !fl && error("Impossible inverse in divexact")
   return q
end

All of these are going to have to be dealt with differently if we want to take advantage of the option to do faster divexact.

[wbhart]

The main disadvantage of check=true is we have to check the code very carefully. For example one doesn't want divides or isprime or some other predicate raising an exception because of a non-exact division (though it would probably indicate a bug anyway).

[wbhart]

I got the tests to pass with check=false. But if @thofma and @fieker don't object, I think I am going to go back and change the default to check=true for root, sqrt and divexact. The new default simply doesn't make any sense.

The downsides will be:

• sqrt will no longer give the truncated sqrt of integers, you will have to switch all such instances to isqrt instead (you have good test code to detect this right?)
• it's a fair bit of work to go back and change this now
• it will be slower until someone writes some faster unchecked divexact functions (especially important for many cases where it is currently implemented in terms of divides)
• it will be slower until someone goes through and sets check=false by hand in all cases where we know for sure the division is exact

[thofma]

I am happy with the default being check = true. I can handle the changes in Hecke and Oscar.

P.S.: At the moment sqrt(::fmpz) is not giving the truncated square root of integers. I probably misunderstood the first comment.

[wbhart]

Possibly the Flint function doesn't do this. I haven't looked at Nemo yet, only AbstractAlgebra. Unfortunately this is a huge job!

[wbhart]

I have finished all the AbstractAlgebra changes. The default is now check=true.

Here are some hints for the Oscar/Hecke work:

• check=false does not necessarily mean no exception will be raised, we just try to remove as many checks as possible
• sometimes a divides function calls a divexact function. If the inputs are in a field, divexact can take check=false since you know the division is possible in a field
• sometimes there is code of the form

  flag, q = divides(a, b)
  !flag && error("Not an exact division")
  

  This can now be replaced with

  q = divexact(a, b; check=check)
  

  assuming of course that this does not occur in a divexact function with a and b as inputs (else a stack overflow will occur)
• Quite a few docstrings will currently say things like "assumes exact division". Those need changing when the check flag is added.

Later we will need to check all the sqrt functions as some of them use divexact. But because I've kept the PR's separate I can't pass the check=check flag to the divexact calls.

Eventually we will also be able to use issquare_with_sqrt in various places (usually in other sqrt code). But this is only possible at the moment for integers, as it isn't implemented for all the other types (except for fallbacks, which will be slow).

I don't think there is a fast divexact for Generic.MPoly implemented at the moment, and there are huge savings possible there.

Documentation for divexact in the ring.md and ring_interface.md files awaits merging of our rearranged docs. Adding them now will simply result in merge conflicts.

[wbhart]

@thofma Ok, looks like sqrt(::fmpz) was already set up to raise an exception if it wasn't a square. That's different to the behaviour we had in AbstractAlgebra. That's good news. It means you can ignore my first comment.

[wbhart]

There are now AbstractAlgebra and Nemo PR's for each of

• iroot (which also adds the check=true flag to root)
• sqrt_functions
• divides_functions

These

• implement the entire interface above for integers
• add the check=true flags for every instance of the functions in the third column in the table above

Many of these functions do not honour the check=true flag, e.g. many Flint functions do not actually check if the division was exact, e.g. the scalar_divexact functions and there is no option to make them do so currently. But everywhere I could, I made the functions honour the flag.

Moreover, there are hundreds or even thousands of invocations of divexact which will need to be inspected by hand to see if they should be checking or not. In many cases we know mathematically that the division is exact and do not need to check. This is obviously a huge project and won't be done in these PR's. I propose we deal with this on the Flint side first (I'll bring it up at the next Flint Friday's meeting) and then deal with it module by module on the Nemo/Oscar side.

[thofma]

I will have a look today and see how much work it requires for the rest of the ecosystem (just getting it to work with the default check = true everywhere).

[wbhart]

Yeah I thought it might work to add generic fallbacks for all the check=true versions that fall back to the versions without. That way when the generic code calls the version with check=true it will still keep working.

When I did this, I easily got stack overflows. I think this is because of generic fallbacks we already have, but it also might just be that I had bugs in the code I wrote.

[wbhart]

I'll try to add generic fallbacks to each of the AbstractAlgebra PR's now and see what happens.

[thofma]

Ah, you mean a fallback like divexact(x::T, y::T; check = true) where {T} = divexact(x, y)?

[wbhart]

yeah, I'd probably restrict to T <: RingElem since I think we definitely get conflicts if not. But yes, basically that.

[wbhart]

The main problem is there are very, very many adhoc divexacts of the form:

divexact(x::RingElem, y::Integer; check::Bool=true)
divexact(x::RingElem, y::Rational; check::Bool=true)
divexact(x::RingElem, y::fmpz; check::Bool=true)
divexact(x::RingElem, y::fmpq; check::Bool=true)
divexact(x::SomeRingElem{T}, y::T; check::Bool=true)

I can't really add fallbacks for all those. But maybe at least the simple version of the fallback already helps, though it could just cover up bugs.

[wbhart]

Actually, on second thoughts, let's not add fallbacks for all these check=true functions. It just covers up cases we missed.

I did add fallbacks for issquare_with_sqrt and other new functions that currently only exist for Integers and a handful of other types. That seems to be less problematic.

I did find some cases where generic functions we added were being called and there was actually test code which was trying to test functionality that hadn't been implemented, and yet the functions actually existed in Flint and were faster. So generics do have the problem that they cover up problems.

[thofma]

Agreed.

[thofma]

How should we move forward with merging the PRs? There are

1. the breaking documentation changes that we would like to get in before next week and
2. the breaking divexact/root/etc changes, which need flint 2.8.

Should we get in 1. today and make a new breaking release? Or should we wait till (tomorrow?) for the next flint release and try to get 1. and 2. in?

[wbhart]

The missing check=true flags are now in. There weren't quite hundreds I'd missed, but still a lot.

I think documentation is a priority for next week so I think we should do that as soon as possible, as we don't want to be trying to sort out a documentation release during the week before anyone can start work.

The other PR's can wait. Some of them touch the docs for various reasons, but it should be a relatively straightforward rebase so long as we don't massively rearrange the doc files again next week.

[wbhart]

Oh apparently I never implemented ispower or ispower_with_root for integer types in either AbstractAlgebra or Nemo I think. That's good because it would have failed until flint-2.9 anyway. I'll open a separate ticket for it just so it is not forgotten.

[wbhart]

@fieker Regarding ispower(a), we can implement this using GMP's mpz_perfect_power_p [1], but the function doesn't tell you what power it is.

We implemented such a function in flint by pulling Jason Moxham's 2008 implementation out of MPIR [2]. However, this has the problem of not guaranteeing to return either the smallest or largest power. Such a thing is really quite a lot harder.

And even if we had it for fmpz, we could not easily implement it for mpz_t or Int. Unfortunately, until someone finds time to add these functions to Flint, we won't get anywhere with it.

Moreover, these really need to be implemented at the mpn level for efficiency, and that is hard or even impossible to do in Flint itself. It would really want to be in GMP, and it isn't straightforward to get something accepted there due for example, to the extreme optimisation level and maturity of the library.

[1] https://gmplib.org/manual/Integer-Roots
[2] https://github.com/wbhart/flint2/blob/trunk/fmpz/is_perfect_power.c

[fieker]

There are different problems.

• in Hecke we have ispower, a homebrew combination using the above sources. There is no choice, it is needed in several algorithms.
• in Oscar/examples there is a ((highly) optimised) implementation using p-adic arithmetic for the n-th roots. And please, before commenting that this is inefficient, slow, ... have a look. For large numbers it is MUCH MUCH faster than the above sources
• as long as the mpn-layer is documented, there is no problem even using it from Julia (I was too lazy)
• it is desirable to have this in flint, however not necessary for having this in Oscar and Hecke (or even Nemo (or AA))
• the fast O(n log n...) method of Bernstein/ Lenstra is based on the fast coprime base - which is also not in flint. This too would be nice to have. It will not replace the one in Hecke as the Hecke one is generic, it works for "all" euclidean rings (or it should; it might be resticted to rings with gcd!, lcm! and a few others. The other foundation would be easy to port to flint, it is a function to compute the exact correct p-th root of n (p prime) and is allowed to return garbage if this is not possible. Testing the root would make this approach O(n^2) - too slow. The flint comment even mentions the p-adic lifting (I actually use 2-adic and 2^64 adic lifting)
• there is the famous quote by Voltaire: “Perfect is the enemy of good.”