New mid-range backend and dispatch code for polynomial multiplication

[fredrik-johansson]

This large PR does a few related things:

• Take the packed-mpn classical+Karatsuba polynomial multiplication code from the mpn_mod module, move it to mpn_extras, optimize it, and generalize it to signed coefficients and unbalanced sizes. This multiplication code is now usable by both mpn_mod_poly and fmpz_poly, and we add fmpz_poly_mul_mpn/fmpz_poly_mulmid_mpn that convert to packed mpn format internally.

• Also adds fmpz_poly_mul(mid)_toom_karatsuba -- recursive Karatsuba with toom_scalar as the basecase.

• Improve and consolidate the fmpz_poly_mul* and sqr* algorithm dispatch. All multiplications (mul, mullow, sqr, sqrlow, mulmid) now go through _fmpz_poly_mulmid. This dispatcher fixes some inaccuracies in the previous tuning and adds use of mpn, toom_scalar, toom_karatsuba and classical_fft_small where these win.

• As a result of the refactoring, _fmpz_poly_mul and _fmpz_poly_mullow now accept their arguments in either order. Several callers have been simplified accordingly.

The bulk of patch is thousands of lines of generated kernel mpn code: m x n limb unsigned and signed dot products and basecase polynomial multiplication routines. The dot products are meant to be reusable for other things. Note that these kernel routines very likely could be optimized by generating x86-64 assembly instead of C.

To compensate a bit for the added code, all the fmpz_poly_mul*_karatsuba functions are removed. These were completely obsolete; they are always slower than toom_karatsuba or mpn.

The cyclotomic multiplication code in the APRCL module is also refactored and optimized.

Developed using several Claude models, mainly Fable 5, Opus 5 and Fable 5.1.

Performance impact

Gives more than a 2x speedup for fmpz_poly arithmetic if you happen to do lots of short to medium length multiplications exactly in one of the sweet spots for the new algorithms (~100-1000 bits for mpn; several thousand bits and up for toom_scalar and classical_fft_small). On realistic, complex tasks one is more likely to see a 5-10% speedup. Operations on very short polynomials with small coefficients and short sparse polynomials can get slightly slower.

Speedup on tests/examples

Very noisy, best of 5 timings, take with a grain of salt.

                                          old      new   speedup
build/examples/huge_expr -ca            1.371    1.316    1.042x
build/examples/huge_expr                2.153    2.045    1.053x
build/fmpz_poly_factor/tests/main       0.353    0.344    1.026x
build/fmpq_poly/tests/main              2.547    2.363    1.078x
build/fmpz_mod_poly/tests/main          1.613    1.561    1.033x
build/fmpz_mod_poly_factor/tests/main   0.706    0.655    1.078x
build/fmpz_poly_mat/tests/main          1.788    1.564    1.143x
build/nf_elem/tests/main                1.103    0.959    1.150x
build/qqbar/tests/main                  2.285    2.205    1.036x
build/fq/tests/main                     0.295    0.311    0.949x
build/fq_poly/tests/main                1.133    1.086    1.043x
build/fq_poly_factor/tests/main         2.018    1.937    1.042x
build/fq_mat/tests/main                 0.808    0.800    1.010x

The fq module test slowdown is real, but I think it's not necessarily representative of real usage (note that fq_poly and fq_mat tests got faster).

Speedup for fmpz_poly_mul

