khipu/arith.hpp

/* Declarations for arithmetic types and functions.

   This file is part of khipu.

   khipu 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 3 of the License, or
   (at your option) any later version.

   This program 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 this program.  If not, see <https://www.gnu.org/licenses/>.  */

#ifndef __KP_ARITH__
#define __KP_ARITH__   1

#include "interp.hpp"

KP_DECLS_BEGIN

/* We'll be using a similar convention to GMP, but with vastly
 * simpler (and slower) algorithms. */

/* Limbs make up long integers and long floats. We settle at 32 bits in
 * order to simplify things, and so that they're easier to convert from
 * regular ints. */
#define KP_LIMB_BITS   32
const int LIMB_BITS = KP_LIMB_BITS;

typedef uint32_t limb_t;
typedef uint64_t twolimb_t;

// Normalize the limbs at [PTR .. PTR + LEN)
template <typename T>
inline void uitrim (limb_t *ptr, T& len)
{
  for (; len > 0 && ptr[len - 1] == 0; --len) ;
}

// Compare XP and YP, up to N limbs.
KP_EXPORT int uicmpn (const limb_t *xp, const limb_t *yp, int n);

// Compare [XP .. XP + XL) with [YP .. YP + YL)
KP_EXPORT int uicmp (const limb_t *xp, int xl, const limb_t *yp, int yl);

// Add VAL to [SRC .. SRC + LEN), placing the result at DST.
KP_EXPORT limb_t uiadd1 (limb_t *dst, const limb_t *src, int len, limb_t val);

// Add LEN limbs from XP and YP, placing the result at DST.
KP_EXPORT int uiaddn (limb_t *dst,
  const limb_t *xp, const limb_t *yp, int len);

// Add [YP .. YP + YL) to [XP .. XP + XL), placing the result at DST.
KP_EXPORT int uiadd (limb_t *dst,
  const limb_t *xp, int xl,
  const limb_t *yp, int yl);

// Subtract VAL from [SRC .. SRC + LEN), placing the result at DST.
KP_EXPORT limb_t uisub1 (limb_t *dst, const limb_t *src, int len, limb_t val);

// Subtract LEN limbs from XP and YP, placing the result at DST. */
KP_EXPORT int uisubn (limb_t *dst,
  const limb_t *xp, const limb_t *yp, int len);

// Subtract [YP .. YP + YL) from [XP .. XP + XL), placing the result at DST.
KP_EXPORT int uisub (limb_t *dst, const limb_t *xp,
  int xl, const limb_t *yp, int yl);

/* Left shift [SRC .. SRC + LEN) by SHIFT bits, placing the result at DST.
 * SHIFT cannot exceed 32 bits. */
KP_EXPORT int uilsh (limb_t *dst, const limb_t *src, int len, int shift);

// Multiply [SRC .. SRC + LEN) by 2^EXPO, and place the result at DST.
KP_EXPORT int uimul2exp (limb_t *dst, const limb_t *src, int len, int expo);

/* Right shift [SRC .. SRC + LEN) by SHIFT bits, placing the result at DST.
 * SHIFT cannot exceed 32 bits. */
KP_EXPORT int uirsh (limb_t *dst, const limb_t *src, int len, int shift);

// Divide [SRC .. SRC + LEN) by 2^EXPO, placing the result at DST.
KP_EXPORT int uidiv2exp (limb_t *dst, const limb_t *src, int len, int expo);

// Multiply [SRC .. SRC + LEN) by VAL, placing the result in DST.
KP_EXPORT int uimul1 (limb_t *dst, const limb_t *src, int len, limb_t val);

// Multiply [XP .. XP + XL) by [YP ..YP + YL), placing the result at DST.
KP_EXPORT limb_t uimul (interpreter *interp, limb_t *dst,
  const limb_t *xp, int xl,
  const limb_t *yp, int yl);

// Square [SRC .. SRC + LEN), placing the results at DST.
KP_EXPORT limb_t uisqr (interpreter *interp,
  limb_t *dst, const limb_t *src, int len);

/* Divide [SRC .. SRC + LEN) by VAL, placing the result at DST, 
 * and returning the remainder. */
