Examples

Four-square identity

The first example shows that the four-square identity holds:

which was originally proved by Euler. The code can also be found in this test of the package.

First, we reset the maximum degree of the polynomial to 4, since the RHS of the equation has a priori terms of fourth order, and define the 8 independent variables.

julia> using TaylorSeries

julia> # Define the variables α₁, ..., α₄, β₁, ..., β₄
       make_variable(name, index::Int) = string(name, TaylorSeries.subscriptify(index))
make_variable (generic function with 1 method)

julia> variable_names = [make_variable("α", i) for i in 1:4]
4-element Array{String,1}:
 "α₁"
 "α₂"
 "α₃"
 "α₄"

julia> append!(variable_names, [make_variable("β", i) for i in 1:4])
8-element Array{String,1}:
 "α₁"
 "α₂"
 "α₃"
 "α₄"
 "β₁"
 "β₂"
 "β₃"
 "β₄"

julia> # Create the TaylorN variables (order=4, numvars=8)
       a1, a2, a3, a4, b1, b2, b3, b4 = set_variables(variable_names, order=4)
8-element Array{TaylorSeries.TaylorN{Float64},1}:
  1.0 α₁ + 𝒪(‖x‖⁵)
  1.0 α₂ + 𝒪(‖x‖⁵)
  1.0 α₃ + 𝒪(‖x‖⁵)
  1.0 α₄ + 𝒪(‖x‖⁵)
  1.0 β₁ + 𝒪(‖x‖⁵)
  1.0 β₂ + 𝒪(‖x‖⁵)
  1.0 β₃ + 𝒪(‖x‖⁵)
  1.0 β₄ + 𝒪(‖x‖⁵)

julia> a1 # variable a1
 1.0 α₁ + 𝒪(‖x‖⁵)

Now we compute each term appearing in Eq. (\ref{eq:Euler})

julia> # left-hand side
       lhs1 = a1^2 + a2^2 + a3^2 + a4^2 ;

julia> lhs2 = b1^2 + b2^2 + b3^2 + b4^2 ;

julia> lhs = lhs1 * lhs2;

julia> # right-hand side
       rhs1 = (a1*b1 - a2*b2 - a3*b3 - a4*b4)^2 ;

julia> rhs2 = (a1*b2 + a2*b1 + a3*b4 - a4*b3)^2 ;

julia> rhs3 = (a1*b3 - a2*b4 + a3*b1 + a4*b2)^2 ;

julia> rhs4 = (a1*b4 + a2*b3 - a3*b2 + a4*b1)^2 ;

julia> rhs = rhs1 + rhs2 + rhs3 + rhs4;

We now compare the two sides of the identity,

julia> lhs == rhs
true

The identity is satisfied. $\square$

Fateman test

Richard J. Fateman, from Berkley, proposed as a stringent test of polynomial multiplication the evaluation of $s*(s+1)$, where $s = (1+x+y+z+w)^{20}$. This is implemented in the function fateman1 below. We shall also consider the form $s^2+s$ in fateman2, which involves fewer operations (and makes a fairer comparison to what Mathematica does).

julia> using TaylorSeries

julia> const order = 20
20

julia> const x, y, z, w = set_variables(Int128, "x", numvars=4, order=2order)
4-element Array{TaylorSeries.TaylorN{Int128},1}:
  1 x₁ + 𝒪(‖x‖⁴¹)
  1 x₂ + 𝒪(‖x‖⁴¹)
  1 x₃ + 𝒪(‖x‖⁴¹)
  1 x₄ + 𝒪(‖x‖⁴¹)

julia> function fateman1(degree::Int)
           T = Int128
           s = one(T) + x + y + z + w
           s = s^degree
           s * ( s + one(T) )
       end
fateman1 (generic function with 1 method)

(In the following lines, which are run when the documentation is built, by some reason the timing appears before the command executed.)

julia> @time fateman1(0);
  0.829796 seconds (158.29 k allocations: 47.186 MiB, 1.30% gc time)

julia> @time f1 = fateman1(20);
  6.129314 seconds (3.79 k allocations: 60.200 MiB, 0.19% gc time)

Another implementation of the same, but exploiting optimizations related to ^2 yields:

julia> function fateman2(degree::Int)
           T = Int128
           s = one(T) + x + y + z + w
           s = s^degree
           s^2 + s
       end
fateman2 (generic function with 1 method)

julia> fateman2(0);

julia> @time f2 = fateman2(20); # the timing appears above
  3.112194 seconds (4.08 k allocations: 64.357 MiB, 0.18% gc time)

We note that the above functions use expansions in Int128. This is actually required, since some coefficients are larger than typemax(Int):

julia> getcoeff(f2, [1,6,7,20]) # coefficient of x y^6 z^7 w^{20}
128358585324486316800

julia> ans > typemax(Int)
true

julia> length(f2)
41

julia> sum(TaylorSeries.size_table)
135751

These examples show that fateman2 is nearly twice as fast as fateman1, and that the series has 135751 monomials in 4 variables.

Bechmarks

The functions described above have been compared against Mathematica v11.1. The relevant files used for benchmarking can be found here. Running on a MacPro with Intel-Xeon processors 2.7GHz, we obtain that Mathematica requires on average (5 runs) 3.075957 seconds for the computation, while for fateman1 and fateman2 above we obtain 2.811391 and 1.490256, respectively.

Then, with the current version of TaylorSeries.jl, our implementation of fateman1 is about 10% faster, and fateman2 is about a factor 2 faster. (The original test by Fateman corresponds to fateman1 above, which avoids some optimizations related to squaring.)