Skip to content

Conversation

@jagot
Copy link
Contributor

@jagot jagot commented Nov 16, 2019

This helper struct simplifies calls to Arb functions that expect vector inputs. There were no tests for the DFT routines, so I am not sure if they have been broken or not, but they should not be.

@jagot
Copy link
Contributor Author

jagot commented Nov 16, 2019

(I wanted this for a wrapper of the acb_hypgeom_coulomb_jet function, which will come in a separate PR).

Copy link
Owner

@JeffreySarnoff JeffreySarnoff left a comment

Choose a reason for hiding this comment

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

We cannot offer AcbVector without also offering ArbVector and ArfVector.
I have not reviewed the code yet.

@jagot
Copy link
Contributor Author

jagot commented Nov 16, 2019

What's the difference between these types?

@JeffreySarnoff
Copy link
Owner

I missed the T parameter -- The name AcbVector implies the vector contains ArbComplex elements.
And given the way other parts of the package interact, making three sorts of vector, one for each sort of element ArbFloat, ArbReal, and ArbComplex would be more flexible (even if they were just different parameterizations of an underlying ArblibVector{T}. Often it is necessary to map args of one sort to another for complete coverage of some functionality.

@jagot
Copy link
Contributor Author

jagot commented Nov 16, 2019

I just realized. There is no arf_vec_{init,clear}, however. ArbVector and AcbVector could be made type aliases of the underlying type, is that enough?

So const ArbVector = AxbVector{ArbReal} and const AcbVector = AxbVector{ArbComplex}.

@JeffreySarnoff
Copy link
Owner

The library offers _acb_vec_init, _arb_vec_init but does not offer _arf_vec_init (similarly for other _vec_ functions). This happens with other functions too. I have adopted the policy of always supporting functionality for each of the three main types and I have done it like this:

function acb_vec_init(..)
   result = zero(ArbComplex)
   ccall(..)
   return result
end

function arb_vec_init(..)
   result = zero(ArbReal)
   ccall(..)
   return result
end

function arf_vec_init(x,..)
   xx = ArbReal(x)
   res = arb_vec_init(xx,..)
   result =  ArbFloat.(res)
   return result
end

see https://github.com/JeffreySarnoff/ArbNumerics.jl/blob/master/src/float/elliptic.jl
and https://github.com/JeffreySarnoff/ArbNumerics.jl/blob/master/src/float/bessel.jl

@jagot
Copy link
Contributor Author

jagot commented Nov 16, 2019

Can you just freely cast a pointer to an ArbReal to an ArbFloat?

@JeffreySarnoff
Copy link
Owner

JeffreySarnoff commented Nov 16, 2019

no -- one must explicitly cast values of each type to the other, target type. Otherwise bad things happen.

@jagot
Copy link
Contributor Author

jagot commented Nov 16, 2019

Hmm, reading the API documentation, there seems to be no functions accepting a vector of ArbFloats, indeed, there is no arf_ptr, whereas both arb_ptr and acb_ptr are defined.

When would a ArblibVector{ArbFloat{P}} be needed?

EDIT: Typo

@JeffreySarnoff
Copy link
Owner

The idea of ArbNumerics is that users have a Julia-centric way of working with arblib features and that clients should be free to program Julia using an Arb extended functional vocabulary. The reason to support Array{ArbFloat{P}, 1} is that people will use it to simplify their own work once any of the three types has vector-of support available. using LinearAlgebra; dot(...) for example.

@jagot
Copy link
Contributor Author

jagot commented Nov 16, 2019

Then I think it's easiest to just add a helper function, that automatically promotes a vector of ArbFloats to a vector of ArbReals, on the Julia side:

ArblibVector(v::AbstractVector{ArbFloat{P}}) where P =
    ArblibVector(ArbReal.(v))

This works like so:

julia> f = rand(ArbFloat{128}, 5)
5-element Array{ArbFloat{128},1}:
 0.5953977471776002041348507759413
 0.4982507430703087645009724322786
 0.7010188632230956411011589105102
 0.5026129464948069049729481221645
 0.4253157939327851802021822646549

julia> ArblibVector(v::AbstractVector{ArbFloat{P}}) where P =
           ArblibVector(ArbReal.(v))
ArblibVector

julia> F = ArblibVector(f)
5-element ArblibVector{ArbReal{128}}:
 0.5953977471776002041348507759413
 0.4982507430703087645009724322786
 0.7010188632230956411011589105102
 0.5026129464948069049729481221645
 0.4253157939327851802021822646549

Converting back to a Vector{ArbFloat{P}} already works:

julia> Vector(F)
5-element Array{ArbReal{128},1}:
 0.5953977471776002041348507759413
 0.4982507430703087645009724322786
 0.7010188632230956411011589105102
 0.5026129464948069049729481221645
 0.4253157939327851802021822646549

julia> Vector{ArbFloat{128}}(F)
5-element Array{ArbFloat{128},1}:
 0.5953977471776002041348507759413
 0.4982507430703087645009724322786
 0.7010188632230956411011589105102
 0.5026129464948069049729481221645
 0.4253157939327851802021822646549

@JeffreySarnoff
Copy link
Owner

I understand your approach. It is concise. And it would be workable if it existed behind the scenes, rather than as the "way to do it", because the way do it must be syntactically similar for each of the types to keep harmony with the way all of the rest of the package works (and exports for use). Also, it does not automagically allow the correct type to flow through .. a specific correction, maybe outside of the logical flow where the initial conversion occurred must be attached. In all other functional use, any type facilitating conversion -> operations -> type recovering conversion flow happens inside of the client's functional call, so no mistyping can occur accidentally. This is both nice as an interface standard and important when working with a C library that manipulates C structs internally in a way that requires Julia follow along properly.

I lived this when it became necessary to provide Matrix types and Julia <-> Arb interconversion. Vectors are easier, but the constraints and considerations are much the same. A look at
https://github.com/JeffreySarnoff/ArbNumerics.jl/blob/master/src/libarb/ArbFloatMatrix.jl and how it relies on https://github.com/JeffreySarnoff/ArbNumerics.jl/blob/master/src/libarb/ArbRealMatrix.jl
shows the excessive detail necessary to keep this all working correctly. That implementation is the result of many near misses where I tried to favor simplicity first.

So your approach will work, and still we need a smoother approach. I will post some scaffolding here later on. And that you are more than welcome to alter or reshape.

@jagot
Copy link
Contributor Author

jagot commented Nov 16, 2019

I'm afraid you've lost me. I never intended for the ArblibVector to be a user-facing type, only to help wrapping the C API, reducing the risk of memory leaks. Here is how I intended to use them for wrapping e.g. the Coulomb functions:

function coulomb!(F::AcbVector{P}, G::AcbVector{P},
                  H⁺::AcbVector{P}, H⁻::AcbVector{P},
                  λ::ArbComplex{P}, η::ArbComplex{P}, z::ArbComplex{P};
                  prec::Int=P) where P
    len = F.len
    @assert G.len  len
    @assert H⁺.len  len
    @assert H⁻.len  len
    ccall(@libarb(acb_hypgeom_coulomb_jet), Cvoid,
          (Ref{ArbComplex{P}}, Ref{ArbComplex{P}}, Ref{ArbComplex{P}}, Ref{ArbComplex{P}},
           Ref{ArbComplex}, Ref{ArbComplex}, Ref{ArbComplex}, Clong, Clong),
          F.ptr, G.ptr, H⁺.ptr, H⁻.ptr, λ, η, z, len, prec)
    F,G,H⁺,H⁻
end

# Taylor series version
function coulomb::ArbComplex{P}, η::ArbComplex{P}, z::ArbComplex{P}, len::Int; kwargs...) where P
    F = ArblibVector(ArbComplex{P}, len)
    G = ArblibVector(ArbComplex{P}, len)
    H⁺ = ArblibVector(ArbComplex{P}, len)
    H⁻ = ArblibVector(ArbComplex{P}, len)

    coulomb!(F,G,H⁺,H⁻,λ,η,z;kwargs...)

    take!(F),take!(G),take!(H⁺),take!(H⁻)
