Skip to content

Changed the function called by sqrt from mpmath to math#25111

Merged
oscarbenjamin merged 4 commits intosympy:masterfrom
haru-44:gmpy_isqrt
May 11, 2023
Merged

Changed the function called by sqrt from mpmath to math#25111
oscarbenjamin merged 4 commits intosympy:masterfrom
haru-44:gmpy_isqrt

Conversation

@haru-44
Copy link
Copy Markdown
Member

@haru-44 haru-44 commented May 7, 2023

sympy.external.gmpy.sqrt used to use mpmath when gmpy was not used, but math.isqrt has been implemented since Python 3.8, so this is now used. Also, sympy.core.power.isqrt now calls sympy.external.gmpy.sqrt.

References to other Issues or PRs

#25067

Brief description of what is fixed or changed

Other comments

Release Notes

NO ENTRY

@haru-44
Changed the function called by sqrt from mpmath to math
0a0a5bc 
`sympy.external.gmpy.sqrt` used to use mpmath when gmpy was not used,
but `math.isqrt` has been implemented since Python 3.8, so this is now used.
Also, `sympy.core.power.isqrt` now calls `sympy.external.gmpy.sqrt`.
@sympy-bot
Copy link
Copy Markdown

sympy-bot commented May 7, 2023

✅

Hi, I am the SymPy bot. I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

  • No release notes entry will be added for this pull request.
Click here to see the pull request description that was parsed. 
`sympy.external.gmpy.sqrt` used to use mpmath when gmpy was not used, but `math.isqrt` has been implemented since Python 3.8, so this is now used. Also, `sympy.core.power.isqrt` now calls `sympy.external.gmpy.sqrt`.

<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->
#25067

#### Brief description of what is fixed or changed


#### Other comments


#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

* physics.units
  * Corrected a semantical error in the conversion between volt and statvolt which
    reported the volt as being larger than the statvolt.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
NO ENTRY
<!-- END RELEASE NOTES -->

@haru-44 haru-44 mentioned this pull request May 7, 2023
Closed
7 tasks
@github-actions
Copy link
Copy Markdown

github-actions bot commented May 7, 2023 •
edited
Loading

Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with + are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with - are speedups.

Significantly changed benchmark results (PR vs master)



Significantly changed benchmark results (master vs previous release)


       before           after         ratio
     [2c3de5f4]       [e0235b4b]
-        83.4±1ms       55.4±0.2ms     0.66  integrate.TimeIntegrationRisch02.time_doit(10)
-        82.5±1ms       54.7±0.2ms     0.66  integrate.TimeIntegrationRisch02.time_doit_risch(10)



Full benchmark results can be found as artifacts in GitHub Actions

(click on checks at the top of the PR).


      		
    

  





        


            

            
            

              
            

        

      


      

            
  
  

    
  
    
    

  

  

  


    












      

  
        





      

  
      



  

  @oscarbenjamin






  



    

      

  

    

      

          
    
    

  


        
            
  
      
              Copy link

  

              
  
      
                Copy Markdown

  

        
      

    


    

        

  
    Collaborator


        

    

  


  

    

      
          
      

      
            oscarbenjamin
  

      

      

      commented


        May 7, 2023


      
    


  





      


        



  
    

      
          
Have you run any timings to see how math.isqrt compares with mpmath's version?


It is not necessarily faster so we should check: #17038 (comment)


It is important to test with both small examples but also extremely large examples and to test both square and non-square numbers.


      		
    

  





        


            

            
            

              
            

        

      


      

            
  
  

    
  
    
    

  

  

  


    












      

  
      



  

  @haru-44






  



    

      

  

    

      

          
    
    

  


        
            
  
      
              Copy link

  

              
  
      
                Copy Markdown

  

        
      

    


    

        

  
  Member


        

  Author


    

  


  

    

      
          
      

      
            haru-44
  

      

      

      commented


        May 7, 2023


      
    


  





      


        



  
    

      
          
The following three functions are compared in time.


from math import sqrt, isqrt
from sympy.core.power import integer_nthroot
import mpmath.libmp as mlib

def isqrt_sympy(n):
    n = int(n)
    if n < 4503599761588224:
        s = int(sqrt(n))
        if 0 <= n - s*s <= 2*s:
            return s
    return integer_nthroot(n, 2)[0]

def isqrt_mpmath(n):
    return mlib.isqrt(n)

def isqrt_math(n):
    return isqrt(n)