Speedup x 100 to multiply two length-n polynomial with randbits(bits) coefficients, e.g. 123 means the new code is 1.23x faster.

 bits \ n     1   2   3   4   5   6   7   8   9  10  11  13  15  18  21  25  30  36  43  51  61  73  87 104 124 148 177 212 
       20   123  94  94  88  98 102 103  94  98 100 102  98 101 104  99 101 103  98 101 101  99 102 100 100 100 100 100  99 
       60   113  99 102 103 106 105 105 106 107 105 104 108 108 106 106 105 105 104 104 109 103 103 110 100 100 100 100 100 
      100   110  85  97  96 113 129 221 222 217 221 207 208 181 169 162 153 145 144 136 130 122 127 100 100 100  99  99 100 
      200    98  84  92  94  96  98 186 148 142 128 124 124 123 122 123 119 120 121 123 103 102  91 100 100 100 100 100 100 
      400    99  87  97  98  99  95 184 239 256 226 207 187 173 177 170 155 139 137 132 121 111 155 100 100 100 100 101 100 
      800   100  95  98  99  98 101 130 121 174 164 158 147 141 132 113 106 102 112 113 118 105 117 100 100 100 100 101 101 
     1000   100 108  95 105  90 107 119 112 155 147 140 131 129 114 102 102  97 110 109  99  97 107 101 100 100 100 100 100 
     1200   101 112 105 120  95 113 111 104 147 138 134 126 124 109 110 110 109 125 107 107 103 119 100 100 100  99 100 100 
     1500   101 119 120 133 118 142 113 105 151 142 137 124 119 116 115 117 102 115 109 112 104 126  99 100 101 101 101 100 
     2000   100 120 130 144 132 158 111 103 156 147 143 128 121 128 129 113 115 120 125 119 118 132 100 100 100 100 100 100 
     3000    99 130 148 163 174 186 120 111 117 117 104 100 100 100 101 100  99 100 100 100 101 100 100 101 100 101 100 101 
     4000   101 129 152 172 187 202 146 124 108 110 100 100 100 100 100 100 100  99 101 101 100 101 101 100 100 100 101 100 
     6000   100 111 162 121 150 137 126 135 124 117 100  99 100 101 100 100 101 100 100 101 100 100 100 101 100 102 101 101 
    10000   100 105 159 120 135 121 115 112 101 100 101 100 100 100 101 100 100 100 100 101 101 101 100 100 101 100 100 100 
    20000   100 103 156 120 112 103 101 100 100 100 100 101 100 100 100 100 101 101 100 101 101 100 102 100 100 100 101  99 
    40000   101 101 155 120 110 103 100  99  99 100 100  99 100 101  98 100 100 100  99 101 100 100 100 101 102 101  99 100 
   100000   100 101 156 119 103 100 110 100  91 101  97 103 107 102 101 100 100 100 100 100 100 100 100 100 100 100 100 100 
   200000    98 101 155 120 106 104 130 114 108 109 112 115 115 100 100 102 100 100 100 100 100 100 100 100 100 100 100 100 
   400000   101 100 154 120 118 123 130 121 128 124 120 121 115 100 100 100 100 100 100 100  97  99 100 100  99 100 100 100 
  1000000   100 102 156 121 135 131 128 122 111 127 127 111 108 100 100 100 100 100 100 101 100 100 100 100 100 100 100 100 

Speedup for fmpz_poly_mullow