end

# Scalar version
function coulomb::ArbComplex{P}, η::ArbComplex{P}, z::ArbComplex{P}; prec::Int=P) where P
    F = ArbComplex{P}()
    G = ArbComplex{P}()
    H⁺ = ArbComplex{P}()
    H⁻ = ArbComplex{P}()
    ccall(@libarb(acb_hypgeom_coulomb), Cvoid,
          (Ref{ArbComplex}, Ref{ArbComplex}, Ref{ArbComplex}, Ref{ArbComplex},
           Ref{ArbComplex}, Ref{ArbComplex}, Ref{ArbComplex}, Clong),
          F, G, H⁺, H⁻, λ, η, z, prec)
    F,G,H⁺,H⁻
end

coulomb(λ, η, z, args...; kwargs...) =
    coulomb(ArbComplex(λ), ArbComplex(η), ArbComplex(z), args...; kwargs...)

I am of course happy to adapt the code to your preferences, but I don't immediately see what's missing.

@JeffreySarnoff
Copy link
Owner

Thank you for the clarification. We have discussing two different things.

@orkolorko
Copy link
Contributor

Just my 2 cents.

While working on the DFT, I saw a couple of things:

  • one of these is that convert(Array{ArbComplex}) does not return a well formed pointer to an array of acb.
  • it seems that letting most of the memory management happen in Julia space is much safer but the issue above makes necessary some workaround, that's why in the older versions I intoduced a ArbComplexVector also