python==3.11.1
mpmath==1.2.1
sympy==1.11.1



import time
import matplotlib.pyplot as plt
import seaborn as sns

sympy_time = []
mpmath_time = []
math_time = []

r = list(range(1, 300))
times = 100_000
micro = 1_000_000

for k in r:
    n = 2**k

    start = time.perf_counter()
    for _ in range(times):
        isqrt_sympy(n)
    sympy_time.append((time.perf_counter() - start) / times * micro)

    start = time.perf_counter()
    for _ in range(times):
        isqrt_mpmath(n)
    mpmath_time.append((time.perf_counter() - start) / times * micro)

    start = time.perf_counter()
    for _ in range(times):
        isqrt_math(n)
    math_time.append((time.perf_counter() - start) / times * micro)

sns.lineplot(x=r, y=sympy_time, label='sympy')
sns.lineplot(x=r, y=mpmath_time, label='mpmath')
sns.lineplot(x=r, y=math_time, label='math')
plt.xlabel('log_2 n')
plt.ylabel('elapsed time [micro sec]')
plt.legend()
plt.savefig('001.png')


This result seems to indicate that math.isqrt is faster, but for larger n, the speed seems to be reversed (approximately n=2**55_000 is the boundary). Further study seems to be needed.


001


      		
    

  





        


            

            
            

              
            

        

      


      

            
  
  

    
  
    
    

  

  

  


    












      

  
      



  

  @oscarbenjamin






  



    

      

  

    

      

          
    
    

  


        
            
  
      
              Copy link

  

              
  
      
                Copy Markdown

  

        
      

    


    

        

  
    Collaborator


        

    

  


  

    

      
          
      

      
            oscarbenjamin
  

      

      

      commented


        May 7, 2023


      
    


  





      


        



  
    

      
          
I would be wary about only timing with powers of 2. Computing the square root of a power of 2 means that you have a number with a binary representation like 100000000000000.... As a big integer it will have many trailing zero limbs each of which could be handled very efficiently by isqrt if special cased.




the speed seems to be reversed (approximately n=2**55_000 is the boundary).




I see the same. I expect it is because the algorithm is different as highlighted by @mdickinson in #17038 (comment):




There's also a good chance that the pure Python variant of isqrt in mpmath is faster than math.isqrt - it's significantly more optimised and complicated, and doesn't depend on slow divisions for large inputs.




For now it looks like we can do:


def isqrt_math_mpmath(n):
    if n.bit_length() < 55_000:
        return isqrt(n)
    else:
        return isqrt_mpmath(n)


Actually though it would be better to add that logic to mpmath. I think that SymPy should use mpmath's isqrt but mpmath should be changed to use math.isqrt for smaller inputs.


      		
    

  





        


            

            
            

              
            

        

      


      

            
  
  

    
  
    
    

  

  

  


    












      

  
      



  

  @mdickinson






  



    

      

  

    

      

          
    
    

  


        
            
  
      
              Copy link

  

              
  
      
                Copy Markdown

  

        
      

    


    

        

  
    Contributor


        

    

  


  

    

      
          
      

      
            mdickinson
  

      

      

      commented


        May 7, 2023


      
    


  





      


        



  
    

      
          
Out of curiosity, what version of Python were the timings run with? It could be interesting to try with Python 3.12, where integer division is no longer asymptotically quadratic.


      		
    

  





        


            

            
            

              
            

        

      


      

            
  
  

    
  
    
    

  

  

  


    












      

  
      



  

  @oscarbenjamin






  



    

      

  

    

      

          
    
    

  


        
            
  
      
              Copy link

  

              
  
      
                Copy Markdown

  

        
      

    


    

        

  
    Collaborator


        

    

  


  

    

      
          
      

      
            oscarbenjamin
  

      

      

      commented


        May 7, 2023


      
    


  





      


        



  
    

      
          


It could be interesting to try with Python 3.12




I tried this with Python 3.12a7 to time sympy, mpmath, and math.isqrt up to a million bits:


from math import sqrt, isqrt
from sympy.core.power import integer_nthroot
import mpmath.libmp as mlib
import random
import time

def isqrt_sympy(n):
    n = int(n)
    if n < 4503599761588224:
        s = int(sqrt(n))
        if 0 <= n - s*s <= 2*s:
            return s
    return integer_nthroot(n, 2)[0]

def isqrt_mpmath(n):
    return mlib.isqrt(n)

def isqrt_math(n):
    return isqrt(n)