Length n inputs, product truncated to length n.

 bits \ n     1   2   3   4   5   6   7   8   9  10  11  13  15  18  21  25  30  36  43  51  61  73  87 104 124 148 177 212 
       20   107  95  93  96 111 106 104  94  96  96  95  98  98  98 100  99 100  99  99  99 100 100 100 100 100  99 100  99 
       60   110 115 116 124 124 121 126 121 127 127 126 127 128 126 127 127 128 127 128 120 127 127 127 128 127 164 171 100 
      100   110  88  94  80 105 110 222 232 229 226 221 207 208 217 208 194 182 168 165 155 143 135  99 100 100 100 101 100 
      200   108  92  97  96  96  99 201 195 184 181 175 160 146 147 147 146 138 132 151 121 124 108  99 100 100 100 100 100 
      400   104  94  96  97  98  92 231 388 423 381 355 317 283 266 233 210 184 197 177 155 143 177 100 101 100 100 101 100 
      800   102  96  97  99  98  96 155 127 273 245 224 195 160 195 172 148 129 133 128 134 120 132 100 102 101 100 100 101 
     1000   101  92  95  97  98  98 141 115 229 206 188 162 139 165 147 126 109 113 125 112 110 120 101 100 101 101 100  99 
     1200   101  92  95  96  97  98 130 107 210 189 173 150 129 157 140 120 103 105 121 120 115 136  98 101 100 101 100  99 
     1500   101  93  97  97  99  99 123  99 204 180 177 149 137 149 129 132 114 131 130 126 116 139 100 101 100  99 100 101 
     2000   100  97  98  99  99  99 115  93 211 189 182 156 148 153 147 127 131 136 139 132 132 149 100 101 100 100 100 101 
     3000   100  99  99  99 100 100 119  97 125 110  99  99 100 100 100 100 100  99 100 100 101 100 100 100 100 100 100 101 
     4000   100  99 100 100  99 100 114  93 114 112  98 100 101 100  98 100 100 100  99 101 101 100 100 101 100 100 100 100 
     6000   100  95 103 111 135 122 114 121 115 106  95 102 101 100 100 100 100 100 100 100 100 100 100 101 100 101  99 101 
    10000   100  98 109 121 103 113 106 115 129 117 119 113 114 110 109 100 100 100 101 100 100 101 100  99 102 100 100 101 
    20000   100 135 112 151 129 129 128 126 133 131 127 132 120 119 116 111 101 107  96  99 100  99 101 101 101 100 101 100 
    40000   101 137 134 165 153 139 122 141 136 140 141 128 127 129 131 121 113 116 107  99  97 102 100  99 102 100 100 101 
   100000   100 136 174 212 156 129 153 152 145 149 142 141 134 134 128 131 124 117 105 103  96  99 100 100 100 100 100 100 
   200000   100 136 141 173 137 144 143 156 156 149 146 156 152 140 139 140 135 124 109 102  93 100 100  99  99  99 100  99 
   400000   100 135 132 175 166 159 178 173 159 158 162 159 160 144 139 132 122 112 102  96  87  99  99 100 100 100 100 100 
  1000000   100 130 150 189 166 176 166 154 161 157 168 146 145 135 131 132 124 116 109 102  94 100  99 101 100 100 100 100 

Speedup for fmpz_poly_sqr

 bits \ n     1   2   3   4   5   6   7   8   9  10  11  13  15  18  21  25  30  36  43  51  61  73  87 104 124 148 177 212 
       20   103  79  92  91  92 105  85  85  85  89  90  93  95  94  97  96  97  98  98  99  99  99 100 100  99 100  99  99 
       60   101 101 105 109 109 110 112 113 116 115 118 120 121 123 122 122 124 127 126 125 127 133 134 114 144 213 101 101 
      100   102  76  81 103 135 150 208 215 211 215 212 200 198 189 190 182 172 162 153 143 134 134  99  99 100 100 100 100 
      200   101  79  88  89  92  93 144 136 126 128 123 141 141 142 138 130 124 133 126 121 115 109 100 100 100  99  99 100 
      400   104  88  94  99 101  96 183 235 262 235 218 194 179 187 169 160 138 144 130 122 108 100 100  99 100 100 100 100 
      800   101  91  95  97  98  88 141 136 194 186 165 149 141 142 125 112 105 100 100 101 100 100 100 101 101 101  99 101 
     1000   102  86  94  97  85  89 132 125 176 171 149 133 128 124 110  98  92 100 100 100 100 100 101 100 100 101 100  99 
     1200    99  98  95  97  92  95 110 107 159 153 133 117 115 113 100  88  84 100 100 100 100 100 100  99 101 100 100 100 
     1500   101  96  91 109  90 106 112 106 139 130 124 113 108 100 100 100 100 100 100 100 100 100 100 100 100 100  99 100 
     2000   100 103  99 110  96 117 111 103 143 134 129 117 111 100 100 100 100 100  99 101 101 100  99 100 101 100 100 100 
     3000   100 109 114 124 128 135 114 106 120 119 105  99  99 100 100 100 100 100 101 101 101 100 100 100 100 100 100 100 
     4000   102 109 116 127 135 144 139 120 108 111 100 100 100 100 100  99 100 100  98 100 100 100 101 100 100 101 100 100 
     6000   101 113 160 117 166 147 133 140 128 119 100  99 100 100 100 100 100 100 101 100 101 100 101 101 100 101 100 102 
    10000   100 106 157 116 147 129 124 116 100 100 100 100 100 100 100 101 101 100 101 100 101 101 100 100 100  99  99 100 
    20000   100 103 155 119 113 104 100 100 100 101 100 100  99 100 100 100 100 101 100 101 100 100 102 100 100 101  99 100 
    40000   101 101 154 105  99 100 100 100 100 100 100 100 100 100 100 100 100 101 100 100 100 100 100 102 101 100 101  99 
   100000   101 102 153 118 102  96  97  98  93  91  90 104 106 101 100 100 101 100 101 100 100 100 100 101 100 100 100 100 
   200000   100 101 154 118 106 101 124  98 101 103  95 108 111 100 100  99 101  99 100 100 100 101 100 100 100 100 100 100 
   400000   100 101 155 118 115 124 120 112 113 105 104 118 119  99 100 100 100 100 100 100 100 100 100 100 100 100 100 100 
  1000000   100 101 153 119 134 134 115 114 114 120 123 112 111  99 100 100 100 100 100  99 100 100 100 100  99 100 100 100 

