/Users/mourrain/Devel/mmx/numerix/include/numerix/tangent.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : tangent.hpp
* DESCRIPTION: Efficient computations in first order jet spaces C[X]/(X^2)
* COPYRIGHT  : (C) 2006  Joris van der Hoeven
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#ifndef __MMX_TANGENT_HPP
#define __MMX_TANGENT_HPP
#include <numerix/rational.hpp>
#include <basix/function.hpp>
namespace mmx {
#define TMPL_DEF template<typename C, typename D=C>
#define TMPL template<typename C,typename D>
#define Tangent tangent<C,D>
#define Abs_tangent tangent<Abs_type(C),Abs_type(D) >
#define Real_tangent tangent<Real_type(C),Real_type(D) >
TMPL class tangent;
TMPL inline C base  (const Tangent& z);
TMPL inline D slope (const Tangent& z);
TMPL inline C& base  (Tangent& z);
TMPL inline D& slope (Tangent& z);
TMPL inline Tangent operator - (const Tangent& z1, const Tangent& z2);

/******************************************************************************
* The tangent type
******************************************************************************/

TMPL_DEF
class tangent {
MMX_ALLOCATORS
private:
  C c;
  D d;

public:
  inline tangent (): c (0), d (0) {}
  inline tangent (const format<C>& fm1, const format<D>& fm2):
    c (promote (0, fm1)), d (promote (0, fm2)) {}
  template<typename T> tangent (const T& c2): c (c2), d (0) {}
  template<typename CT,typename DT> inline tangent (const tangent<CT,DT>& z):
    c (as<C> (base (z))), d (as<D> (slope (z))) {}
  inline tangent (const C& z): c (copy (z)), d (0) {}
  inline tangent (const C& c2, const D& d2): c (c2), d (d2) {}
  friend C base LESSGTR (const Tangent& z);
  friend D slope LESSGTR (const Tangent& z);
  friend C& base LESSGTR (Tangent& z);
  friend D& slope LESSGTR (Tangent& z);
};
DEFINE_BINARY_FORMAT_2 (tangent)

TMPL struct rounding_helper<Tangent >: public rounding_helper<C> {};

STYPE_TO_TYPE(TMPL,scalar_type,Tangent,C);
UNARY_RETURN_TYPE(TMPL,Re_op,Tangent,Real_tangent);
UNARY_RETURN_TYPE(TMPL,abs_op,Tangent,Abs_tangent);

/******************************************************************************
* Basic routines
******************************************************************************/

TMPL inline void scalar_tangent (Tangent& r, const C& c) {
  r= Tangent (c, promote (0, slope (r))); }
TMPL inline C base (const Tangent& z) { return z.c; }
TMPL inline D slope (const Tangent& z) { return z.d; }
TMPL inline C& base (Tangent& z) { return z.c; }
TMPL inline D& slope (Tangent& z) { return z.d; }
TMPL inline format<C> CF1 (const Tangent& z) { return get_format (base (z)); }
TMPL inline format<D> CF2 (const Tangent& z) { return get_format (slope (z)); }

TMPL inline Tangent copy (const Tangent& z) {
  return Tangent (copy (base (z)), copy (slope (z))); }
TMPL inline Tangent duplicate (const Tangent& z) {
  return Tangent (duplicate (base (z)), duplicate (slope (z))); }

template<typename Op, typename C, typename D> nat
unary_hash (const Tangent& x) {
  nat h= Op::op (base (x));
  return (h<<3) ^ (h<<11) ^ (h>>29) ^ Op::op (slope (x));
}

template<typename Op, typename C, typename D> bool
binary_test (const Tangent& x1, const Tangent& x2) {
  return Op::op (base (x1), base (x2)) &&
         Op::op (slope (x1), slope (x2));
}

template<typename CS, typename CD, typename DS, typename DD> tangent<CD,DD>
map (const function_1<CD,Argument(CS) >& funC,
     const function_1<DD,Argument(DS) >& funD,
     const tangent<CS,DS>& z,
     const format<CD>& fmC,
     const format<DD>& fmD) {
  return tangent<CD,DD> (funC (base (z)), funD (slope (z)));
}

TRUE_IDENTITY_OP_SUGAR(TMPL,Tangent)
EXACT_IDENTITY_OP_SUGAR(TMPL,Tangent)
HARD_IDENTITY_OP_SUGAR(TMPL,Tangent)

TMPL syntactic
flatten (const Tangent& z) {
  return apply ("tangent", flatten (base (z)), flatten (slope (z)));
}

TMPL
struct binary_helper<Tangent >: public void_binary_helper<Tangent > {
  static inline string short_type_name () {
    return "Ta" * Short_type_name (C) * Short_type_name (D); }
  static inline generic full_type_name () {
    return gen ("Tangent", Full_type_name (C)); }
  static inline generic disassemble (const Tangent& z) {
    return gen_vec (as<generic> (base (z)), as<generic> (slope (z))); }
  static inline Tangent assemble (const generic& v) {
    return Tangent (as<C> (vector_access (v, 0)),
                    as<D> (vector_access (v, 1))); }
  static inline void write (const port& out, const Tangent& z) {
    binary_write<C> (out, base (z));
    binary_write<D> (out, slope (z)); }
  static inline Tangent read (const port& in) {
    C b= binary_read<C> (in);
    D s= binary_read<D> (in);
    return Tangent (b, s); }
};

/******************************************************************************
* Constants
******************************************************************************/

TMPL inline void set_nan (Tangent& z) {
  set_nan (base (z)); set_nan (slope (z)); }
TMPL inline void set_infinity (Tangent& z) {
  set_infinity (base (z)); set_zero (slope (z)); }
TMPL inline void set_fuzz (Tangent& z) {
  set_fuzz (base (z)); set_fuzz (slope (z)); }
TMPL inline void set_smallest (Tangent& z) {
  set_smallest (base (z)); set_zero (slope (z)); }
TMPL inline void set_largest (Tangent& z) {
  set_largest (base (z)); set_zero (slope (z)); }
TMPL inline void set_accuracy (Tangent& z) {
  set_accuracy (base (z)); set_zero (slope (z)); }
TMPL inline void set_log2 (Tangent& z) {
  set_log2 (base (z)); set_zero (slope (z)); }
TMPL inline void set_pi (Tangent& z) {
  set_pi (base (z)); set_zero (slope (z)); }
TMPL inline void set_euler (Tangent& z) {
  set_euler (base (z)); set_zero (slope (z)); }
TMPL inline void set_imaginary (Tangent& z) {
  set_imaginary (base (z)); set_zero (slope (z)); }
TMPL inline Tangent times_infinity (const Tangent& x) {
  return x * infinity_cst<Tangent > (); }

/******************************************************************************
* Basic arithmetic
******************************************************************************/

TMPL inline Tangent
operator - (const Tangent& z) {
  return Tangent (-base (z), -slope (z));
}

TMPL inline Tangent
operator + (const Tangent& z1, const Tangent& z2) {
  return Tangent (base (z1) + base (z2), slope (z1) + slope (z2));
}

TMPL inline Tangent
operator - (const Tangent& z1, const Tangent& z2) {
  return Tangent (base (z1) - base (z2), slope (z1) - slope (z2));
}

TMPL inline Tangent
operator * (const Tangent& z1, const Tangent& z2) {
  return Tangent (base (z1) * base (z2),
                  base (z1) * slope (z2) + slope (z1) * base (z2));
}

TMPL inline Tangent
operator * (const C& z1, const Tangent& z2) {
  return Tangent (z1 * base (z2), z1 * slope (z2));
}

TMPL inline Tangent
operator * (const Tangent& z1, const C& z2) {
  return Tangent (base (z1) * z2, slope (z1) * z2);
}

TMPL inline Tangent
square (const Tangent& z) {
  return Tangent (square (base (z)), (base (z) + base (z)) * slope (z));
}

TMPL inline Tangent
invert (const Tangent& z) {
  C inv= invert (base (z));
  return Tangent (inv, - square (inv) * slope (z));
}

TMPL inline Tangent
operator / (const Tangent& z1, const Tangent& z2) {
  C inv= invert (base (z2));
  return Tangent (inv * base (z1),
                  inv * (slope (z1) - inv * base (z1) * slope (z2)));
}

TMPL inline Tangent
operator / (const C& z1, const Tangent& z2) {
  C inv= invert (base (z2));
  return Tangent (inv * z1, - square (inv) * z1 * slope (z2));
}

TMPL inline Tangent
operator / (const Tangent& z1, const C& z2) {
  C inv= invert (z2);
  return Tangent (inv * base (z1), inv * slope (z1));
}

ARITH_INT_SUGAR(TMPL,Tangent);

/******************************************************************************
* Elementary functions
******************************************************************************/

TMPL inline Tangent
sqrt (const Tangent& z) {
  C s= sqrt (base (z));
  return Tangent (s, invert (s + s) * slope (z));
}

TMPL inline Tangent
exp (const Tangent& z) {
  C e= exp (base (z));
  return Tangent (e, e * slope (z));
}

TMPL inline Tangent
log (const Tangent& z) {
  return Tangent (log (base (z)), invert (base (z)) * slope (z));
}

template<typename C, typename D, typename K> inline Tangent
pow (const Tangent& z, const K& y) {
  return exp (Tangent (y) * log (z));
}

TMPL inline Tangent
cos (const Tangent& z) {
  return Tangent (cos (base (z)), -sin (base (z)) * slope (z));
}

TMPL inline Tangent
sin (const Tangent& z) {
  return Tangent (sin (base (z)), cos (base (z)) * slope (z));
}

TMPL inline Tangent
tan (const Tangent& z) {
  return sin (z) / cos (z);
}

TMPL inline Tangent
cosh (const Tangent& z) {
  return Tangent (cosh (base (z)), sinh (base (z)) * slope (z));
}

TMPL inline Tangent
sinh (const Tangent& z) {
  return Tangent (sinh (base (z)), cosh (base (z)) * slope (z));
}

TMPL Tangent
tanh (const Tangent& z) {
  return sinh (z) / cosh (z);
}

TMPL Tangent
acos (const Tangent& z) {
  C df= -invert (sqrt (promote (1, base (z)) - square (base (z))));
  return Tangent (acos (base (z)), df * slope (z));
}

TMPL Tangent
asin (const Tangent& z) {
  C df= invert (sqrt (promote (1, base (z)) - square (base (z))));
  return Tangent (asin (base (z)), df * slope (z));
}

TMPL Tangent
atan (const Tangent& z) {
  C df= invert (promote (1, base (z)) + square (base (z)));
  return Tangent (asin (base (z)), df * slope (z));
}

TMPL Tangent
atan2 (const Tangent& y, const Tangent& x) {
  ERROR ("not yet implemented");
  return Tangent (atan2 (base (y), base (x)), 0);
}

HYPOT_SUGAR(TMPL,Tangent)
INV_TRIGO_SUGAR(TMPL,Tangent)
INV_HYPER_SUGAR(TMPL,Tangent)
ARG_HYPER_SUGAR(TMPL,Tangent)

/******************************************************************************
* Comparisons and other order related functions
******************************************************************************/

TMPL inline bool
operator <= (const Tangent& z1, const Tangent& z2) {
  return base (z1) < base (z2) || z1 == z2;
}

TMPL inline bool
operator < (const Tangent& z1, const Tangent& z2) {
  return base (z1) < base (z2);
}

TMPL inline bool
operator >= (const Tangent& z1, const Tangent& z2) {
  return base (z1) > base (z2) || z1 == z2;
}

TMPL inline bool
operator > (const Tangent& z1, const Tangent& z2) {
  return base (z1) > base (z2);
}

TMPL inline Tangent
min (const Tangent& z1, const Tangent& z2) {
  if (z1 <= z2) return z1; else return z2;
}

TMPL Tangent
max (const Tangent& z1, const Tangent& z2) {
  if (z1 >= z2) return z1; else return z2;
}

EQUAL_INT_SUGAR(TMPL,Tangent);
COMPARE_INT_SUGAR(TMPL,Tangent);

/******************************************************************************
* Floating point related functions
******************************************************************************/

TMPL inline bool is_finite (const Tangent& z) {
  return is_finite (base (z)) && is_finite (slope (z)); }
TMPL inline bool is_nan (const Tangent& z) {
  return is_nan (base (z)) || is_nan (slope (z)); }
TMPL inline bool is_infinite (const Tangent& z) {
  return !is_nan (z) && (is_infinite (base (z)) || is_infinite (slope (z))); }
TMPL inline bool is_fuzz (const Tangent& z) {
  return !is_nan (z) && (is_fuzz (base (z)) || is_fuzz (slope (z))); }
TMPL inline bool is_reliable (const Tangent& z) {
  return is_reliable (base (z)); }

TMPL inline Tangent change_precision (const Tangent& z, xnat prec) {
  return Tangent (change_precision (base (z), prec),
                  change_precision (slope (z), prec)); }
TMPL inline xnat precision (const Tangent& z) {
  return precision (base (z)); }
TMPL inline Tangent additive_error (const Tangent& z) {
  return Tangent (additive_error (base (z)), additive_error (slope (z))); }

TMPL inline double magnitude (const Tangent& z) {
  return magnitude (base (z)); }
TMPL inline xint exponent (const Tangent& z) {
  return exponent (base (z)); }
TMPL inline Tangent operator << (const Tangent& z, const xint& shift) {
  return Tangent (incexp2 (base (z), shift), incexp2 (slope (z), shift)); }
TMPL inline Tangent operator >> (const Tangent& z, const xint& shift) {
  return Tangent (decexp2 (base (z), shift), decexp2 (slope (z), shift)); }

TMPL inline Tangent sharpen (const Tangent& z) {
  return Tangent (sharpen (base (z)), sharpen (slope (z))); }
template<typename C, typename D, typename C2, typename D2> inline Tangent
blur (const Tangent& z, const tangent<C2,D2>& r) {
  return Tangent (blur (base (z), base (r)), blur (slope (z), slope (r))); }

/******************************************************************************
* Complex tangent numbers
******************************************************************************/

TMPL inline Abs_tangent
abs (const Tangent& z) {
  // FIXME: this is not really nice
  return Abs_tangent (abs (base (z)), abs (slope (z)));
}

TMPL inline Real_tangent
Re (const Tangent& z) {
  return Real_tangent (Re (base (z)), Re (slope (z)));
}

TMPL inline Real_tangent
Im (const Tangent& z) {
  return Real_tangent (Im (base (z)), Im (slope (z)));
}

#undef TMPL_DEF
#undef TMPL
#undef Tangent
#undef Abs_tangent
#undef Real_tangent
} // namespace mmx
#endif // __MMX_TANGENT_HPP