sympy_time = []
mpmath_time = []
math_time = []

r = sorted(random.sample(range(1, 1_000_000), 200))

times = 10
micro = 1_000_000

for k in r:
    ns = [random.randint(0, 2**k) for _ in range(times)]

    start = time.perf_counter()
    for n in ns:
        isqrt_sympy(n)
    sympy_time.append((time.perf_counter() - start) / times * micro)

    start = time.perf_counter()
    for n in ns:
        isqrt_mpmath(n)
    mpmath_time.append((time.perf_counter() - start) / times * micro)

    start = time.perf_counter()
    for n in ns:
        isqrt_math(n)
    math_time.append((time.perf_counter() - start) / times * micro)

print('r = ', r)
print('sympy_time = ', sympy_time)
print('mpmath_time = ', mpmath_time)
print('math_time = ', math_time)


Then switching back to Python 3.11 with matplotlib and gmpy2 to time gmpy2 (with the same bitsizes but not the same numbers) and plot the results:


import matplotlib.pyplot as plt
from gmpy2 import isqrt
import random
import time

def isqrt_gmpy2(n):
    return isqrt(n)

# Output from previous script in 3.12:
r =  [5676, 13862, 24534, 27734, 38915, 41670, 43728, 46537, 47615, 49432, 58199, 64848, 73910, 77901, 83966, 85451, 92757, 94695, 103240, 114262, 119448, 121780, 125952, 145384, 148254, 153651, 157144, 157187, 162258, 182326, 182945, 201953, 209356, 211665, 215130, 216537, 218991, 228517, 229803, 229916, 236122, 241872, 247667, 250866, 255019, 258817, 264498, 266655, 266731, 266781, 268135, 273375, 276701, 277938, 286679, 290686, 292456, 296262, 299173, 300593, 304171, 306946, 307144, 308352, 309559, 309666, 314840, 329643, 331021, 339015, 347558, 348377, 357612, 357770, 358304, 364235, 365660, 372699, 376861, 377172, 379422, 380414, 381063, 382590, 385746, 396540, 399069, 405709, 411611, 418945, 421853, 430425, 441141, 447515, 450898, 464118, 487991, 488742, 495884, 497115, 497632, 510090, 511857, 515564, 516082, 525101, 538990, 543763, 545064, 551397, 552869, 553380, 555205, 558100, 574808, 578718, 578748, 584399, 587409, 592792, 621746, 623590, 631034, 643440, 648118, 657576, 658709, 658783, 660505, 666066, 666802, 666871, 667279, 668663, 672109, 677728, 683136, 699154, 700497, 701733, 706240, 710378, 715499, 720039, 725988, 732666, 735130, 736628, 739047, 739373, 741365, 744504, 750703, 752238, 757187, 764388, 765109, 774230, 787616, 794845, 798436, 798888, 800614, 809784, 814593, 827739, 832164, 837520, 838006, 838490, 840884, 849899, 859944, 860778, 864389, 869347, 884833, 885800, 892924, 893673, 900690, 903683, 906446, 907052, 917592, 926022, 941967, 946440, 948254, 951815, 955220, 955865, 959230, 959727, 962781, 971734, 972831, 979391, 983269, 991004]
sympy_time =  [76.43100007044268, 268.65529998758575, 652.4931000058132, 789.333499960776, 1318.4137000280316, 1479.1998000873718, 1591.510400066909, 1768.6190999484097, 1823.4338999718602, 1941.6392000493943, 2573.2961999892723, 3154.4184999802383, 3678.174099968601, 3952.380000009725, 4507.442799967976, 4592.549599965423, 5246.924800030683, 5420.893799964688, 6238.697100070567, 7451.059899995016, 7947.937799963256, 8341.434500016476, 8694.5284999274, 10480.679899956158, 10971.514999982901, 11359.07430007137, 11891.0393999613, 12162.807899949257, 12790.080200011289, 14952.526399974886, 15156.803300033062, 18240.495299960457, 19195.857499926205, 19768.139400002838, 19910.09869998379, 20193.666300019686, 20920.574400042824, 23801.170500064472, 26446.027400015737, 29325.380800037237, 26663.973400081886, 26273.58189993174, 26083.15040006346, 26415.26170000361, 27459.476300009555, 43625.96169994504, 32039.410100060195, 32140.72780001516, 32675.480100078854, 32152.487100029248, 31749.814899922054, 29303.61589997119, 35958.82330000677, 30208.00379999855, 51466.301799973735, 43467.53390000231, 34286.949199940864, 35376.342899962765, 35643.53669999036, 35508.10540000384, 38576.64050001404, 37902.75240007759, 37502.0241999664, 35630.39810005648, 40036.19839995736, 37369.223200039414, 39613.117199951375, 40356.555600010324, 39519.60209997196, 44447.70079999216, 45534.68100002647, 45525.63799998097, 46019.31669994883, 45090.59199999683, 46623.14290007998, 48431.83600005432, 49000.079199959146, 49191.82499997987, 49950.5402000068, 51893.74170004157, 59561.918700001115, 56069.21690005038, 49945.44700002734, 50930.080600028305, 50303.090999932465, 52938.470499975665, 53159.256399976584, 55954.70670004943, 56924.4286999492, 58326.38589999988, 59747.173600044334, 61994.41720000323, 64336.53120002418, 65670.10070002652, 65303.536600004016, 69396.73249999032, 76553.1146000285, 77519.06150006107, 77695.14760002494, 78838.44759999192, 81035.92409997873, 81588.73130005304, 83833.95440005188, 82268.30500007054, 81755.60350000524, 87818.3720999914, 84587.53120003166, 80158.57420004977, 81274.4343000304, 83365.86109999189, 84499.35400003596, 85162.34920007264, 83045.73909999817, 83512.83829997556, 95625.66399999923, 94686.28829999943, 94684.01160002031, 98372.88839999019, 98456.24630006569, 106691.34400004623, 111243.6446999709, 110749.88620002841, 108681.77490001472, 115302.93269997856, 115010.89559997126, 118945.7100999789, 120016.8078000388, 118990.2684000117, 120213.74720006861, 119528.035500025, 121517.04650004831, 121244.13689998619, 117929.56479994245, 121266.23429994652, 123122.5538000217, 123659.87589992073, 124263.29749996512, 131035.22379997229, 128809.58530004136, 132069.7038000617, 134614.64060001163, 134046.43849999047, 132174.92020003192, 135859.10679994413, 139947.82760000817, 136026.67249997467, 141891.59119996475, 142411.1676999928, 142460.68369993736, 144299.2403000062, 148975.89689999222, 223166.36750001635, 243599.11989995453, 146130.69540000652, 148247.379500026, 151499.1648999967, 154984.88310004177, 159141.90100002088, 152913.34480007208, 163063.76379998255, 160725.37009995358, 168197.61440001457, 172434.79649996516, 168330.58709999023, 170267.02060002208, 174346.0065000363, 170924.40499991426, 174000.5156000734, 174136.39690003038, 172363.8456000117, 173591.1759999908, 180051.9145999715, 188018.34060004694, 189726.3675999966, 185704.70029999342, 187436.18580001566, 190457.3949999758, 193131.42669998342, 213800.21599998145, 196155.05679994385, 194888.3475999537, 201789.17869998259, 200857.6331999393, 200579.9837000268, 199483.52329993213, 207831.66129995152, 212728.27319999124, 215464.99579999363, 217928.85360000582, 222085.31289998064, 222129.58700001764, 219316.60430000193, 216856.84759995638, 222010.25490003303, 223598.4508999536, 227914.7601000659, 227841.90880001916, 225763.38629996826, 233641.03670000986, 236984.2143999449]
mpmath_time =  [69.43509997654473, 262.72040004187147, 638.6766999639804, 791.2547999694652, 1312.9036000464112, 1599.8642000340624, 1572.8258000308415, 1756.1177000061434, 1816.8060999414593, 1936.2775000445254, 2698.116600004141, 3111.579200049164, 3652.9293999592483, 3951.7639000223426, 4423.9586999538005, 4625.098700034869, 5280.245700032538, 5432.665700027428, 6323.520299974916, 7475.593100025435, 7949.640500009991, 8471.96039994742, 8714.551599950937, 10560.999999961496, 10887.730499962345, 11334.643200007122, 11750.38469991705, 12109.38129997885, 12701.351499981683, 15412.28880005292, 15134.33210002404, 18087.84889999515, 19047.313099963503, 19973.879699955432, 19881.04640004167, 20595.830000002024, 20984.545400006027, 23737.67900007806, 22822.612799973285, 31387.005800024777, 24581.918099920586, 24729.634699997405, 25994.114900004206, 26286.68409997772, 27318.65169998855, 28582.7061999953, 30970.112700015306, 32108.66380004518, 30986.860199936928, 33324.97439996587, 31859.522700051457, 28823.99030004308, 31011.48599998851, 29583.588300010888, 45225.833899985446, 39910.662100010086, 33657.58790005202, 36124.28189999264, 35493.59000007826, 37014.13359995058, 37141.45379999536, 40471.402399998624, 36709.07209998404, 36291.76939994068, 37207.71789994615, 36470.46960004445, 39099.61419994943, 39973.64300003028, 39374.090199999046, 45879.938599955494, 43589.12349998718, 45908.25879995464, 45254.05719996343, 43920.19690003508, 46072.854599970015, 48756.89049995344, 47881.10569998025, 52902.16980001787, 51624.61279996933, 54077.357100049994, 53836.58379996632, 52870.35709998236, 49888.208100037446, 51740.93689993242, 50640.97029999175, 52451.271599966276, 52749.176099951, 55534.397100018396, 57542.10469995087, 58229.27440003696, 59538.68690003219, 64173.8393000196, 64525.73990000018, 66096.50070004136, 65631.43060002403, 68916.51400001138, 76280.06830000231, 76855.92830002861, 77575.55550006145, 81532.80229998927, 79422.43200004668, 81468.20379997735, 83388.12550000512, 82567.12940001307, 81620.24529992777, 87597.79360007087, 84885.11339992328, 80221.48580002977, 80712.13450002688, 83409.23480000129, 84135.821900054, 82603.6651999857, 83970.09649997926, 84509.48370000333, 94988.27329998676, 95020.25390002018, 94640.27539997915, 99176.72960000345, 97944.9678000492, 102907.21229994233, 119387.5549999575, 107161.22929998164, 108982.98069996599, 113916.22060000373, 116187.99360003322, 118919.32499993345, 118814.42559997595, 119001.80640004692, 120573.21859992953, 119081.51389998238, 121115.14889993487, 120749.08139993568, 119590.18509996895, 121530.07870001602, 122503.472399967, 123883.77130000663, 124472.95940000913, 131513.29190004617, 130972.71510005157, 132000.09960000898, 131190.57090007118, 136362.41539998082, 132702.52680003978, 136066.8333999456, 140286.17370004213, 136560.79480006156, 142452.76559995546, 142282.64949997538, 144270.80300001762, 147472.67229995487, 173051.17939995398, 232039.0290999967, 210736.96520006706, 145721.05819997887, 148356.87139993752, 151942.26019993948, 153275.4572000158, 156068.23419993816, 156268.1500000508, 163361.89710000326, 160561.4512000102, 168275.754699971, 171761.7997999696, 167021.59289998235, 169492.54219998693, 174291.6250000235, 171143.9681000229, 175228.1271999891, 174463.74829996785, 172775.3147999465, 173398.51970000382, 178362.30670000077, 188424.95359995155, 185429.36170006215, 184249.23760003367, 187320.9977999977, 190298.76490003517, 193664.07019997496, 195147.90049997828, 195701.29109997652, 194307.17499999446, 201178.79179997544, 200779.64580004846, 200467.64999997322, 199279.81539995017, 210367.31479998707, 214825.0459999872, 216047.1394000524, 218293.06829995403, 221774.4446999859, 221527.01859995432, 218525.1096999309, 220408.07609992044, 223787.08779997396, 226744.82369993711, 227726.4371000456, 227351.30519995437, 225372.05549997452, 236307.22779998905, 234974.3572999614]
math_time =  [28.32870004567667, 120.63430003763642, 309.2103000199131, 387.43449995308765, 769.5067999520688, 840.6287000070733, 869.4179000485747, 967.0434999861754, 993.7395000633842, 1058.5403000732185, 1483.7442000498413, 1766.6712000391271, 2100.9709000281873, 2276.2783999496605, 2552.276099959272, 2666.295499966509, 3004.3145000490767, 3111.7676000576466, 3633.1768999843916, 4244.2048999873805, 4570.793099992443, 4776.273900006345, 5176.04790002224, 6316.2188999740465, 6455.522200030828, 6857.673000013165, 7068.838300074276, 7218.572800047696, 7572.051100032695, 8959.394699922996, 9203.830000024027, 10555.238200049644, 11227.720099941507, 11669.757399977243, 11893.722800050455, 11905.372700039152, 13082.415299959393, 13655.178000044543, 13124.880899977143, 15464.625999993586, 13834.353099991858, 16019.606300051237, 15018.909399987024, 15274.212199983594, 18655.048399978114, 18381.727999985742, 28807.106600015686, 18679.366200012737, 18983.472000036272, 18894.491700029903, 18632.475799950043, 20392.36290002009, 18311.582500064105, 17959.79200005604, 32525.638800052548, 25809.523699990677, 21876.109200002247, 20515.083000009326, 21677.733800061105, 25240.82999998427, 21523.8994000174, 23210.67820003009, 24375.40250002712, 23413.185800018255, 23328.026900071563, 21708.64330000768, 22428.153900000325, 23905.12289994149, 26285.841899971274, 24888.019500031078, 26364.76039997433, 25524.21630007302, 27175.639500001125, 30003.797999961535, 28970.09190000972, 29584.82570002161, 28580.61449996967, 29461.653100042895, 29894.856499959133, 29633.978299989394, 29400.0822000271, 31413.52000002371, 31232.662000002165, 29643.81580004556, 30036.81210002469, 32606.050400045202, 31919.720000041707, 32682.47159994644, 33554.080900012195, 34789.88059996482, 34687.721699992835, 36503.6537999913, 37760.85259996762, 38600.19570001896, 38949.27199999074, 40874.94469995363, 45575.08730003974, 45378.17769996764, 46894.20340000652, 48128.93759999497, 47378.1302000134, 48905.090399966866, 48947.18280002053, 50074.72770003005, 50166.57989999658, 52002.072499999485, 51208.463800048776, 48841.52020003967, 51827.3837000379, 50203.87210006447, 51023.96319998661, 50712.576799924136, 50418.33900004349, 50892.18749999418, 56534.82309999163, 57219.19759998855, 57419.33330000393, 58473.1897999518, 61175.69170000934, 60378.20689998625, 64303.81399995895, 64624.430899948486, 65654.56509997603, 67777.4271999624, 68301.22459996346, 71192.23569998212, 70440.25519999195, 70991.72530006399, 71608.15340002955, 71272.23460001915, 75627.90460005999, 71078.53659999819, 88812.60629996177, 71810.89400000928, 72556.58840003889, 74052.0485000161, 74217.71669996815, 77904.38639995045, 78124.91630002114, 77944.20820000596, 78815.49959993208, 79286.71970003052, 81471.00770001998, 81946.11230001101, 82783.19729997747, 83024.89069992589, 84379.20890000896, 84554.7080999495, 85426.3223000089, 85843.80380007133, 131179.65640003604, 129167.96310000791, 123574.5921000671, 87320.58349996805, 90817.4563999637, 90887.90900004824, 90383.8072000326, 94335.03919999566, 94089.21390004252, 95773.22169998297, 96105.29869996753, 100062.4042999334, 96422.04649999258, 103640.46459999372, 99791.46870000477, 102807.94820000665, 102525.8195000788, 103610.08069994568, 104963.44390003287, 106158.3475999214, 104625.9894999821, 107156.53870001915, 109741.35550004576, 109984.6441000409, 110038.10150004938, 110588.4696999965, 114403.75310003219, 114318.6844000411, 115799.18890001863, 115485.60239998551, 117009.55749993227, 119570.7171000322, 120778.69240001746, 120989.3109999939, 121655.54659995905, 123491.11929997889, 126834.22410000276, 127929.6465000698, 128843.69380008137, 130512.51420001792, 129819.35389998398, 130560.53199998132, 131543.0951000053, 130310.47859994943, 131317.31909998052, 134311.42710005588, 134495.13460000162, 136669.98930002592, 139281.7663000642, 140939.7929000079]