Speedup for fmpz_poly_sqrlow

 bits \ n     1   2   3   4   5   6   7   8   9  10  11  13  15  18  21  25  30  36  43  51  61  73  87 104 124 148 177 212 
       20   114  86  92  82  86  86  87  83  89  88  94  92  95  97  95  97  98  97  97  98 100 100 101 101 101  99 100  99 
       60   115 113 132 134 130 130 140 139 129 140 132 134 136 135 132 131 130 135 133 133 130 132 134 134 132 194 236 256 
      100   113  61  84  96 113 119 238 239 241 260 247 248 250 246 256 252 242 233 224 216 209 197  99 100 100 100 100  99 
      200   108  72  92  94  95  97 212 201 201 194 191 179 175 222 213 208 208 204 197 185 176 151 100 100 100 100 100 100 
      400   105  83  93  92  96  89 291 409 458 406 385 350 311 313 285 253 213 215 190 171 151 100 100 100 100 100 100 100 
      800   103  88  98  97  98  95 186 155 322 290 270 238 203 204 176 152 126 139 100 100 100 100 100 100 100 100 100 100 
     1000   102  90  91  94  96  96 172 141 296 263 243 207 182 174 149 129 108 120 100 100 100 100 100 100 100 100 100 100 
     1200   101  90  94  95  96  97 148 121 250 223 208 180 161 154 135 116  97 107 101 100 100 100 101 100 100 100 101 100 
     1500   102  93  96  96  97  98 140 117 239 217 198 170 148 100 100 100 100 100 100 100 100 100 100  99 100 100 100 100 
     2000   101  94  95  97  98  97 134 111 230 201 187 162 142 100 100 100 100 100 100 100  99 100 100 100 100 100 100 100 
     3000   100  98 100  99 100 100 139 116 162 149 132 115 115 100 100 100 100 101  99 100  99 101 100 100 100 100 100 100 
     4000   100  99  99  99  99  99 130 110 132 123 109 105  91 100 100 100 100 101  99 100 100 100 100 100 101 100 100 100 
     6000   101 101 100 100 162 145 135 141 131 123 127 118 121 100 100 100 101 100 101  99 100 100 101  99 102 100 100  99 
    10000   100  99 100 100 119 131 127 123 131 132 119 126 121 126 120 106 108 107 103  91  86 100 101 100 101  99 100 100 
    20000   100 102  83  91 124 125 131 121 125 135 137 134 125 130 127 119 109 119 113 107  96 100 100 100 100 100  99 100 
    40000   100 104 100 111 147 155 147 126 151 148 131 141 121 147 143 123 124 156 130 124 130 125 105 100  99 101  99 101 
   100000   100 102 133 147 154 140 158 162 164 152 142 175 184 153 157 162 167 136 131 137 138 112 104 100 100 100 100 100 
   200000    99 104 131 146 178 179 210 189 179 162 178 184 195 168 171 170 178 158 143 142 142 116 106 101 101 100 100 100 
   400000   100 100 130 146 182 196 192 186 180 184 179 181 203 183 170 167 172 153 139 137 129 111 101 100 100 100  99 100 
  1000000   100 105 123 129 205 221 183 184 197 185 199 180 166 166 166 166 155 152 148 133 124 119 119 100 100 100 100 100 

