chualynn / theta.jl Goto Github PK

Hello,

I was very excited to find this package since I need to compute Riemann Theta functions in Julia. However, the theta function exported by this package seems to disagree with the SiegelTheta function in Mathematica on a simple test case. The Mathematica code
SiegelTheta[{{I/2, -1/2}, {-1/2, I/2}}, {0.5, 0.5}] gives 2.21457*10^-11 - 9.61931*10^-24 I. This makes sense, as one can prove this value is exactly 0. However, I find that theta([0.5, 0.5], RiemannMatrix([[1im/2, -1/2] [-1/2, 1im/2]])) evaluates to 0.9925441784910575 - 5.823949418449035e-35im in Julia 1.7. Am I misunderstanding something? Thanks very much for your help!

I get (version julia-1.9.0-win64) :

julia> theta([0], RiemannMatrix([im]))
ERROR: MethodError: no method matching +(::Float64, ::Vector{Float64})
For element-wise addition, use broadcasting with dot syntax: scalar .+ array

Closest candidates are:
+(::Any, ::Any, ::Any, ::Any...)
@ Base operators.jl:578
+(::T, ::T) where T<:Union{Float16, Float32, Float64}
@ Base float.jl:408
+(::Union{Float16, Float32, Float64}, ::BigFloat)
@ Base mpfr.jl:423
...

Stacktrace:
[1] (::Theta.var"#9#12"{RiemannMatrix, Vector{Float64}, Vector{Any}, Matrix{Float64}, Matrix{Float64}, Int64})(p::Vector{Float64})
@ Theta .\none:0
[2] iterate
@ .\generator.jl:47 [inlined]
[3] collect(itr::Base.Generator{Vector{Vector{Float64}}, Theta.var"#9#12"{RiemannMatrix, Vector{Float64}, Vector{Any}, Matrix{Float64}, Matrix{Float64}, Int64}})
@ Base .\array.jl:782
[4] oscillatory_part(R::RiemannMatrix, x::Vector{Float64}, y0::Matrix{Float64}, shift_n::Vector{Float64}, derivs::Vector{Any})
@ Theta xxx.julia\packages\Theta\2Upgk\src\eval_theta.jl:61
[5] theta(z::Vector{Int64}, R::RiemannMatrix; char::Vector{Any}, derivs::Vector{Any}, derivs_t::Vector{Any})
@ Theta xxx.julia\packages\Theta\2Upgk\src\eval_theta.jl:46
[6] theta(z::Vector{Int64}, R::RiemannMatrix)
@ Theta xxx.julia\packages\Theta\2Upgk\src\eval_theta.jl:20
[7] top-level scope
@ REPL[2]:1

Wrong :

julia> theta([0,0], RiemannMatrix([im/5 0;0 im/5]))
1.000000602807001 + 0.0im

Right :

julia> theta([0,0], RiemannMatrix([5im 0;0 5im]))
1.000000602807001 + 0.0im

My julia Version : julia-1.9.0-win64

Let me say first that your package saved my life 😉 I might have a suggestion for an enhancement though.

I needed the Jacobi theta function, of which Riemann theta is a generalization. While wrapping variables into arrays works perfectly fine, it would be nice to have a dedicated method for scalar inputs in order to avoid allocating temporary 1-element arrays.

I can imagine it could be too difficult due to how your implementation relies on inputs being arrays or it could bring negligible benefits in comparison with evaluation of the function itself. In such a case, a simple wrapper would be still very convenient.

EDIT

I was a bit fast with this one... I don't know what exactly happened, but a combination of kernel restarts and updates broke something and now I can't wrap scalars any more. Boiling your example down to single numbers results in:

M = hcat(0.794612+1.9986im)
R = RiemannMatrix(M)
z = [0.30657351+0.34017115im]
theta(z, RiemannMatrix(M))

MethodError: no method matching +(::Float64, ::Array{Float64,1})
For element-wise addition, use broadcasting with dot syntax: scalar .+ array
Closest candidates are:
  +(::Any, ::Any, ::Any, ::Any...) at operators.jl:538
  +(::Float64, ::Float64) at float.jl:401
  +(::AbstractFloat, ::Bool) at bool.jl:106
  ...

Here's my environment:

Pkg.status()

Status `~/.julia/environments/v1.5/Project.toml`
  [7073ff75] IJulia v1.23.2
  [76087f3c] NLopt v0.5.1
  [91a5bcdd] Plots v1.10.6
  [d330b81b] PyPlot v2.9.0
  [42bdb5c4] Theta v0.1.2

Any help appreciated!