gmpy2_time = []

times = 10
micro = 1_000_000

for k in r:
    ns = [random.randint(0, 2**k) for _ in range(times)]

    start = time.perf_counter()
    for n in ns:
        isqrt_gmpy2(n)
    gmpy2_time.append((time.perf_counter() - start) / times * micro)

print('r = ', r)
print('sympy_time = ', sympy_time)
print('mpmath_time = ', mpmath_time)
print('math_time = ', math_time)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(r, sympy_time, label='sympy')
ax.plot(r, mpmath_time, label='mpmath')
ax.plot(r, math_time, label='math')
ax.plot(r, gmpy2_time, label='gmpy2')
ax.set_yscale('log')
ax.set_xlabel('bits')
ax.set_ylabel('time (microseconds)')
ax.legend()
plt.show()
fig.savefig('isqrt.png')


isqrt

So it looks like math.isqrt is about 2x faster than mpmath but still about 20x slower than gmpy2's version. For the tested range there doesn't seem to be an obvious asymptotic difference between the three but I guess we're probably not up to FFT-sized inputs.


I guess that means three things:


    

  1. If gmpy2 is installed it should always be used.
    2. 

  2. For 3.12 upwards math.isqrt is best if gmpy2 is not installed.
    3. 

  3. Both sympy and mpmath should be changed to use either gmpy2 or math.isqrt for version 3.12 upwards.
    4. 



