Improve sine and cosine of power series

[fredrik-johansson]

Until now, we computed sine and cosine of power series of high order via the tangent, which in turn uses Newton iteration to invert the arctangent. It is faster to use a direct simultaneous Newton iteration for the sine and cosine, which can be derived from that for $\exp(i h)$ (without actually using complex arithmetic). Detailed PR content:

• Add _gr_poly_sin_cos_series_newton / gr_poly_sin_cos_series_newton (implemented using Claude, first by asking it to adapt _gr_poly_exp_series for sine and cosine computation, then making various improvements to the initial working version).

• Add the following functions (and their underscore versions), with automatic selection between basecase and Newton algorithms, plus test code:

  • gr_poly_sin_cos_pi_series
  • gr_poly_sin_pi_series
  • gr_poly_cos_pi_series
  • gr_poly_sin_cos_series
  • gr_poly_sin_series
  • gr_poly_cos_series

• Add their wrappers for power series and truncated power series in the gr_series module:

  • gr_series_mod_sin
  • gr_series_mod_cos
  • gr_series_mod_sin_pi
  • gr_series_mod_cos_pi
  • gr_series_mod_sin_cos
  • gr_series_mod_sin_cos_pi
  • gr_series_sin
  • gr_series_cos
  • gr_series_sin_pi
  • gr_series_cos_pi
  • gr_series_sin_cos
  • gr_series_sin_cos_pi

• Use the new GR algorithm in nmod_poly_sin_series, nmod_poly_cos_series, arb_poly_sin_series, arb_poly_cos_series, arb_poly_sin_pi_series, arb_poly_cos_pi_series, arb_poly_sin_cos_series, arb_poly_sin_cos_pi_series, acb_poly_sin_series, acb_poly_cos_series, acb_poly_sin_pi_series, acb_poly_cos_pi_series, acb_poly_sin_cos_series, acb_poly_sin_cos_pi_series, improving performance of these methods

• Major bug fix: _acb_poly_sin_cos_series always computed with 53 bit precision (for n > 2).

• Small optimization (reciprocal precomputation) of _gr_poly_sin_cos_series_basecase when used in finite characteristic

• Other minor cleanup and fixes

Still missing is an fmpq_poly version of the Newton algorithm.

Timings for nmod_poly_sin_series, 64-bit modulus:

      n          old          new  speedup
      4      2.95e-07    1.01e-07  (2.921)   
      8      6.65e-07    2.53e-07  (2.628)   
     16      1.51e-06    6.43e-07  (2.348)   
     32      4.31e-06    1.58e-06  (2.728)   
     64      1.14e-05    4.23e-06  (2.695)   
    128      3.11e-05    1.25e-05  (2.488)   
    256      8.63e-05    4.12e-05  (2.095)   
    512      0.000201    0.000114  (1.763)   
   1024      0.000443    0.000269  (1.647)   
   2048       0.00103    0.000619  (1.664)   
   4096       0.00226     0.00138  (1.638)   
   8192       0.00486     0.00297  (1.636)   
  16384        0.0104      0.0063  (1.651)   
  32768        0.0218      0.0132  (1.652)   
  65536        0.0452      0.0275  (1.644)   

Timings for arb_poly_sin_series, 1024-bit precision, input $\log(1+x)$:

      n           old         new  speedup     [x^(n-1)]  old         [x^(n-1)]  new
      4      7.63e-07    7.43e-07  (1.027)   0.167 +/- 4.1720e-309    0.167 +/- 4.1720e-309    
      8       2.9e-06     2.9e-06  (1.000)   -0.145 +/- 7.6239e-309    -0.145 +/- 7.6239e-309    
     16      1.62e-05    1.61e-05  (1.006)   -0.122 +/- 8.8582e-309    -0.122 +/- 8.8582e-309    
     32      8.54e-05    8.53e-05  (1.001)   -0.0510 +/- 9.2007e-309    -0.0510 +/- 9.2007e-309    
     64      0.000398    0.000399  (0.997)   -0.00821 +/- 9.3724e-309    -0.00821 +/- 9.3724e-309    
    128       0.00175     0.00161  (1.087)   0.00626 +/- 9.5871e-309    0.00626 +/- 1.0855e-306    
    256       0.00737     0.00442  (1.667)   0.00678 +/- 9.6380e-309    0.00678 +/- 1.0523e-304    
    512        0.0165      0.0103  (1.602)   0.00363 +/- 2.5300e-280    0.00363 +/- 1.4184e-302    
   1024        0.0361      0.0223  (1.619)   0.00109 +/- 4.2209e-263    0.00109 +/- 2.3708e-300    
   2048        0.0842      0.0519  (1.622)   -6.69e-5 +/- 2.2312e-234    -6.69e-5 +/- 5.0254e-298    
   4096         0.232       0.104  (2.231)   -0.000324 +/- 2.8829e-184    -0.000324 +/- 1.3357e-295    
   8192         0.608       0.218  (2.789)   -0.000233 +/- 4.3381e-93    -0.000233 +/- 4.5419e-293    
  16384         1.622       0.468  (3.466)   7.56e+7 +/- 1.5416e+79    -9.78e-5 +/- 1.8537e-290    

You may note that the new algorithm is a lot more numerically stable for large n, though this may vary with the input.

Timings for arb_poly_sin_cos_series with the same input:

      n           old         new  speedup     [x^(n-1)]  old         [x^(n-1)]  new
      4      7.16e-07    7.49e-07  (0.956)   0.167 +/- 4.1720e-309    0.167 +/- 4.1720e-309    
      8      2.77e-06    2.79e-06  (0.993)   -0.145 +/- 7.6239e-309    -0.145 +/- 7.6239e-309    
     16      1.55e-05    1.54e-05  (1.006)   -0.122 +/- 8.8582e-309    -0.122 +/- 8.8582e-309    
     32      8.19e-05     8.2e-05  (0.999)   -0.0510 +/- 9.2007e-309    -0.0510 +/- 9.2007e-309    
     64      0.000388     0.00039  (0.995)   -0.00821 +/- 9.3724e-309    -0.00821 +/- 9.3724e-309    
    128       0.00171     0.00171  (1.000)   0.00626 +/- 9.5871e-309    0.00626 +/- 2.8561e-306    
    256        0.0072     0.00474  (1.519)   0.00678 +/- 9.6380e-309    0.00678 +/- 2.6019e-304    
    512        0.0162      0.0104  (1.558)   0.00363 +/- 2.5300e-280    0.00363 +/- 3.3135e-302    
   1024        0.0354      0.0232  (1.526)   0.00109 +/- 4.2209e-263    0.00109 +/- 5.7644e-300    
   2048        0.0806      0.0512  (1.574)   -6.69e-5 +/- 2.2312e-234    -6.69e-5 +/- 1.2752e-297    
   4096         0.227       0.106  (2.142)   -0.000324 +/- 2.8829e-184    -0.000324 +/- 3.3368e-295    
   8192           0.6       0.222  (2.703)   -0.000233 +/- 4.3381e-93    -0.000233 +/- 1.0825e-292    
  16384         1.556       0.489  (3.182)   7.56e+7 +/- 1.5416e+79    -9.78e-5 +/- 4.3923e-290