CAS Best Practices

Specific Domains Machine Learning
47

Hi Petr @PetrKryslUCSD

Ugh :cold_face:, I really sympathize w you about " … by the morning the next day it choked :upside_down_face: and stopped with an :nauseated_face: error."

And you have always been uber helpful :innocent:, so I wanted to throw out there for your consideration “… arbitrary-precision interval arithmetic using the midpoint-radius representation, also known as ball arithmetic.” using Nemo and/or IntervalArithmetic.jl described better by @dpsanders here @@ Calculation Pi^Pi^Pi in hexadecimal in High precision - #60 by dpsanders

As for my take on arbitrary-precision interval arithmetic (so far)
is that it is deterministic (no “next day bombs” - UGH!),
AND rigorous (?almost?) to the extent of a mathematical proof of correctness. :innocent:

So I’ve been using Nemo instead of SpecialFunctions for a while when using ApproxFun,
but my particular use of Nemo code might be only obliquely related to your needs,
so I whipped up some quick code using Nemo and using Hecke for you below
that just might be better suited to your purposes (of iteration), memory-conservative,
and fast symbolic code ?

NOTE WARNING and TIP >> To perform the SAME task using floating point arithmetic would require an a priori error analysis. … In contrast Arb’s ball arithmetic allows the error analysis to be carried along with the computation. … << NOTE WARNING and TIP

NOTE TIP >> using Nemo allows for fast development, while at the same time maintaining Guaranteed results. << NOTE TIP

####
####

# Example Code using Nemo using Hecke Computer Algebra / CAS and Number Theory Packages for Julia - 1705.06134.pdf @@ https://arxiv.org/pdf/1705.06134.pdf
#
# NOTE Don't whack yourself lol, Buy some insurance by First Creating a New Experimental TEST Environment
# (v1.4-EXP3) pkg> activate v1.4-EXP3
# add Nemo     # Nemo  v0.18.4
# add Hecke    # Hecke v0.8.6
using Nemo
using Hecke

## Marc adding ...
# Nemo/Hecke: >> Computer Algebra (CAS) and Number Theory Packages for the Julia Programming Language <<
# Source: https://arxiv.org/pdf/1705.06134.pdf
#
# 5.2  Arb: arbitrary precision ball arithmetic Computing over R and C requires using some approximate model for these fields. The most common model is floating-point arithmetic.
# NOTE >> However, for many computer algebra / CAS algorithms,
# the error analysis necessary to guarantee correct results with floating-point arithmetic becomes impractical. Interval arithmetic solves this problem by effectively making error analysis automatic << NOTE
#
# Many problems can be solved using lazy evaluation: the user can try a computation with some tentative precision P and restart with precision 2*P if that fails. The precision can be set optimally when a good estimate for the minimal required P available; that is, the intervals can be used as if they were plain floating-point numbers,and NOTE >> the Automatic error bounds simply provide a Certificate. << NOTE
#
# 6.2  The use of interval arithmetic @@ https://arxiv.org/pdf/1705.06134.pdf
# Finally we check whether diam(B(i)α) ≤2−p.
# NOTE If not, we Increase the Working Precision P and repeat.
# NOTE >> When finished, the midpoints of the ballsB(i)α are approximationsˆσi(α)
# with the desired property. << NOTE
# The following Hecke code illustrates this.
#
QQx, x = PolynomialRing(QQ, "x") # NIX FlintQQ
# function NumberField(f::fmpq_poly, s::AbstractString; cached::Bool = true, check::Bool = true)
flintfun_f = x^3 + 3*x + 1
K, gen_po_a = Nemo.NumberField(flintfun_f, "pgen_a",cached=true,check=true)
alpha = gen_po_a^2 + gen_po_a + 1
# WARNING Original paper seems to have a Code BUG ish ? IOW appeared initially IMPOSSIBLE w global int_target_precision_p = 128
#  per julia> log(1/typemax(Int64))/log(2) == -63.0
#      julia> 2^64 == 0
# BUG w Original -int_target_precision_p = 128
global int_target_precision_p = 16
global int_initial_precision_p_prime = 4
println("")
println("Increase the Working Precision P and repeat UNTIL all([ Hecke.radiuslttwopower(u, -int_target_precision_p) for u in z])")
while true
 CCy, y = PolynomialRing(AcbField(int_initial_precision_p_prime), "y")
 g = CCy(x^3 + 3*x + 1)
 rts = roots(g) # Attempts to isolate the complex roots of the complex polynomial xx by iteratively refining balls in which they lie. .. Q: ? Related to https://en.wikipedia.org/wiki/Fundamental_theorem_of_algebra , https://en.wikipedia.org/wiki/Weierstrass_factorization_theorem ?
 z = map(x^2 + x + 1, rts)
 # 3-element Array{acb,1}: # .. # [-1.8908090 +/- 3.32e-8] + i*[-2.3196168 +/- 2.12e-8]
 #
 # https://fredrikj.net/arb/mag.html#c.mag_cmp_2exp_si
   #  int mag_cmp_2exp_si(const mag_t x, slong y)
   #   Returns negative, zero, or positive, depending on whether x is smaller, equal, or larger than 2y
 if all( [ Hecke.radiuslttwopower(u, -1*int_target_precision_p) for u in z ] )
   # Possible BUG w Original code where -int_target_precision_p = 128 likely due to 2^128 Julia overflow/underflow from 64Bit numbers per julia> 2^128 == 0
   # TODO Try again with ALL numbers as arbitrary high precision Arbs ?
   # TODO IOW ALL Arbs, and NO 32Bit or 64Bit numbers
   break
 else
  global int_initial_precision_p_prime = 2*int_initial_precision_p_prime
  println("Increase/DOUBLE the Working Precision P int_initial_precision_p_prime:=",int_initial_precision_p_prime," and repeat UNTIL all([ Hecke.radiuslttwopower(u, -int_target_precision_p) for u in z])")
 end
end

# NOTE WARNING and TIP >> To perform the SAME task using floating point arithmetic would require an a priori error analysis. .. In contrast Arb’s ball arithmetic allows the error analysis to be carried along with the computation. .. << NOTE WARNING and TIP
# NOTE TIP >> using Nemo allows for fast development, while at the same time maintaining Guaranteed results. << NOTE TIP 

println("##125 !! DONE EOJ Increased the Working Precision P_prime and repeated UNTIL all([ Hecke.radiuslttwopower(u, -int_target_precision_p) for u in z]) is >=", 2.0^(-1*int_target_precision_p))

# === Original source of my motivation / inspiraton
# arbitrary-precision interval arithmetic using the midpoint-radius representation, also known as ball arithmetic.
# https://discourse.julialang.org/t/calculation-pi-pi-pi-in-hexadecimal-in-high-precision/51785/60?u=marc.cox
# dpsanders ** # Committer
# @Marc.Cox I believe it’s using Fredrik Johansson’s arb software. He has a nice paper describing how that works.

## File: CITATION.bib
## Citing Hecke
# If your research depends on computations done with Hecke, please consider giving us a formal citation:
# - Claus Fieker, William Hart, Tommy Hofmann and Fredrik Johansson, [Nemo/Hecke: Computer Algebra and Number Theory Packages
#   for the Julia Programming Language](http://doi.acm.org/10.1145/3087604.3087611). In: Proceedings of ISSAC '17, pages 157–164, New York, NY, USA, 2017. ACM.
#
## NOTE Documentation Yay !
# Documentation for Nemo.jl can be found here: https://nemocas.github.io/Nemo.jl/latest/index.html
# Documentation for Hecke.jl can be found here: https://www.thofma.com/Hecke.jl/latest/
1 Like