Test of expand and substitute in Symbolics.jl, SymEngine.jl and Mathematica

General Usage Performance
1

My computer is Mac Pro M1, 2020.
General Comments: for expand, SymEngine is nearly as good as Mathematica, Symbolics‘s “expand” function still needs to be greatly improved. substitute are similar for the three CAS in performance.

  1. expand,
    Symbolics (code+time)
 f=(x+big(2)y)^6000|>expand;  # 64.2s 
 f=(x+big(2)y)^100000|>expand;  # looks infinity

SymEngine

 f=(x+2y)^6000|>expand; #0.5s
 f=(x+2y)^100000|>expand;  # 10.6s

Mathematica takes 0.01s and 3.1s respectively.

  1. substitute
    Symbolics takes 0.9s
 substitute(f,Dict(x=>big(1),y=>big(2)));

SymEngine takes 0.4s

subs(f,Dict(x=>1,y=>2));

Mathematica takes 0.06s

1 Like
2

Yes, it does need to be improved. https://github.com/JuliaSymbolics/SymbolicUtils.jl/pull/432 should help a lot but is not released right now to Symbolics.

3

great news! Looking forward to that.

4

Now I tested again. Symbolics is much faster than SymEngine now

5 
julia> using Symbolics

julia> @variables x, y
julia> @time f = (x + y)^6000 |>expand;
 14.399519 seconds (30.27 M allocations: 3.254 GiB, 3.97% gc time, 25.09% compilation time)

julia> using SymEngine
julia> @vars x, y
julia> f = (x + y)^6000 |> expand;
julia> @time f = (x + y)^6000 |>expand;
  0.012988 seconds (36.10 k allocations: 12.602 MiB)

In my Mac M1, SymEngine is faster than Symbolics.

6

Is that using the unityper form on master?

7 
using Symbolics
@variables x, y
@time f = (x + 2y)^50000 |>expand;
0.026855 seconds (70.92 k allocations: 6.684 MiB, 30.24% gc time, 9.12% compilation time)

Strange! looks too good to be true! Maybe I did a wrong test.

8 
julia> @time f = (x + 2y)^50000 |>expand;
  0.022756 seconds (51.71 k allocations: 4.784 MiB)

julia> @time f = (x + y)^6000 |>expand;
 10.740392 seconds (15.51 M allocations: 2.443 GiB, 4.71% gc time)

It’s so strange.

9

That’s because of the type-instability that unityper fixes.

10

why is (x + 2y)^50000 |>expand is faster than f = (x + y)^50000?

11 
julia> @time f = (x + 2y)^50000 |>expand;
  0.023723 seconds (51.76 k allocations: 4.787 MiB)

julia> @time f = (x + y)^50000 |>expand;
12 
julia> @time f = (x + 4y)^50000 |>expand;
  0.005870 seconds (16.61 k allocations: 1.465 MiB)

julia> @time f = (x + 6y)^50000 |>expand;
  0.044559 seconds (109.41 k allocations: 8.064 MiB, 56.32% compilation time)

julia> @time f = (x + 8y)^50000 |>expand;
  0.051483 seconds (40.72 k allocations: 2.619 MiB, 92.70% compilation time)

julia> @time f = (x + 10y)^50000 |>expand;
  0.028289 seconds (51.71 k allocations: 4.784 MiB, 31.20% gc time)
13

You may have hit a compilation phase?

14

It seems that when the coefficient of y is even, it will be very fast.
But when the coefficient is odd, it will run for a long long long time.

15 
julia> for i in 1:10
           @time f = (x + i * y)^5000 |>expand
       end
  5.883341 seconds (10.36 M allocations: 1.695 GiB, 6.25% gc time, 0.07% compilation time)
  0.008549 seconds (39.72 k allocations: 3.583 MiB)
  5.866237 seconds (10.36 M allocations: 1.694 GiB, 5.82% gc time)
  0.002496 seconds (13.64 k allocations: 1.154 MiB)
  6.197537 seconds (10.36 M allocations: 1.694 GiB, 9.01% gc time)
  0.008952 seconds (39.72 k allocations: 3.582 MiB)
  5.930975 seconds (10.36 M allocations: 1.694 GiB, 5.82% gc time)
  0.001464 seconds (7.69 k allocations: 619.297 KiB)
  6.237053 seconds (10.36 M allocations: 1.694 GiB, 10.19% gc time)
  0.011223 seconds (39.72 k allocations: 3.582 MiB)
