// file kernel/x/c/div.c: division of extensible integers
/*-----------------------------------------------------------------------+
 |  Copyright 2005-2006, Michel Quercia (michel.quercia@prepas.org)      |
 |                                                                       |
 |  This file is part of Numerix. Numerix is free software; you can      |
 |  redistribute it and/or modify it under the terms of the GNU Lesser   |
 |  General Public License as published by the Free Software Foundation; |
 |  either version 2.1 of the License, or (at your option) any later     |
 |  version.                                                             |
 |                                                                       |
 |  The Numerix Library is distributed in the hope that it will be       |
 |  useful, but WITHOUT ANY WARRANTY; without even the implied warranty  |
 |  of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU  |
 |  Lesser General Public License for more details.                      |
 |                                                                       |
 |  You should have received a copy of the GNU Lesser General Public     |
 |  License along with the GNU MP Library; see the file COPYING. If not, |
 |  write to the Free Software Foundation, Inc., 59 Temple Place -       |
 |  Suite 330, Boston, MA 02111-1307, USA.                               |
 +-----------------------------------------------------------------------+
 |                                                                       |
 |                                Division                               |
 |                                                                       |
 +-----------------------------------------------------------------------*/

                         /* +---------------------+
                            |  Division gnrale  |
                            +---------------------+ */