KP_EXPORT limb_t uidivrem1 (limb_t *dst,
  const limb_t *src, int len, limb_t val);

/* Divide [XP .. XP + XL) by [YP .. YP + YL). If REM_P is true,
 * compute the remainder as well. The quotient is stored at 
 * [DST + YL .. DST + XL - YL), and the remainder at [DST .. DST + YL). */
KP_EXPORT int uidiv (interpreter *interp,
  const limb_t *xp, int xl,
  const limb_t *yp, int yl,
  limb_t *dst, bool rem_p);

/* For LEN limbs, the most significant being MSD, compute the needed number
 * characters to produce a string in radix RADIX. */
KP_EXPORT int uibsize (int radix, limb_t msd, int len);

/* Given a string of length LEN in radix RADIX, compute the number of limbs
 * needed to represent the integer number. */
KP_EXPORT int invbsize (int radix, int len);

/* Convert [SRC .. SRC + LEN) into a string of size CAP in *DST,
 * using RADIX as the radix. Returns the resulting string's length. */
KP_EXPORT int uitostr (interpreter *interp, char **dst,
  int cap, const limb_t *src, int len, int radix);

/* Convert the string at [SRC .. SRC + LEN) into a limb array, placing
 * the result at DST, and using RADIX as the radix. */
KP_EXPORT int strtoui (limb_t *dst,
  const char *src, int len, int radix);

// Same as above, for single limbs.
KP_EXPORT char* uitostr1 (char *dst, limb_t val, int radix);

KP_EXPORT intptr_t strtoui1 (const char *src, int len, int radix);

#ifdef __GNUC__

  inline unsigned int uiclz (limb_t val)
    {
      return (__builtin_clz (val));
    }

#elif defined (KP_PLATFORM_UNIX)
  #include <strings.h>

  inline unsigned int uiclz (limb_t val)
    {
      return (LIMB_BITS - 1 - ffsl (val));
    }
#else

  extern int uiclz (limb_t val);

#endif

// Convert a floating point value VAL into a limb array.
KP_EXPORT int dbltoui (double val, limb_t *dst, int& len);

// Convert [SRC .. SRC + LEN) with exponent EXPO to a floating point value.
KP_EXPORT double uitodbl (const limb_t *src, int len, int expo);

/* Convert a long float indicated by [SRC .. SRC + LEN), with exponent *EXPP
 * into a string, placing the result a DST. Produce no more than DIGS digits,
 * and use RADIX as the radix. */
KP_EXPORT int lftostr (interpreter *interp, char *dst,
  int *expp, const limb_t *src, int len, int digs, int radix);

struct num_info
{
  int type;
  int sign;
  int radix;
  int dec_start;
  int dec_end;
  int frac_end;
  int got_dot;
  int expo_sign;
  int expo_start;
};

/* Parse the string [SRC .. SRC + LEN), and determine the numeric 
 * type it evaluates to. Store the information at INFO. */
KP_EXPORT int parse_num (interpreter *interp, const char *src,
  int len, num_info& info);

/* Transform [SRC .. SRC + SLEN) into a series of limbs, and multiply
 * it by radix^EXPO, taking the radix, as well as any necessary numeric data
 * from INFO. Place the result at OUTP. */
KP_EXPORT int strtolf (interpreter *interp, const char *src,
  int slen, const num_info& info, int& expo, limb_t *outp);

union uival
{
  twolimb_t qv;
  struct
    {
#ifdef KP_LITTLE_ENDIAN
      limb_t lo, hi;
#else
      limb_t hi, lo;
#endif
    } limbs;

  void cuthi ()
    {
#ifdef KP_ARCH_WIDE
      this->qv >>= LIMB_BITS;
#else
      this->limbs.lo = this->limbs.hi;
      this->limbs.hi = 0;
#endif
    }

  // Compute M1 * M2 and store the result in *this.
  void qmul (limb_t m1, limb_t m2)
    {
      this->qv = (twolimb_t)m1 * m2;
    }

  // Compute M1 * M2 + ADD and store the result in *this.
  void qmuladd (limb_t m1, limb_t m2, limb_t add)
    {
      this->qmul (m1, m2);
      this->qv += add;
    }
};

KP_DECLS_END

#endif