1  /* mpq_div -- divide two rational numbers.
       2  
       3  Copyright 1991, 1994-1996, 2000, 2001, 2015, 2018 Free Software
       4  Foundation, Inc.
       5  
       6  This file is part of the GNU MP Library.
       7  
       8  The GNU MP Library is free software; you can redistribute it and/or modify
       9  it under the terms of either:
      10  
      11    * the GNU Lesser General Public License as published by the Free
      12      Software Foundation; either version 3 of the License, or (at your
      13      option) any later version.
      14  
      15  or
      16  
      17    * the GNU General Public License as published by the Free Software
      18      Foundation; either version 2 of the License, or (at your option) any
      19      later version.
      20  
      21  or both in parallel, as here.
      22  
      23  The GNU MP Library is distributed in the hope that it will be useful, but
      24  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
      25  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
      26  for more details.
      27  
      28  You should have received copies of the GNU General Public License and the
      29  GNU Lesser General Public License along with the GNU MP Library.  If not,
      30  see https://www.gnu.org/licenses/.  */
      31  
      32  #include "gmp-impl.h"
      33  
      34  
      35  void
      36  mpq_div (mpq_ptr quot, mpq_srcptr op1, mpq_srcptr op2)
      37  {
      38    mpz_t gcd1, gcd2;
      39    mpz_t tmp1, tmp2;
      40    mp_size_t op1_size;
      41    mp_size_t op2_size;
      42    mp_size_t alloc;
      43    TMP_DECL;
      44  
      45    op2_size = SIZ(NUM(op2));
      46  
      47    if (UNLIKELY (op2_size == 0))
      48      DIVIDE_BY_ZERO;
      49  
      50    if (UNLIKELY (quot == op2))
      51      {
      52        if (UNLIKELY (op1 == op2))
      53  	{
      54  	  mpq_set_ui (quot, 1, 1);
      55  	  return;
      56  	}
      57  
      58        /* We checked for op1 == op2: we are not in the x=x/x case.
      59  	 We compute x=y/x by computing x=inv(x)*y */
      60        MPN_PTR_SWAP (PTR(NUM(quot)), ALLOC(NUM(quot)),
      61  		    PTR(DEN(quot)), ALLOC(DEN(quot)));
      62        if (op2_size > 0)
      63  	{
      64  	  SIZ(NUM(quot)) = SIZ(DEN(quot));
      65  	  SIZ(DEN(quot)) = op2_size;
      66  	}
      67        else
      68  	{
      69  	  SIZ(NUM(quot)) = - SIZ(DEN(quot));
      70  	  SIZ(DEN(quot)) = - op2_size;
      71  	}
      72        mpq_mul (quot, quot, op1);
      73        return;
      74      }
      75  
      76    op1_size = ABSIZ(NUM(op1));
      77  
      78    if (op1_size == 0)
      79      {
      80        /* We special case this to simplify allocation logic; gcd(0,x) = x
      81  	 is a singular case for the allocations.  */
      82        SIZ(NUM(quot)) = 0;
      83        MPZ_NEWALLOC (DEN(quot), 1)[0] = 1;
      84        SIZ(DEN(quot)) = 1;
      85        return;
      86      }
      87  
      88    op2_size = ABS(op2_size);
      89  
      90    TMP_MARK;
      91  
      92    alloc = MIN (op1_size, op2_size);
      93    MPZ_TMP_INIT (gcd1, alloc);
      94  
      95    alloc = MAX (op1_size, op2_size);
      96    MPZ_TMP_INIT (tmp1, alloc);
      97  
      98    op2_size = SIZ(DEN(op2));
      99    op1_size = SIZ(DEN(op1));
     100  
     101    alloc = MIN (op1_size, op2_size);
     102    MPZ_TMP_INIT (gcd2, alloc);
     103  
     104    alloc = MAX (op1_size, op2_size);
     105    MPZ_TMP_INIT (tmp2, alloc);
     106  
     107    /* QUOT might be identical to OP1, so don't store the result there
     108       until we are finished with the input operand.  We can overwrite
     109       the numerator of QUOT when we are finished with the numerator of
     110       OP1. */
     111  
     112    mpz_gcd (gcd1, NUM(op1), NUM(op2));
     113    mpz_gcd (gcd2, DEN(op2), DEN(op1));
     114  
     115    mpz_divexact_gcd (tmp1, NUM(op1), gcd1);
     116    mpz_divexact_gcd (tmp2, DEN(op2), gcd2);
     117  
     118    mpz_mul (NUM(quot), tmp1, tmp2);
     119  
     120    mpz_divexact_gcd (tmp1, NUM(op2), gcd1);
     121    mpz_divexact_gcd (tmp2, DEN(op1), gcd2);
     122  
     123    mpz_mul (DEN(quot), tmp1, tmp2);
     124  
     125    /* Keep the denominator positive.  */
     126    if (SIZ(DEN(quot)) < 0)
     127      {
     128        SIZ(DEN(quot)) = -SIZ(DEN(quot));
     129        SIZ(NUM(quot)) = -SIZ(NUM(quot));
     130      }
     131  
     132    TMP_FREE;
     133  }