/*
  entre :
  a,b = entiers extensibles
  _c  = NULL ou adresse d'un entier extensible
  _d  = NULL ou adresse d'un entier extensible
  mode = long

  contraintes :
  b != 0
  si _c et _d ne sont pas gaux  NULL alors ils sont distincts

  sortie :
  si mode & 3 = 0 : c <- floor(a/b),     d <- a - b*c (division tronque    )
  si mode & 3 = 1 : c <- floor(a/b+1/2), d <- a - b*c (division centre up  )
  si mode & 3 = 2 : c <- ceil(a/b),      d <- a - b*c (division majore     )
  si mode & 3 = 3 : c <- ceil(a/b-1/2),  d <- a - b*c (division centre down)
  si _c != NULL : *_c <- c
  si _d != NULL : *_d <- d

  valeur de retour  : c_api   ml_api
  si mode & 12 =  0 :  NULL     unit
  si mode & 12 =  4 :     c        c
  si mode & 12 =  8 :     d        d
  si mode & 12 = 12 :  NULL    (c,d)

  erreur :
  ZERO_DIVISOR si b = 0
  MULTIPLE_RESULT si _c == _d != NULL
*/
#if defined(c_api)
xint
#elif defined(caml_api) || defined(ocaml_api)
value
#endif
xx(private_quomod)(xint *_c, xint *_d, xint a, xint b, long mode) {
    long la = xx_lg(a), sa = xx_sgn(a);
    long lb = xx_lg(b), sb = xx_sgn(b);
    long lc,ld,sc,sd;
    int  want_c, want_d, inc_c;
    xx_push_roots_42(a,b,_c,_d,c,d);
#ifdef caml_api
#define  a __lr.a
#define  b __lr.b
#define  c __lr.c
#define  d __lr.d
#define _c __lr._c
#define _d __lr._d
#endif

    /* rsultats  calculer */
    want_c = ((_c != xx_null) || (mode & 4));
    want_d = ((_d != xx_null) || (mode & 8));
#ifdef c_api
    if (!(want_c | want_d)) return NULL;
#endif

    /* contrle */
    if (lb ==  0) xx(failwith)(ZERO_DIVISOR);
    if ((_c == _d) && (_d != xx_null)) xx(failwith)(MULTIPLE_RESULT);

    /* Algorithme pour la division avec reste :
       1. on divise les valeurs absolues : |a| = |b|*q + r

       2. correction pour la division tronque
       a >= 0, b > 0 (a =  bq+r)              c <-  q,   d <- r
       a <  0, b > 0 (a = -bq-r) r > 0        c <- -q-1, d <- |b|-r
                                 r = 0        c <- -q,   d <- r
       a >= 0, b < 0 (a = -bq+r) r > 0        c <- -q-1, d <- -(|b|-r)
                                 r = 0        c <- -q,   d <- -r
       a < 0,  b < 0 (a =  bq-r)              c <-  q,   d <- -r

       3. correction pour la division centre up
       a >= 0, b > 0 (a =  bq+r) 2r <  |b|    c <-  q,   d <- r
                                 2r >= |b|    c <-  q+1, d <- -(|b|-r)
       a <  0, b > 0 (a = -bq-r) 2r <= |b|    c <- -q,   d <- -r
                                 2r >  |b|    c <- -q-1, d <- |b|-r
       a >= 0, b < 0 (a = -bq+r) 2r <= |b|    c <- -q,   d <- r
                                 2r >  |b[    c <- -q-1, d <- -(|b|-r)
       a < 0,  b < 0 (a =  bq-r) 2r <  |b|    c <-  q,   d <- -r
                                 2r >= |b|    c <-  q+1, d <- |b|-r

       3. correction pour la division majore
       a >= 0, b > 0 (a =  bq+r) r > 0        c <-  q+1, d <- -(|b|-r)
                                 r = 0        c <-  q,   d <- -r
       a <  0, b > 0 (a = -bq-r)              c <- -q,   d <- -r
       a >= 0, b < 0 (a = -bq+r)              c <- -q,   d <- r
       a < 0,  b < 0 (a =  bq-r) r > 0        c <-  q+1, d <- |b|-r
                                 r = 0        c <-  q,   d <- r

       3. correction pour la division centre down
       a >= 0, b > 0 (a =  bq+r) 2r <= |b|    c <-  q,   d <- r
                                 2r >  |b|    c <-  q+1, d <- -(|b|-r)
       a <  0, b > 0 (a = -bq-r) 2r <  |b|    c <- -q,   d <- -r
                                 2r >= |b|    c <- -q-1, d <- |b|-r
       a >= 0, b < 0 (a = -bq+r) 2r <  |b|    c <- -q,   d <- r
                                 2r >= |b[    c <- -q-1, d <- -(|b|-r)
       a < 0,  b < 0 (a =  bq-r) 2r <= |b|    c <-  q,   d <- -r
                                 2r >  |b|    c <-  q+1, d <- |b|-r

       4. moralit :
       le signe de c est sa ^ sb
       on incrmente q et on complmente r lorsque :
       mode = 0, sa != sb, r > 0
       mode = 1, 2r > |b|
       mode = 1, sa = sb, 2r = |b|
       mode = 2, sa = sb, r > 0
       mode = 3, 2r > |b|
       mode = 3, sa != sb, 2r = |b|
       le signe de d est
       celui de  b en mode 0
       celui de  a en mode 1 si q n'est pas incrment
       celui de -a en mode 1 si q est incrment
       celui de -b en mode 2
       celui de  a en mode 3 si q n'est pas incrment
       celui de -a en mode 3 si q est incrment

       l'incrmentation de q et la fixation des signes seront effectues 
       la fin de la fonction. La complmentation de r est effectue juste
       aprs chaque division dans N, pendant qu'on dispose encore de |b|.

       Algorithme pour la division sans reste :
       1. on utilise l'algorithme prcdent lorsque :
       |b|   est petit (lb <= moddiv_lim)
       |a/b| est petit (la - lb + 3 <= div_small_c_lim)

       2. sinon, on calcule q = approx(x/|b|) avec :
       sa =  sb, mode = 0 -> x = |a|*BASE
       sa =  sb, mode = 1 -> x = (|a| + |b|/2)*BASE
       sa =  sb, mode = 2 -> x = (|a| + |b|)*BASE   - 1
       sa =  sb, mode = 3 -> x = (|a| + |b|/2)*BASE - 1
       sa != sb, mode = 0 -> x = (|a| + |b|)*BASE   - 1
       sa != sb, mode = 1 -> x = (|a| + |b|/2)*BASE - 1
       sa != sb, mode = 2 -> x = |a|*BASE
       sa != sb, mode = 3 -> x = (|a| + |b|/2)*BASE

       3. c <- floor(q/base), sc <- sa^sb

    */

    /* ------------------------------ division  un chiffre */
    if (lb <= chiffres_per_long) {
        unsigned long bb,dd;

        lc = la; ld = chiffres_per_long;

        /* division dans N */
#if chiffres_per_long == 1
        bb = b->val[0];
#else
        bb = (lb > 1) ? (unsigned long)b->val[0] + ((unsigned long)b->val[1] << HW)
                      : (unsigned long)b->val[0];
#endif
        if (want_c) {
             c  = xx(enlarge)(_c,lc+1);
             dd = xn(div_1)(a->val,la,bb,c->val);
        }
        else dd = xn(mod_1)(a->val,la,bb);

        /* correction */
        switch(mode & 3) {
            case 0: inc_c = ((sa != sb) && (dd));                      break;
            case 1: inc_c = ((dd>bb-dd) || ((dd==bb-dd) && (sa==sb))); break;
            case 2: inc_c = ((sa == sb) && (dd));                      break;
            default:inc_c = ((dd>bb-dd) || ((dd==bb-dd) && (sa!=sb))); break;
        }
        if (want_d) {
            d = xx(enlarge)(_d,ld);
#if chiffres_per_long == 1
            d->val[0] = (inc_c) ? bb-dd : dd;
#else
            if (inc_c) dd = bb-dd;
            d->val[0] = dd;
            d->val[1] = dd >> HW;
#endif
        }
    }

    /* ------------------------------ |dividende| < |diviseur| */
    else if (xn(cmp)(a->val,la,b->val,lb) < 0) {

        /* correction */
        switch(mode & 3) {
            case 0: inc_c = ((sa != sb) && (la));                             break;
            case 1: inc_c = ((sa == sb) + xn(cmp2)(a->val,la,b->val,lb) > 0); break;
            case 2: inc_c = ((sa == sb) && (la));                             break;
            default:inc_c = ((sa != sb) + xn(cmp2)(a->val,la,b->val,lb) > 0); break;
        }

        /* copie le reste */
        lc = 0;
        ld = (inc_c) ? lb : la;
        if (want_c) c = xx(enlarge)(_c,1);
        if (want_d) {
            d = xx(enlarge)(_d,ld);
            if (inc_c) {xn(sub)(b->val,lb,a->val,la,d->val);}
            else       {if (a != d) xn(move)(a->val,la,d->val);}
        }
    }

    /* ------------------------------ division avec reste  plusieurs chiffres */
    else if ((want_d) || (lb <= moddiv_lim) || (la - lb + 3 <= div_small_c_lim)) {
        chiffre *aa,*bb,*cc,*x,*y,r;
        long n,lx;

        /* taille des rsultats */
        lc = la - lb + 1; if (want_c) c = xx(enlarge)(_c,lc+1);
        ld = lb;          if (want_d) d = xx(enlarge)(_d,ld);

        /* dcalage  appliquer pour avoir msb(b) = 1 */
        for (r=b->val[lb-1], n=0; (r & (BASE_2)) == 0; r <<=1, n++);

        /* copie des oprandes :
           il faut copier b s'il va tre cras ou si n > 0
           il faut copier a si d n'est pas assez grand pour recevoir a*2^n
           augment d'un bit nul
        */
        lx = (want_c) ? 0 : lc;
        if ((n) || (b==c) || (b==d))             lx += lb;
        if ((!want_d) || (xx_capacity(d) <= la)) lx += la+1;
        x = y = (lx) ? xn(alloc)(lx) : NULL;

        if (!want_c) {
            cc = y; y += lc;
        } else  cc = c->val;

        if ((n) || (b==c) || (b==d)) {
            bb = y; y += lb;
            xn(shift_up)(b->val,lb,bb,n);
        } else bb = b->val;

        aa = ((!want_d) || (xx_capacity(d) <= la)) ? y : d->val;
        aa[la] = ((a->val != aa) || (n)) ? xn(shift_up)(a->val,la,aa,n): 0;

        /* Division dans N :
           on bascule directement vers div_n2 si b ou c est assez petit,
           sinon on utilise karpdiv qui renverra vers les autres algorithmes
           le cas chant.
        */
        ((lb <= burnidiv_lim) || (lc <= div_small_c_lim)) ?
            xn(div_n2) (aa,lc,bb,lb,cc) :
            xn(karpdiv)(aa,lc,bb,lb,cc,1);

        /* correction */
        switch(mode & 3) {
            case 0: inc_c = ((sa != sb) && (xn(cmp)(aa,lb,aa,0)));    break;
            case 1: inc_c = ((sa == sb) + xn(cmp2)(aa,lb,bb,lb) > 0); break;
            case 2: inc_c = ((sa == sb) && (xn(cmp)(aa,lb,aa,0)));    break;
            default:inc_c = ((sa != sb) + xn(cmp2)(aa,lb,bb,lb) > 0); break;
        }

        /* recopie le reste dans d avec dcalage s'il y a lieu */
        if (want_d) {
            if (inc_c) {
                xn(sub)(bb,lb,aa,lb,d->val);
                if (n) xn(shift_down)(d->val,lb,d->val,n);
            }
            else if ((d->val != aa) || (n)) xn(shift_down)(aa,lb,d->val,n);
        }

        /* libre la mmoire temporaire */
        xn(free)(x); /* on peut avoir x = NULL ici */
    }

    /* ------------------------------ division sans reste  plusieurs chiffres */
    else {
        chiffre *aa,*bb,r;
        long n;

        /* taille des rsultats */
        lc = la - lb + 2;
        ld = lb;
        c = xx(enlarge)(_c,lc);

        /* dcalage  appliquer pour avoir msb(b) = 1 */
        for (r=b->val[lb-1], n=0; (r & (BASE_2)) == 0; r <<=1, n++);

        /* copie des oprandes
           il faut copier b s'il va tre cras ou si n > 0
           il faut copier a puisqu'on ne calcule pas d
        */
        if ((n) || (b==c)) {aa=xn(alloc)(la+lb+3); bb = aa+la+3;}
        else               {aa=xn(alloc)(la+3);    bb = b->val;   }

        aa[0] = 0; aa[la+1] = xn(shift_up)(a->val,la,aa+1,n); aa[la+2] = 0;
        if (b->val != bb)     xn(shift_up)(b->val,lb,bb,n);

        /* arrondit le numrateur en fonction des signes de a et b et du mode */
        switch(mode & 3) {

            case 0: if (sa != sb) {
                      xn(inc)(aa+1,la+2,bb,lb);
                      xn(dec1)(aa,la+3);
                    }
                    break;

            case 2: if (sa == sb) {
                      xn(inc)(aa+1,la+2,bb,lb);
                      xn(dec1)(aa,la+3);
                    }
                    break;

            default:xn(shift_up)(aa+1,la+2,aa+1,1);
                    xn(inc)(aa+1,la+2,bb,lb);
                    xn(shift_down)(aa,la+3,aa,1);
                    if ((sa == sb)  != ((mode & 2) == 0)) xn(dec1)(aa,la+3);
                    break;

        }

        /* division dans N avec reste optionnel */
        xn(karpdiv)(aa,lc+1,bb,lb,c->val-1,2);

        /* libre la mmoire temporaire */
        xn(free)(aa);
        inc_c = 0;
    }

    /* ------------------------------ longueur et signe des rsultats */
    if (want_c) {
        sc = sa ^ sb;
        if ((inc_c) && (xn(inc1)(c->val,lc))) c->val[lc++] = 1;
        xx(make_head)(c,lc,sc);
        xx(update)(_c,c);
    }

    if (want_d) {
        switch(mode & 3) {
            case 0: sd = sb;                           break;
            case 2: sd = sb ^ SIGN_m;                  break;
            default:sd = ((inc_c) ? sa ^ SIGN_m : sa); break;
        }
        xx(make_head)(d,ld,sd);
        xx(update)(_d,d);
    }

#if defined(c_api)
    switch(mode & 12) {
        case  0: return NULL;
        case  4: return c;
        case  8: return d;
        default: return NULL;
    }
#elif defined(caml_api) || defined(ocaml_api)
    switch(mode & 12) {
        case  0: xx_pop_roots(); return Val_unit;
        case  4: xx_pop_roots(); return (value)c;
        case  8: xx_pop_roots(); return (value)d;
        default: {
            xint *r = (xint *)alloc_tuple(2);
            r[0] = c;
            r[1] = d;
            xx_pop_roots();
            return (value)r;
        }
    }
#endif /* api */

#undef  a
#undef  b
#undef  c
#undef  d
#undef _c
#undef _d
}

