Ciao,

Il Dom, 18 Febbraio 2018 9:51 am, Niels Möller ha scritto:
> I was thinking that if many of the mpq functions are written with simple
> code on top of mpz_

> But I might be mistaken, e.g, mpq/aors.c look a lot more more complex
> than in your sketch of mini-mpq.c.

A new version is attached. It does not support any string conversion, no
{set,get,inp,out}_str, as the previous one. But I added some support for
the type double.

Does the following loop make sense to "shift" the float x until it
represents an integer?

      B = 2.0 * (double) GMP_LIMB_HIGHBIT;
      for (e = 0; x != x + 0.5; e += GMP_LIMB_BITS)
        x *= B;

I began writing some tests, so that I can also claim that the code should
be correct. It should be usable now...

Comments will be welcome!

Ĝis,
m

-- 
http://bodrato.it/papers/
/* mini-mpq, a minimalistic implementation of a GNU GMP subset.

Copyright 1991-1997, 1999-2018 Free Software Foundation, Inc.
Copyright 2018 Marco BODRATO.

This file is intendet for the GNU MP Library.

This is free software; you can redistribute it and/or modify it under
the terms of the GNU Affero General Public License as published by
the Free Software Foundation; either version 3 of the License, or (at
your option) any later version.

This software 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 Affero
General Public License for more details.

You should have received copies of the GNU Affero General Public
License with the software.
If not, see https://www.gnu.org/licenses/.  */

/* Header */

#ifndef __MINI_MPQ_H__
#define __MINI_MPQ_H__

#include "mini-gmp.h"

#if defined (__cplusplus)
extern "C" {
#endif

typedef struct
{
  __mpz_struct _mp_num;
  __mpz_struct _mp_den;
} __mpq_struct;

typedef __mpq_struct mpq_t[1];

typedef const __mpq_struct *mpq_srcptr;
typedef __mpq_struct *mpq_ptr;

#define mpq_numref(Q) (&((Q)->_mp_num))
#define mpq_denref(Q) (&((Q)->_mp_den))

void mpq_abs (mpq_t, const mpq_t);
void mpq_add (mpq_t, const mpq_t, const mpq_t);
void mpq_canonicalize (mpq_t);
void mpq_clear (mpq_t);
int mpq_cmp (const mpq_t, const mpq_t);
int mpq_cmp_si (const mpq_t, signed long, unsigned long);
int mpq_cmp_ui (const mpq_t, unsigned long, unsigned long);
int mpq_cmp_z (const mpq_t, const mpz_t);
void mpq_div (mpq_t, const mpq_t, const mpq_t);
void mpq_div_2exp (mpq_t, const mpq_t, mp_bitcnt_t);
int mpq_equal (const mpq_t, const mpq_t);
double mpq_get_d (const mpq_t);
void mpq_get_den (mpz_t, const mpq_t);
void mpq_get_num (mpz_t, const mpq_t);
void mpq_init (mpq_t);
void mpq_inv (mpq_t, const mpq_t);
void mpq_mul (mpq_t, const mpq_t, const mpq_t);
void mpq_mul_2exp (mpq_t, const mpq_t, mp_bitcnt_t);
void mpq_neg (mpq_t, const mpq_t);
void mpq_set (mpq_t, const mpq_t);
void mpq_set_d (mpq_t, double);
void mpq_set_den (mpq_t, const mpz_t);
void mpq_set_num (mpq_t, const mpz_t);
void mpq_set_si (mpq_t, signed long, unsigned long);
void mpq_set_ui (mpq_t, unsigned long, unsigned long);
void mpq_set_z (mpq_t, const mpz_t);
int mpq_sgn (const mpq_t);
void mpq_sub (mpq_t, const mpq_t, const mpq_t);
void mpq_swap (mpq_t, mpq_t);

void mpz_set_q (mpz_t, const mpq_t);

#if defined (__cplusplus)
}
#endif
#endif /* __MINI_MPQ_H__ */
/* mini-mpq, a minimalistic implementation of a GNU GMP subset.

Copyright 1991-1997, 1999-2018 Free Software Foundation, Inc.
Copyright 2018 Marco BODRATO.

This file is intendet for the GNU MP Library.

This is free software; you can redistribute it and/or modify it under
the terms of the GNU Affero General Public License as published by
the Free Software Foundation; either version 3 of the License, or (at
your option) any later version.

This software 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 Affero
General Public License for more details.

You should have received copies of the GNU Affero General Public
License with the software.
If not, see https://www.gnu.org/licenses/.  */

#include "mini-mpq.h"


/* Macros */
#define NUM(x) mpq_numref(x)
#define DEN(x) mpq_denref(x)