It would be good to test even larger like 10**9 bits but it takes a long time on this computer.


The other interesting comparison is always very small numbers since those are probably most common anyway but I'm pretty sure that the conclusions there are the same: use gmpy2 if possible or otherwise math.isqrt. The only difference with small numbers is that math.isqrt should always be used for small numbers regardless of Python version.


      		
    

  





        


            

            
            

              
            

        

      


      

            
  
  

    
  
    
    

  

  

  


    












      

  
      



  

  @mdickinson






  



    

      

  

    

      

          
    
    

  


        
            
  
      
              Copy link

  

              
  
      
                Copy Markdown

  

        
      

    


    

        

  
    Contributor


        

    

  


  

    

      
          
      

      
            mdickinson
  

      

      

      commented


        May 8, 2023


      
    


  





      


        



  
    

      
          
@oscarbenjamin Very interesting; thank you for doing the testing. I've been meaning to retest the performance of math.isqrt ever since python/cpython#96673 got merged, but never got around to it.


That gmpy2 massively outperforms math.isqrt is, of course, expected.


      		
    

  





        


            

            
            

              
            

        

      


      

            
  
  

    
  
    
    

  

  

  


    












      

  
      



  

  @oscarbenjamin






  



    

      

  

    

      

          
    
    

  


        
            
  
      
              Copy link

  

              
  
      
                Copy Markdown

  

        
      

    


    

        

  
    Collaborator


        

    

  


  

    

      
          
      

      
            oscarbenjamin
  

      

      

      commented


        May 8, 2023


      
    


  





      


        



  
    

      
          