#if defined(caml_api) || defined(ocaml_api)

/* versions spcialises */
value xx(gquomod)  (value mode, xint *_c, xint *_d, xint a, xint b) {return xx(private_quomod)(_c,     _d,      a,b,  0|(Round_val(mode)));}
value xx(gquo)     (value mode, xint *_c,           xint a, xint b) {return xx(private_quomod)(_c,     xx_null, a,b,  0|(Round_val(mode)));}
value xx(gmod)     (value mode,           xint *_d, xint a, xint b) {return xx(private_quomod)(xx_null,_d,      a,b,  0|(Round_val(mode)));}

value xx(f_gquomod)(value mode,                     xint a, xint b) {return xx(private_quomod)(xx_null,xx_null, a,b, 12|(Round_val(mode)));}
value xx(f_gquo)   (value mode,                     xint a, xint b) {return xx(private_quomod)(xx_null,xx_null, a,b,  4|(Round_val(mode)));}
value xx(f_gmod)   (value mode,                     xint a, xint b) {return xx(private_quomod)(xx_null,xx_null, a,b,  8|(Round_val(mode)));}

value xx(quomod)   (            xint *_c, xint *_d, xint a, xint b) {return xx(private_quomod)(_c,     _d,      a,b,  0|0);}
value xx(quo)      (            xint *_c,           xint a, xint b) {return xx(private_quomod)(_c,     xx_null, a,b,  0|0);}
value xx(mod)      (                      xint *_d, xint a, xint b) {return xx(private_quomod)(xx_null,_d,      a,b,  0|0);}