16 
julia> for i in 1:10
           @time f = (x + i * y)^10000 |>expand
       end
 38.717680 seconds (39.79 M allocations: 5.797 GiB, 26.24% gc time, 0.01% compilation time)
  0.009729 seconds (42.83 k allocations: 3.898 MiB)
 31.223451 seconds (39.79 M allocations: 5.797 GiB, 11.21% gc time)
  0.002769 seconds (14.51 k allocations: 1.247 MiB)
 32.939406 seconds (39.79 M allocations: 5.797 GiB, 13.77% gc time)
  0.009396 seconds (42.83 k allocations: 3.898 MiB)
 31.781937 seconds (39.79 M allocations: 5.797 GiB, 10.80% gc time)
  0.043063 seconds (37.99 k allocations: 2.418 MiB, 95.46% compilation time)
 32.121851 seconds (39.79 M allocations: 5.797 GiB, 10.63% gc time)
  0.009469 seconds (42.83 k allocations: 3.898 MiB)
17

The pattern goes away if you use floating point numbers. This looks like integer overflow to me.

julia> for i in 1.0:10.0
       @time f=(x+i*y)^5000 |> expand
       end
  8.599273 seconds (10.36 M allocations: 1.696 GiB, 8.55% gc time)
  8.491495 seconds (10.36 M allocations: 1.696 GiB, 9.21% gc time)
  8.400282 seconds (10.36 M allocations: 1.696 GiB, 8.55% gc time)
  8.377576 seconds (10.36 M allocations: 1.696 GiB, 8.16% gc time)
  8.375058 seconds (10.37 M allocations: 1.696 GiB, 8.39% gc time, 0.09% compilation time)
  8.483856 seconds (10.36 M allocations: 1.696 GiB, 8.57% gc time)
  8.438680 seconds (10.36 M allocations: 1.696 GiB, 8.32% gc time)
  8.378851 seconds (10.36 M allocations: 1.696 GiB, 8.54% gc time)
  8.499839 seconds (10.36 M allocations: 1.696 GiB, 9.27% gc time)
  8.608331 seconds (10.37 M allocations: 1.696 GiB, 8.55% gc time)

Interestingly,

julia> (2*y)^5000
ERROR: OverflowError: 4294967296 * 4294967296 overflowed for type Int64

but

julia> (x+2*y)^5000
(x + 2y)^5000
2 Likes
18 
julia> (x + y)^66 |> expand
x^66 + y^66 + 66y*(x^65) + 7219428434016265740(x^33)*(y^33) + 4472995859186094240(x^37)*(y^29) + 6406484391866534976(x^31)*(y^35) + 3413602103063071920(x^28)*(y^38) + 4922879481520(x^12)*(y^54) + 2515943049305400(x^49)*(y^17) + 77413632286320(x^52)*(y^14) + 17302625882942400(x^19)*(y^47) + 3413602103063071920(x^38)*(y^28) + 17302625882942400(x^47)*(y^19) + 37014131440(x^9)*(y^57) + 6406484391866534976(x^35)*(y^31) + 4922879481520(x^54)*(y^12) + 2450791253481179840(x^39)*(y^27) + 348524321356411200(x^23)*(y^43) + 77413632286320(x^14)*(y^52) + 40661170824914640(x^20)*(y^46) + 268367258592576(x^15)*(y^51) + 5743572120(x^8)*(y^58) + 20448884000160(x^53)*(y^13) + 45760(x^63)*(y^3) + 855420636763836(x^16)*(y^50) + 855420636763836(x^50)*(y^16) + 2145(x^64)*(y^2) + 210980549208(x^56)*(y^10) + 1049058207282797712(x^41)*(y^25) + 2145(x^2)*(y^64) + 1074082795968(x^55)*(y^11) + 5516694892996182896(x^30)*(y^36) + 90858768(x^6)*(y^60) + 1074082795968(x^11)*(y^55) + 210980549208(x^10)*(y^56) + 90858768(x^60)*(y^6) + 1049058207282797712(x^25)*(y^41) + 7007092303604022630(x^34)*(y^32) + 720720(x^62)*(y^4) + 6848956078664700(x^18)*(y^48) + 268367258592576(x^51)*(y^15) + 348524321356411200(x^43)*(y^23) + 778789440(x^59)*(y^7) + 7007092303604022630(x^32)*(y^34) + 4472995859186094240(x^29)*(y^37) + 37014131440(x^57)*(y^9) + 2515943049305400(x^17)*(y^49) + 45760(x^3)*(y^63) + 1654284096099796392(x^26)*(y^40) + 89067326568860640(x^21)*(y^45) + 182183167981760400(x^22)*(y^44) + 5743572120(x^58)*(y^8) + 89067326568860640(x^45)*(y^21) + 778789440(x^7)*(y^59) + 5516694892996182896(x^36)*(y^30) + 8936928(x^5)*(y^61) + 20448884000160(x^13)*(y^53) + 6848956078664700(x^48)*(y^18) + 720720(x^4)*(y^62) + 182183167981760400(x^44)*(y^22) + 2450791253481179840(x^27)*(y^39) + 66x*(y^65) + 40661170824914640(x^46)*(y^20) + 624439409096903400(x^24)*(y^42) + 624439409096903400(x^42)*(y^24) + 8936928(x^61)*(y^5) + 1654284096099796392(x^40)*(y^26)