That gmpy2 massively outperforms math.isqrt is, of course, expected.




They look asymptotically similar though which is good. Which operations exactly have changed in 3.12?


I can't see the _pylong.py module in cpython main. Did it get moved somewhere?


      		
    

  





        


            

            
            

              
            

        

      


      

            
  
  

    
  
    
    

  

  

  


    












      

  
      



  

  @haru-44






  



    

      

  

    

      

          
    
    

  


        
            
  
      
              Copy link

  

              
  
      
                Copy Markdown

  

        
      

    


    

        

  
  Member


        

  Author


    

  


  

    

      
          
      

      
            haru-44
  

      

      

      commented


        May 8, 2023


      
    


  





      


        



  
    

      
          
I'm trying to revert gmpy.py to an previous version, do I need to cast mlib.isqrt(x) to int?


      		
    

  





        


            

            
            

              
            

        

      


      

            
  
  

    
  
    
    

  

  

  


    












      

  
      



  

  @mdickinson






  



    

      

  

    

      

          
    
    

  


        
            
  
      
              Copy link

  

              
  
      
                Copy Markdown

  

        
      

    


    

        

  
    Contributor


        

    

  


  

    

      
          
      

      
            mdickinson
  

      

      

      commented


        May 8, 2023


      
    


  





      


        



  
    

      
          