value xx(f_quomod) (                                xint a, xint b) {return xx(private_quomod)(xx_null,xx_null, a,b, 12|0);}
value xx(f_quo)    (                                xint a, xint b) {return xx(private_quomod)(xx_null,xx_null, a,b,  4|0);}
value xx(f_mod)    (                                xint a, xint b) {return xx(private_quomod)(xx_null,xx_null, a,b,  8|0);}

#endif /* api */

                      /* +--------------------------+
                         |  Division par un entier  |
                         +--------------------------+ */

/*
   entre :
   a  = entier extensible
   b  = long ou Caml/Ocaml int
   _c = NULL ou adresse d'un entier extensible
   mode = long

   sortie :
   si mode & 3 = 0 : c <- floor(a/b),     d <- a - b*c (division tronque)
   si mode & 3 = 1 : c <- floor(a/b+1/2), d <- a - b*c (division centre up)
   si mode & 3 = 2 : c <- ceil(a/b),      d <- a - b*c (division majore)
   si mode & 3 = 3 : c <- ceil(a/b-1/2),  d <- a - b*c (division centre down)
   si _c != NULL : *_c <- c

   valeur de retour  : c_api   ml_api
   si mode & 12 =  0 :     d     unit
   si mode & 12 =  4 :     c        c
   si mode & 12 =  8 :     d        d
   si mode & 12 = 12 :     c    (c,d)

   erreur :
   ZERO_DIVISOR si b = 0
*/
#if defined(c_api)
long
#elif defined(caml_api) || defined(ocaml_api)
value
#endif
xx(private_quomod_1)(xint *_c, xint a, long b, long mode) {
    long la = xx_lg(a), sa = xx_sgn(a);
    long sb = b & SIGN_m;
    long d;
    long lc,sc,sd;
    int  inc_c, want_c;
    xx_push_roots_21(a,_c,c);
#ifdef caml_api
#define  a __lr.a
#define  c __lr.c
#define _c __lr._c
#endif

    /* calculer le quotient ? */
    want_c = ((_c != xx_null) || (mode & 4));

    /* b <- |b| */
#ifdef c_api
    if (sb) b = -b;
#else
    b = (sb) ? -Long_val(b) : Long_val(b);
    if (b ==  0) failwith(ZERO_DIVISOR);
#endif

    /* division dans N */
    lc = la;
    if (want_c) {
        c = xx(enlarge)(_c,lc+1);
        d = xn(div_1)(a->val,la,b,c->val);
    }
    else d = xn(mod_1)(a->val,la,b);

    /* correction */
    switch(mode & 3) {
        case 0: inc_c = ((sa != sb) && (d)); sd = sb; break;
        case 1: if ((d>b-d) || ((d==b-d) && (sa==sb))) {inc_c = 1; sd = sa^SIGN_m;}
                else                                   {inc_c = 0; sd = sa;}
                break;
        case 2: inc_c = ((sa == sb) && (d)); sd = sb^SIGN_m; break;
        default:if ((d>b-d) || ((d==b-d) && (sa!=sb))) {inc_c = 1; sd = sa^SIGN_m;}
                else                                   {inc_c = 0; sd = sa;}
                break;
    }
    if (inc_c) d = b-d;
    if (sd)    d = -d;

    /* -------------------- longueur et signe du quotient */
    if (want_c) {
        sc = sa ^ sb;
        if ((inc_c) && (xn(inc1)(c->val,lc))) c->val[lc++] = 1;
        xx(make_head)(c,lc,sc);
        xx(update)(_c,c);
    }

#if defined(c_api)
    switch(mode & 4) {
        case  0: return d;
        default: return (long)c;
    }
#elif defined(caml_api) || defined(ocaml_api)
    switch(mode & 12) {
        case  0: xx_pop_roots(); return Val_unit;
        case  4: xx_pop_roots(); return (value)c;
        case  8: xx_pop_roots(); return Val_long(d);
        default: {
            value *r = (value *)alloc_tuple(2);
            r[0] = (value)c;
            r[1] = Val_long(d);
            xx_pop_roots();
            return (value)r;
        }
    }
#endif /* api */

#undef  a
#undef  c
#undef _c
}


