performance improvement to tdiv_qr.c

[Joe Weening]

mpn/generic/tdiv_qr.c contains the following comment:

    /* In theory, we could fall out to the cute code below
       since we now have exactly the situation that code
       is designed to handle.  We botch this badly and call
       the basic mpn_sb_divrem_mn!  */

followed by a call to mpn_sb_divrem_mn (or mpn_divrem_2), when for
large operands it should be using a faster method, as the comment
says.  This has a major impact on performance for large operands.

The following change fixes this problem.  Instead of "falling into the
cute code below", which does not work given the structure of the code
at that point, it performs a recursive call to mpn_tdiv_qr, so that
the appropriate method is chosen for any size operands.

The recursive call overwrites the high-order limb of its quotient,
whose value in qp[nn-dn] needs to be preserved, so it is saved and
restored around this call.  The other tricky thing is proving there is
enough space in the buffer at n2p for the remainder computed by the
recursive call.  I've convinced myself this is true, and can write it
up in detail if you'd like.

Joe Weening
IDA Center for Communications Research


*** tdiv_qr.c.orig	2002-11-11 15:34:43.000000000 -0800
--- tdiv_qr.c	2004-10-07 08:35:50.615063000 -0700
***************
*** 161,172 ****
  		  {
  		    n2p -= nn - dn;
! 		    /* In theory, we could fall out to the cute code below
! 		       since we now have exactly the situation that code
! 		       is designed to handle.  We botch this badly and call
! 		       the basic mpn_sb_divrem_mn!  */
! 		    if (dn == 2)
! 		      mpn_divrem_2 (qp, 0L, n2p, nn, d2p);
! 		    else
! 		      mpn_sb_divrem_mn (qp, n2p, nn, d2p, dn);
  		  }
  	      }
--- 161,170 ----
  		  {
  		    n2p -= nn - dn;
! 		    {
! 		      /* Preserve quotient limb overwritten by recursive call */
! 		      mp_limb_t qs = qp[nn-dn];
! 		      mpn_tdiv_qr (qp, n2p, 0L, n2p, nn, d2p, dn);
! 		      qp[nn-dn] = ds;
! 		    }
  		  }
  	      }