Which operations exactly have changed in 3.12?




A Python (yes, really!) fast path was added for int<->decimal string and int<->Decimal conversions, as well as divmod.




I can't see the _pylong.py module in cpython main. Did it get moved somewhere?




It's here. I don't think it got moved.


      		
    

  





        


            

            
            

              
            

        

      


      

            
  
  

    
  
    
    

  

  

  


    












      

  
      



  

  @oscarbenjamin






  



    

      

  

    

      

          
    
    

  


        
            
  
      
              Copy link

  

              
  
      
                Copy Markdown

  

        
      

    


    

        

  
    Collaborator


        

    

  


  

    

      
          
      

      
            oscarbenjamin
  

      

      

      commented


        May 8, 2023


      
    


  





      


        



  
    

      
          




Which operations exactly have changed in 3.12?




A Python (yes, really!) fast path was added for int<->decimal string and int<->Decimal conversions, as well as divmod.




Okay so no FFT multiplication yet? I'm sure I saw an implementation of that somewhere.


Presumably then isqrt became faster because it uses divmod. Is it fair to say that the mpmath implementation is largely a workaround for poor asymptotics of divmod? Or is it something that could also be improved by taking advantage of faster divmod (I haven't looked at how it actually works)?


      		
    

  





        


            

            
            

              
            

        

      


      

            
  
  

    
  
    
    

  

  

  


    












      

  
      



  

  @oscarbenjamin






  



    

      

  

    

      

          
    
    

  


        
            
  
      
              Copy link

  

              
  
      
                Copy Markdown

  

        
      

    


    

        

  
    Collaborator


        

    

  


  

    

      
          
      

      
            oscarbenjamin
  

      

      

      commented


        May 8, 2023


      
    


  





      


        



  
    

      
          


I'm trying to revert gmpy.py to an previous version, do I need to cast mlib.isqrt(x) to int?




We would only use it if gmpy2 is not installed in which case mpmath would be using plain int as well. Although I suppose it is also possible that someone could set SYMPY_GROUND_TYPES to disable use of gmpy2 in sympy but mpmath would still be using it.


      		
    

  





        


            

            
            

              
            

        

      


      

            
  
  

    
  
    
    

  

  

  


    












      

  
      



  

  @haru-44






  



    

      

  

    

      

          
    
    

  


        
            
  
      
              Copy link

  

              
  
      
                Copy Markdown

  

        
      

    


    

        

  
  Member


        

  Author


    

  


  

    

      
          
      

      
            haru-44
  

      

      

      commented


        May 8, 2023


      
    


  





      


        



  
    

      
          
Okay, the cast to int will remain.


      		
    

  





        


            

            
            

              
            

        

      


      

            
  
  

    
  
    
    

  

  

  


    












      

  
        





      

  
          

  

    
  


  

      


    @haru-44
haru-44



    mentioned this pull request
    
      May 8, 2023
    





  


  

      
   Merged

  







  









      

  
      



  

  @mdickinson






  



    

      

  

    

      

          
    
    

  


        
            
  
      
              Copy link

  

              
  
      
                Copy Markdown

  

        
      

    


    

        

  
    Contributor


        

    

  


  

    

      
          
      

      
            mdickinson
  

      

      

      commented


        May 8, 2023


      
    


  





      


        



  
    

      
          
@oscarbenjamin




Okay so no FFT multiplication yet?




Right. Maybe someday.




Presumably then isqrt became faster because it uses divmod.




Yes - for large inputs, the math.isqrt running time is (overwhelmingly) dominated by the final division, which is a call to PyNumber_FloorDivide.


I haven't looked at the mpmath implementation, so I don't have a basis for comparison.


      		
    

  





        


            

            
            

              
            

        

      


      

            
  
  

    
  
    
    

  

  

  


    












      

  
      



  

  @oscarbenjamin






  



    

      

  

    

      

          
    
    

  


        
            
  
      
              Copy link

  

              
  
      
                Copy Markdown

  

        
      

    


    

        

  
    Collaborator


        

    

  


  

    

      
          
      

      
            oscarbenjamin
  

      

      

      commented


        May 11, 2023


      
    


  





      


        



  
    

      
          
Looks good, merging.


      		
    

  





        


            

            
            

              
            

        

      


      

            
  
  

    
  
    
    

  

  

  


    












      

  
          
 
  

    
  

  



        @oscarbenjamin
  oscarbenjamin



        merged commit 26ad10b
      into

      
  
    sympy:master
  


    May 11, 2023

  







  

  
  

  
      @haru-44
  haru-44


    
    deleted the
    
      
        gmpy_isqrt
    
    branch

    May 14, 2023 14:24    
    













  
  

    

      

  










  








  
      

  
    Sign up for free
    to join this conversation on GitHub.
    Already have an account?
    Sign in to comment


  










      

  

    
    

    

  

    Reviewers
  


    

    No reviews






    
    
    

    

  


      
  

    Assignees
  



        


    No one assigned







      


  


  

    Labels
  



    

    None yet







      

  

    

        

    Projects
  


        




    None yet




  



      


  

    
  

    Milestone
  


      No milestone





    
      
          


  

    

      
        

          
  

    Development
  



            


Successfully merging this pull request may close these issues.





  
  



      
    

  




      

    

  
        

  

    

      4 participants
    

    

        
          @haru-44 
        
          @sympy-bot 
        
          @oscarbenjamin 
        
          @mdickinson