#if defined(caml_api) || defined(ocaml_api)

/* versions spcialises */
value xx(gquomod_1)  (value mode, xint *_c, xint a, long b) {return xx(private_quomod_1)(_c,      a,b,  8|(Round_val(mode)));}
value xx(gquo_1)     (value mode, xint *_c, xint a, long b) {return xx(private_quomod_1)(_c,      a,b,  0|(Round_val(mode)));}

value xx(f_gquomod_1)(value mode,           xint a, long b) {return xx(private_quomod_1)(xx_null, a,b, 12|(Round_val(mode)));}
value xx(f_gquo_1)   (value mode,           xint a, long b) {return xx(private_quomod_1)(xx_null, a,b,  4|(Round_val(mode)));}
value xx(f_gmod_1)   (value mode,           xint a, long b) {return xx(private_quomod_1)(xx_null, a,b,  8|(Round_val(mode)));}

value xx(quomod_1)   (            xint *_c, xint a, long b) {return xx(private_quomod_1)(_c,      a,b,  8|0);}
value xx(quo_1)      (            xint *_c, xint a, long b) {return xx(private_quomod_1)(_c,      a,b,  0|0);}

value xx(f_quomod_1) (                      xint a, long b) {return xx(private_quomod_1)(xx_null, a,b, 12|0);}
value xx(f_quo_1)    (                      xint a, long b) {return xx(private_quomod_1)(xx_null, a,b,  4|0);}
value xx(f_mod_1)    (                      xint a, long b) {return xx(private_quomod_1)(xx_null, a,b,  8|0);}

#endif /* api */