Speedup for fmpz_poly_mulmid

Balanced mulmid window: 2n x n -> [n, 2n)

 bits \ n     1   2   3   4   5   6   7   8   9  10  11  13  15  18  21  25  30  36  43  51  61  73  87 104 124 148 177 212 
       20    86 103 101  99 119 111 117  96  99  98  99 100 101 101 100 100 100  99 101 100 100  99 100  99 100  99 100 100 
       60    92 126 123 129 134 134 130 131 136 133 132 130 131 130 127 129 130 130 129 124 129 128 126 128 129 166 171  99 
      100    86  94  95  88 102 115 246 234 219 224 230 226 208 212 202 193 184 169 166 154 143 136  98 100 100 100 100 100 
      200    82  93  97  98  97  99 201 185 182 180 169 159 146 147 146 145 136 135 149 121 123 108 100 100 101 100 100 100 
      400    89  94  97  99 100  92 236 377 429 385 357 322 276 274 234 212 187 198 179 156 142 178 100 100 100 101  99 101 
      800    96  95  96  98  98  96 154 128 272 243 226 195 162 194 172 147 128 133 129 134 119 133 100 101 101 100 100 101 
     1000    97  89  94  96  97  97 142 114 232 203 187 162 142 166 148 126 108 112 125 112 110 119 101 100 100 100 100 100 
     1200   100  92  95  97  97  98 132 108 213 190 174 151 130 156 140 120 104 105 121 120 115 136 100 100 100 100 101 100 
     1500    99  94  95  97  96  98 120  97 200 180 175 149 138 148 130 132 114 130 129 127 116 140 102 100 100 100 101 100 
     2000    99  95  97  98  99  99 114  92 211 188 181 154 145 153 147 128 130 136 140 133 133 152 100 100 100 100 100 100 
     3000   100  98  99  99  99 100 119  96 125 110  98 100 100 100 100 100 100  99 101 101 100 100  99 101 100 100 100 100 
     4000    99  98  98  99 100  99 114  93 114 112  98 100 100  99 100 101 101 101 101 101 101 100 100 100 100 101 101 101 
     6000   100  97 106 111 134 121 113 121 114 107  96 102 100  99 100  99 100 100 100 101 100 100 101 100 100 100 100 100 
    10000   100  97 109 121 119 125 103 110 123 118 115 106 111 108 100 100 100 100 100 100 100 100 101 100 101 100 102 100 
    20000   100 135 120 137 121 118 134 123 123 130 127 124 114 125 114 109 101 104  94 101 100 101  99 101 100 101  99 101 
    40000   100 132 127 156 143 148 147 119 136 135 143 132 125 132 124 119 113 120 104  99  95 100 100 100 100 101 101 100 
   100000   100 134 170 205 155 115 145 138 146 140 133 141 145 134 129 128 128 113 110 103  96 101 100 100 100 100 100 100 
   200000   100 135 139 152 157 156 163 144 142 141 145 146 152 133 144 142 134 122 113 104  92 100 101 100 100 101 100 100 
   400000   100 131 142 177 139 171 155 173 167 156 157 160 163 148 137 139 128 111 101  94  86 100 100 100 100 100 100 100 
  1000000   100 139 165 171 183 165 180 169 158 161 158 156 145 136 127 131 122 115 109 102  93 100  98 100 100 100 101  99 

Speedup for APRCL

Time to call fmpz_is_prime on nextprime(2^(bits+1)/3):

    bits    old        new       speedup
     256    0.00797    0.00733   1.087x
     512    0.0727     0.0678    1.072x
     768    0.383      0.349     1.097x
    1024    1.249      1.126     1.109x
    1280    2.613      2.283     1.145x
    1536    4.581      3.986     1.149x
    1792    8.435      7.437     1.134x
    2048    17.318     16.594    1.044x
    2304    22.036     20.727    1.063x
    2560    37.815     34.15     1.107x
    2816    50.005     45.159    1.107x
    3072    107.49     99.106    1.085x
    3328    96.905     89.124    1.087x
    3584    118.673    110.511   1.074x