avr/libf7/libf7.c

/* Copyright (C) 2019-2021 Free Software Foundation, Inc.

   This file is part of LIBF7, which is part of GCC.

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

   GCC 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 General Public License
   for more details.

   Under Section 7 of GPL version 3, you are granted additional
   permissions described in the GCC Runtime Library Exception, version
   3.1, as published by the Free Software Foundation.

   You should have received a copy of the GNU General Public License and
   a copy of the GCC Runtime Library Exception along with this program;
   see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
   <http://www.gnu.org/licenses/>.  */

#include "libf7.h"

#ifndef __AVR_TINY__

#define ALIAS(X, Y) \
  F7_WEAK __attribute__((__alias__(F7_STRINGY(X)))) __typeof__(X) Y;

#define DALIAS(...) // empty
#define LALIAS(...) // empty

#ifndef IN_LIBGCC2

#include <stdio.h>
#include <assert.h>

#define in_libgcc false

_Static_assert (sizeof (f7_t) == 10 && F7_MANT_BYTES == 7,
		"libf7 will only work with 7-byte mantissa.");
#else

#define in_libgcc true

#if __SIZEOF_DOUBLE__ == 8
#undef  DALIAS
#define DALIAS(X,Y) \
  F7_WEAK __attribute__((__alias__(F7_STRINGY(X)))) __typeof__(X) Y;
#endif

#if __SIZEOF_LONG_DOUBLE__ == 8
#undef  LALIAS
#define LALIAS(X,Y) \
  F7_WEAK __attribute__((__alias__(F7_STRINGY(X)))) __typeof__(X) Y;
#endif

#endif // in libgcc

static F7_INLINE
void f7_assert (bool x)
{
  if (!in_libgcc && !x)
    __builtin_abort();
}

static F7_INLINE
int16_t abs_ssat16 (int16_t a)
{
  _Sat _Fract sa = __builtin_avr_rbits (a);
  return __builtin_avr_bitsr (__builtin_avr_absr (sa));
}

static F7_INLINE
int16_t add_ssat16 (int16_t a, int16_t b)
{
  _Sat _Fract sa = __builtin_avr_rbits (a);
  _Sat _Fract sb = __builtin_avr_rbits (b);
  return __builtin_avr_bitsr (sa + sb);
}

static F7_INLINE
int16_t sub_ssat16 (int16_t a, int16_t b)
{
  _Sat _Fract sa = __builtin_avr_rbits (a);
  _Sat _Fract sb = __builtin_avr_rbits (b);
  return __builtin_avr_bitsr (sa - sb);
}

static F7_INLINE
int8_t ssat8_range (int16_t a, int8_t range)
{
  if (a >= range)
    return range;
  if (a <= -range)
    return -range;
  return a;
}