@JeffreySarnoff
Copy link
Owner

You now understand this quite well -- if you like this as is, I will merge it,
if there is more to be done for it to work as you intend, just let me know.

@JeffreySarnoff
Copy link
Owner

It has been a while -- "one of these is that convert(Array{ArbComplex}) does not return a well formed pointer to an array of acb", what does it return (and I assume you are referring to a Matrix)?
The fact that Julia is column-oriented and C is row-oriented added much complexity to getting anything to happen with Matricies without segfaults or lesser errors.

@jagot
Copy link
Contributor Author

jagot commented Nov 18, 2019

@orkolorko Could you provide a unit test for the DFT and its inverse? Such that we know that it works as it should?

I'll make a PR (as soon as this PR is merged) for the Coulomb routines I posted above, along with unit tests, and then we would have two different cases of real usage of this PR.

@JeffreySarnoff
Copy link
Owner

good work, nice to see dot in the mix.

All functions from LinearAlgebra that had been avoided were omitted after the inappropriateness of some fallthroughs when used with Arb types became clear. There has been improvement in LinearAlgebra, which along with some alternative versions of very low level float functionality implemented in this package, suggests revisiting additional vector-related ops. The Arb polynomial routines are not available because providing a Julian approach to each next Arb structured type requires a good deal of time/testing/automatic-interconversion-provision.

@orkolorko
Copy link
Contributor

@orkolorko Could you provide a unit test for the DFT and its inverse? Such that we know that it works as it should?

I'll make a PR (as soon as this PR is merged) for the Coulomb routines I posted above, along with unit tests, and then we would have two different cases of real usage of this PR.

Yes, I will provide a unit test

Hmm, reading the API documentation, there seems to be no functions accepting a vector of ArbFloats, indeed, there is no arf_ptr, whereas both arb_ptr and acb_ptr are defined.

When would a ArblibVector{ArbFloat{P}} be needed?

EDIT: Typo

