Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add support for ArbPoly #32

Merged
merged 26 commits into from
Sep 8, 2020
Merged

Add support for ArbPoly #32

merged 26 commits into from
Sep 8, 2020

Conversation

Joel-Dahne
Copy link
Collaborator

So far I have just added one commit with

  • ArbPoly and arb_poly_struct types
  • Support for parsing arb_poly_t
  • Parsed methods from arb_poly.h

One of the main things to discuss is what interface to use for ArbPoly. There is Polynomials.jl but I think that might only work for polynomials with a parametrized type. We could maybe still take inspiration from it and possibly overload some of their methods for ArbPoly. Do you have any other ideas?

Another thing is that I want to add a dedicated ArbSeries type, basically an ArbPoly with a specified order. Something like

struct ArbSeries
    arb_poly::arb_poly_struct
    order::Int
    prec::Int

    ArbPoly(order::Int; prec::Integer = DEFAULT_PRECISION[]) = new(arb_poly_struct(), order, prec)
end

It would work identically to ArbPoly when it comes to interfacing directly with functions from Arb. The difference would be that we overload different methods for it. For example * would for ArbPoly use arb_poly_mul whereas for ArbSeries it would use arb_poly_mullow, for ArbPoly the length would depend on the highest non-zero term whereas for ArbSeries it would be fixed to the given degree. There are a number of differences like that. Here I'm also not sure what interface to use, we could maybe take inspiration from TaylorSeries.jl.

@kalmarek
Copy link
Owner

I'm pretty amazed by how fast we can bootstrap a new Arb type with almost all Arb functionality wrapped!

deps/build.jl Outdated
@@ -0,0 +1 @@

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

don't re-add build.jl ;)

arbcall"slong arb_poly_allocated_bytes(const arb_poly_t x)"

### Basic manipulation
arbcall"slong arb_poly_length(const arb_poly_t poly)"
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This needs an import Base: length at the top

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There will be a similar function in acb_poly.jl as well. I'm still thinking about how to handle this in a good way.

arbcall"void arb_poly_evaluate2(arb_t y, arb_t z, const arb_poly_t f, const arb_t x, slong prec)"
arbcall"void _arb_poly_evaluate2_acb_horner(acb_t y, acb_t z, arb_srcptr f, slong len, const acb_t x, slong prec)"
arbcall"void arb_poly_evaluate2_acb_horner(acb_t y, acb_t z, const arb_poly_t f, const acb_t x, slong prec)"
arbcall"void _arb_poly_evaluate2_acb_rectangular(acb_t y, acb_t z, arb_srcptr f, slong len, const acb_t x, slong prec)"
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here we should probably change the code gen to produce a function which also accepts an ArbVector for the two arguments arb_srcptr f, slong len

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My plan here was to later implement methods which are not in-place for which we add more conveniences, like automatically picking up the length from AcbVector and such. It would be great if we could change the code-gen to automatically make length a keyword argument which takes the length of AcbVector as default value, but I'm worried this would be hard to do. The length parameter is sometimes called len, sometimes flen or glen, so it's hard to automatically identify it.

Copy link
Collaborator

@saschatimme saschatimme Aug 13, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Makes sense, but note that I added some code exactly doing this

It would be great if we could change the code-gen to automatically make length a keyword argument which takes the length of AcbVector as default value

in #27 for the vector types 🙃

@Joel-Dahne Joel-Dahne marked this pull request as ready for review September 2, 2020 09:41
@Joel-Dahne
Copy link
Collaborator Author

Now I have finally taken the time to finish this pull request. I have added support for ArbPoly, ArbSeries, AcbPoly and AcbSeries. So far the implemented methods are quite basic and there is a lot of convenience methods to add later on, but I think this is a good start. I have added tests for everything I have added and also cleaned up some of the other tests a little bit.

Later on I might do some rearranging of code, moving things from poly.jl to series.jl, but that can wait for now.