/* MPQ helper functions */
void
mpq_init (mpq_t x)
{
  mpz_init (NUM(x));
  mpz_init_set_ui (DEN(x), 1);
}

void
mpq_clear (mpq_t x)
{
  mpz_clear (NUM(x));
  mpz_clear (DEN(x));
}

static void
__mpq_canonical_sign (mpq_t r)
{
  int cmp = DEN(r)->_mp_size;
  if (cmp <= 0)
    {
      if (cmp == 0)
	gmp_die("mpq: Fraction with zero denominator.");
      mpz_neg (DEN(r), DEN(r));
      mpz_neg (NUM(r), NUM(r));
    }
}

static void
__mpq_helper_canonicalize (mpq_t r, mpz_t g)
{
  if (NUM(r)->_mp_size == 0)
    mpz_set_ui (DEN(r), 1);
  else
    {
      __mpq_canonical_sign (r);
      mpz_gcd (g, NUM(r), DEN(r));
      mpz_div_qr (NUM(r), NULL, NUM(r), g, GMP_DIV_TRUNC);
      mpz_div_qr (DEN(r), NULL, DEN(r), g, GMP_DIV_TRUNC);
    }
}

void
mpq_canonicalize (mpq_t r)
{
  mpz_t t;

  mpz_init (t);
  __mpq_helper_canonicalize (r, t);
  mpz_clear (t);
}

void
mpq_swap (mpq_t a, mpq_t b)
{
  mpz_swap (NUM(a), NUM(b));
  mpz_swap (DEN(a), DEN(b));
}


/* MPQ assignment and conversions. */
void
mpz_set_q (mpz_t r, const mpq_t q)
{
  mpz_div_qr (r, NULL, NUM(q), DEN(q), GMP_DIV_TRUNC);
}

void
mpq_set (mpq_t r, const mpq_t q)
{
  mpz_set (NUM(r), NUM(q));
  mpz_set (DEN(r), DEN(q));
}

void
mpq_set_ui (mpq_t r, unsigned long n, unsigned long d)
{
  assert (d != 0);

  mpz_set_ui (NUM(r), n);
  mpz_set_ui (DEN(r), d);
}

void
mpq_set_si (mpq_t r, signed long n, unsigned long d)
{
  assert (d != 0);

  mpz_set_si (NUM(r), n);
  mpz_set_ui (DEN(r), d);
}

void
mpq_set_z (mpq_t r, const mpz_t n)
{
  mpz_set_ui (DEN(r), 1);
  mpz_set (NUM(r), n);
}

void
mpq_set_num (mpq_t r, const mpz_t z)
{
  mpz_set (NUM(r), z);
}

void
mpq_set_den (mpq_t r, const mpz_t z)
{
  assert (z->_mp_size != 0);
  mpz_set (DEN(r), z);
}

void
mpq_get_num (mpz_t r, const mpq_t q)
{
  mpz_set (r, NUM(q));
}

void
mpq_get_den (mpz_t r, const mpq_t q)
{
  mpz_set (r, DEN(q));
}


/* MPQ comparisons and the like. */
int
mpq_cmp (const mpq_t a, const mpq_t b)
{
  mpz_t t1, t2;
  int res;

  mpz_init (t1);
  mpz_init (t2);
  mpz_mul (t1, NUM(a), DEN(b));
  mpz_mul (t2, NUM(b), DEN(a));
  res = mpz_cmp (t1, t2);
  mpz_clear (t1);
  mpz_clear (t2);

  return res;
}

int
mpq_cmp_z (const mpq_t a, const mpz_t b)
{
  mpz_t t;
  int res;

  mpz_init (t);
  mpz_mul (t, b, DEN(a));
  res = mpz_cmp (NUM(a), t);
  mpz_clear (t);

  return res;
}

int
mpq_equal (const mpq_t a, const mpq_t b)
{
  return (mpz_cmp (NUM(a), NUM(b)) == 0) &&
    (mpz_cmp (DEN(a), DEN(b)) == 0);
}

int
mpq_cmp_ui (const mpq_t q, unsigned long n, unsigned long d)
{
  mpq_t t;
  unsigned long l_n = n;
  unsigned long l_d = d;
  
  assert (d != 0);

  mpz_roinit_n (NUM (t), &l_n, 1);
  mpz_roinit_n (DEN (t), &l_d, 1);
  return mpq_cmp (q, t);
}