I tink the reason why is that in ARB you're not really supposed to use ArbFloats directly: ARB is a library made for rigorous and efficient numerical types, with ball representation.
Quoting ARB`s homepage: " Arb is a C library for rigorous real and complex arithmetic with arbitrary precision. Arb tracks numerical errors automatically using ball arithmetic, a form of interval arithmetic based on a midpoint-radius representation."

@JeffreySarnoff
Copy link
Owner

For the Julian perspective, we manually create the functional equivalents that allow the user to work as if the Arb C library better supported the arf type.

@JeffreySarnoff
Copy link
Owner

is this ready to merge?

@jagot
Copy link
Contributor Author

jagot commented Dec 4, 2019

From my point-of-view, yes, barring the fact that there still are no unit tests ensuring that the (I)DFT routines work as they should.

@JeffreySarnoff
Copy link
Owner

JeffreySarnoff commented Dec 4, 2019

Please jot down two tests each for dft and idft here -- I need to give examples of use.

@jagot
Copy link
Contributor Author

jagot commented Dec 5, 2019

Please jot down two tests each for dft and idft here -- I need to give examples of use.

Well, that would be up to @orkolorko. Although I've used DFTs a long time ago, I have no example at hand presently.

@JeffreySarnoff
Copy link
Owner

@orkolorko I need to provide others with a few how to use this examples.

@orkolorko
Copy link
Contributor

orkolorko commented Dec 5, 2019

I'm sorry it is taking me some time, it is the end of the semester and I'm in grading swamp. How should I do, clone your branch @jagot and then work on it and then send you a pull-request?
When I fork I get @JeffreySarnoff's repository, not yours

What is the right process to do this? Should I add the unit test to my branch and then send a pull request?

@jagot
Copy link
Contributor Author

jagot commented Dec 5, 2019

I think the absolute easiest is to just take an example you already have and paste it here in the PR conversation, and then me or @JeffreySarnoff can turn it into a test. E.g. something like this:

x = ...
y = ArbNumerics.dft(x)
@test y  reference_solution
x̃ = ArbNumerics.inverse_dft(y)
@test x

The last step assumes that ArbNumerics.inverse_dft(ArbNumerics.dft(x)) ≈ x, which could be wrong if there is some normalization factor (like √(2π)). The test should be repeated for the various possible input types for dft and inverse_dft.

@orkolorko
Copy link
Contributor

orkolorko commented Dec 5, 2019

These were some tests I checked by hand

dft([ArbComplex(1.0), ArbComplex(2.0), ArbComplex(3.0), ArbComplex(4.0)])== [ArbComplex(10.0), ArbComplex(-2.0,+2.0), ArbComplex(-2.0),ArbComplex(-2.0,-2.0)]
i.e. FFT(1,2,3,4)=(10,-2+2im, -2, -2-2im)

inverse_dft(dft([ArbComplex(1.0), ArbComplex(2.0), ArbComplex(3.0), ArbComplex(4.0)]))== [ArbComplex(1.0), ArbComplex(2.0), ArbComplex(3.0),ArbComplex(4.0)]
i.e. inverse(FFT(1,2,3,4))==(1,2,3,4)

dft([ArbComplex(1.0), ArbComplex(0), ArbComplex(0), ArbComplex(0),ArbComplex(0), ArbComplex(0), ArbComplex(0), ArbComplex(0) ])==[ArbComplex(1.0), ArbComplex(1.0), ArbComplex(1.0), ArbComplex(1.0), ArbComplex(1.0), ArbComplex(1.0), ArbComplex(1.0), ArbComplex(1.0)]
i.e. FFT(1,0,0,0,0,0,0,0)=(1,1,1,1,1,1,1,1)

dft([ArbComplex(0,1), ArbComplex(0), ArbComplex(0), ArbComplex(0),ArbComplex(0), ArbComplex(0), ArbComplex(0), ArbComplex(0) ])==[ArbComplex(0,1.0), ArbComplex(0,1.0), ArbComplex(0,1.0), ArbComplex(0,1.0), ArbComplex(0,1.0), ArbComplex(0,1.0), ArbComplex(0,1.0), ArbComplex(0,1.0)]
i.e. FFT(1im, 0,0,0,0,0,0,0)=(1,1,1,1,1,1,1,1)*im

I will prepare some examples with the dft of samples of trigonometric functions to show the use case.

@JeffreySarnoff
Copy link
Owner

much appreciated -- what you prepare, I will incorporate into the docs

@jagot
Copy link
Contributor Author

jagot commented Dec 6, 2019

Thanks @orkolorko, those tests enabled me to find some bugs in the DFT wrappers.

@JeffreySarnoff
Copy link
Owner

I am merging this so adding to the documentation will be easier. When that is done, a new release.

@JeffreySarnoff JeffreySarnoff merged commit 565bba2 into JeffreySarnoff:master Dec 6, 2019
@JeffreySarnoff
Copy link
Owner

@orkolorko @jagot Something is not correct with the merged #master:

julia> using ArbNumerics
julia> using Test
julia> using LinearAlgebra
julia> import ArbNumerics: ArblibVector, take!, free!

julia> T = [ArbReal{128},ArbComplex{128}]
2-element Array{DataType,1}:
 ArbReal{128}
 ArbComplex{128}

julia> v = rand(T, 5)
5-element Array{DataType,1}:
 ArbComplex{128}
 ArbComplex{128}
 ArbReal{128}
 ArbComplex{128}
 ArbReal{128}

julia>  vv = ArblibVector(v)
ERROR: MethodError: no method matching init(::Type{DataType}, ::Int64)
Closest candidates are:
  init(::Type{T}, ::Integer) where T<:ArbReal at C:\Users\jas\.julia\packages\ArbNumerics\Yqi2b\src\support\ArblibVector.jl:4
  init(::Type{T}, ::Integer) where T<:ArbComplex at C:\Users\jas\.julia\packages\ArbNumerics\Yqi2b\src\support\ArblibVector.jl:4
Stacktrace:
 [1] ArblibVector(::Type{DataType}, ::Int64) at C:\Users\jas\.julia\packages\ArbNumerics\Yqi2b\src\support\ArblibVector.jl:27
 [2] ArblibVector(::Array{DataType,1}) at C:\Users\jas\.julia\packages\ArbNumerics\Yqi2b\src\support\ArblibVector.jl:48
 [3] top-level scope at REPL[7]:1

@jagot
Copy link
Contributor Author

jagot commented Dec 6, 2019

@JeffreySarnoff Well, there is no C library function for initializing vectors of mixed real/complex numbers. As you see in your REPL, eltype(v)=DataType, so there is no way Arb could know how to handle that.

As I said above, I don't think ArblibVector is useful for users of ArbNumerics.jl, rather for implementing wrappers arounds Arb functions, which will only ever accept vectors of either real or complex elements (at least, I think so).

EDIT: I read too quickly, it seems you're trying to create an ArblibVector of data types, whereas it is supposed to be a container for numbers of one data type. Example usage:

julia> using ArbNumerics

julia> using LinearAlgebra

julia> import ArbNumerics: ArblibVector, take!, free!

julia> v = ones(ArbReal{128}, 5)
5-element Array{ArbReal{128},1}:
 1.0
 1.0
 1.0
 1.0
 1.0

julia> vv = ArblibVector(v)
5-element ArblibVector{ArbReal{128}}:
 1.0
 1.0
 1.0
 1.0
 1.0

# Be sure to free! the allocated memory, either by free!, or take! which also returns the contents:
julia> free!(vv)
0

julia> vv
0-element ArblibVector{ArbReal{128}}

julia> vv = ArblibVector(v)
5-element ArblibVector{ArbReal{128}}:
 1.0
 1.0
 1.0
 1.0
 1.0

julia> take!(vv)
5-element Array{ArbReal{128},1}:
 1.0
 1.0
 1.0
 1.0
 1.0

julia> vv
0-element ArblibVector{ArbReal{128}}

@orkolorko
Copy link
Contributor

orkolorko commented Dec 6, 2019

I think something like this could enter the documentation of the function

The DFT permits is often used as a tool to analyze a collection of samples as a sum of periodic signals DFT on Wikipedia.
ArbNumerics implements a rigorous DFT, this means that, given a set of samples as complex balls, the DFT of these samples is enclosed in the result of Arb DFT (each component has a midpoint and a radius).

To see how it works, and how it can be used, we can look at its behavior on a sample from a periodic signal

julia>n = 8;
julia>v = ArbComplex.([im*2*j*pi/n for j in 0:n-1]);
julia>dft(exp.(v))
8-element Array{ArbComplex{128},1}:
   -2.956558911618339712192e-16 + 1.2246467991473535393e-16im
                          8.0 - 8.572527594031472240582e-16im
  2.956558911618339756119e-16 + 1.224646799147352815152e-16im
   1.224646799147353289708e-16 + 1.22464679914735302725e-16im
   5.07265313323633507716e-17 + 1.224646799147353115104e-16im
                  1.31229e-32 + 1.224646799147353177226e-16im
  -5.07265313323633251691e-17 + 1.224646799147353239348e-16im
 -1.224646799147353064744e-16 + 1.224646799147353327202e-16im

The result of the DFT (ignoring the small terms) is the vector with value 8 (which divided by n is 1) at the first component, this tells us that the original signal is a periodic signal generated by $e^{2\pi i x}$.
Please remember that $e^{2\pi i k x}=\sin(2\pi k x)+i \cos(2\pi k x)$, so we are looking at a periodic function.

In the same way

julia>dft(exp.(3*v)) 
8-element Array{ArbComplex{128},1}:
      -1.521795939970898986998e-16 + 3.673940397442060090777e-16im
       -3.67394039744205851934e-16 + 3.673940397442060881462e-16im
      -8.869676734855019004796e-16 + 3.673940397442062790345e-16im
 7.999999999999999999999999999999 - 2.5717582782094416721747e-15im
       8.869676734855019400138e-16 + 3.673940397442056273012e-16im
       3.673940397442060544016e-16 + 3.673940397442058181894e-16im
 \      1.521795939970901291223e-16 + 3.673940397442058972579e-16im
                      1.181061e-31 + 3.673940397442059531678e-16im

tells us that the sample was generated by $e^{3*2\pi i x}$.

The output of the DFT is a vector whose components contain, by index:
index 1: the zero frequency, i.e., the sum of the samples
index 2 to N/2: the positive frequencies
index N/2+1: the Nyquist frequency
index N/2+2 to N: the negative frequencies in decreasing order with respect to the modulus

Remember that the result is multiplied by a factor N.

julia>dft(exp.(-1*v))
8-element Array{ArbComplex{128},1}:
   -2.956558911618339712192e-16 - 1.2246467991473535393e-16im
 -1.224646799147353064744e-16 - 1.224646799147353327202e-16im
  -5.07265313323633251691e-17 - 1.224646799147353239348e-16im
                  1.31229e-32 - 1.224646799147353177226e-16im
   5.07265313323633507716e-17 - 1.224646799147353115104e-16im
   1.224646799147353289708e-16 - 1.22464679914735302725e-16im
  2.956558911618339756119e-16 - 1.224646799147352815152e-16im
                          8.0 + 8.572527594031472240582e-16im
julia> dft(exp.(-3*v))
8-element Array{ArbComplex{128},1}:
      -1.521795939970898986998e-16 - 3.673940397442060090777e-16im
                      1.181061e-31 - 3.673940397442059531678e-16im
       1.521795939970901291223e-16 - 3.673940397442058972579e-16im
       3.673940397442060544016e-16 - 3.673940397442058181894e-16im
       8.869676734855019400138e-16 - 3.673940397442056273012e-16im
 7.999999999999999999999999999999 + 2.5717582782094416721747e-15im
      -8.869676734855019004796e-16 - 3.673940397442062790345e-16im
       -3.67394039744205851934e-16 - 3.673940397442060881462e-16im

The DFT is a linear operator

julia> dft(exp.(-3*v)+5*exp.(2*v))
8-element Array{ArbComplex{128},1}:
         -1.376826393144442850962e-15 + 8.572527594031474681e-16im
       -2.956558911618339572123e-15 + 8.57252759403147948206e-16im
 40.00000000000000015217959399709 - 8.9399216337756781378404e-15im
        3.323952951362545832484e-15 + 8.57252759403146634889e-16im
        2.111614472632855342204e-15 + 8.57252759403147249973e-16im
  8.000000000000000507265313323633 + 3.796405077356794725157e-15im
       -8.86967673485501638022e-16 + 8.572527594031468981916e-16im
        -8.74659353067838975613e-16 + 8.57252759403147213324e-16im

Resuming, the idea is that the DFT permits us to see the frequency composition of a signal

julia> n = 4;
julia> v = ArbComplex.(rand(Float64, 4)+im*rand(Float64,4))
4-element Array{ArbComplex{128},1}:
  0.9469540787992749564239147730405 + 0.9011349970035253953426490625134im
 0.9407191613523724349477106443373 + 0.09744835485762926019503993302351im
  0.6974901676912010906050909397891 + 0.6783474788880952122127609982272im
  0.370033099250 dft(v)/4
julia>dft(v)/4
4-element Array{ArbComplex{128},1}:
    0.7387991267733211020463102158828 + 0.5290776432152716757251198487211im
 -0.02311686903653342728404140871135 - 0.08697463599662658140232451842166im
   0.08342299647191692146819264053192 + 0.2606635947305386280525851816492im
     0.147848824590570360193453325337 + 0.1983683950543416729672685505648im

means that the original signal is approximately $(0.74+0.53im)+(-0.023-0.087im)e^{i 2\pi x}+(0.15+0.2im)e^{-i 2\pi x}+(0.083+0.27im)e^{2 i 2\pi x}$.

@JeffreySarnoff
Copy link
Owner

When I test ArbNumerics#master on Windows it crashes because the dot product tested in test/ArblibVector.jl line 27 when run causes Julia v1.3 to crash (with T=ArbComplex{128} and with T=ArbReal{128}). The dot product tested in line 56 works.

julia> using ArbNumerics, Test
julia> import ArbNumerics: ArblibVector, take!, free!
julia> using LinearAlgebra

julia> T=ArbReal{128}
ArbReal{128}

julia>     o = ones(T, 5)
5-element Array{ArbReal{128},1}:
 1.0
 1.0
 1.0
 1.0
 1.0

julia>        d=dot(o,o)      # repl crashes, window vanishes

@JeffreySarnoff
Copy link
Owner

@orkolorko @jagot I am working to integrate the examples into the docs. Help me with the dot problem above.

@jagot
Copy link
Contributor Author

jagot commented Dec 11, 2019

@JeffreySarnoff I'm afraid I can't reproduce; I am normally only using Mac & Linux, but I tried the following in a VirtualBox instance of Win10:

               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.3.0 (2019-11-26)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> cd(mkpath("test"))

(v1.3) pkg> activate .
Activating new environment at `C:\Users\Jagot\AppData\Local\Julia-1.3.0\test\Project.toml`

(test) pkg> add LinearAlgebra, ArbNumerics#master
  Updating registry at `C:\Users\Jagot\.julia\registries\General`
  Updating git-repo `https://github.com/JuliaRegistries/General.git`
  Updating git-repo `https://github.com/JeffreySarnoff/ArbNumerics.jl.git`
 Resolving package versions...
  Updating `C:\Users\Jagot\AppData\Local\Julia-1.3.0\test\Project.toml`
  [7e558dbc] + ArbNumerics v1.0.1 #master (https://github.com/JeffreySarnoff/ArbNumerics.jl.git)
  [37e2e46d] + LinearAlgebra
  Updating `C:\Users\Jagot\AppData\Local\Julia-1.3.0\test\Manifest.toml`
  [7e558dbc] + ArbNumerics v1.0.1 #master (https://github.com/JeffreySarnoff/ArbNumerics.jl.git)
  [b99e7846] + BinaryProvider v0.5.8
  [01680d73] + GenericSVD v0.2.2
  [c145ed77] + GenericSchur v0.3.0
  [efe28fd5] + OpenSpecFun_jll v0.5.3+1
  [0d4725de] + Readables v0.3.3
  [276daf66] + SpecialFunctions v0.9.0
  [2a0f44e3] + Base64
  [ade2ca70] + Dates
  [8ba89e20] + Distributed
  [b77e0a4c] + InteractiveUtils
  [76f85450] + LibGit2
  [8f399da3] + Libdl
  [37e2e46d] + LinearAlgebra
  [56ddb016] + Logging
  [d6f4376e] + Markdown
  [44cfe95a] + Pkg
  [de0858da] + Printf
  [3fa0cd96] + REPL
  [9a3f8284] + Random
  [ea8e919c] + SHA
  [9e88b42a] + Serialization
  [6462fe0b] + Sockets
  [8dfed614] + Test
  [cf7118a7] + UUIDs
  [4ec0a83e] + Unicode