@saschatimme saschatimme mentioned this pull request Sep 4, 2020
6 tasks
@Joel-Dahne Joel-Dahne mentioned this pull request Sep 4, 2020
@Joel-Dahne
Copy link
Collaborator Author

Now I have merged this with the updated arbcall. The merge went very smooth, I really like the changes you did to arbcall. The code in the pull request is mostly self contained so it should not be to dangerous to merge it. I have added tests for all new functionality as well as some old functionality. Once this is merged I'll look into parsing more of the Arb header files.

One thing I noted is that now that the length of the vectors become key-word arguments we have some clashes. Specifically void _arb_vec_add(arb_ptr C, arb_srcptr A, arb_srcptr B, slong n, slong prec) and void _arb_poly_add(arb_ptr C, arb_srcptr A, slong lenA, arb_srcptr B, slong lenB, slong prec) parses to the same method. My solution now was to comment out the latter. If you do want to add them as polynomials you can always do that by using the length of the shorter one, so I don't see it as a problem for now.

src/poly.jl Outdated

Base.@propagate_inbounds function Base.setindex!(
poly::Union{ArbPoly,ArbSeries},
x::Arb,
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
x::Arb,
x::ArbLike,

src/poly.jl Outdated

Base.@propagate_inbounds function Base.setindex!(
poly::Union{AcbPoly,AcbSeries},
x::Acb,
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
x::Acb,
x::AcbLike,

Copy link
Collaborator

@saschatimme saschatimme left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good! Some comments are nitpicky, but I think we should definitely widen the setindex! methods and add a ref version for getindex.

src/poly.jl Outdated
end

# Constructors
for (TPoly, T) in [(:ArbPoly, :Arb), (:AcbPoly, :Acb)]
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
for (TPoly, T) in [(:ArbPoly, :Arb), (:AcbPoly, :Acb)]
for (TPoly, T) in [(:ArbPoly, :ArbLike), (:AcbPoly, :AcbLike)]

src/poly.jl Outdated
end

@eval function $TPoly(
coeffs::AbstractVector{$T};
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
coeffs::AbstractVector{$T};
coeffs::AbstractVector{<:$T};

end
end

for (TSeries, T) in [(:ArbSeries, :Arb), (:AcbSeries, :Acb)]
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Prob. the same changes as for the poly case above

res = eltype(T)(prec = precision(poly))
get_coeff!(res, poly, i)
return res
end
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think should also add a ref(poly:T, i) method which returns an ArbRef/AcbRef

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This was my initial plan, but it's slightly complicated. Arb allows getting the value of any index i >= 0, for i higher than any set index it will just give the value 0. This means that while poly[i] will be defined for any i >= 0 it might not correspond to an actual arb_struct in memory.

In practice my experience is that you don't change polynomials element-wise very often, or at least it's definitely not where most of the computational cost lies, so it's not such a big problem.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok that's a fair point :)

src/poly.jl Outdated
Comment on lines 109 to 113
function Base.one(::Type{T}) where {T<:Union{ArbPoly,AcbPoly}}
res = T()
one!(res)
return res
end
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
function Base.one(::Type{T}) where {T<:Union{ArbPoly,AcbPoly}}
res = T()
one!(res)
return res
end
Base.one(::Type{T}) where {T<:Union{ArbPoly,AcbPoly}} = one!(T())

We can do these way shorter by now :)

@Joel-Dahne
Copy link
Collaborator Author

I'll wait for #62 to be merged with master and then start updating things!

@Joel-Dahne
Copy link
Collaborator Author

I merged it with the updates to master and took care of your comments! The only thing I have not taken care of is adding a ref method, but as I explained this is not really applicable.

@saschatimme saschatimme merged commit 7a57f32 into master Sep 8, 2020
@kalmarek kalmarek deleted the poly branch July 21, 2022 15:01
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

3 participants