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

gangchern · March 6, 2022, 10:18am

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.

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.

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

ChrisRackauckas · March 6, 2022, 12:42pm

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.

gangchern · March 9, 2022, 5:29am

great news! Looking forward to that.

gangchern · June 5, 2022, 5:41am

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

Raymond · June 5, 2022, 7:21pm

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.

ChrisRackauckas · June 5, 2022, 10:56pm

Is that using the unityper form on master?

gangchern · June 6, 2022, 3:49am

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.

Raymond · June 8, 2022, 8:31pm

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.

ChrisRackauckas · June 8, 2022, 8:36pm

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

Raymond · June 8, 2022, 8:38pm

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

Raymond · June 8, 2022, 8:41pm

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

julia> @time f = (x + y)^50000 |>expand;

Raymond · June 8, 2022, 8:51pm

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)

ChrisRackauckas · June 8, 2022, 8:52pm

You may have hit a compilation phase?

Raymond · June 8, 2022, 8:54pm

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.

Raymond · June 8, 2022, 8:57pm

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)

Raymond · June 8, 2022, 9:05pm

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)

contradict · June 8, 2022, 11:26pm

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

Raymond · June 9, 2022, 12:13am

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.

ChrisRackauckas · June 9, 2022, 2:38am

Yeah it’s overflow. Change the integer type