raw
smg_comms_c_wrappers    1 /* mpi-div.c  -  MPI functions
smg_comms_c_wrappers    2  * Modified by No Such Labs. (C) 2015. See README.
smg_comms_c_wrappers    3  *
smg_comms_c_wrappers    4  * This file was originally part of Gnu Privacy Guard (GPG), ver. 1.4.10,
smg_comms_c_wrappers    5  * SHA256(gnupg-1.4.10.tar.gz):
smg_comms_c_wrappers    6  *        0bfd74660a2f6cedcf7d8256db4a63c996ffebbcdc2cf54397bfb72878c5a85a
smg_comms_c_wrappers    7  * (C) 1994-2005 Free Software Foundation, Inc.
smg_comms_c_wrappers    8  *
smg_comms_c_wrappers    9  * This program is free software: you can redistribute it and/or modify
smg_comms_c_wrappers   10  * it under the terms of the GNU General Public License as published by
smg_comms_c_wrappers   11  * the Free Software Foundation, either version 3 of the License, or
smg_comms_c_wrappers   12  * (at your option) any later version.
smg_comms_c_wrappers   13  *
smg_comms_c_wrappers   14  * This program is distributed in the hope that it will be useful,
smg_comms_c_wrappers   15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
smg_comms_c_wrappers   16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
smg_comms_c_wrappers   17  * GNU General Public License for more details.
smg_comms_c_wrappers   18  *
smg_comms_c_wrappers   19  * You should have received a copy of the GNU General Public License
smg_comms_c_wrappers   20  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
smg_comms_c_wrappers   21  */
smg_comms_c_wrappers   22 
smg_comms_c_wrappers   23 #include <stdio.h>
smg_comms_c_wrappers   24 #include <stdlib.h>
smg_comms_c_wrappers   25 
smg_comms_c_wrappers   26 #include "knobs.h"
smg_comms_c_wrappers   27 #include "mpi-internal.h"
smg_comms_c_wrappers   28 #include "longlong.h"
smg_comms_c_wrappers   29 
smg_comms_c_wrappers   30 
smg_comms_c_wrappers   31 
smg_comms_c_wrappers   32 void
smg_comms_c_wrappers   33 mpi_fdiv_r( MPI rem, MPI dividend, MPI divisor )
smg_comms_c_wrappers   34 {
smg_comms_c_wrappers   35     int divisor_sign = divisor->sign;
smg_comms_c_wrappers   36     MPI temp_divisor = NULL;
smg_comms_c_wrappers   37 
smg_comms_c_wrappers   38     /* We need the original value of the divisor after the remainder has been
smg_comms_c_wrappers   39      * preliminary calculated.	We have to copy it to temporary space if it's
smg_comms_c_wrappers   40      * the same variable as REM.  */
smg_comms_c_wrappers   41     if( rem == divisor ) {
smg_comms_c_wrappers   42 	temp_divisor = mpi_copy( divisor );
smg_comms_c_wrappers   43 	divisor = temp_divisor;
smg_comms_c_wrappers   44     }
smg_comms_c_wrappers   45 
smg_comms_c_wrappers   46     mpi_tdiv_r( rem, dividend, divisor );
smg_comms_c_wrappers   47 
smg_comms_c_wrappers   48     if( ((divisor_sign?1:0) ^ (dividend->sign?1:0)) && rem->nlimbs )
smg_comms_c_wrappers   49 	mpi_add( rem, rem, divisor);
smg_comms_c_wrappers   50 
smg_comms_c_wrappers   51     if( temp_divisor )
smg_comms_c_wrappers   52 	mpi_free(temp_divisor);
smg_comms_c_wrappers   53 }
smg_comms_c_wrappers   54 
smg_comms_c_wrappers   55 
smg_comms_c_wrappers   56 
smg_comms_c_wrappers   57 /****************
smg_comms_c_wrappers   58  * Division rounding the quotient towards -infinity.
smg_comms_c_wrappers   59  * The remainder gets the same sign as the denominator.
smg_comms_c_wrappers   60  * rem is optional
smg_comms_c_wrappers   61  */
smg_comms_c_wrappers   62 
smg_comms_c_wrappers   63 ulong
smg_comms_c_wrappers   64 mpi_fdiv_r_ui( MPI rem, MPI dividend, ulong divisor )
smg_comms_c_wrappers   65 {
smg_comms_c_wrappers   66     mpi_limb_t rlimb;
smg_comms_c_wrappers   67 
smg_comms_c_wrappers   68     rlimb = mpihelp_mod_1( dividend->d, dividend->nlimbs, divisor );
smg_comms_c_wrappers   69     if( rlimb && dividend->sign )
smg_comms_c_wrappers   70 	rlimb = divisor - rlimb;
smg_comms_c_wrappers   71 
smg_comms_c_wrappers   72     if( rem ) {
smg_comms_c_wrappers   73 	rem->d[0] = rlimb;
smg_comms_c_wrappers   74 	rem->nlimbs = rlimb? 1:0;
smg_comms_c_wrappers   75     }
smg_comms_c_wrappers   76     return rlimb;
smg_comms_c_wrappers   77 }
smg_comms_c_wrappers   78 
smg_comms_c_wrappers   79 
smg_comms_c_wrappers   80 void
smg_comms_c_wrappers   81 mpi_fdiv_q( MPI quot, MPI dividend, MPI divisor )
smg_comms_c_wrappers   82 {
smg_comms_c_wrappers   83     MPI tmp = mpi_alloc( mpi_get_nlimbs(quot) );
smg_comms_c_wrappers   84     mpi_fdiv_qr( quot, tmp, dividend, divisor);
smg_comms_c_wrappers   85     mpi_free(tmp);
smg_comms_c_wrappers   86 }
smg_comms_c_wrappers   87 
smg_comms_c_wrappers   88 void
smg_comms_c_wrappers   89 mpi_fdiv_qr( MPI quot, MPI rem, MPI dividend, MPI divisor )
smg_comms_c_wrappers   90 {
smg_comms_c_wrappers   91     int divisor_sign = divisor->sign;
smg_comms_c_wrappers   92     MPI temp_divisor = NULL;
smg_comms_c_wrappers   93 
smg_comms_c_wrappers   94     if( quot == divisor || rem == divisor ) {
smg_comms_c_wrappers   95 	temp_divisor = mpi_copy( divisor );
smg_comms_c_wrappers   96 	divisor = temp_divisor;
smg_comms_c_wrappers   97     }
smg_comms_c_wrappers   98 
smg_comms_c_wrappers   99     mpi_tdiv_qr( quot, rem, dividend, divisor );
smg_comms_c_wrappers  100 
smg_comms_c_wrappers  101     if( (divisor_sign ^ dividend->sign) && rem->nlimbs ) {
smg_comms_c_wrappers  102 	mpi_sub_ui( quot, quot, 1 );
smg_comms_c_wrappers  103 	mpi_add( rem, rem, divisor);
smg_comms_c_wrappers  104     }
smg_comms_c_wrappers  105 
smg_comms_c_wrappers  106     if( temp_divisor )
smg_comms_c_wrappers  107 	mpi_free(temp_divisor);
smg_comms_c_wrappers  108 }
smg_comms_c_wrappers  109 
smg_comms_c_wrappers  110 
smg_comms_c_wrappers  111 /* If den == quot, den needs temporary storage.
smg_comms_c_wrappers  112  * If den == rem, den needs temporary storage.
smg_comms_c_wrappers  113  * If num == quot, num needs temporary storage.
smg_comms_c_wrappers  114  * If den has temporary storage, it can be normalized while being copied,
smg_comms_c_wrappers  115  *   i.e no extra storage should be allocated.
smg_comms_c_wrappers  116  */
smg_comms_c_wrappers  117 
smg_comms_c_wrappers  118 void
smg_comms_c_wrappers  119 mpi_tdiv_r( MPI rem, MPI num, MPI den)
smg_comms_c_wrappers  120 {
smg_comms_c_wrappers  121     mpi_tdiv_qr(NULL, rem, num, den );
smg_comms_c_wrappers  122 }
smg_comms_c_wrappers  123 
smg_comms_c_wrappers  124 void
smg_comms_c_wrappers  125 mpi_tdiv_qr( MPI quot, MPI rem, MPI num, MPI den)
smg_comms_c_wrappers  126 {
smg_comms_c_wrappers  127     mpi_ptr_t np, dp;
smg_comms_c_wrappers  128     mpi_ptr_t qp, rp;
smg_comms_c_wrappers  129     mpi_size_t nsize = num->nlimbs;
smg_comms_c_wrappers  130     mpi_size_t dsize = den->nlimbs;
smg_comms_c_wrappers  131     mpi_size_t qsize, rsize;
smg_comms_c_wrappers  132     mpi_size_t sign_remainder = num->sign;
smg_comms_c_wrappers  133     mpi_size_t sign_quotient = num->sign ^ den->sign;
smg_comms_c_wrappers  134     unsigned normalization_steps;
smg_comms_c_wrappers  135     mpi_limb_t q_limb;
smg_comms_c_wrappers  136     mpi_ptr_t marker[5];
smg_comms_c_wrappers  137     int markidx=0;
smg_comms_c_wrappers  138 
smg_comms_c_wrappers  139     /* Ensure space is enough for quotient and remainder.
smg_comms_c_wrappers  140      * We need space for an extra limb in the remainder, because it's
smg_comms_c_wrappers  141      * up-shifted (normalized) below.  */
smg_comms_c_wrappers  142     rsize = nsize + 1;
smg_comms_c_wrappers  143     mpi_resize( rem, rsize);
smg_comms_c_wrappers  144 
smg_comms_c_wrappers  145     qsize = rsize - dsize;	  /* qsize cannot be bigger than this.	*/
smg_comms_c_wrappers  146     if( qsize <= 0 ) {
smg_comms_c_wrappers  147 	if( num != rem ) {
smg_comms_c_wrappers  148 	    rem->nlimbs = num->nlimbs;
smg_comms_c_wrappers  149 	    rem->sign = num->sign;
smg_comms_c_wrappers  150 	    MPN_COPY(rem->d, num->d, nsize);
smg_comms_c_wrappers  151 	}
smg_comms_c_wrappers  152 	if( quot ) {
smg_comms_c_wrappers  153 	    /* This needs to follow the assignment to rem, in case the
smg_comms_c_wrappers  154 	     * numerator and quotient are the same.  */
smg_comms_c_wrappers  155 	    quot->nlimbs = 0;
smg_comms_c_wrappers  156 	    quot->sign = 0;
smg_comms_c_wrappers  157 	}
smg_comms_c_wrappers  158 	return;
smg_comms_c_wrappers  159     }
smg_comms_c_wrappers  160 
smg_comms_c_wrappers  161     if( quot )
smg_comms_c_wrappers  162 	mpi_resize( quot, qsize);
smg_comms_c_wrappers  163 
smg_comms_c_wrappers  164     /* Read pointers here, when reallocation is finished.  */
smg_comms_c_wrappers  165     np = num->d;
smg_comms_c_wrappers  166     dp = den->d;
smg_comms_c_wrappers  167     rp = rem->d;
smg_comms_c_wrappers  168 
smg_comms_c_wrappers  169     /* Optimize division by a single-limb divisor.  */
smg_comms_c_wrappers  170     if( dsize == 1 ) {
smg_comms_c_wrappers  171 	mpi_limb_t rlimb;
smg_comms_c_wrappers  172 	if( quot ) {
smg_comms_c_wrappers  173 	    qp = quot->d;
smg_comms_c_wrappers  174 	    rlimb = mpihelp_divmod_1( qp, np, nsize, dp[0] );
smg_comms_c_wrappers  175 	    qsize -= qp[qsize - 1] == 0;
smg_comms_c_wrappers  176 	    quot->nlimbs = qsize;
smg_comms_c_wrappers  177 	    quot->sign = sign_quotient;
smg_comms_c_wrappers  178 	}
smg_comms_c_wrappers  179 	else
smg_comms_c_wrappers  180 	    rlimb = mpihelp_mod_1( np, nsize, dp[0] );
smg_comms_c_wrappers  181 	rp[0] = rlimb;
smg_comms_c_wrappers  182 	rsize = rlimb != 0?1:0;
smg_comms_c_wrappers  183 	rem->nlimbs = rsize;
smg_comms_c_wrappers  184 	rem->sign = sign_remainder;
smg_comms_c_wrappers  185 	return;
smg_comms_c_wrappers  186     }
smg_comms_c_wrappers  187 
smg_comms_c_wrappers  188 
smg_comms_c_wrappers  189     if( quot ) {
smg_comms_c_wrappers  190 	qp = quot->d;
smg_comms_c_wrappers  191 	/* Make sure QP and NP point to different objects.  Otherwise the
smg_comms_c_wrappers  192 	 * numerator would be gradually overwritten by the quotient limbs.  */
smg_comms_c_wrappers  193 	if(qp == np) { /* Copy NP object to temporary space.  */
smg_comms_c_wrappers  194 	    np = marker[markidx++] = mpi_alloc_limb_space(nsize,
smg_comms_c_wrappers  195 							  mpi_is_secure(quot));
smg_comms_c_wrappers  196 	    MPN_COPY(np, qp, nsize);
smg_comms_c_wrappers  197 	}
smg_comms_c_wrappers  198     }
smg_comms_c_wrappers  199     else /* Put quotient at top of remainder. */
smg_comms_c_wrappers  200 	qp = rp + dsize;
smg_comms_c_wrappers  201 
smg_comms_c_wrappers  202     count_leading_zeros( normalization_steps, dp[dsize - 1] );
smg_comms_c_wrappers  203 
smg_comms_c_wrappers  204     /* Normalize the denominator, i.e. make its most significant bit set by
smg_comms_c_wrappers  205      * shifting it NORMALIZATION_STEPS bits to the left.  Also shift the
smg_comms_c_wrappers  206      * numerator the same number of steps (to keep the quotient the same!).
smg_comms_c_wrappers  207      */
smg_comms_c_wrappers  208     if( normalization_steps ) {
smg_comms_c_wrappers  209 	mpi_ptr_t tp;
smg_comms_c_wrappers  210 	mpi_limb_t nlimb;
smg_comms_c_wrappers  211 
smg_comms_c_wrappers  212 	/* Shift up the denominator setting the most significant bit of
smg_comms_c_wrappers  213 	 * the most significant word.  Use temporary storage not to clobber
smg_comms_c_wrappers  214 	 * the original contents of the denominator.  */
smg_comms_c_wrappers  215 	tp = marker[markidx++] = mpi_alloc_limb_space(dsize,mpi_is_secure(den));
smg_comms_c_wrappers  216 	mpihelp_lshift( tp, dp, dsize, normalization_steps );
smg_comms_c_wrappers  217 	dp = tp;
smg_comms_c_wrappers  218 
smg_comms_c_wrappers  219 	/* Shift up the numerator, possibly introducing a new most
smg_comms_c_wrappers  220 	 * significant word.  Move the shifted numerator in the remainder
smg_comms_c_wrappers  221 	 * meanwhile.  */
smg_comms_c_wrappers  222 	nlimb = mpihelp_lshift(rp, np, nsize, normalization_steps);
smg_comms_c_wrappers  223 	if( nlimb ) {
smg_comms_c_wrappers  224 	    rp[nsize] = nlimb;
smg_comms_c_wrappers  225 	    rsize = nsize + 1;
smg_comms_c_wrappers  226 	}
smg_comms_c_wrappers  227 	else
smg_comms_c_wrappers  228 	    rsize = nsize;
smg_comms_c_wrappers  229     }
smg_comms_c_wrappers  230     else {
smg_comms_c_wrappers  231 	/* The denominator is already normalized, as required.	Copy it to
smg_comms_c_wrappers  232 	 * temporary space if it overlaps with the quotient or remainder.  */
smg_comms_c_wrappers  233 	if( dp == rp || (quot && (dp == qp))) {
smg_comms_c_wrappers  234 	    mpi_ptr_t tp;
smg_comms_c_wrappers  235 
smg_comms_c_wrappers  236 	    tp = marker[markidx++] = mpi_alloc_limb_space(dsize, mpi_is_secure(den));
smg_comms_c_wrappers  237 	    MPN_COPY( tp, dp, dsize );
smg_comms_c_wrappers  238 	    dp = tp;
smg_comms_c_wrappers  239 	}
smg_comms_c_wrappers  240 
smg_comms_c_wrappers  241 	/* Move the numerator to the remainder.  */
smg_comms_c_wrappers  242 	if( rp != np )
smg_comms_c_wrappers  243 	    MPN_COPY(rp, np, nsize);
smg_comms_c_wrappers  244 
smg_comms_c_wrappers  245 	rsize = nsize;
smg_comms_c_wrappers  246     }
smg_comms_c_wrappers  247 
smg_comms_c_wrappers  248     q_limb = mpihelp_divrem( qp, 0, rp, rsize, dp, dsize );
smg_comms_c_wrappers  249 
smg_comms_c_wrappers  250     if( quot ) {
smg_comms_c_wrappers  251 	qsize = rsize - dsize;
smg_comms_c_wrappers  252 	if(q_limb) {
smg_comms_c_wrappers  253 	    qp[qsize] = q_limb;
smg_comms_c_wrappers  254 	    qsize += 1;
smg_comms_c_wrappers  255 	}
smg_comms_c_wrappers  256 
smg_comms_c_wrappers  257 	quot->nlimbs = qsize;
smg_comms_c_wrappers  258 	quot->sign = sign_quotient;
smg_comms_c_wrappers  259     }
smg_comms_c_wrappers  260 
smg_comms_c_wrappers  261     rsize = dsize;
smg_comms_c_wrappers  262     MPN_NORMALIZE (rp, rsize);
smg_comms_c_wrappers  263 
smg_comms_c_wrappers  264     if( normalization_steps && rsize ) {
smg_comms_c_wrappers  265 	mpihelp_rshift(rp, rp, rsize, normalization_steps);
smg_comms_c_wrappers  266 	rsize -= rp[rsize - 1] == 0?1:0;
smg_comms_c_wrappers  267     }
smg_comms_c_wrappers  268 
smg_comms_c_wrappers  269     rem->nlimbs = rsize;
smg_comms_c_wrappers  270     rem->sign	= sign_remainder;
smg_comms_c_wrappers  271     while( markidx )
smg_comms_c_wrappers  272 	mpi_free_limb_space(marker[--markidx]);
smg_comms_c_wrappers  273 }
smg_comms_c_wrappers  274 
smg_comms_c_wrappers  275 void
smg_comms_c_wrappers  276 mpi_tdiv_q_2exp( MPI w, MPI u, unsigned count )
smg_comms_c_wrappers  277 {
smg_comms_c_wrappers  278     mpi_size_t usize, wsize;
smg_comms_c_wrappers  279     mpi_size_t limb_cnt;
smg_comms_c_wrappers  280 
smg_comms_c_wrappers  281     usize = u->nlimbs;
smg_comms_c_wrappers  282     limb_cnt = count / BITS_PER_MPI_LIMB;
smg_comms_c_wrappers  283     wsize = usize - limb_cnt;
smg_comms_c_wrappers  284     if( limb_cnt >= usize )
smg_comms_c_wrappers  285 	w->nlimbs = 0;
smg_comms_c_wrappers  286     else {
smg_comms_c_wrappers  287 	mpi_ptr_t wp;
smg_comms_c_wrappers  288 	mpi_ptr_t up;
smg_comms_c_wrappers  289 
smg_comms_c_wrappers  290 	RESIZE_IF_NEEDED( w, wsize );
smg_comms_c_wrappers  291 	wp = w->d;
smg_comms_c_wrappers  292 	up = u->d;
smg_comms_c_wrappers  293 
smg_comms_c_wrappers  294 	count %= BITS_PER_MPI_LIMB;
smg_comms_c_wrappers  295 	if( count ) {
smg_comms_c_wrappers  296 	    mpihelp_rshift( wp, up + limb_cnt, wsize, count );
smg_comms_c_wrappers  297 	    wsize -= !wp[wsize - 1];
smg_comms_c_wrappers  298 	}
smg_comms_c_wrappers  299 	else {
smg_comms_c_wrappers  300 	    MPN_COPY_INCR( wp, up + limb_cnt, wsize);
smg_comms_c_wrappers  301 	}
smg_comms_c_wrappers  302 
smg_comms_c_wrappers  303 	w->nlimbs = wsize;
smg_comms_c_wrappers  304     }
smg_comms_c_wrappers  305 }
smg_comms_c_wrappers  306 
smg_comms_c_wrappers  307 /****************
smg_comms_c_wrappers  308  * Check whether dividend is divisible by divisor
smg_comms_c_wrappers  309  * (note: divisor must fit into a limb)
smg_comms_c_wrappers  310  */
smg_comms_c_wrappers  311 int
smg_comms_c_wrappers  312 mpi_divisible_ui(MPI dividend, ulong divisor )
smg_comms_c_wrappers  313 {
smg_comms_c_wrappers  314     return !mpihelp_mod_1( dividend->d, dividend->nlimbs, divisor );
smg_comms_c_wrappers  315 }
smg_comms_c_wrappers  316