julia> using ArbNumerics

julia> using LinearAlgebra

julia>

julia> T=ArbComplex{128}
ArbComplex{128}

julia> o = ones(T, 5)
5-element Array{ArbComplex{128},1}:
 1.0 + 0im
 1.0 + 0im
 1.0 + 0im
 1.0 + 0im
 1.0 + 0im

julia> d=dot(o,o)
5.0 + 0im

julia> T=ArbReal{128}
ArbReal{128}

julia> o = ones(T, 5)
5-element Array{ArbReal{128},1}:
 1.0
 1.0
 1.0
 1.0
 1.0

julia> d = dot(o,o)
5.0

julia> versioninfo()
Julia Version 1.3.0
Commit 46ce4d7933 (2019-11-26 06:09 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-7920HQ CPU @ 3.10GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)

julia>

@JeffreySarnoff
Copy link
Owner

That was very helpful. The problem only exists in 1.3.0-rc not in 1.3.0 -- thanks, I will clean up doc and rerelease

@JeffreySarnoff
Copy link
Owner

it has been registered -- I do not know how to explain [let alone simply] that inverse_dft(dft(x)) may not equal x and how to apply the length of x as a scale factor to have that become equal to x.

@JeffreySarnoff JeffreySarnoff mentioned this pull request Dec 28, 2019
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.

3 participants