int
mpq_cmp_si (const mpq_t q, signed long n, unsigned long d)
{
  assert (d != 0);

  if (n >= 0)
    return mpq_cmp_ui (q, n, d);
  else
    {
      mpq_t t;
      unsigned long l_n = GMP_NEG_CAST (unsigned long, n);
      unsigned long l_d = d;

      mpz_roinit_n (NUM (t), &l_n, -1);
      mpz_roinit_n (DEN (t), &l_d, 1);
      return mpq_cmp (q, t);
    }
}

int
mpq_sgn (const mpq_t a)
{
  return mpz_sgn (NUM (a));
}


/* MPQ arithmetic. */
void
mpq_abs (mpq_t r, const mpq_t q)
{
  mpz_abs (NUM(r), NUM(q));
  mpz_set (DEN(r), DEN(q));
}

void
mpq_neg (mpq_t r, const mpq_t q)
{
  mpz_neg (NUM(r), NUM(q));
  mpz_set (DEN(r), DEN(q));
}

void
mpq_add (mpq_t r, const mpq_t a, const mpq_t b)
{
  mpz_t t;

  mpz_init (t);
  mpz_mul (t, NUM(a), DEN(b));
  mpz_addmul (t, NUM(b), DEN(a));
  mpz_mul (DEN(r), DEN(a), DEN(b));
  mpz_swap (NUM(r), t);
  __mpq_helper_canonicalize (r, t);
  mpz_clear (t);
}

void
mpq_sub (mpq_t r, const mpq_t a, const mpq_t b)
{
  mpz_t t;

  mpz_init (t);
  mpz_mul (t, NUM(a), DEN(b));
  mpz_submul (t, NUM(b), DEN(a));
  mpz_mul (DEN(r), DEN(a), DEN(b));
  mpz_swap (NUM(r), t);
  __mpq_helper_canonicalize (r, t);
  mpz_clear (t);
}

void
mpq_div (mpq_t r, const mpq_t a, const mpq_t b)
{
  mpz_t t;

  mpz_init (t);
  mpz_mul (t, NUM(a), DEN(b));
  mpz_mul (DEN(r), DEN(a), NUM(b));
  mpz_swap (NUM(r), t);
  __mpq_helper_canonicalize (r, t);
  mpz_clear (t);
}

void
mpq_mul (mpq_t r, const mpq_t a, const mpq_t b)
{
  mpz_mul (NUM(r), NUM(a), NUM(b));
  mpz_mul (DEN(r), DEN(a), DEN(b));
  mpq_canonicalize (r);
}

void
mpq_div_2exp (mpq_t r, const mpq_t q, mp_bitcnt_t e)
{
  mpz_mul_2exp (DEN (r), DEN(q), e);
  mpz_set (NUM (r), NUM(q));
  mpq_canonicalize (r);
}

void
mpq_mul_2exp (mpq_t r, const mpq_t q, mp_bitcnt_t e)
{
  mpz_mul_2exp (NUM (r), NUM(q), e);
  mpz_set (DEN (r), DEN(q));
  mpq_canonicalize (r);
}

void
mpq_inv (mpq_t r, const mpq_t q)
{
  mpq_set (r, q);
  mpz_swap (DEN(r), NUM(r));
  __mpq_canonical_sign (r);
}


/* MPQ to/from double. */
void
mpq_set_d (mpq_t r, double x)
{
  mpz_set_ui (DEN (r), 1);

  /* x != x is true when x is a NaN, and x == x * 0.5 is true when x is
     zero or infinity. */
  if (x != x || x == x * 0.5)
    NUM (r)->_mp_size = 0;
  else
    {
      double B;
      mp_bitcnt_t e;

      B = 2.0 * (double) GMP_LIMB_HIGHBIT;
      for (e = 0; x != x + 0.5; e += GMP_LIMB_BITS)
	x *= B;

      mpz_set_d (NUM (r), x);
      mpq_div_2exp (r, r, e);
    }
}

double
mpq_get_d (const mpq_t u)
{
  mp_bitcnt_t ne, de, ee;
  mpz_t z;
  double B, ret;

  ne = mpz_sizeinbase (NUM (u), 2);
  de = mpz_sizeinbase (DEN (u), 2);

  ee = 8 * sizeof (double);
  if (de == 1 || ne > de + ee)
    ee = 0;
  else
    ee = (ee + de - ne) / GMP_LIMB_BITS + 1;

  mpz_init (z);
  mpz_mul_2exp (z, NUM (u), ee * GMP_LIMB_BITS);
  mpz_div_qr (z, NULL, z, DEN(u), GMP_DIV_TRUNC);
  ret = mpz_get_d (z);
  mpz_clear (z);

  B = 2.0 * (double) GMP_LIMB_HIGHBIT;
  for (B = 1 / B; ee != 0; --ee)
    ret *= B;

  return ret;
}
_______________________________________________
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel

Reply via email to