#define IN_LIBF7_H
  #define F7_CONST_DEF(NAME, FLAGS, M6, M5, M4, M3, M2, M1, M0, EXPO) \
    F7_UNUSED static const uint8_t F7_(const_##NAME##_msb)  = M6;     \
    F7_UNUSED static const int16_t F7_(const_##NAME##_expo) = EXPO;
  #include "libf7-const.def"
  #undef F7_CONST_DEF
#undef IN_LIBF7_H


/*
  libgcc naming converntions for conversions:

   __float<fmode><fmode>  : Convert float modes.
   __floatun<imode><fmode>: Convert unsigned integral to float.
   __fix<fmode><imode>    : Convert float to signed integral.
   __fixuns<fmode><imode> : Convert float to unsigned integral.
*/


#ifdef F7MOD_floatundidf_
F7_WEAK
f7_double_t __floatundidf (uint64_t x)
{
  f7_t xx;
  f7_set_u64 (&xx, x);
  return f7_get_double (&xx);
}
#endif // F7MOD_floatundidf_


#ifdef F7MOD_floatdidf_
F7_WEAK
f7_double_t __floatdidf (int64_t x)
{
  f7_t xx;
  f7_set_s64 (&xx, x);
  return f7_get_double (&xx);
}
#endif // F7MOD_floatdidf_


#ifdef F7MOD_init_
f7_t* f7_init_impl (uint64_t mant, uint8_t flags, f7_t *cc, int16_t expo)
{
  flags &= F7_FLAGS;
  if (f7_class_number (flags))
    {
      uint8_t msb;
      while ((__builtin_memcpy (&msb, (uint8_t*) &mant + 7, 1), msb))
	{
	  mant >>= 1;
	  expo = add_ssat16 (expo, 1);
	}
      *(uint64_t*) cc->mant = mant;
      cc->expo = add_ssat16 (expo, F7_MANT_BITS-1);

      cc = f7_normalize_asm (cc);
    }

  cc->flags = flags;

  return cc;
}
#endif // F7MOD_init_


#ifdef F7MOD_set_s16_
f7_t* f7_set_s16_impl (f7_t *cc, int16_t i16)
{
  uint16_t u16 = (uint16_t) i16;
  uint8_t flags = 0;
  if (i16 < 0)
    {
      u16 = -u16;
      flags = F7_FLAG_sign;
    }
  f7_set_u16_impl (cc, u16);
  cc->flags = flags;
  return cc;
}
#endif // F7MOD_set_s16_


#ifdef F7MOD_set_u16_
f7_t* f7_set_u16_impl (f7_t *cc, uint16_t u16)
{
  f7_clr (cc);
  F7_MANT_HI2 (cc) = u16;
  cc->expo = 15;
  return f7_normalize_asm (cc);
}
#endif // F7MOD_set_u16_


#ifdef F7MOD_set_s32_
f7_t* f7_set_s32 (f7_t *cc, int32_t i32)
{
  uint32_t u32 = (uint32_t) i32;
  uint8_t flags = 0;
  if (i32 < 0)
    {
      u32 = -u32;
      flags = F7_FLAG_sign;
    }
  cc = f7_set_u32 (cc, u32);
  cc->flags = flags;
  return cc;
}
ALIAS (f7_set_s32, f7_floatsidf)
#endif // F7MOD_set_s32_


#ifdef F7MOD_set_u32_
f7_t* f7_set_u32 (f7_t *cc, uint32_t u32)
{
  f7_clr (cc);
  F7_MANT_HI4 (cc) = u32;
  cc->expo = 31;
  return f7_normalize_asm (cc);
}
ALIAS (f7_set_u32, f7_floatunsidf)
#endif // F7MOD_set_u32_


// IEEE 754 single
// float =  s  bbbbbbbb mmmmmmmmmmmmmmmmmmmmmmm
//	   31
// s = sign
// b = biased exponent, bias = 127
// m = mantissa

// +0	      =	 0 0 0
// -0	      =	 1 0 0
// Inf	      =	 S B 0	=  S * Inf, B = 0xff
// NaN	      =	 S B M,	   B = 0xff, M != 0
// Sub-normal =	 S 0 M	=  S * 0.M * 2^{1 - bias}, M != 0
// Normal     =  S B M  =  S * 1.M * 2^{B - bias}, B = 1 ... 0xfe

#define FLT_DIG_EXP   8
#define FLT_DIG_MANT  (31 - FLT_DIG_EXP)
#define FLT_MAX_EXP   ((1 << FLT_DIG_EXP) - 1)
#define FLT_EXP_BIAS  (FLT_MAX_EXP >> 1)

#ifdef F7MOD_set_float_
F7_WEAK
void f7_set_float (f7_t *cc, float f)
{
  uint32_t val32;

  _Static_assert (__SIZEOF_FLOAT__ == 4, "");
  _Static_assert (__FLT_MANT_DIG__ == 24, "");
  __builtin_memcpy (&val32, &f, __SIZEOF_FLOAT__);

  uint16_t val16 = val32 >> 16;
  val16 >>= FLT_DIG_MANT - 16;

  uint8_t expo_biased = val16 & FLT_MAX_EXP;
  bool sign = val16 & (1u << FLT_DIG_EXP);

  f7_clr (cc);

  uint32_t mant = val32 & ((1ul << FLT_DIG_MANT) -1);

  if (mant == 0)
    {
      if (expo_biased == 0)
	return;
      if (expo_biased >= FLT_MAX_EXP)
	return f7_set_inf (cc, sign);
    }

  if (expo_biased == 0)
    expo_biased = 1;   // Sub-normal: biased expo of 1 was encoded as 0.
  else if (expo_biased < FLT_MAX_EXP)
    mant |= (1ul << FLT_DIG_MANT);
  else
    return f7_set_nan (cc);

  __builtin_memcpy (& F7_MANT_HI4 (cc), &mant, 4);

  cc->expo = expo_biased - FLT_EXP_BIAS + 31 - FLT_DIG_MANT;
  f7_normalize_asm (cc);
  f7_set_sign (cc, sign);
}
ALIAS (f7_set_float, f7_extendsfdf2)
#endif // F7MOD_set_float_


#ifdef F7MOD_get_float_
static F7_INLINE
float make_float (uint32_t x)
{
  float ff;
  __builtin_memcpy (&ff, &x, 4);
  return ff;

}

F7_WEAK
float f7_get_float (const f7_t *aa)
{
  uint8_t a_class = f7_classify (aa);

  if (f7_class_nan (a_class))
    return make_float (0xffc00000 /* NaN: Biased expo = 0xff, mant != 0 */);

  uint32_t mant;
  __builtin_memcpy (&mant, &F7_MANT_CONST_HI4 (aa), 4);

  uint8_t expo8 = 0;
  uint8_t mant_offset = FLT_DIG_EXP;
  int16_t c_expo = add_ssat16 (aa->expo, FLT_EXP_BIAS);

  if (f7_class_zero (a_class) || c_expo <= -FLT_DIG_MANT)
    {
      // Zero or tiny.
      return 0.0f;
    }
  else if (c_expo >= FLT_MAX_EXP || f7_class_inf (a_class))
    {
      // Inf or overflow.
      expo8 = FLT_MAX_EXP;
      mant = 0;
    }
  else if (c_expo > 0)
    {
      // Normal.
      expo8 = c_expo;
    }
  else
    {
      // Sub-normal:  -DIG_MANT < c_expo <= 0.
      // Encoding of 0 represents a biased exponent of 1.
      // mant_offset in 9...31.
      expo8 = 0;
      mant_offset += 1 - c_expo;
    }

  uint16_t expo16 = expo8 << (FLT_DIG_MANT - 16);

  if (f7_class_sign (a_class))
    expo16 |= 1u << (FLT_DIG_EXP + FLT_DIG_MANT - 16);

  mant >>= mant_offset;

  __asm ("cbr %T0%t2, 1 << (7 & %2)"  "\n\t"
	 "or  %C0, %A1"		      "\n\t"
	 "or  %D0, %B1"
	 : "+d" (mant)
	 : "r" (expo16), "n" (FLT_DIG_MANT));

  return make_float (mant);
}
F7_PURE ALIAS (f7_get_float, f7_truncdfsf2)
#endif // F7MOD_get_float_

#define DBL_DIG_EXP   11
#define DBL_DIG_MANT  (63 - DBL_DIG_EXP)
#define DBL_MAX_EXP   ((1 << DBL_DIG_EXP) - 1)
#define DBL_EXP_BIAS  (DBL_MAX_EXP >> 1)


#ifdef F7MOD_set_double_
void f7_set_double_impl (f7_double_t val64, f7_t *cc)
{
  f7_clr (cc);
  register uint64_t mant __asm ("r18") = val64 & ((1ull << DBL_DIG_MANT) -1);

  uint16_t val16 = 3[(uint16_t*) & val64];
  val16 >>= DBL_DIG_MANT - 48;

  uint16_t expo_biased = val16 & DBL_MAX_EXP;
  bool sign = val16 & (1u << DBL_DIG_EXP);

  if (mant == 0)
    {
      if (expo_biased == 0)
	return;
      if (expo_biased >= DBL_MAX_EXP)
	return f7_set_inf (cc, sign);
    }
  __asm ("" : "+r" (mant));

  if (expo_biased == 0)
    expo_biased = 1;   // Sub-normal: biased expo of 1 was encoded as 0.
  else if (expo_biased < DBL_MAX_EXP)
    mant |= (1ull << DBL_DIG_MANT);
  else
    return f7_set_nan (cc);

  *(uint64_t*) & cc->mant = mant;

  cc->expo = expo_biased - DBL_EXP_BIAS + 63 - DBL_DIG_MANT - 8;
  f7_normalize_asm (cc);
  f7_set_sign (cc, sign);
}
#endif // F7MOD_set_double_


#ifdef F7MOD_set_pdouble_
void f7_set_pdouble (f7_t *cc, const f7_double_t *val64)
{
  f7_set_double (cc, *val64);
}
#endif // F7MOD_set_pdouble_


#ifdef F7MOD_get_double_
static F7_INLINE
uint64_t clr_r18 (void)
{
  extern void __clr_8 (void);
  register uint64_t r18 __asm ("r18");
  __asm ("%~call %x[f]" : "=r" (r18) : [f] "i" (__clr_8));
  return r18;
}

static F7_INLINE
f7_double_t make_double (uint64_t x)
{
  register f7_double_t r18 __asm ("r18") = x;
  __asm ("" : "+r" (r18));
  return r18;
}

F7_WEAK
f7_double_t f7_get_double (const f7_t *aa)
{
  uint8_t a_class = f7_classify (aa);

  if (f7_class_nan (a_class))
    {
      uint64_t nan = clr_r18() | (0x7fffull << 48);
      return make_double (nan);
    }

  uint64_t mant;
  __builtin_memcpy (&mant, & aa->mant, 8);

  mant &= 0x00ffffffffffffff;

  // FIXME: For subnormals, rounding is premature and should be
  //	    done *after* the mantissa has been shifted into place
  //	    (or the round value be shifted left accordingly).
  // Round.
  mant += 1u << (F7_MANT_BITS - (1 + DBL_DIG_MANT) - 1);

  uint8_t dex;
  register uint64_t r18 __asm ("r18") = mant;
  // dex = Overflow ? 1 : 0.
  __asm ("bst %T[mant]%T[bitno]"  "\n\t"
	 "clr %0"		  "\n\t"
	 "bld %0,0"
	 : "=r" (dex), [mant] "+r" (r18)
	 : [bitno] "n" (64 - 8));

  mant = r18 >> dex;

  uint16_t expo16 = 0;
  uint16_t mant_offset = DBL_DIG_EXP - 8;
  int16_t c_expo = add_ssat16 (aa->expo, DBL_EXP_BIAS + dex);

  if (f7_class_zero (a_class) || c_expo <= -DBL_DIG_MANT)
    {
      // Zero or tiny.
      return make_double (clr_r18());
    }
  else if (c_expo >= DBL_MAX_EXP || f7_class_inf (a_class))
    {
      // Inf or overflow.
      expo16 = DBL_MAX_EXP;
      mant = clr_r18();
    }
  else if (c_expo > 0)
    {
      // Normal.
      expo16 = c_expo;
    }
  else
    {
      // Sub-normal:  -DIG_MANT < c_expo <= 0.
      // Encoding expo of 0 represents a biased exponent of 1.
      // mant_offset in 5...55 = 63-8.
      mant_offset += 1 - c_expo;
    }

  expo16 <<= (DBL_DIG_MANT - 48);

  if (f7_class_sign (a_class))
    expo16 |= 1u << (DBL_DIG_EXP + DBL_DIG_MANT - 48);

  // mant >>= mant_offset;
  mant = f7_lshrdi3 (mant, mant_offset);

  r18 = mant;
  __asm ("cbr %T0%t2, 1 << (7 & %2)"  "\n\t"
	 "or  %r0+6, %A1"	      "\n\t"
	 "or  %r0+7, %B1"
	 : "+r" (r18)
	 : "r" (expo16), "n" (DBL_DIG_MANT));

  return make_double (r18);
}
#endif // F7MOD_get_double_


#ifdef F7MOD_fabs_
F7_WEAK
void f7_fabs (f7_t *cc, const f7_t *aa)
{
  f7_abs (cc, aa);
}
#endif // F7MOD_fabs_


#ifdef F7MOD_neg_
F7_WEAK
f7_t* f7_neg (f7_t *cc, const f7_t *aa)
{
  f7_copy (cc, aa);

  uint8_t c_class = f7_classify (cc);

  if (! f7_class_zero (c_class))
    cc->sign = ! f7_class_sign (c_class);

  return cc;
}
#endif // F7MOD_neg_


#ifdef F7MOD_frexp_
F7_WEAK
void f7_frexp (f7_t *cc, const f7_t *aa, int *expo)
{
  uint8_t a_class = f7_classify (aa);

  if (f7_class_nan (a_class))
    return f7_set_nan (cc);

  if (f7_class_inf (a_class) || aa->expo == INT16_MAX)
    return f7_set_inf (cc, f7_class_sign (a_class));

  if (! f7_msbit (aa))
    {
      *expo = 0;
      return f7_clr (cc);
    }

  *expo = 1 + aa->expo;
  cc->flags = a_class & F7_FLAG_sign;
  cc->expo = -1;
  f7_copy_mant (cc, aa);
}
#endif // F7MOD_frexp_

#ifdef F7MOD_get_s16_
F7_WEAK
int16_t f7_get_s16 (const f7_t *aa)
{
  extern int16_t to_s16 (const f7_t*, uint8_t) F7ASM(f7_to_integer_asm);
  return to_s16 (aa, 0xf);
}
#endif // F7MOD_get_s16_


#ifdef F7MOD_get_s32_
F7_WEAK
int32_t f7_get_s32 (const f7_t *aa)
{
  extern int32_t to_s32 (const f7_t*, uint8_t) F7ASM(f7_to_integer_asm);
  return to_s32 (aa, 0x1f);
}
F7_PURE ALIAS (f7_get_s32, f7_fixdfsi)
#endif // F7MOD_get_s32_


#ifdef F7MOD_get_s64_
  F7_WEAK
  int64_t f7_get_s64 (const f7_t *aa)
{
  extern int64_t to_s64 (const f7_t*, uint8_t) F7ASM(f7_to_integer_asm);
  return to_s64 (aa, 0x3f);
}
F7_PURE ALIAS (f7_get_s64, f7_fixdfdi)
#endif // F7MOD_get_s64_

#ifdef F7MOD_get_u16_
  F7_WEAK
  uint16_t f7_get_u16 (const f7_t *aa)
{
  extern uint16_t to_u16 (const f7_t*, uint8_t) F7ASM(f7_to_unsigned_asm);
  return to_u16 (aa, 0xf);
}
#endif // F7MOD_get_u16_


#ifdef F7MOD_get_u32_
F7_WEAK
uint32_t f7_get_u32 (const f7_t *aa)
{
  extern uint32_t to_u32 (const f7_t*, uint8_t) F7ASM(f7_to_unsigned_asm);
  return to_u32 (aa, 0x1f);
}
F7_PURE ALIAS (f7_get_u32, f7_fixunsdfsi)
#endif // F7MOD_get_u32_


#ifdef F7MOD_get_u64_
F7_WEAK
uint64_t f7_get_u64 (const f7_t *aa)
{
  extern int64_t to_u64 (const f7_t*, uint8_t) F7ASM(f7_to_unsigned_asm);
  return to_u64 (aa, 0x3f);
}
F7_PURE ALIAS (f7_get_u64, f7_fixunsdfdi)
#endif // F7MOD_get_u64_


#ifdef F7MOD_cmp_unordered_
F7_NOINLINE
static int8_t cmp_u8 (uint8_t a_class, uint8_t b_class, bool sign_a);

F7_WEAK
int8_t f7_cmp_unordered (const f7_t *aa, const f7_t *bb, bool with_sign)
{
  uint8_t a_class = f7_classify (aa);
  uint8_t b_class = f7_classify (bb);

  uint8_t a_sign = f7_class_sign (a_class) & with_sign;
  uint8_t b_sign = f7_class_sign (b_class) & with_sign;
  uint8_t ab_class = a_class | b_class;
  ab_class &= with_sign - 2;

  if (f7_class_nan (ab_class))
    return INT8_MIN;

  if (a_sign != b_sign)
    return b_sign - a_sign;

  if (f7_class_inf (ab_class))
    return cmp_u8 (a_class, b_class, a_sign);

  if (f7_class_zero (ab_class))
    return cmp_u8 (b_class, a_class, a_sign);

  if (aa->expo < bb->expo)
    return a_sign ? 1 : -1;

  if (aa->expo > bb->expo)
    return a_sign ? -1 : 1;

  return cmp_u8 (1 + f7_cmp_mant (aa, bb), 1, a_sign);
}


int8_t cmp_u8 (uint8_t a_class, uint8_t b_class, bool sign_a)
{
  int8_t c;
  __asm ("sub  %[a], %[b]"    "\n\t"
	 "breq 1f"	      "\n\t"
	 "sbc  %[c], %[c]"    "\n\t"
	 "sbci %[c], -1"      "\n\t"
	 "sbrc %[s], 0"	      "\n\t"
	 "neg  %[c]"	      "\n\t"
	 "1:"
	 : [c] "=d" (c)
	 : [a] "0" (a_class), [b] "r" (b_class), [s] "r" (sign_a));
  return c;
}
#endif // F7MOD_cmp_unordered_


#ifdef F7MOD_cmp_abs_
F7_WEAK
int8_t f7_cmp_abs (const f7_t *aa, const f7_t *bb)
{
  return f7_cmp_unordered (aa, bb, false /* no signs */);
}
#endif // F7MOD_cmp_abs_


#ifdef F7MOD_cmp_
F7_WEAK
int8_t f7_cmp (const f7_t *aa, const f7_t *bb)
{
  return f7_cmp_unordered (aa, bb, true /* with signs */);
}
#endif // F7MOD_cmp_


#ifdef F7MOD_abscmp_msb_ge_
// Compare absolute value of Number aa against a f7_t represented
// by msb and expo.
F7_WEAK
bool f7_abscmp_msb_ge (const f7_t *aa, uint8_t msb, int16_t expo)
{
  uint8_t a_msb = aa->mant[F7_MANT_BYTES - 1];

  if (0 == (0x80 & a_msb))
    // 0 or subnormal.
    return false;

  return aa->expo == expo
    ? a_msb >= msb
    : aa->expo > expo;
}
#endif // F7MOD_abscmp_msb_ge_

#ifdef F7MOD_lt_
F7_WEAK
bool f7_lt_impl (const f7_t *aa, const f7_t *bb)
{
  return f7_lt (aa, bb);
}
#endif // F7MOD_lt_

#ifdef F7MOD_le_
F7_WEAK
bool f7_le_impl (const f7_t *aa, const f7_t *bb)
{
  return f7_le (aa, bb);
}
#endif // F7MOD_le_

#ifdef F7MOD_gt_
F7_WEAK
bool f7_gt_impl (const f7_t *aa, const f7_t *bb)
{
  return f7_gt (aa, bb);
}
#endif // F7MOD_gt_

#ifdef F7MOD_ge_
F7_WEAK
bool f7_ge_impl (const f7_t *aa, const f7_t *bb)
{
  return f7_ge (aa, bb);
}
#endif // F7MOD_ge_

#ifdef F7MOD_ne_
F7_WEAK
bool f7_ne_impl (const f7_t *aa, const f7_t *bb)
{
  return f7_ne (aa, bb);
}
#endif // F7MOD_ne_

#ifdef F7MOD_eq_
F7_WEAK
bool f7_eq_impl (const f7_t *aa, const f7_t *bb)
{
  return f7_eq (aa, bb);
}
#endif // F7MOD_eq_


#ifdef F7MOD_unord_
F7_WEAK
bool f7_unord_impl (const f7_t *aa, const f7_t *bb)
{
  return f7_unordered (aa, bb);
}
#endif // F7MOD_unord_


#ifdef F7MOD_minmax_
F7_WEAK
f7_t* f7_minmax (f7_t *cc, const f7_t *aa, const f7_t *bb, bool do_min)
{
  int8_t cmp = f7_cmp_unordered (aa, bb, true /* with signs */);

  if (cmp == INT8_MIN)
    return (f7_set_nan (cc), cc);

  if (do_min)
    cmp = -cmp;

  return f7_copy (cc, cmp >= 0 ? aa : bb);
}
#endif // F7MOD_minmax_


#ifdef F7MOD_fmax_
F7_WEAK
f7_t* f7_fmax (f7_t *cc, const f7_t *aa, const f7_t *bb)
{
  return f7_minmax (cc, aa, bb, false);
}
ALIAS (f7_fmax, f7_max)
#endif // F7MOD_fmax_


#ifdef F7MOD_fmin_
F7_WEAK
f7_t* f7_fmin (f7_t *cc, const f7_t *aa, const f7_t *bb)
{
  return f7_minmax (cc, aa, bb, true);
}
ALIAS (f7_fmin, f7_min)
#endif // F7MOD_fmin_


#ifdef F7MOD_mulx_
F7_WEAK
uint8_t f7_mulx (f7_t *cc, const f7_t *aa, const f7_t *bb, bool no_rounding)
{
  uint8_t a_class = f7_classify (aa);
  uint8_t b_class = f7_classify (bb);
  // From this point on, no more access aa->flags or bb->flags
  // to avoid early-clobber when writing cc->flags.

  uint8_t ab_class = a_class | b_class;
  // If either value is NaN, return NaN.
  if (f7_class_nan (ab_class)
      // Any combination of Inf and 0.
      || (f7_class_zero (ab_class) && f7_class_inf (ab_class)))
    {
      cc->flags = F7_FLAG_nan;
      return 0;
    }
  // If either value is 0.0, return 0.0.
  if (f7_class_zero (ab_class))
    {
      f7_clr (cc);
      return 0;
    }
  // We have 2 non-zero numbers-or-INF.

  uint8_t c_sign = (a_class ^ b_class) & F7_FLAG_sign;
  uint8_t c_inf  = ab_class & F7_FLAG_inf;
  cc->flags = c_sign | c_inf;
  if (c_inf)
    return 0;

  int16_t expo = add_ssat16 (aa->expo, bb->expo);
  // Store expo and handle expo = INT16_MIN  and INT16_MAX.
  if (f7_store_expo (cc, expo))
    return 0;

  return f7_mul_mant_asm (cc, aa, bb, no_rounding);
}
#endif // F7MOD_mulx_


#ifdef F7MOD_square_
F7_WEAK
void f7_square (f7_t *cc, const f7_t *aa)
{
  f7_mul (cc, aa, aa);
}
#endif // F7MOD_square_


#ifdef F7MOD_mul_
F7_WEAK
void f7_mul (f7_t *cc, const f7_t *aa, const f7_t *bb)
{
  f7_mulx (cc, aa, bb, false);
}
#endif // F7MOD_mul_


#ifdef F7MOD_Iadd_
F7_WEAK void f7_Iadd (f7_t *cc, const f7_t *aa) { f7_add (cc, cc, aa); }
#endif

#ifdef F7MOD_Isub_
F7_WEAK void f7_Isub (f7_t *cc, const f7_t *aa) { f7_sub (cc, cc, aa); }
#endif

#ifdef F7MOD_Imul_
F7_WEAK void f7_Imul (f7_t *cc, const f7_t *aa) { f7_mul (cc, cc, aa); }
#endif

#ifdef F7MOD_Idiv_
F7_WEAK void f7_Idiv (f7_t *cc, const f7_t *aa) { f7_div (cc, cc, aa); }
#endif

#ifdef F7MOD_IRsub_
F7_WEAK void f7_IRsub (f7_t *cc, const f7_t *aa) { f7_sub (cc, aa, cc); }
#endif

#ifdef F7MOD_Ineg_
F7_WEAK void f7_Ineg (f7_t *cc) { f7_neg (cc, cc); }
#endif

#ifdef F7MOD_Isqrt_
F7_WEAK void f7_Isqrt (f7_t *cc) { f7_sqrt (cc, cc); }
#endif

#ifdef F7MOD_Isquare_
F7_WEAK void f7_Isquare (f7_t *cc) { f7_square (cc, cc); }
#endif

#ifdef F7MOD_Ildexp_
F7_WEAK f7_t* f7_Ildexp (f7_t *cc, int ex) { return f7_ldexp (cc, cc, ex); }
#endif


#ifdef F7MOD_add_
F7_WEAK
void f7_add (f7_t *cc, const f7_t *aa, const f7_t *bb)
{
  f7_addsub (cc, aa, bb, false);
}
#endif // F7MOD_add_


#ifdef F7MOD_sub_
F7_WEAK
void f7_sub (f7_t *cc, const f7_t *aa, const f7_t *bb)
{
  f7_addsub (cc, aa, bb, true);
}
#endif // F7MOD_sub_


#ifdef F7MOD_addsub_
static void return_with_sign (f7_t *cc, const f7_t *aa, int8_t c_sign)
{
  __asm (";;; return with sign");
  f7_copy (cc, aa);
  if (c_sign != -1)
    f7_set_sign (cc, c_sign);
}

F7_WEAK
void f7_addsub (f7_t *cc, const f7_t *aa, const f7_t *bb, bool neg_b)
{
  uint8_t a_class = f7_classify (aa);
  uint8_t b_class = f7_classify (bb);
  // From this point on, no more access aa->flags or bb->flags
  // to avoid early-clobber when writing cc->flags.

  // Hande NaNs.
  if (f7_class_nan (a_class | b_class))
    return f7_set_nan (cc);

  bool a_sign = f7_class_sign (a_class);
  bool b_sign = f7_class_sign (b_class) ^ neg_b;

  // Add the mantissae?
  bool do_add = a_sign == b_sign;

  // Handle +Infs and -Infs.
  bool a_inf = f7_class_inf (a_class);
  bool b_inf = f7_class_inf (b_class);

  if (a_inf && b_inf)
    {
      if (do_add)
	return f7_set_inf (cc, a_sign);
      else
	return f7_set_nan (cc);
    }
  else if (a_inf)
    return f7_set_inf (cc, a_sign);
  else if (b_inf)
    return f7_set_inf (cc, b_sign);

  int16_t shift16 = sub_ssat16 (aa->expo, bb->expo);

  // aa + 0 = aa.
  // Also check MSBit to get rid of Subnormals and 0.
  if (shift16 > F7_MANT_BITS || f7_is0 (bb))
    return return_with_sign (cc, aa, -1);

  // 0 + bb = bb.
  // 0 - bb = -bb.
  // Also check MSBit to get rid of Subnormals and 0.
  if (shift16 < -F7_MANT_BITS || f7_is0 (aa))
    return return_with_sign (cc, bb, b_sign);

  // Now aa and bb are non-zero, non-NaN, non-Inf.
  // shift > 0 ==> |a| > |b|
  // shift < 0 ==> |a| < |b|
  int8_t shift = (int8_t) shift16;
  bool c_sign = a_sign;
  if (shift < 0
      || (shift == 0 && f7_cmp_mant (aa, bb) < 0))
    {
      const f7_t *p = aa; aa = bb; bb = p;
      c_sign = b_sign;
      shift = -shift;
    }
  uint8_t shift2 = (uint8_t) (shift << 1);

  cc->expo = aa->expo;
  // From this point on, no more access aa->expo or bb->expo
  // to avoid early-clobber when writing cc->expo.

  cc->flags = c_sign;  _Static_assert (F7_FLAGNO_sign == 0, "");

  // This function uses neither .expo nor .flags from either aa or bb,
  // hence there is early-clobber for cc->expo and cc->flags.
  f7_addsub_mant_scaled_asm (cc, aa, bb, shift2 | do_add);
}
#endif // F7MOD_addsub_


#ifdef F7MOD_madd_msub_
F7_WEAK
void f7_madd_msub (f7_t *cc, const f7_t *aa, const f7_t *bb, const f7_t *dd,
                   bool neg_d)
{
  f7_t xx7, *xx = &xx7;
  uint8_t x_lsb = f7_mulx (xx, aa, bb, true /* no rounding */);
  uint8_t x_sign = f7_signbit (xx);
  int16_t x_expo = xx->expo;
  f7_addsub (xx, xx, dd, neg_d);
  // Now add LSB.  If cancellation occured in the add / sub, then we have the
  // chance of extra 8 bits of precision.  Turn LSByte into f7_t.
  f7_clr (cc);
  cc->expo = sub_ssat16 (x_expo, F7_MANT_BITS);
  cc->mant[F7_MANT_BYTES - 1] = x_lsb;
  cc = f7_normalize_asm (cc);
  cc->flags = x_sign;
  f7_Iadd (cc, xx);
}
#endif // F7MOD_madd_msub_

#ifdef F7MOD_madd_
F7_WEAK
void f7_madd (f7_t *cc, const f7_t *aa, const f7_t *bb, const f7_t *dd)
{
  f7_madd_msub (cc, aa, bb, dd, false);
}
#endif // F7MOD_madd_

#ifdef F7MOD_msub_
F7_WEAK
void f7_msub (f7_t *cc, const f7_t *aa, const f7_t *bb, const f7_t *dd)
{
  f7_madd_msub (cc, aa, bb, dd, true);
}
#endif // F7MOD_msub_


#ifdef F7MOD_ldexp_
F7_WEAK
f7_t* f7_ldexp (f7_t *cc, const f7_t *aa, int delta)
{
  uint8_t a_class = f7_classify (aa);

  cc->flags = a_class;

  // Inf and NaN.
  if (! f7_class_number (a_class))
    return cc;

  if (f7_msbit (aa) == 0)
    return (f7_clr (cc), cc);

  int16_t expo = add_ssat16 (delta, aa->expo);
  // Store expo and handle expo = INT16_MIN  and INT16_MAX.
  if (! f7_store_expo (cc, expo))
    f7_copy_mant (cc, aa);

  return cc;
}
#endif // F7MOD_ldexp_


#if USE_LPM
#elif USE_LD
#else
#error need include "asm-defs.h"
#endif // USE_LPM

/*
  Handling constants:

  F7_PCONST (PVAR, X)

      Set  f7_t [const] *PVAR  to an LD address for one
      of the  f7_const_X[_P]  constants.
      PVAR might be set to point to a local auto that serves
      as temporary storage for f7_const_X_P.  PVAR is only
      valid in the current block.

  const f7_t* F7_PCONST_U16 (PVAR, <ident> X)       // USE_LD
  f7_t*       F7_PCONST_U16 (PVAR, uint16_t X)      // USE_LPM

      Set  f7_t [const] *PVAR  to an LD address for one of the
      f7_const_X[_P]  constants.  PVAR might be set to point to a
      local auto that serves as temporary storage for X.  PVAR is
      only valid in the current block.

  F7_PCONST_VAR (PVAR, VAR)

      VAR is a pointer variable holding the address of some f7_const_X[_P]
      constant.  Set  [const] f7_t *PVAR  to a respective LD address.
      PVAR might be set to point to a local auto that serves
      as temporary storage for f7_const_X_P.  PVAR is only
      valid in the current block.

  F7_CONST_ADDR (<ident> CST, f7_t* PTMP)

      Return an LD address to for some f7_const_X[_P] constant.
      *PTMP might be needed to hold a copy of f7_const_X_P in RAM.

  f7_t*       F7_U16_ADDR (uint16_t     X, f7_t* PTMP)   // USE_LPM
  const f7_t* F7_U16_ADDR (<cst-ident>  X, <unused>)     // USE_LD

      Return an LD address to compile-time constant  uint16_t X  which is
      also known as  f7_const_X[_P].  *PTMP might be set to  (f7_t) X.

  f7_t* f7_const (f7_t *PVAR, <cst-ident> X)

      Copy  f7_const_X[_P]  to *PVAR.

  f7_t* f7_copy_flash (f7_t *DST, const f7_t *SRC)

      Copy to *DST with LD (from .rodata in flash) if the address
      space is linear, or with  LPM (from .progmem.data) if the
      address space is not linear.

  f7_t* f7_copy (f7_t *DST, const f7_t* SRC)

      Copy to RAM using LD.

  f7_t* f7_copy_P (f7_t *DST, const f7_t *SRC)

      Copy to RAM using LPM.
*/

#if USE_LPM
  #define F7_RAW_CONST_ADDR(CST) \
      & F7_(const_##CST##_P)

  #define F7_PCONST(PVAR, CST)				    \
      f7_t _var_for_##CST;				    \
      f7_copy_P (& _var_for_##CST, & F7_(const_##CST##_P)); \
      PVAR = & _var_for_##CST

  #define F7_PCONST_U16(PVAR, CST)		\
      f7_t _var_for_##CST;			\
      PVAR = f7_set_u16 (& _var_for_##CST, CST)

  #define F7_PCONST_VAR(PVAR, VAR)		\
      f7_t _var_for_##VAR;			\
      f7_copy_P (& _var_for_##VAR, VAR);	\
      PVAR = & _var_for_##VAR

  #define MAYBE_const // empty

  #define F7_CONST_ADDR(CST, PTMP) \
      f7_copy_P ((PTMP), & F7_(const_##CST##_P))

  #define F7_U16_ADDR(CST, PTMP)   \
      f7_set_u16 ((PTMP), CST)

#elif USE_LD
  #define F7_RAW_CONST_ADDR(CST)   \
      & F7_(const_##CST)

  #define F7_PCONST(PVAR, CST)	   \
      PVAR = & F7_(const_##CST)

  #define F7_PCONST_U16(PVAR, CST) \
      PVAR = & F7_(const_##CST)

  #define F7_PCONST_VAR(PVAR, VAR) \
      PVAR = (VAR)

  #define F7_CONST_ADDR(CST, PTMP) \
      (& F7_(const_##CST))

  #define F7_U16_ADDR(CST, PTMP)   \
      (& F7_(const_##CST))

  #define MAYBE_const const
#endif


#define DD(str, X)		\
  do {				\
    LOG_PSTR (PSTR (str));	\
    f7_dump (X);		\
  } while (0)

#undef DD
#define DD(...) (void) 0


#ifdef F7MOD_sqrt_
static void sqrt_worker (f7_t *cc, const f7_t *rr)
{
  f7_t tmp7, *tmp = &tmp7;
  f7_t aa7, *aa = &aa7;

  // aa in  [1/2, 2)  =>  aa->expo in { -1, 0 }.
  int16_t a_expo = -(rr->expo & 1);
  int16_t c_expo = (rr->expo - a_expo) >> 1;  // FIXME: r_expo = INT_MAX???

  __asm ("" : "+r" (aa));

  f7_copy (aa, rr);
  aa->expo = a_expo;

  // No use of rr or *cc past this point:  We may use cc as temporary.
  // Approximate square-root of  A  by  X <-- (X + A / X) / 2.

  f7_sqrt_approx_asm (cc, aa);

  // Iterate  X <-- (X + A / X) / 2.
  // 3 Iterations with 16, 32, 58 bits of precision for the quotient.

  for (uint8_t prec = 16; (prec & 0x80) == 0; prec <<= 1)
    {
      f7_divx (tmp, aa, cc, (prec & 64) ? 2 + F7_MANT_BITS : prec);
      f7_Iadd (cc, tmp);
      // This will never underflow because |c_expo| is small.
      cc->expo--;
    }

  // Similar: |c_expo| is small, hence no ldexp needed.
  cc->expo += c_expo;
}

F7_WEAK
void f7_sqrt (f7_t *cc, const f7_t *aa)
{
  uint8_t a_class = f7_classify (aa);

  if (f7_class_nan (a_class) || f7_class_sign (a_class))
    return f7_set_nan (cc);

  if (f7_class_inf (a_class))
    return f7_set_inf (cc, 0);

  if (f7_class_zero (a_class))
    return f7_clr (cc);

  sqrt_worker (cc, aa);
}
#endif // F7MOD_sqrt_


#ifdef F7MOD_hypot_
F7_WEAK
void f7_hypot (f7_t *cc, const f7_t *aa, const f7_t *bb)
{
  f7_t xx7, *xx = &xx7;

  f7_square (xx, aa);
  f7_square (cc, bb);
  f7_Iadd (cc, xx);
  f7_Isqrt (cc);
}
#endif // F7MOD_hypot_


#ifdef F7MOD_const_m1_
#include "libf7-constdef.h"
#endif // -1

#ifdef F7MOD_const_1_2_
#include "libf7-constdef.h"
#endif // 1/2

#ifdef F7MOD_const_1_3_
#include "libf7-constdef.h"
#endif // 1/3

#ifdef F7MOD_const_ln2_
#include "libf7-constdef.h"
#endif // ln2

#ifdef F7MOD_const_1_ln2_
#include "libf7-constdef.h"
#endif // 1_ln2

#ifdef F7MOD_const_ln10_
#include "libf7-constdef.h"
#endif // ln10

#ifdef F7MOD_const_1_ln10_
#include "libf7-constdef.h"
#endif // 1_ln10

#ifdef F7MOD_const_1_
#include "libf7-constdef.h"
#endif // 1

#ifdef F7MOD_const_sqrt2_
#include "libf7-constdef.h"
#endif // sqrt2

#ifdef F7MOD_const_2_
#include "libf7-constdef.h"
#endif // 2

#ifdef F7MOD_const_pi_
#include "libf7-constdef.h"
#endif // pi


#ifdef F7MOD_divx_

// C /= A
extern void f7_div_asm (f7_t*, const f7_t*, uint8_t);

F7_WEAK
void f7_divx (f7_t *cc, const f7_t *aa, const f7_t *bb, uint8_t quot_bits)
{
  uint8_t a_class = f7_classify (aa);
  uint8_t b_class = f7_classify (bb);
  // From this point on, no more access aa->flags or bb->flags
  // to avoid early-clobber when writing cc->flags.

  // If either value is NaN, return NaN.
  if (f7_class_nan (a_class | b_class)
      // If both values are Inf or both are 0, return NaN.
      || f7_class_zero (a_class & b_class)
      || f7_class_inf (a_class & b_class)
      // Inf / 0 = NaN.
      || (f7_class_inf (a_class) && f7_class_zero (b_class)))
    {
      return f7_set_nan (cc);
    }

  // 0 / B   = 0  for non-zero, non-NaN B.
  // A / Inf = 0  for non-zero numbers A.
  if (f7_class_zero (a_class) || f7_class_inf (b_class))
    return f7_clr (cc);

  uint8_t c_sign = (a_class ^ b_class) & F7_FLAG_sign;

  if (f7_class_inf (a_class) || f7_class_zero (b_class))
    return f7_set_inf (cc, c_sign);

  cc->flags = c_sign;     _Static_assert (F7_FLAGNO_sign == 0, "");
  int16_t expo = sub_ssat16 (aa->expo, bb->expo);

  // Store expo and handle expo = INT16_MIN  and INT16_MAX.
  if (f7_store_expo (cc, expo))
    return;

  f7_t ss7, *ss = &ss7;
  ss->flags = cc->flags;
  ss->expo  = cc->expo;

  f7_copy_mant (ss, aa);
  f7_div_asm (ss, bb, quot_bits);
  f7_copy (cc, ss);
}
#endif // F7MOD_divx_


#ifdef F7MOD_div_
F7_WEAK
void f7_div (f7_t *cc, const f7_t *aa, const f7_t *bb)
{
  /* When f7_divx calls f7_div_asm, dividend and divisor are valid
     mantissae, i.e. their MSBit is set.  Therefore, the quotient will
     be in  [0x0.ff..., 0x0.40...]  and to adjust it, at most 1 left-shift
     is needed.  Compute F7_MANT_BITS + 2 bits of the quotient:
     One bit is used for rounding, and one bit might be consumed by the
     mentioned left-shift.  */

  f7_divx (cc, aa, bb, 2 + F7_MANT_BITS);
}
#endif // F7MOD_div_


#ifdef F7MOD_div1_
F7_WEAK
void f7_div1 (f7_t *cc, const f7_t *aa)
{
  F7_PCONST_U16 (const f7_t *one, 1);
  f7_div (cc, one, aa);
}
#endif // F7MOD_div_


#ifdef F7MOD_fmod_
F7_WEAK
void f7_fmod (f7_t *cc, const f7_t *aa, const f7_t *bb)
{
  uint8_t a_class = f7_classify (aa);
  uint8_t b_class = f7_classify (bb);

  if (! f7_class_number (a_class)
      || f7_class_nan (b_class)
      || f7_class_zero (b_class))
    {
      return f7_set_nan (cc);
    }

  // A == 0 and B != 0  =>  0.
  if (f7_class_zero (a_class))
    return f7_clr (cc);

  f7_t zz7, *zz = & zz7;

  f7_div (zz, aa, bb);

  // Z in Z,  |Z| <= |A/B|.
  f7_trunc (zz, zz);

  // C = A - Z * B.
  f7_msub (cc, zz, bb, aa);
  cc->flags ^= F7_FLAG_sign;
}
#endif // F7MOD_fmod_


#ifdef F7MOD_truncx_
F7_WEAK
f7_t* f7_truncx (f7_t *cc, const f7_t *aa, bool do_floor)
{
  uint8_t a_class = f7_classify (aa);

  if (! f7_class_nonzero (a_class))
    return f7_copy (cc, aa);

  bool sign = f7_class_sign (a_class);

  int16_t a_expo = aa->expo;

  if (a_expo < 0)
    {
      // |A| < 1.
      if (sign & do_floor)
	return f7_set_s16 (cc, -1);

      f7_clr (cc);
      return cc;
    }
  else if (a_expo >= F7_MANT_BITS - 1)
    // A == floor (A).
    return f7_copy (cc, aa);

  f7_t tmp7, *tmp = &tmp7;

  // Needed if aa === cc.
  f7_copy (tmp, aa);

  cc->flags = sign;
  cc->expo = a_expo;
  f7_clr_mant_lsbs (cc, aa, F7_MANT_BITS - 1 - a_expo);

  if (do_floor && cc->sign && f7_cmp_mant (cc, tmp) != 0)
    {
      F7_PCONST_U16 (const f7_t *one, 1);
      f7_Isub (cc, one);
    }

  return cc;
}
#endif // F7MOD_truncx_


#ifdef F7MOD_floor_
F7_WEAK
f7_t* f7_floor (f7_t *cc, const f7_t *aa)
{
  return f7_truncx (cc, aa, true);
}
#endif // F7MOD_floor_


#ifdef F7MOD_trunc_
F7_WEAK
f7_t* f7_trunc (f7_t *cc, const f7_t *aa)
{
  return f7_truncx (cc, aa, false);
}
#endif // F7MOD_trunc_


#ifdef F7MOD_ceil_
F7_WEAK
void f7_ceil (f7_t *cc, const f7_t *aa)
{
  cc = f7_copy (cc, aa);
  cc->flags ^= F7_FLAG_sign;
  cc = f7_floor (cc, cc);
  f7_Ineg (cc);
}
#endif // F7MOD_ceil_


#ifdef F7MOD_round_
F7_WEAK
void f7_round (f7_t *cc, const f7_t *aa)
{
  f7_t tmp;
  (void) tmp;
  const f7_t *half = F7_CONST_ADDR (1_2, &tmp);

  f7_addsub   (cc, aa, half, f7_signbit (aa));
  f7_trunc (cc, cc);
}
#endif // F7MOD_round_


#ifdef F7MOD_horner_

// Assertion when using this function is that either cc != xx,
// or if cc == xx, then tmp1 must be non-NULL and tmp1 != xx.
// In General, the calling functions have a spare f7_t object available
// and can pass it down to save some stack.
// Moreover, the power series must have degree 1 at least.

F7_WEAK
void f7_horner (f7_t *cc, const f7_t *xx, uint8_t n_coeff, const f7_t *coeff,
                f7_t *tmp1)
{
  f7_assert (n_coeff > 1);

  if (cc != xx)
    tmp1 = cc;
  else
    f7_assert (tmp1 != NULL && tmp1 != xx);

  f7_t *yy = tmp1;
  f7_t tmp27, *tmp2 = &tmp27;

  n_coeff--;
  const f7_t *pcoeff = coeff + n_coeff;

  f7_copy_flash (yy, pcoeff);

  while (1)
    {
      --pcoeff;
#if 1
      f7_Imul (yy, xx);
      const f7_t *cst = USE_LD ? pcoeff : f7_copy_P (tmp2, pcoeff);
      if (coeff == pcoeff)
	return f7_add (cc, yy, cst);
      f7_Iadd (yy, cst);
#else
      const f7_t *cst = USE_LD ? pcoeff : f7_copy_P (tmp2, pcoeff);
      f7_madd (yy, yy, xx, cst);
      if (coeff == pcoeff)
        {
	  f7_copy (cc, yy);
	  return;
        }
#endif
    }

  __builtin_unreachable();
}
#endif // F7MOD_horner_


#ifdef F7MOD_log_
F7_WEAK
void f7_log (f7_t *cc, const f7_t *aa)
{
    f7_logx (cc, aa, NULL);
}
#endif // F7MOD_log_


#ifdef F7MOD_log2_
F7_WEAK
void f7_log2 (f7_t *cc, const f7_t *aa)
{
  f7_logx (cc, aa, F7_RAW_CONST_ADDR (1_ln2));
}
#endif // F7MOD_log2_


#ifdef F7MOD_log10_
F7_WEAK
void f7_log10 (f7_t *cc, const f7_t *aa)
{
  f7_logx (cc, aa, F7_RAW_CONST_ADDR (1_ln10));
}
#endif // F7MOD_log10_


#ifdef F7MOD_logx_

#define ARRAY_NAME coeff_artanh
#include "libf7-array.def"
#undef ARRAY_NAME

// Compute P * ln(A)  if P != NULL and ln(A), otherwise.
// P is a LD-address if USE_LD and a LPM-address if USE_LPM.
// Assumption is that P > 0.

F7_WEAK
void f7_logx (f7_t *cc, const f7_t *aa, const f7_t *p)
{
  uint8_t a_class = f7_classify (aa);

  if (f7_class_nan (a_class) || f7_class_sign (a_class))
    return f7_set_nan (cc);

  if (f7_class_inf (a_class))
    return f7_set_inf (cc, 0);

  if (f7_class_zero (a_class))
    return f7_set_inf (cc, 1);

  f7_t *yy = cc;
  f7_t xx7, *xx = &xx7;
  f7_t tmp7, *tmp = &tmp7;

  // Y in [1, 2]  =  A * 2 ^ (-a_expo).
  int16_t a_expo = aa->expo;
  f7_copy (yy, aa);
  yy->expo = 0;

  // Y in [1 / sqrt2, sqrt2].

  if (f7_abscmp_msb_ge (yy, F7_(const_sqrt2_msb), F7_(const_sqrt2_expo)))
    {
      yy->expo = -1;
      a_expo = add_ssat16 (a_expo, 1);
    }

  const f7_t *one = F7_U16_ADDR (1, & tmp7);

  // X := (Y - 1) / (Y + 1),  |X| <= (sqrt2 - 1) / (sqrt2 + 1)  ~  0.172.
  f7_sub (xx, yy, one);
  f7_Iadd (yy, one);
  f7_Idiv (xx, yy);

  // Y := X^2,  |Y| < 0.03.
  f7_square (yy, xx);

  // Y := artanh (X^2) / X
  f7_horner (yy, yy, n_coeff_artanh, coeff_artanh, tmp);

  // C = X * Y = ln A - a_expo * ln2.
  f7_mul (cc, xx, yy);

  // X := a_expo * ln2.
  f7_set_s16 (xx, a_expo);
  f7_Imul (xx, F7_CONST_ADDR (ln2, & tmp7));

  // C = ln A.
  f7_Iadd (cc, xx);

  if (p && USE_LPM)
    f7_Imul (cc, f7_copy_P (tmp, p));
  if (p && USE_LD)
    f7_Imul (cc, p);
}
#endif // F7MOD_logx_


#ifdef F7MOD_exp_

#define ARRAY_NAME coeff_exp
#include "libf7-array.def"
#undef ARRAY_NAME

#define STATIC static
#include "libf7-constdef.h" // ln2_low
#undef STATIC

F7_WEAK
void f7_exp (f7_t *cc, const f7_t *aa)
{
  uint8_t a_class = f7_classify (aa);

  if (f7_class_nan (a_class))
    return f7_set_nan (cc);

  /* The maximal exponent of 2 for a double is 1023, hence we may limit
     to  |A| < 1023 * ln2 ~ 709.  We limit to  1024 ~ 1.99 * 2^9  */

  if (f7_class_inf (a_class)
      || (f7_class_nonzero (a_class) && aa->expo >= 9))
    {
      if (f7_class_sign (a_class))
	return f7_clr (cc);
      else
	return f7_set_inf (cc, 0);
    }

  f7_t const *cst;
  f7_t qq7, *qq = &qq7;

  F7_PCONST (cst, ln2);

  // We limited |A| to 1024 and are now dividing by ln2, hence Q will
  // be at most 1024 / ln2 ~ 1477 and fit into 11 bits.  We will
  // round Q anyway, hence only request 11 bits from the division and
  // one additional bit that might be needed to normalize the quotient.
  f7_divx (qq, aa, cst, 1 + 11);

  // Use the smallest (by absolute value) remainder system.
  f7_round (qq, qq);
  int16_t q = f7_get_s16 (qq);

  // Reducing A mod ln2 gives |C| <= ln2 / 2,  C = -A mod ln2.
  f7_msub (cc, qq, cst, aa);

  // Corrigendum:  We added Q * ln2; now add Q times the low part of ln2
  // for better precision.  Due to |C| < |A| this is not a no-op in general.
  const f7_t *yy = F7_CONST_ADDR (ln2_low, &_var_for_ln2);
  f7_madd (cc, qq, yy, cc);

  // Because we computed C = -A mod ...
  cc->flags ^= F7_FLAG_sign;

  // Reduce further to |C| < ln2 / 8 which is the range of our MiniMax poly.
  const uint8_t MAX_LN2_RED = 3;
  int8_t scal2 = 0;

  while (f7_abscmp_msb_ge (cc, F7_(const_ln2_msb),
			   F7_(const_ln2_expo) - MAX_LN2_RED))
    {
      scal2++;
      cc->expo--;
    }

  f7_horner (cc, cc, n_coeff_exp, coeff_exp, qq);

  while (--scal2 >= 0)
    f7_Isquare (cc);

  f7_Ildexp (cc, q);
}
#endif // F7MOD_exp_


#ifdef F7MOD_pow10_
F7_WEAK
void f7_pow10 (f7_t *cc, const f7_t *aa)
{
  const f7_t *p_ln10;
  F7_PCONST (p_ln10, ln10);
  f7_mul (cc, aa, p_ln10);
  f7_exp (cc, cc);
}
ALIAS (f7_pow10, f7_exp10)
#endif // F7MOD_pow10_


#ifdef F7MOD_cbrt_
F7_WEAK
void f7_cbrt (f7_t *cc, const f7_t *aa)
{
  f7_copy (cc, aa);
  const f7_t *p_1_3;
  uint8_t c_flags = cc->flags;
  cc->flags &= ~F7_FLAG_sign;
  f7_log (cc, cc);
  F7_PCONST (p_1_3, 1_3);
  f7_Imul (cc, p_1_3);
  f7_exp (cc, cc);

  if (c_flags & F7_FLAG_sign)
    cc->flags |= F7_FLAG_sign;
}
#endif // F7MOD_cbrt_


#ifdef F7MOD_pow_
F7_WEAK
void f7_pow (f7_t *cc, const f7_t *aa, const f7_t *bb)
{
#if 0
  f7_t slots[cc == bb];
  f7_t *yy = cc == bb ? slots : cc;
#else
  f7_t yy7, *yy = &yy7;
#endif
  f7_log (yy, aa);
  f7_Imul (yy, bb);
  f7_exp (cc, yy);
}
#endif // F7MOD_pow_


#ifdef F7MOD_powi_
F7_WEAK
void f7_powi (f7_t *cc, const f7_t *aa, int ii)
{
  uint16_t u16 = ii;
  f7_t xx27, *xx2 = &xx27;

  if (ii < 0)
    u16 = -u16;

  f7_copy (xx2, aa);

  f7_set_u16 (cc, 1);

  while (1)
    {
      if (u16 & 1)
	f7_Imul (cc, xx2);

      if (! f7_is_nonzero (cc))
	break;

      u16 >>= 1;
      if (u16 == 0)
	break;
      f7_Isquare (xx2);
    }

  if (ii < 0)
    f7_div1 (xx2, aa);
}
#endif // F7MOD_powi_


#ifdef F7MOD_sincos_

#define ARRAY_NAME coeff_sin
  #define FOR_SIN
  #include "libf7-array.def"
  #undef  FOR_SIN
#undef ARRAY_NAME

#define ARRAY_NAME coeff_cos
  #define FOR_COS
  #include "libf7-array.def"
  #undef  FOR_COS
#undef ARRAY_NAME

#define STATIC static
#include "libf7-constdef.h" // pi_low
#undef STATIC

typedef union
{
  struct
  {
    bool    neg_sin : 1; // Must be bit F7_FLAGNO_sign.
    bool    neg_cos : 1;
    bool    do_sin: 1;
    bool    do_cos: 1;
    bool    swap_sincos : 1;
    uint8_t res : 3;
  };
  uint8_t bits;
} sincos_t;


F7_WEAK
void f7_sincos (f7_t *ss, f7_t *cc, const f7_t *aa)
{
  uint8_t a_class = f7_classify (aa);

  sincos_t sc = { .bits = a_class & F7_FLAG_sign };
  if (ss != NULL) sc.do_sin = 1;
  if (cc != NULL) sc.do_cos = 1;

  if (f7_class_nan (a_class) || f7_class_inf (a_class))
    {
      if (sc.do_sin)  f7_set_nan (ss);
      if (sc.do_cos)  f7_set_nan (cc);
      return;
    }

  f7_t pi7, *pi = &pi7;
  f7_t xx7, *xx = &xx7;
  f7_t yy7, *yy = &yy7;
  f7_t *hh = sc.do_sin ? ss : cc;

  // X = |A|
  f7_copy (xx, aa);
  xx->flags = 0;

  // Y is how often we subtract PI from X.
  f7_clr (yy);
  f7_const (pi, pi);

  if (f7_abscmp_msb_ge (xx, F7_(const_pi_msb), F7_(const_pi_expo) + 1))
    {
      pi->expo = 1 + F7_(const_pi_expo);  // 2*pi

      // Y = X / 2pi.
      f7_div (yy, xx, pi);

      // The integral part of |A| / pi mod 2 is bit 55 - x_expo.
      if (yy->expo >= F7_MANT_BITS && !f7_is_zero (yy))
        {
	  // Too big for sensible calculation:  Should this be NaN instead?
	  if (sc.do_sin)  f7_clr (ss);
	  if (sc.do_cos)  f7_clr (cc);
	  return;
        }

      // X -= 2pi * [ X / 2pi ]
      f7_floor (yy, yy);

      f7_msub (xx, yy, pi, xx);
      xx->flags ^= F7_FLAG_sign;

      // We divided by 2pi, but Y should count times we subtracted pi.
      yy->expo++;
    }

  pi->expo = F7_(const_pi_expo); // pi
  f7_sub (hh, xx, pi);
  if (!f7_signbit (hh))
    {
      // H = X - pi >= 0  =>  X >= pi
      // sin(x) = -sin(x - pi)
      // cos(x) = -cos(x - pi)
      f7_copy (xx, hh);
      // Y: We subtracted pi one more time.
      f7_Iadd (yy, f7_set_u16 (hh, 1));
      sc.neg_sin ^= 1;
      sc.neg_cos ^= 1;
    }

  pi->expo = F7_(const_pi_expo) - 1; // pi/2
  if (f7_gt (xx, pi))
    {
      // x > pi/2
      // sin(x) =  sin(pi - x)
      // cos(x) = -cos(pi - x)
      pi->expo = F7_(const_pi_expo); // pi
      f7_IRsub (xx, pi);
      // Y: We subtracted pi one more time (and then negated).
      f7_Iadd (yy, f7_set_u16 (hh, 1));
      yy->flags ^= F7_FLAG_sign;
      sc.neg_cos ^= 1;
    }

  pi->expo = F7_(const_pi_expo) - 2; // pi/4
  if (f7_gt (xx, pi))
    {
      // x > pi/4
      // sin(x) = cos(pi/2 - x)
      // cos(x) = sin(pi/2 - x)
      pi->expo = F7_(const_pi_expo) - 1; // pi/2
      f7_IRsub (xx, pi);
      // Y: We subtracted pi/2 one more time (and then negated).
      f7_Iadd (yy, f7_set_1pow2 (hh, -1, 0));
      yy->flags ^= F7_FLAG_sign;
      sc.swap_sincos = 1;
    }

  if (!f7_is0 (yy))
    {
      // Y counts how often we subtracted pi from X in order to
      // get 0 <= X < pi/4 as small as possible (Y is 0 mod 0.5).
      // Now also subtract the low part of pi:
      // f7_const_pi_low = pi - f7_const_pi  in order to get more precise
      // results in the cases where the final result is close to 0.
      const f7_t *pi_low = F7_CONST_ADDR (pi_low, pi);
      //f7_const (pi, pi_low);
      f7_Imul (yy, pi_low);
      f7_Isub (xx, yy);
    }

  // X   in [0, pi/4].
  // X^2 in [0, pi^2/16] ~ [0, 0.6169]

  f7_square (yy, xx);

  f7_t *x_sin = xx;
  f7_t *x_cos = yy;

  if ((sc.do_sin && !sc.swap_sincos)
      || (sc.do_cos && sc.swap_sincos))
    {
      f7_horner (hh, yy, n_coeff_sin, coeff_sin, NULL);
      f7_mul (x_sin, hh, xx);
    }

  if ((sc.do_cos && !sc.swap_sincos)
      || (sc.do_sin && sc.swap_sincos))
    {
      f7_horner (x_cos, yy, n_coeff_cos, coeff_cos, hh);
    }

  if (sc.swap_sincos)
    {
      x_sin = yy;
      x_cos = xx;
    }

  if (sc.do_sin)
    {
      x_sin->flags ^= sc.bits;
      x_sin->flags &= F7_FLAG_sign;
      f7_copy (ss, x_sin);
    }

  if (sc.do_cos)
    {
      x_cos->flags = sc.neg_cos;
      f7_copy (cc, x_cos);
    }
}
#endif // F7MOD_sincos_

#ifdef F7MOD_sin_
F7_WEAK
void f7_sin (f7_t *ss, const f7_t *aa)
{
  f7_sincos (ss, NULL, aa);
}
#endif // F7MOD_sin_

#ifdef F7MOD_cos_
F7_WEAK
void f7_cos (f7_t *cc, const f7_t *aa)
{
  f7_sincos (NULL, cc, aa);
}
#endif // F7MOD_cos_


#ifdef F7MOD_tan_
F7_WEAK
void f7_tan (f7_t *tt, const f7_t *aa)
{
  f7_t xcos;
  f7_sincos (tt, & xcos, aa);
  f7_Idiv (tt, & xcos);
}
#endif // F7MOD_tan_


#ifdef F7MOD_cotan_
F7_WEAK
void f7_cotan (f7_t *ct, const f7_t *aa)
{
  f7_t xcos;
  f7_sincos (ct, & xcos, aa);
  f7_div (ct, & xcos, ct);
}
#endif // F7MOD_cotan_


#ifdef F7MOD_sinhcosh_
F7_WEAK
void f7_sinhcosh (f7_t *cc, const f7_t *aa, bool do_sinh)
{
  f7_t xx7, *xx = &xx7;
  // C = exp(A)
  f7_exp (cc, aa);
  // X = exp(-A)
  f7_div (xx, f7_set_u16 (xx, 1), cc);
  // sinh(A) = (exp(A) - exp(-A)) / 2
  // cosh(A) = (exp(A) + exp(-A)) / 2
  f7_addsub (cc, cc, xx, do_sinh);
  cc->expo--;
}
#endif // F7MOD_sinhcosh_


#ifdef F7MOD_sinh_
F7_WEAK
void f7_sinh (f7_t *cc, const f7_t *aa)
{
  f7_sinhcosh (cc, aa, true);
}
#endif // F7MOD_sinh_


#ifdef F7MOD_cosh_
F7_WEAK
void f7_cosh (f7_t *cc, const f7_t *aa)
{
  f7_sinhcosh (cc, aa, false);
}
#endif // F7MOD_cosh_


#ifdef F7MOD_tanh_
F7_WEAK
void f7_tanh (f7_t *cc, const f7_t *aa)
{
  // tanh(A) = (exp(2A) - 1) / (exp(2A) + 1)
  f7_t xx7, *xx = &xx7;
  F7_PCONST_U16 (const f7_t *one, 1);
  // C = 2A
  f7_copy (cc, aa);
  cc->expo++;
  // C = exp(2A)
  f7_exp (cc, cc);
  // X = exp(2A) + 1
  f7_add (xx, cc, one);
  // C = exp(2A) - 1
  f7_Isub (cc, one);
  // C = tanh(A)
  f7_Idiv (cc, xx);
}
#endif // F7MOD_tanh_


#ifdef F7MOD_atan_

#define MINIMAX_6_6_IN_0_1

#define ARRAY_NAME coeff_atan_zahler
#define FOR_NUMERATOR
#include "libf7-array.def"
#undef FOR_NUMERATOR
#undef ARRAY_NAME

#define ARRAY_NAME coeff_atan_nenner
#define FOR_DENOMINATOR
#include "libf7-array.def"
#undef FOR_DENOMINATOR
#undef ARRAY_NAME

#include "libf7-constdef.h"

F7_WEAK
void f7_atan (f7_t *cc, const f7_t *aa)
{
  uint8_t a_class = f7_classify (aa);
  uint8_t flags = a_class & F7_FLAG_sign;

  if (f7_class_nan (a_class))
    return f7_set_nan (cc);

  f7_t yy7, *yy = &yy7;
  f7_t zz7, *zz = &zz7;

  if (f7_class_inf (a_class))
    {
      f7_set_u16 (cc, 0);
      goto do_Inf;
    }

  // C = |A|
  f7_copy (cc, aa);
  cc->flags = 0;

  if (!f7_class_zero (a_class) && cc->expo >= 0)
    {
      // C >= 1:  use  atan (x) + atan (1/x) = pi / 2  to reduce to [0, 1].
      flags |= 1 << 1;
      f7_div (cc, f7_set_u16 (yy, 1), cc);
    }
#if !defined (MINIMAX_6_6_IN_0_1)
  const uint8_t const_a_msb = 0x89;
  const int16_t const_a_expo = -2;
  if (f7_abscmp_msb_ge (cc, const_a_msb, const_a_expo))
    {
      // We have C in [0, 1] and we want to use argument reduction by means
      // of addition theorem  atan(x) - atan(y) = atan((x - y) / (1 + xy)).
      // We want to split [0, 1] into  [0, a] u [a, 1]  in such a way that
      // the upper interval will be mapped to  [-a, a].  The system is easy
      // to solve and yiels
      //    y = 1 / sqrt (3)       ~  0.57735     atan(y) = pi / 6
      //    a = (1 - y) / (1 + y)  ~  0.26795  ~  0x0.8930A2F * 2^-1.
      flags |= 1 << 2;
      // C <- (C - Y) / (1 + C * Y)  in  [-a, a]
      const f7_t *cst = F7_CONST_ADDR (1_sqrt3, zz);
      f7_mul (yy, cc, cst);
      f7_Isub (cc, cst);
      f7_Iadd (yy, F7_U16_ADDR (1, zz));
      f7_Idiv (cc, yy);
    }
#endif
  // C <- C * p(C^2) / q(C^2)
  f7_square (yy, cc);
  f7_horner (zz, yy, n_coeff_atan_zahler, coeff_atan_zahler, NULL);
  f7_Imul (zz, cc);
  f7_horner (cc, yy, n_coeff_atan_nenner, coeff_atan_nenner, NULL);
  f7_div (cc, zz, cc);

#if !defined (MINIMAX_6_6_IN_0_1)
  if (flags & (1 << 2))
    f7_Iadd (cc, F7_CONST_ADDR (pi_6, yy));
#endif

  if (flags & (1 << 1))
    {
    do_Inf:;
      // Y = pi / 2
      f7_const (yy, pi);
      yy->expo = F7_(const_pi_expo) - 1;
      f7_IRsub (cc, yy);
    }

  cc->flags = a_class & F7_FLAG_sign;
}
#undef MINIMAX_6_6_IN_0_1
#endif // F7MOD_atan_


#ifdef F7MOD_asinacos_

#define ARRAY_NAME coeff_func_a_zahler
#define FOR_NUMERATOR
#include "libf7-array.def"
#undef  FOR_NUMERATOR
#undef  ARRAY_NAME

#define ARRAY_NAME coeff_func_a_nenner
#define FOR_DENOMINATOR
#include "libf7-array.def"
#undef  FOR_DENOMINATOR
#undef  ARRAY_NAME

typedef union
{
  struct
  {
    bool    sign : 1;       // Must be bit F7_FLAGNO_sign.
    bool    do_acos : 1;    // From caller.
    bool    have_acos : 1;  // What we compute from rational approx p/q.
    uint8_t res : 5;
  };
  uint8_t bits;
} asinacos_t;

F7_WEAK
void f7_asinacos (f7_t *cc, const f7_t *aa, uint8_t what)
{
  f7_t xx7, *xx = &xx7;
  f7_t xx27, *xx2 = &xx27;

  asinacos_t flags = { .bits = what | f7_signbit (aa) };

  f7_abs (xx, aa);

  int8_t cmp = f7_cmp (xx, f7_set_u16 (cc, 1));

  if (cmp == INT8_MIN
      || cmp > 0)
    {
      return f7_set_nan (cc);
    }

  if (xx->expo <= -2 || f7_is_zero (xx))
    {
      // |A| < 1/2:  asin(x) = x * a(2*x^2)
      f7_square (xx2, xx);
      xx2->expo ++;
    }
  else
    {
      // |A| > 1/2: acos (1-x) = sqrt(2*x) * a(x)
      // C is 1 from above.
      f7_IRsub (xx, cc);
      f7_copy (xx2, xx);
      flags.have_acos = 1;
    }

  // MiniMax [5/4] numerator.
  f7_horner (cc, xx2, n_coeff_func_a_zahler, coeff_func_a_zahler, NULL);

  if (flags.have_acos)
    {
      xx->expo ++;
      f7_Isqrt (xx);
    }
  f7_Imul (cc, xx);

  // MiniMax [5/4] denominator.
  f7_horner (xx, xx2, n_coeff_func_a_nenner, coeff_func_a_nenner, NULL);

  f7_Idiv (cc, xx);

  /*
      With the current value of C, we have:

		|	     |	      do_asin	    |	    do_acos
		|     C	     |	A <= 0	 |  A >= 0  |  A <= 0  |  A >= 0
      ----------+------------+-----------+----------+----------+----------
      have_asin | asin (|A|) |	  -C	 |    C	    | pi/2 + C | pi/2 - C
      have_acos | acos (|A|) | -pi/2 + C | pi/2 - C |  pi - C  |     C

      Result = n_pi2 * pi/2 + C * (c_sign ? -1 : 1)
      Result (A, do_asin) = asin (A)
      Result (A, do_acos) = acos (A)

      with
	  c_sign = do_acos ^ have_acos ^ a_sign
	  n_pi2	 = do_acos + have_acos * (a_sign ^ do_acos) ? (-1 : 1)
	  n_pi2 in { -1, 0, 1, 2 }
  */

  // All what matters for c_sign is bit 0.
  uint8_t c_sign = flags.bits;
  int8_t n_pi2 = flags.do_acos;
  c_sign ^= flags.do_acos;
  if (flags.have_acos)
    {
      n_pi2++;
      __asm ("" : "+r" (n_pi2));
      if (c_sign & 1)  // c_sign & 1  =  a_sign ^ do_acos
	n_pi2 -= 2;
      c_sign++;
    }

  cc->flags = c_sign & F7_FLAG_sign;

  if (n_pi2)
    {
      f7_const (xx, pi);
      if (n_pi2 < 0)
	xx->sign = 1;
      if (n_pi2 != 2)
	xx->expo = F7_(const_pi_expo) - 1;

      f7_Iadd (cc, xx);
    }
}
#endif // F7MOD_asinacos_


#ifdef F7MOD_asin_
F7_WEAK
void f7_asin (f7_t *cc, const f7_t *aa)
{
  f7_asinacos (cc, aa, 0);
}
#endif // F7MOD_asin_


#ifdef F7MOD_acos_
F7_WEAK
void f7_acos (f7_t *cc, const f7_t *aa)
{
  f7_asinacos (cc, aa, 1 << 1);
}
#endif // F7MOD_acos_


#ifndef IN_LIBGCC2

#ifdef F7MOD_put_C_

#include <stdio.h>
#include <avr/pgmspace.h>

static F7_INLINE
uint8_t f7_hex_digit (uint8_t nibble)
{
  nibble = (uint8_t) (nibble + '0');
  if (nibble > '9')
    nibble = (uint8_t) (nibble + ('a' - '0' - 10));
  return nibble;
}

static void f7_put_hex2 (uint8_t x, FILE *stream)
{
  putc ('0', stream);
  if (x)
    {
      putc ('x', stream);
      putc (f7_hex_digit (x >> 4), stream);
      putc (f7_hex_digit (x & 0xf), stream);
    }
}

#define XPUT(str) \
  fputs_P (PSTR (str), stream)

// Write to STREAM a line that is appropriate for usage in libf7-const.def.

F7_WEAK
void f7_put_CDEF (const char *name, const f7_t *aa, FILE *stream)
{
  char buf[7];
  XPUT ("F7_CONST_DEF (");
  fputs (name, stream);
  XPUT (",\t");
  uint8_t a_class = f7_classify (aa);
  if (! f7_class_nonzero (a_class))
    {
      f7_put_hex2 (a_class & F7_FLAGS, stream);
      XPUT (",\t0,0,0,0,0,0,0,\t0)");
      return;
    }
  putc ('0' + (a_class & F7_FLAGS), stream);
  XPUT (",\t");

  for (uint8_t i = 0; i < F7_MANT_BYTES; i++)
    {
      f7_put_hex2 (aa->mant[F7_MANT_BYTES-1 - i], stream);
      putc (',', stream);
    }
  putc ('\t', stream);

  itoa (aa->expo, buf, 10);
  fputs (buf, stream);
  XPUT (")");
}

void f7_put_C (const f7_t *aa, FILE *stream)
{
  char buf[7];

  uint8_t a_class = f7_classify (aa);
  if (f7_class_nan (a_class))
    {
      XPUT ("{ .is_nan = 1 }");
      return;
    }
  bool sign = a_class & F7_FLAG_sign;

  if (f7_class_inf (a_class))
    {
      XPUT ("{ .is_nan = 1, .sign = ");
      putc ('0' + sign, stream);
      XPUT (" }");
      return;
    }

  XPUT ("{ .sign = ");
  putc ('0' + sign, stream);

  XPUT (", .mant = { ");
  for (uint8_t i = 0; i < F7_MANT_BYTES; i++)
    {
      f7_put_hex2 (aa->mant[F7_MANT_BYTES-1 - i], stream);
      if (i != F7_MANT_BYTES - 1)
	putc (',', stream);
    }

  XPUT (" }, .expo = ");
  itoa (aa->expo, buf, 10);
  fputs (buf, stream);
  XPUT (" }");
}
#endif //F7MOD_put_C_


#ifdef F7MOD_dump_

#include <avr/pgmspace.h>

#ifndef AVRTEST_H

#include <stdio.h>

static void LOG_PSTR (const char *str)
{
  fputs_P (str, stdout);
}

static void LOG_PFMT_U16 (const char *fmt, uint16_t x)
{
  printf_P (fmt, x);
}

static void LOG_PFMT_FLOAT (const char *fmt, float x)
{
  printf_P (fmt, x);
}

#define LOG_X8(X)               LOG_PFMT_U16 (PSTR (" 0x%02x "), (uint8_t)(X))
#define LOG_PFMT_S16(FMT, X)    LOG_PFMT_U16 (FMT, (unsigned)(X))
#define LOG_PFMT_ADDR(FMT, X)   LOG_PFMT_U16 (FMT, (unsigned)(X))

#endif // AVRTEST_H

static void dump_byte (uint8_t b)
{
  LOG_PSTR (PSTR (" "));
  for (uint8_t i = 0; i < 8; i++)
    {
      LOG_PSTR ((b & 0x80) ? PSTR ("1") : PSTR ("0"));
      b = (uint8_t) (b << 1);
    }
}

void f7_dump_mant (const f7_t *aa)
{
  LOG_PSTR (PSTR ("\tmant   ="));
  for (int i = F7_MANT_BYTES - 1; i >= 0; i--)
    LOG_X8 (aa->mant[i]);
  LOG_PSTR (PSTR ("\n\t       ="));

  for (int i = F7_MANT_BYTES - 1; i >= 0; i--)
    dump_byte (aa->mant[i]);
  LOG_PSTR (PSTR ("\n"));
}

void f7_dump (const f7_t *aa)
{
  LOG_PFMT_ADDR (PSTR ("\n0x%04x\tflags  = "), aa);
  dump_byte (aa->flags);
  uint8_t a_class = f7_classify (aa);
  LOG_PSTR (PSTR ("  = "));
  LOG_PSTR (f7_class_sign (a_class) ? PSTR ("-") : PSTR ("+"));
  if (f7_class_inf (a_class))    LOG_PSTR (PSTR ("Inf "));
  if (f7_class_nan (a_class))    LOG_PSTR (PSTR ("NaN "));
  if (f7_class_zero (a_class))   LOG_PSTR (PSTR ("0 "));
  if (f7_class_number (a_class)) LOG_PSTR (PSTR ("Number "));

  LOG_PFMT_FLOAT (PSTR (" = %.10g\n"), f7_get_float (aa));
  LOG_PFMT_S16 (PSTR ("\texpo   = %d\n"), aa->expo);

  f7_dump_mant (aa);
}
#endif // F7MOD_dump_

#endif // ! libgcc

#endif // !AVR_TINY