julia> (x + y)^67 |> expand
x^67 + y^67 + 67y*(x^66) + 766480(x^4)*(y^63) + 766480(x^63)*(y^4) + 7886597962249166160(x^29)*(y^38) + 972963730453314600(x^43)*(y^24) + 3371363686069236(x^50)*(y^17) + 3371363686069236(x^17)*(y^50) + 129728497393775280(x^46)*(y^21) + 99795696(x^6)*(y^61) + 4105075349580976232(x^40)*(y^27) + 972963730453314600(x^24)*(y^43) + 67x*(y^66) + 2703342303382594104(x^41)*(y^26) + 1673497616379701112(x^42)*(y^25) + 5996962277488(x^12)*(y^55) + 99795696(x^61)*(y^6) + 2703342303382594104(x^26)*(y^41) + 1673497616379701112(x^25)*(y^42) + 869648208(x^7)*(y^60) + 2211(x^65)*(y^2) + 9364899127970100(x^49)*(y^18) + 345780890878896(x^52)*(y^15) + 247994680648(x^57)*(y^10) + 47905(x^3)*(y^64) + 9657648(x^62)*(y^5) + 271250494550621040(x^22)*(y^45) + 42757703560(x^9)*(y^58) + 2211(x^2)*(y^65) + 247994680648(x^10)*(y^57) + 42757703560(x^58)*(y^9) + 4105075349580976232(x^27)*(y^40) + 1285063345176(x^56)*(y^11) + 47905(x^64)*(y^3) + 97862516286480(x^14)*(y^53) + 57963796707857040(x^47)*(y^20) + 97862516286480(x^53)*(y^14) + 1123787895356412(x^16)*(y^51) + 24151581961607100(x^48)*(y^19) + 57963796707857040(x^20)*(y^47) + 1285063345176(x^11)*(y^56) + 5996962277488(x^55)*(y^12) + 530707489338171600(x^23)*(y^44) + 5864393356544251760(x^39)*(y^28) + 271250494550621040(x^45)*(y^22) + 869648208(x^60)*(y^7) + 6522361560(x^8)*(y^59) + 5864393356544251760(x^28)*(y^39) + 9364899127970100(x^18)*(y^49) + 25371763481680(x^54)*(y^13) + 25371763481680(x^13)*(y^54) + 129728497393775280(x^21)*(y^46) + 9657648(x^5)*(y^62) + 24151581961607100(x^19)*(y^48) + 345780890878896(x^15)*(y^52) + 530707489338171600(x^44)*(y^23) + 1123787895356412(x^51)*(y^16) + 7886597962249166160(x^38)*(y^29) + 6522361560(x^59)*(y^8) - 4220223336089263246(x^33)*(y^34) - 5033167378238994010(x^32)*(y^35) - 6523564788846833744(x^31)*(y^36) - 6523564788846833744(x^36)*(y^31) - 8457053321527274480(x^37)*(y^30) - 4220223336089263246(x^34)*(y^33) - 8457053321527274480(x^30)*(y^37) - 5033167378238994010(x^35)*(y^32)

it seems that when the exponent is larger than 66, the result is wrong.

19

Yeah it’s overflow. Change the integer type