diff --git a/src/base/atomicops_internals_mips64_gcc.h b/src/base/atomicops_internals_mips64_gcc.h
index 85b4e46..cf2e194 100644
--- a/src/base/atomicops_internals_mips64_gcc.h
+++ b/src/base/atomicops_internals_mips64_gcc.h
@@ -91,18 +91,19 @@
                                           Atomic32 increment) {
   Atomic32 temp, temp2;
 
-  __asm__ __volatile__(".set push\n"
-                       ".set noreorder\n"
-                       "1:\n"
-                       "ll %0, %2\n"  // temp = *ptr
-                       "addu %1, %0, %3\n"  // temp2 = temp + increment
-                       "sc %1, %2\n"  // *ptr = temp2 (with atomic check)
-                       "beqz %1, 1b\n"  // start again on atomic error
-                       "addu %1, %0, %3\n"  // temp2 = temp + increment
-                       ".set pop\n"
-                       : "=&r" (temp), "=&r" (temp2), "=m" (*ptr)
-                       : "Ir" (increment), "m" (*ptr)
-                       : "memory");
+  __asm__ __volatile__(
+      ".set push\n"
+      ".set noreorder\n"
+      "1:\n"
+      "ll %0, %2\n"        // temp = *ptr
+      "addu %1, %0, %3\n"  // temp2 = temp + increment
+      "sc %1, %2\n"        // *ptr = temp2 (with atomic check)
+      "beqz %1, 1b\n"      // start again on atomic error
+      "addu %1, %0, %3\n"  // temp2 = temp + increment
+      ".set pop\n"
+      : "=&r"(temp), "=&r"(temp2), "=ZC"(*ptr)
+      : "Ir"(increment), "m"(*ptr)
+      : "memory");
   // temp2 now holds the final value.
   return temp2;
 }
@@ -228,18 +229,19 @@
                                           Atomic64 increment) {
   Atomic64 temp, temp2;
 
-  __asm__ __volatile__(".set push\n"
-                       ".set noreorder\n"
-                       "1:\n"
-                       "lld %0, %2\n"  // temp = *ptr
-                       "daddu %1, %0, %3\n"  // temp2 = temp + increment
-                       "scd %1, %2\n"  // *ptr = temp2 (with atomic check)
-                       "beqz %1, 1b\n"  // start again on atomic error
-                       "daddu %1, %0, %3\n"  // temp2 = temp + increment
-                       ".set pop\n"
-                       : "=&r" (temp), "=&r" (temp2), "=m" (*ptr)
-                       : "Ir" (increment), "m" (*ptr)
-                       : "memory");
+  __asm__ __volatile__(
+      ".set push\n"
+      ".set noreorder\n"
+      "1:\n"
+      "lld %0, %2\n"        // temp = *ptr
+      "daddu %1, %0, %3\n"  // temp2 = temp + increment
+      "scd %1, %2\n"        // *ptr = temp2 (with atomic check)
+      "beqz %1, 1b\n"       // start again on atomic error
+      "daddu %1, %0, %3\n"  // temp2 = temp + increment
+      ".set pop\n"
+      : "=&r"(temp), "=&r"(temp2), "=ZC"(*ptr)
+      : "Ir"(increment), "m"(*ptr)
+      : "memory");
   // temp2 now holds the final value.
   return temp2;
 }
diff --git a/src/base/bits.h b/src/base/bits.h
index 2e6527b..1bb3a0f 100644
--- a/src/base/bits.h
+++ b/src/base/bits.h
@@ -111,7 +111,6 @@
   return result;
 }
 
-
 // CountTrailingZeros32(value) returns the number of zero bits preceding the
 // least significant 1 bit in |value| if |value| is non-zero, otherwise it
 // returns 32.
@@ -147,6 +146,14 @@
 #endif
 }
 
+// Overloaded versions of CountTrailingZeros32/64.
+inline unsigned CountTrailingZeros(uint32_t value) {
+  return CountTrailingZeros32(value);
+}
+
+inline unsigned CountTrailingZeros(uint64_t value) {
+  return CountTrailingZeros64(value);
+}
 
 // Returns true iff |value| is a power of 2.
 inline bool IsPowerOfTwo32(uint32_t value) {
diff --git a/src/base/cpu.cc b/src/base/cpu.cc
index 12a3881..16eb7c9 100644
--- a/src/base/cpu.cc
+++ b/src/base/cpu.cc
@@ -338,7 +338,8 @@
       has_vfp_(false),
       has_vfp3_(false),
       has_vfp3_d32_(false),
-      is_fp64_mode_(false) {
+      is_fp64_mode_(false),
+      has_non_stop_time_stamp_counter_(false) {
   memcpy(vendor_, "Unknown", 8);
 #if V8_OS_NACL
 // Portable host shouldn't do feature detection.
@@ -419,6 +420,13 @@
     has_sahf_ = (cpu_info[2] & 0x00000001) != 0;
   }
 
+  // Check if CPU has non stoppable time stamp counter.
+  const int parameter_containing_non_stop_time_stamp_counter = 0x80000007;
+  if (num_ext_ids >= parameter_containing_non_stop_time_stamp_counter) {
+    __cpuid(cpu_info, parameter_containing_non_stop_time_stamp_counter);
+    has_non_stop_time_stamp_counter_ = (cpu_info[3] & (1 << 8)) != 0;
+  }
+
 #elif V8_HOST_ARCH_ARM
 
 #if V8_OS_LINUX
diff --git a/src/base/cpu.h b/src/base/cpu.h
index 3778d27..19d4102 100644
--- a/src/base/cpu.h
+++ b/src/base/cpu.h
@@ -97,6 +97,9 @@
   bool has_lzcnt() const { return has_lzcnt_; }
   bool has_popcnt() const { return has_popcnt_; }
   bool is_atom() const { return is_atom_; }
+  bool has_non_stop_time_stamp_counter() const {
+    return has_non_stop_time_stamp_counter_;
+  }
 
   // arm features
   bool has_idiva() const { return has_idiva_; }
@@ -148,6 +151,7 @@
   bool has_vfp3_;
   bool has_vfp3_d32_;
   bool is_fp64_mode_;
+  bool has_non_stop_time_stamp_counter_;
 };
 
 }  // namespace base
diff --git a/src/base/file-utils.cc b/src/base/file-utils.cc
new file mode 100644
index 0000000..2262df9
--- /dev/null
+++ b/src/base/file-utils.cc
@@ -0,0 +1,36 @@
+// Copyright 2016 the V8 project authors. All rights reserved.
+// Use of this source code is governed by a BSD-style license that can be
+// found in the LICENSE file.
+
+#include "src/base/file-utils.h"
+
+#include <stdlib.h>
+#include <string.h>
+
+#include "src/base/platform/platform.h"
+
+namespace v8 {
+namespace internal {
+
+char* RelativePath(char** buffer, const char* exec_path, const char* name) {
+  DCHECK(exec_path);
+  int path_separator = static_cast<int>(strlen(exec_path)) - 1;
+  while (path_separator >= 0 &&
+         !base::OS::isDirectorySeparator(exec_path[path_separator])) {
+    path_separator--;
+  }
+  if (path_separator >= 0) {
+    int name_length = static_cast<int>(strlen(name));
+    *buffer =
+        reinterpret_cast<char*>(calloc(path_separator + name_length + 2, 1));
+    *buffer[0] = '\0';
+    strncat(*buffer, exec_path, path_separator + 1);
+    strncat(*buffer, name, name_length);
+  } else {
+    *buffer = strdup(name);
+  }
+  return *buffer;
+}
+
+}  // namespace internal
+}  // namespace v8
diff --git a/src/base/file-utils.h b/src/base/file-utils.h
new file mode 100644
index 0000000..ce9e9a1
--- /dev/null
+++ b/src/base/file-utils.h
@@ -0,0 +1,18 @@
+// Copyright 2016 the V8 project authors. All rights reserved.
+// Use of this source code is governed by a BSD-style license that can be
+// found in the LICENSE file.
+
+#ifndef V8_FILE_UTILS_H_
+#define V8_FILE_UTILS_H_
+
+namespace v8 {
+namespace internal {
+
+// Helper functions to manipulate file paths.
+
+char* RelativePath(char** buffer, const char* exec_path, const char* name);
+
+}  // namespace internal
+}  // namespace v8
+
+#endif  // V8_FILE_UTILS_H_
diff --git a/src/base/hashmap.h b/src/base/hashmap.h
new file mode 100644
index 0000000..efb5dc8
--- /dev/null
+++ b/src/base/hashmap.h
@@ -0,0 +1,352 @@
+// Copyright 2012 the V8 project authors. All rights reserved.
+// Use of this source code is governed by a BSD-style license that can be
+// found in the LICENSE file.
+
+// The reason we write our own hash map instead of using unordered_map in STL,
+// is that STL containers use a mutex pool on debug build, which will lead to
+// deadlock when we are using async signal handler.
+
+#ifndef V8_BASE_HASHMAP_H_
+#define V8_BASE_HASHMAP_H_
+
+#include <stdlib.h>
+
+#include "src/base/bits.h"
+#include "src/base/logging.h"
+
+namespace v8 {
+namespace base {
+
+class DefaultAllocationPolicy {
+ public:
+  V8_INLINE void* New(size_t size) { return malloc(size); }
+  V8_INLINE static void Delete(void* p) { free(p); }
+};
+
+template <class AllocationPolicy>
+class TemplateHashMapImpl {
+ public:
+  typedef bool (*MatchFun)(void* key1, void* key2);
+
+  // The default capacity.  This is used by the call sites which want
+  // to pass in a non-default AllocationPolicy but want to use the
+  // default value of capacity specified by the implementation.
+  static const uint32_t kDefaultHashMapCapacity = 8;
+
+  // initial_capacity is the size of the initial hash map;
+  // it must be a power of 2 (and thus must not be 0).
+  TemplateHashMapImpl(MatchFun match,
+                      uint32_t capacity = kDefaultHashMapCapacity,
+                      AllocationPolicy allocator = AllocationPolicy());
+
+  ~TemplateHashMapImpl();
+
+  // HashMap entries are (key, value, hash) triplets.
+  // Some clients may not need to use the value slot
+  // (e.g. implementers of sets, where the key is the value).
+  struct Entry {
+    void* key;
+    void* value;
+    uint32_t hash;  // The full hash value for key
+    int order;      // If you never remove entries this is the insertion order.
+  };
+
+  // If an entry with matching key is found, returns that entry.
+  // Otherwise, NULL is returned.
+  Entry* Lookup(void* key, uint32_t hash) const;
+
+  // If an entry with matching key is found, returns that entry.
+  // If no matching entry is found, a new entry is inserted with
+  // corresponding key, key hash, and NULL value.
+  Entry* LookupOrInsert(void* key, uint32_t hash,
+                        AllocationPolicy allocator = AllocationPolicy());
+
+  // Removes the entry with matching key.
+  // It returns the value of the deleted entry
+  // or null if there is no value for such key.
+  void* Remove(void* key, uint32_t hash);
+
+  // Empties the hash map (occupancy() == 0).
+  void Clear();
+
+  // The number of (non-empty) entries in the table.
+  uint32_t occupancy() const { return occupancy_; }
+
+  // The capacity of the table. The implementation
+  // makes sure that occupancy is at most 80% of
+  // the table capacity.
+  uint32_t capacity() const { return capacity_; }
+
+  // Iteration
+  //
+  // for (Entry* p = map.Start(); p != NULL; p = map.Next(p)) {
+  //   ...
+  // }
+  //
+  // If entries are inserted during iteration, the effect of
+  // calling Next() is undefined.
+  Entry* Start() const;
+  Entry* Next(Entry* p) const;
+
+  // Some match functions defined for convenience.
+  static bool PointersMatch(void* key1, void* key2) { return key1 == key2; }
+
+ private:
+  MatchFun match_;
+  Entry* map_;
+  uint32_t capacity_;
+  uint32_t occupancy_;
+
+  Entry* map_end() const { return map_ + capacity_; }
+  Entry* Probe(void* key, uint32_t hash) const;
+  void Initialize(uint32_t capacity, AllocationPolicy allocator);
+  void Resize(AllocationPolicy allocator);
+};
+
+typedef TemplateHashMapImpl<DefaultAllocationPolicy> HashMap;
+
+template <class AllocationPolicy>
+TemplateHashMapImpl<AllocationPolicy>::TemplateHashMapImpl(
+    MatchFun match, uint32_t initial_capacity, AllocationPolicy allocator) {
+  match_ = match;
+  Initialize(initial_capacity, allocator);
+}
+
+template <class AllocationPolicy>
+TemplateHashMapImpl<AllocationPolicy>::~TemplateHashMapImpl() {
+  AllocationPolicy::Delete(map_);
+}
+
+template <class AllocationPolicy>
+typename TemplateHashMapImpl<AllocationPolicy>::Entry*
+TemplateHashMapImpl<AllocationPolicy>::Lookup(void* key, uint32_t hash) const {
+  Entry* p = Probe(key, hash);
+  return p->key != NULL ? p : NULL;
+}
+
+template <class AllocationPolicy>
+typename TemplateHashMapImpl<AllocationPolicy>::Entry*
+TemplateHashMapImpl<AllocationPolicy>::LookupOrInsert(
+    void* key, uint32_t hash, AllocationPolicy allocator) {
+  // Find a matching entry.
+  Entry* p = Probe(key, hash);
+  if (p->key != NULL) {
+    return p;
+  }
+
+  // No entry found; insert one.
+  p->key = key;
+  p->value = NULL;
+  p->hash = hash;
+  p->order = occupancy_;
+  occupancy_++;
+
+  // Grow the map if we reached >= 80% occupancy.
+  if (occupancy_ + occupancy_ / 4 >= capacity_) {
+    Resize(allocator);
+    p = Probe(key, hash);
+  }
+
+  return p;
+}
+
+template <class AllocationPolicy>
+void* TemplateHashMapImpl<AllocationPolicy>::Remove(void* key, uint32_t hash) {
+  // Lookup the entry for the key to remove.
+  Entry* p = Probe(key, hash);
+  if (p->key == NULL) {
+    // Key not found nothing to remove.
+    return NULL;
+  }
+
+  void* value = p->value;
+  // To remove an entry we need to ensure that it does not create an empty
+  // entry that will cause the search for another entry to stop too soon. If all
+  // the entries between the entry to remove and the next empty slot have their
+  // initial position inside this interval, clearing the entry to remove will
+  // not break the search. If, while searching for the next empty entry, an
+  // entry is encountered which does not have its initial position between the
+  // entry to remove and the position looked at, then this entry can be moved to
+  // the place of the entry to remove without breaking the search for it. The
+  // entry made vacant by this move is now the entry to remove and the process
+  // starts over.
+  // Algorithm from http://en.wikipedia.org/wiki/Open_addressing.
+
+  // This guarantees loop termination as there is at least one empty entry so
+  // eventually the removed entry will have an empty entry after it.
+  DCHECK(occupancy_ < capacity_);
+
+  // p is the candidate entry to clear. q is used to scan forwards.
+  Entry* q = p;  // Start at the entry to remove.
+  while (true) {
+    // Move q to the next entry.
+    q = q + 1;
+    if (q == map_end()) {
+      q = map_;
+    }
+
+    // All entries between p and q have their initial position between p and q
+    // and the entry p can be cleared without breaking the search for these
+    // entries.
+    if (q->key == NULL) {
+      break;
+    }
+
+    // Find the initial position for the entry at position q.
+    Entry* r = map_ + (q->hash & (capacity_ - 1));
+
+    // If the entry at position q has its initial position outside the range
+    // between p and q it can be moved forward to position p and will still be
+    // found. There is now a new candidate entry for clearing.
+    if ((q > p && (r <= p || r > q)) || (q < p && (r <= p && r > q))) {
+      *p = *q;
+      p = q;
+    }
+  }
+
+  // Clear the entry which is allowed to en emptied.
+  p->key = NULL;
+  occupancy_--;
+  return value;
+}
+
+template <class AllocationPolicy>
+void TemplateHashMapImpl<AllocationPolicy>::Clear() {
+  // Mark all entries as empty.
+  const Entry* end = map_end();
+  for (Entry* p = map_; p < end; p++) {
+    p->key = NULL;
+  }
+  occupancy_ = 0;
+}
+
+template <class AllocationPolicy>
+typename TemplateHashMapImpl<AllocationPolicy>::Entry*
+TemplateHashMapImpl<AllocationPolicy>::Start() const {
+  return Next(map_ - 1);
+}
+
+template <class AllocationPolicy>
+typename TemplateHashMapImpl<AllocationPolicy>::Entry*
+TemplateHashMapImpl<AllocationPolicy>::Next(Entry* p) const {
+  const Entry* end = map_end();
+  DCHECK(map_ - 1 <= p && p < end);
+  for (p++; p < end; p++) {
+    if (p->key != NULL) {
+      return p;
+    }
+  }
+  return NULL;
+}
+
+template <class AllocationPolicy>
+typename TemplateHashMapImpl<AllocationPolicy>::Entry*
+TemplateHashMapImpl<AllocationPolicy>::Probe(void* key, uint32_t hash) const {
+  DCHECK(key != NULL);
+
+  DCHECK(base::bits::IsPowerOfTwo32(capacity_));
+  Entry* p = map_ + (hash & (capacity_ - 1));
+  const Entry* end = map_end();
+  DCHECK(map_ <= p && p < end);
+
+  DCHECK(occupancy_ < capacity_);  // Guarantees loop termination.
+  while (p->key != NULL && (hash != p->hash || !match_(key, p->key))) {
+    p++;
+    if (p >= end) {
+      p = map_;
+    }
+  }
+
+  return p;
+}
+
+template <class AllocationPolicy>
+void TemplateHashMapImpl<AllocationPolicy>::Initialize(
+    uint32_t capacity, AllocationPolicy allocator) {
+  DCHECK(base::bits::IsPowerOfTwo32(capacity));
+  map_ = reinterpret_cast<Entry*>(allocator.New(capacity * sizeof(Entry)));
+  if (map_ == NULL) {
+    FATAL("Out of memory: HashMap::Initialize");
+    return;
+  }
+  capacity_ = capacity;
+  Clear();
+}
+
+template <class AllocationPolicy>
+void TemplateHashMapImpl<AllocationPolicy>::Resize(AllocationPolicy allocator) {
+  Entry* map = map_;
+  uint32_t n = occupancy_;
+
+  // Allocate larger map.
+  Initialize(capacity_ * 2, allocator);
+
+  // Rehash all current entries.
+  for (Entry* p = map; n > 0; p++) {
+    if (p->key != NULL) {
+      Entry* entry = LookupOrInsert(p->key, p->hash, allocator);
+      entry->value = p->value;
+      entry->order = p->order;
+      n--;
+    }
+  }
+
+  // Delete old map.
+  AllocationPolicy::Delete(map);
+}
+
+// A hash map for pointer keys and values with an STL-like interface.
+template <class Key, class Value, class AllocationPolicy>
+class TemplateHashMap : private TemplateHashMapImpl<AllocationPolicy> {
+ public:
+  STATIC_ASSERT(sizeof(Key*) == sizeof(void*));    // NOLINT
+  STATIC_ASSERT(sizeof(Value*) == sizeof(void*));  // NOLINT
+  struct value_type {
+    Key* first;
+    Value* second;
+  };
+
+  class Iterator {
+   public:
+    Iterator& operator++() {
+      entry_ = map_->Next(entry_);
+      return *this;
+    }
+
+    value_type* operator->() { return reinterpret_cast<value_type*>(entry_); }
+    bool operator!=(const Iterator& other) { return entry_ != other.entry_; }
+
+   private:
+    Iterator(const TemplateHashMapImpl<AllocationPolicy>* map,
+             typename TemplateHashMapImpl<AllocationPolicy>::Entry* entry)
+        : map_(map), entry_(entry) {}
+
+    const TemplateHashMapImpl<AllocationPolicy>* map_;
+    typename TemplateHashMapImpl<AllocationPolicy>::Entry* entry_;
+
+    friend class TemplateHashMap;
+  };
+
+  TemplateHashMap(
+      typename TemplateHashMapImpl<AllocationPolicy>::MatchFun match,
+      AllocationPolicy allocator = AllocationPolicy())
+      : TemplateHashMapImpl<AllocationPolicy>(
+            match,
+            TemplateHashMapImpl<AllocationPolicy>::kDefaultHashMapCapacity,
+            allocator) {}
+
+  Iterator begin() const { return Iterator(this, this->Start()); }
+  Iterator end() const { return Iterator(this, NULL); }
+  Iterator find(Key* key, bool insert = false,
+                AllocationPolicy allocator = AllocationPolicy()) {
+    if (insert) {
+      return Iterator(this, this->LookupOrInsert(key, key->Hash(), allocator));
+    }
+    return Iterator(this, this->Lookup(key, key->Hash()));
+  }
+};
+
+}  // namespace base
+}  // namespace v8
+
+#endif  // V8_BASE_HASHMAP_H_
diff --git a/src/base/ieee754.cc b/src/base/ieee754.cc
new file mode 100644
index 0000000..b7178ca
--- /dev/null
+++ b/src/base/ieee754.cc
@@ -0,0 +1,2313 @@
+// The following is adapted from fdlibm (http://www.netlib.org/fdlibm).
+//
+// ====================================================
+// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+//
+// Developed at SunSoft, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+//
+// The original source code covered by the above license above has been
+// modified significantly by Google Inc.
+// Copyright 2016 the V8 project authors. All rights reserved.
+
+#include "src/base/ieee754.h"
+
+#include <cmath>
+#include <limits>
+
+#include "src/base/build_config.h"
+#include "src/base/macros.h"
+
+namespace v8 {
+namespace base {
+namespace ieee754 {
+
+namespace {
+
+/* Disable "potential divide by 0" warning in Visual Studio compiler. */
+
+#if V8_CC_MSVC
+
+#pragma warning(disable : 4723)
+
+#endif
+
+/*
+ * The original fdlibm code used statements like:
+ *  n0 = ((*(int*)&one)>>29)^1;   * index of high word *
+ *  ix0 = *(n0+(int*)&x);     * high word of x *
+ *  ix1 = *((1-n0)+(int*)&x);   * low word of x *
+ * to dig two 32 bit words out of the 64 bit IEEE floating point
+ * value.  That is non-ANSI, and, moreover, the gcc instruction
+ * scheduler gets it wrong.  We instead use the following macros.
+ * Unlike the original code, we determine the endianness at compile
+ * time, not at run time; I don't see much benefit to selecting
+ * endianness at run time.
+ */
+
+/*
+ * A union which permits us to convert between a double and two 32 bit
+ * ints.
+ */
+
+#if V8_TARGET_LITTLE_ENDIAN
+
+typedef union {
+  double value;
+  struct {
+    uint32_t lsw;
+    uint32_t msw;
+  } parts;
+  struct {
+    uint64_t w;
+  } xparts;
+} ieee_double_shape_type;
+
+#else
+
+typedef union {
+  double value;
+  struct {
+    uint32_t msw;
+    uint32_t lsw;
+  } parts;
+  struct {
+    uint64_t w;
+  } xparts;
+} ieee_double_shape_type;
+
+#endif
+
+/* Get two 32 bit ints from a double.  */
+
+#define EXTRACT_WORDS(ix0, ix1, d) \
+  do {                             \
+    ieee_double_shape_type ew_u;   \
+    ew_u.value = (d);              \
+    (ix0) = ew_u.parts.msw;        \
+    (ix1) = ew_u.parts.lsw;        \
+  } while (0)
+
+/* Get a 64-bit int from a double. */
+#define EXTRACT_WORD64(ix, d)    \
+  do {                           \
+    ieee_double_shape_type ew_u; \
+    ew_u.value = (d);            \
+    (ix) = ew_u.xparts.w;        \
+  } while (0)
+
+/* Get the more significant 32 bit int from a double.  */
+
+#define GET_HIGH_WORD(i, d)      \
+  do {                           \
+    ieee_double_shape_type gh_u; \
+    gh_u.value = (d);            \
+    (i) = gh_u.parts.msw;        \
+  } while (0)
+
+/* Get the less significant 32 bit int from a double.  */
+
+#define GET_LOW_WORD(i, d)       \
+  do {                           \
+    ieee_double_shape_type gl_u; \
+    gl_u.value = (d);            \
+    (i) = gl_u.parts.lsw;        \
+  } while (0)
+
+/* Set a double from two 32 bit ints.  */
+
+#define INSERT_WORDS(d, ix0, ix1) \
+  do {                            \
+    ieee_double_shape_type iw_u;  \
+    iw_u.parts.msw = (ix0);       \
+    iw_u.parts.lsw = (ix1);       \
+    (d) = iw_u.value;             \
+  } while (0)
+
+/* Set a double from a 64-bit int. */
+#define INSERT_WORD64(d, ix)     \
+  do {                           \
+    ieee_double_shape_type iw_u; \
+    iw_u.xparts.w = (ix);        \
+    (d) = iw_u.value;            \
+  } while (0)
+
+/* Set the more significant 32 bits of a double from an int.  */
+
+#define SET_HIGH_WORD(d, v)      \
+  do {                           \
+    ieee_double_shape_type sh_u; \
+    sh_u.value = (d);            \
+    sh_u.parts.msw = (v);        \
+    (d) = sh_u.value;            \
+  } while (0)
+
+/* Set the less significant 32 bits of a double from an int.  */
+
+#define SET_LOW_WORD(d, v)       \
+  do {                           \
+    ieee_double_shape_type sl_u; \
+    sl_u.value = (d);            \
+    sl_u.parts.lsw = (v);        \
+    (d) = sl_u.value;            \
+  } while (0)
+
+/* Support macro. */
+
+#define STRICT_ASSIGN(type, lval, rval) ((lval) = (rval))
+
+int32_t __ieee754_rem_pio2(double x, double *y) WARN_UNUSED_RESULT;
+double __kernel_cos(double x, double y) WARN_UNUSED_RESULT;
+int __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec,
+                      const int32_t *ipio2) WARN_UNUSED_RESULT;
+double __kernel_sin(double x, double y, int iy) WARN_UNUSED_RESULT;
+
+/* __ieee754_rem_pio2(x,y)
+ *
+ * return the remainder of x rem pi/2 in y[0]+y[1]
+ * use __kernel_rem_pio2()
+ */
+int32_t __ieee754_rem_pio2(double x, double *y) {
+  /*
+   * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
+   */
+  static const int32_t two_over_pi[] = {
+      0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 0x95993C,
+      0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 0x424DD2, 0xE00649,
+      0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 0xA73EE8, 0x8235F5, 0x2EBB44,
+      0x84E99C, 0x7026B4, 0x5F7E41, 0x3991D6, 0x398353, 0x39F49C, 0x845F8B,
+      0xBDF928, 0x3B1FF8, 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D,
+      0x367ECF, 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
+      0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 0x560330,
+      0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 0x91615E, 0xE61B08,
+      0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 0x4D7327, 0x310606, 0x1556CA,
+      0x73A8C9, 0x60E27B, 0xC08C6B,
+  };
+
+  static const int32_t npio2_hw[] = {
+      0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C,
+      0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C,
+      0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A,
+      0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C,
+      0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB,
+      0x404858EB, 0x404921FB,
+  };
+
+  /*
+   * invpio2:  53 bits of 2/pi
+   * pio2_1:   first  33 bit of pi/2
+   * pio2_1t:  pi/2 - pio2_1
+   * pio2_2:   second 33 bit of pi/2
+   * pio2_2t:  pi/2 - (pio2_1+pio2_2)
+   * pio2_3:   third  33 bit of pi/2
+   * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
+   */
+
+  static const double
+      zero = 0.00000000000000000000e+00,    /* 0x00000000, 0x00000000 */
+      half = 5.00000000000000000000e-01,    /* 0x3FE00000, 0x00000000 */
+      two24 = 1.67772160000000000000e+07,   /* 0x41700000, 0x00000000 */
+      invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
+      pio2_1 = 1.57079632673412561417e+00,  /* 0x3FF921FB, 0x54400000 */
+      pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
+      pio2_2 = 6.07710050630396597660e-11,  /* 0x3DD0B461, 0x1A600000 */
+      pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
+      pio2_3 = 2.02226624871116645580e-21,  /* 0x3BA3198A, 0x2E000000 */
+      pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
+
+  double z, w, t, r, fn;
+  double tx[3];
+  int32_t e0, i, j, nx, n, ix, hx;
+  uint32_t low;
+
+  z = 0;
+  GET_HIGH_WORD(hx, x); /* high word of x */
+  ix = hx & 0x7fffffff;
+  if (ix <= 0x3fe921fb) { /* |x| ~<= pi/4 , no need for reduction */
+    y[0] = x;
+    y[1] = 0;
+    return 0;
+  }
+  if (ix < 0x4002d97c) { /* |x| < 3pi/4, special case with n=+-1 */
+    if (hx > 0) {
+      z = x - pio2_1;
+      if (ix != 0x3ff921fb) { /* 33+53 bit pi is good enough */
+        y[0] = z - pio2_1t;
+        y[1] = (z - y[0]) - pio2_1t;
+      } else { /* near pi/2, use 33+33+53 bit pi */
+        z -= pio2_2;
+        y[0] = z - pio2_2t;
+        y[1] = (z - y[0]) - pio2_2t;
+      }
+      return 1;
+    } else { /* negative x */
+      z = x + pio2_1;
+      if (ix != 0x3ff921fb) { /* 33+53 bit pi is good enough */
+        y[0] = z + pio2_1t;
+        y[1] = (z - y[0]) + pio2_1t;
+      } else { /* near pi/2, use 33+33+53 bit pi */
+        z += pio2_2;
+        y[0] = z + pio2_2t;
+        y[1] = (z - y[0]) + pio2_2t;
+      }
+      return -1;
+    }
+  }
+  if (ix <= 0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */
+    t = fabs(x);
+    n = static_cast<int32_t>(t * invpio2 + half);
+    fn = static_cast<double>(n);
+    r = t - fn * pio2_1;
+    w = fn * pio2_1t; /* 1st round good to 85 bit */
+    if (n < 32 && ix != npio2_hw[n - 1]) {
+      y[0] = r - w; /* quick check no cancellation */
+    } else {
+      uint32_t high;
+      j = ix >> 20;
+      y[0] = r - w;
+      GET_HIGH_WORD(high, y[0]);
+      i = j - ((high >> 20) & 0x7ff);
+      if (i > 16) { /* 2nd iteration needed, good to 118 */
+        t = r;
+        w = fn * pio2_2;
+        r = t - w;
+        w = fn * pio2_2t - ((t - r) - w);
+        y[0] = r - w;
+        GET_HIGH_WORD(high, y[0]);
+        i = j - ((high >> 20) & 0x7ff);
+        if (i > 49) { /* 3rd iteration need, 151 bits acc */
+          t = r;      /* will cover all possible cases */
+          w = fn * pio2_3;
+          r = t - w;
+          w = fn * pio2_3t - ((t - r) - w);
+          y[0] = r - w;
+        }
+      }
+    }
+    y[1] = (r - y[0]) - w;
+    if (hx < 0) {
+      y[0] = -y[0];
+      y[1] = -y[1];
+      return -n;
+    } else {
+      return n;
+    }
+  }
+  /*
+   * all other (large) arguments
+   */
+  if (ix >= 0x7ff00000) { /* x is inf or NaN */
+    y[0] = y[1] = x - x;
+    return 0;
+  }
+  /* set z = scalbn(|x|,ilogb(x)-23) */
+  GET_LOW_WORD(low, x);
+  SET_LOW_WORD(z, low);
+  e0 = (ix >> 20) - 1046; /* e0 = ilogb(z)-23; */
+  SET_HIGH_WORD(z, ix - static_cast<int32_t>(e0 << 20));
+  for (i = 0; i < 2; i++) {
+    tx[i] = static_cast<double>(static_cast<int32_t>(z));
+    z = (z - tx[i]) * two24;
+  }
+  tx[2] = z;
+  nx = 3;
+  while (tx[nx - 1] == zero) nx--; /* skip zero term */
+  n = __kernel_rem_pio2(tx, y, e0, nx, 2, two_over_pi);
+  if (hx < 0) {
+    y[0] = -y[0];
+    y[1] = -y[1];
+    return -n;
+  }
+  return n;
+}
+
+/* __kernel_cos( x,  y )
+ * kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
+ * Input x is assumed to be bounded by ~pi/4 in magnitude.
+ * Input y is the tail of x.
+ *
+ * Algorithm
+ *      1. Since cos(-x) = cos(x), we need only to consider positive x.
+ *      2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0.
+ *      3. cos(x) is approximated by a polynomial of degree 14 on
+ *         [0,pi/4]
+ *                                       4            14
+ *              cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
+ *         where the remez error is
+ *
+ *      |              2     4     6     8     10    12     14 |     -58
+ *      |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  )| <= 2
+ *      |                                                      |
+ *
+ *                     4     6     8     10    12     14
+ *      4. let r = C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  , then
+ *             cos(x) = 1 - x*x/2 + r
+ *         since cos(x+y) ~ cos(x) - sin(x)*y
+ *                        ~ cos(x) - x*y,
+ *         a correction term is necessary in cos(x) and hence
+ *              cos(x+y) = 1 - (x*x/2 - (r - x*y))
+ *         For better accuracy when x > 0.3, let qx = |x|/4 with
+ *         the last 32 bits mask off, and if x > 0.78125, let qx = 0.28125.
+ *         Then
+ *              cos(x+y) = (1-qx) - ((x*x/2-qx) - (r-x*y)).
+ *         Note that 1-qx and (x*x/2-qx) is EXACT here, and the
+ *         magnitude of the latter is at least a quarter of x*x/2,
+ *         thus, reducing the rounding error in the subtraction.
+ */
+V8_INLINE double __kernel_cos(double x, double y) {
+  static const double
+      one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
+      C1 = 4.16666666666666019037e-02,  /* 0x3FA55555, 0x5555554C */
+      C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
+      C3 = 2.48015872894767294178e-05,  /* 0x3EFA01A0, 0x19CB1590 */
+      C4 = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */
+      C5 = 2.08757232129817482790e-09,  /* 0x3E21EE9E, 0xBDB4B1C4 */
+      C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
+
+  double a, iz, z, r, qx;
+  int32_t ix;
+  GET_HIGH_WORD(ix, x);
+  ix &= 0x7fffffff;                           /* ix = |x|'s high word*/
+  if (ix < 0x3e400000) {                      /* if x < 2**27 */
+    if (static_cast<int>(x) == 0) return one; /* generate inexact */
+  }
+  z = x * x;
+  r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * C6)))));
+  if (ix < 0x3FD33333) { /* if |x| < 0.3 */
+    return one - (0.5 * z - (z * r - x * y));
+  } else {
+    if (ix > 0x3fe90000) { /* x > 0.78125 */
+      qx = 0.28125;
+    } else {
+      INSERT_WORDS(qx, ix - 0x00200000, 0); /* x/4 */
+    }
+    iz = 0.5 * z - qx;
+    a = one - qx;
+    return a - (iz - (z * r - x * y));
+  }
+}
+
+/* __kernel_rem_pio2(x,y,e0,nx,prec,ipio2)
+ * double x[],y[]; int e0,nx,prec; int ipio2[];
+ *
+ * __kernel_rem_pio2 return the last three digits of N with
+ *              y = x - N*pi/2
+ * so that |y| < pi/2.
+ *
+ * The method is to compute the integer (mod 8) and fraction parts of
+ * (2/pi)*x without doing the full multiplication. In general we
+ * skip the part of the product that are known to be a huge integer (
+ * more accurately, = 0 mod 8 ). Thus the number of operations are
+ * independent of the exponent of the input.
+ *
+ * (2/pi) is represented by an array of 24-bit integers in ipio2[].
+ *
+ * Input parameters:
+ *      x[]     The input value (must be positive) is broken into nx
+ *              pieces of 24-bit integers in double precision format.
+ *              x[i] will be the i-th 24 bit of x. The scaled exponent
+ *              of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
+ *              match x's up to 24 bits.
+ *
+ *              Example of breaking a double positive z into x[0]+x[1]+x[2]:
+ *                      e0 = ilogb(z)-23
+ *                      z  = scalbn(z,-e0)
+ *              for i = 0,1,2
+ *                      x[i] = floor(z)
+ *                      z    = (z-x[i])*2**24
+ *
+ *
+ *      y[]     output result in an array of double precision numbers.
+ *              The dimension of y[] is:
+ *                      24-bit  precision       1
+ *                      53-bit  precision       2
+ *                      64-bit  precision       2
+ *                      113-bit precision       3
+ *              The actual value is the sum of them. Thus for 113-bit
+ *              precison, one may have to do something like:
+ *
+ *              long double t,w,r_head, r_tail;
+ *              t = (long double)y[2] + (long double)y[1];
+ *              w = (long double)y[0];
+ *              r_head = t+w;
+ *              r_tail = w - (r_head - t);
+ *
+ *      e0      The exponent of x[0]
+ *
+ *      nx      dimension of x[]
+ *
+ *      prec    an integer indicating the precision:
+ *                      0       24  bits (single)
+ *                      1       53  bits (double)
+ *                      2       64  bits (extended)
+ *                      3       113 bits (quad)
+ *
+ *      ipio2[]
+ *              integer array, contains the (24*i)-th to (24*i+23)-th
+ *              bit of 2/pi after binary point. The corresponding
+ *              floating value is
+ *
+ *                      ipio2[i] * 2^(-24(i+1)).
+ *
+ * External function:
+ *      double scalbn(), floor();
+ *
+ *
+ * Here is the description of some local variables:
+ *
+ *      jk      jk+1 is the initial number of terms of ipio2[] needed
+ *              in the computation. The recommended value is 2,3,4,
+ *              6 for single, double, extended,and quad.
+ *
+ *      jz      local integer variable indicating the number of
+ *              terms of ipio2[] used.
+ *
+ *      jx      nx - 1
+ *
+ *      jv      index for pointing to the suitable ipio2[] for the
+ *              computation. In general, we want
+ *                      ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
+ *              is an integer. Thus
+ *                      e0-3-24*jv >= 0 or (e0-3)/24 >= jv
+ *              Hence jv = max(0,(e0-3)/24).
+ *
+ *      jp      jp+1 is the number of terms in PIo2[] needed, jp = jk.
+ *
+ *      q[]     double array with integral value, representing the
+ *              24-bits chunk of the product of x and 2/pi.
+ *
+ *      q0      the corresponding exponent of q[0]. Note that the
+ *              exponent for q[i] would be q0-24*i.
+ *
+ *      PIo2[]  double precision array, obtained by cutting pi/2
+ *              into 24 bits chunks.
+ *
+ *      f[]     ipio2[] in floating point
+ *
+ *      iq[]    integer array by breaking up q[] in 24-bits chunk.
+ *
+ *      fq[]    final product of x*(2/pi) in fq[0],..,fq[jk]
+ *
+ *      ih      integer. If >0 it indicates q[] is >= 0.5, hence
+ *              it also indicates the *sign* of the result.
+ *
+ */
+int __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec,
+                      const int32_t *ipio2) {
+  /* Constants:
+   * The hexadecimal values are the intended ones for the following
+   * constants. The decimal values may be used, provided that the
+   * compiler will convert from decimal to binary accurately enough
+   * to produce the hexadecimal values shown.
+   */
+  static const int init_jk[] = {2, 3, 4, 6}; /* initial value for jk */
+
+  static const double PIo2[] = {
+      1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
+      7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
+      5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
+      3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */
+      1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */
+      1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
+      2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
+      2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
+  };
+
+  static const double
+      zero = 0.0,
+      one = 1.0,
+      two24 = 1.67772160000000000000e+07,  /* 0x41700000, 0x00000000 */
+      twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
+
+  int32_t jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih;
+  double z, fw, f[20], fq[20], q[20];
+
+  /* initialize jk*/
+  jk = init_jk[prec];
+  jp = jk;
+
+  /* determine jx,jv,q0, note that 3>q0 */
+  jx = nx - 1;
+  jv = (e0 - 3) / 24;
+  if (jv < 0) jv = 0;
+  q0 = e0 - 24 * (jv + 1);
+
+  /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
+  j = jv - jx;
+  m = jx + jk;
+  for (i = 0; i <= m; i++, j++) {
+    f[i] = (j < 0) ? zero : static_cast<double>(ipio2[j]);
+  }
+
+  /* compute q[0],q[1],...q[jk] */
+  for (i = 0; i <= jk; i++) {
+    for (j = 0, fw = 0.0; j <= jx; j++) fw += x[j] * f[jx + i - j];
+    q[i] = fw;
+  }
+
+  jz = jk;
+recompute:
+  /* distill q[] into iq[] reversingly */
+  for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--) {
+    fw = static_cast<double>(static_cast<int32_t>(twon24 * z));
+    iq[i] = static_cast<int32_t>(z - two24 * fw);
+    z = q[j - 1] + fw;
+  }
+
+  /* compute n */
+  z = scalbn(z, q0);           /* actual value of z */
+  z -= 8.0 * floor(z * 0.125); /* trim off integer >= 8 */
+  n = static_cast<int32_t>(z);
+  z -= static_cast<double>(n);
+  ih = 0;
+  if (q0 > 0) { /* need iq[jz-1] to determine n */
+    i = (iq[jz - 1] >> (24 - q0));
+    n += i;
+    iq[jz - 1] -= i << (24 - q0);
+    ih = iq[jz - 1] >> (23 - q0);
+  } else if (q0 == 0) {
+    ih = iq[jz - 1] >> 23;
+  } else if (z >= 0.5) {
+    ih = 2;
+  }
+
+  if (ih > 0) { /* q > 0.5 */
+    n += 1;
+    carry = 0;
+    for (i = 0; i < jz; i++) { /* compute 1-q */
+      j = iq[i];
+      if (carry == 0) {
+        if (j != 0) {
+          carry = 1;
+          iq[i] = 0x1000000 - j;
+        }
+      } else {
+        iq[i] = 0xffffff - j;
+      }
+    }
+    if (q0 > 0) { /* rare case: chance is 1 in 12 */
+      switch (q0) {
+        case 1:
+          iq[jz - 1] &= 0x7fffff;
+          break;
+        case 2:
+          iq[jz - 1] &= 0x3fffff;
+          break;
+      }
+    }
+    if (ih == 2) {
+      z = one - z;
+      if (carry != 0) z -= scalbn(one, q0);
+    }
+  }
+
+  /* check if recomputation is needed */
+  if (z == zero) {
+    j = 0;
+    for (i = jz - 1; i >= jk; i--) j |= iq[i];
+    if (j == 0) { /* need recomputation */
+      for (k = 1; jk >= k && iq[jk - k] == 0; k++) {
+        /* k = no. of terms needed */
+      }
+
+      for (i = jz + 1; i <= jz + k; i++) { /* add q[jz+1] to q[jz+k] */
+        f[jx + i] = ipio2[jv + i];
+        for (j = 0, fw = 0.0; j <= jx; j++) fw += x[j] * f[jx + i - j];
+        q[i] = fw;
+      }
+      jz += k;
+      goto recompute;
+    }
+  }
+
+  /* chop off zero terms */
+  if (z == 0.0) {
+    jz -= 1;
+    q0 -= 24;
+    while (iq[jz] == 0) {
+      jz--;
+      q0 -= 24;
+    }
+  } else { /* break z into 24-bit if necessary */
+    z = scalbn(z, -q0);
+    if (z >= two24) {
+      fw = static_cast<double>(static_cast<int32_t>(twon24 * z));
+      iq[jz] = z - two24 * fw;
+      jz += 1;
+      q0 += 24;
+      iq[jz] = fw;
+    } else {
+      iq[jz] = z;
+    }
+  }
+
+  /* convert integer "bit" chunk to floating-point value */
+  fw = scalbn(one, q0);
+  for (i = jz; i >= 0; i--) {
+    q[i] = fw * iq[i];
+    fw *= twon24;
+  }
+
+  /* compute PIo2[0,...,jp]*q[jz,...,0] */
+  for (i = jz; i >= 0; i--) {
+    for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++) fw += PIo2[k] * q[i + k];
+    fq[jz - i] = fw;
+  }
+
+  /* compress fq[] into y[] */
+  switch (prec) {
+    case 0:
+      fw = 0.0;
+      for (i = jz; i >= 0; i--) fw += fq[i];
+      y[0] = (ih == 0) ? fw : -fw;
+      break;
+    case 1:
+    case 2:
+      fw = 0.0;
+      for (i = jz; i >= 0; i--) fw += fq[i];
+      y[0] = (ih == 0) ? fw : -fw;
+      fw = fq[0] - fw;
+      for (i = 1; i <= jz; i++) fw += fq[i];
+      y[1] = (ih == 0) ? fw : -fw;
+      break;
+    case 3: /* painful */
+      for (i = jz; i > 0; i--) {
+        fw = fq[i - 1] + fq[i];
+        fq[i] += fq[i - 1] - fw;
+        fq[i - 1] = fw;
+      }
+      for (i = jz; i > 1; i--) {
+        fw = fq[i - 1] + fq[i];
+        fq[i] += fq[i - 1] - fw;
+        fq[i - 1] = fw;
+      }
+      for (fw = 0.0, i = jz; i >= 2; i--) fw += fq[i];
+      if (ih == 0) {
+        y[0] = fq[0];
+        y[1] = fq[1];
+        y[2] = fw;
+      } else {
+        y[0] = -fq[0];
+        y[1] = -fq[1];
+        y[2] = -fw;
+      }
+  }
+  return n & 7;
+}
+
+/* __kernel_sin( x, y, iy)
+ * kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854
+ * Input x is assumed to be bounded by ~pi/4 in magnitude.
+ * Input y is the tail of x.
+ * Input iy indicates whether y is 0. (if iy=0, y assume to be 0).
+ *
+ * Algorithm
+ *      1. Since sin(-x) = -sin(x), we need only to consider positive x.
+ *      2. if x < 2^-27 (hx<0x3e400000 0), return x with inexact if x!=0.
+ *      3. sin(x) is approximated by a polynomial of degree 13 on
+ *         [0,pi/4]
+ *                               3            13
+ *              sin(x) ~ x + S1*x + ... + S6*x
+ *         where
+ *
+ *      |sin(x)         2     4     6     8     10     12  |     -58
+ *      |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x  +S6*x   )| <= 2
+ *      |  x                                               |
+ *
+ *      4. sin(x+y) = sin(x) + sin'(x')*y
+ *                  ~ sin(x) + (1-x*x/2)*y
+ *         For better accuracy, let
+ *                   3      2      2      2      2
+ *              r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6))))
+ *         then                   3    2
+ *              sin(x) = x + (S1*x + (x *(r-y/2)+y))
+ */
+V8_INLINE double __kernel_sin(double x, double y, int iy) {
+  static const double
+      half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
+      S1 = -1.66666666666666324348e-01,  /* 0xBFC55555, 0x55555549 */
+      S2 = 8.33333333332248946124e-03,   /* 0x3F811111, 0x1110F8A6 */
+      S3 = -1.98412698298579493134e-04,  /* 0xBF2A01A0, 0x19C161D5 */
+      S4 = 2.75573137070700676789e-06,   /* 0x3EC71DE3, 0x57B1FE7D */
+      S5 = -2.50507602534068634195e-08,  /* 0xBE5AE5E6, 0x8A2B9CEB */
+      S6 = 1.58969099521155010221e-10;   /* 0x3DE5D93A, 0x5ACFD57C */
+
+  double z, r, v;
+  int32_t ix;
+  GET_HIGH_WORD(ix, x);
+  ix &= 0x7fffffff;      /* high word of x */
+  if (ix < 0x3e400000) { /* |x| < 2**-27 */
+    if (static_cast<int>(x) == 0) return x;
+  } /* generate inexact */
+  z = x * x;
+  v = z * x;
+  r = S2 + z * (S3 + z * (S4 + z * (S5 + z * S6)));
+  if (iy == 0) {
+    return x + v * (S1 + z * r);
+  } else {
+    return x - ((z * (half * y - v * r) - y) - v * S1);
+  }
+}
+
+/* __kernel_tan( x, y, k )
+ * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
+ * Input x is assumed to be bounded by ~pi/4 in magnitude.
+ * Input y is the tail of x.
+ * Input k indicates whether tan (if k=1) or
+ * -1/tan (if k= -1) is returned.
+ *
+ * Algorithm
+ *      1. Since tan(-x) = -tan(x), we need only to consider positive x.
+ *      2. if x < 2^-28 (hx<0x3e300000 0), return x with inexact if x!=0.
+ *      3. tan(x) is approximated by a odd polynomial of degree 27 on
+ *         [0,0.67434]
+ *                               3             27
+ *              tan(x) ~ x + T1*x + ... + T13*x
+ *         where
+ *
+ *              |tan(x)         2     4            26   |     -59.2
+ *              |----- - (1+T1*x +T2*x +.... +T13*x    )| <= 2
+ *              |  x                                    |
+ *
+ *         Note: tan(x+y) = tan(x) + tan'(x)*y
+ *                        ~ tan(x) + (1+x*x)*y
+ *         Therefore, for better accuracy in computing tan(x+y), let
+ *                   3      2      2       2       2
+ *              r = x *(T2+x *(T3+x *(...+x *(T12+x *T13))))
+ *         then
+ *                                  3    2
+ *              tan(x+y) = x + (T1*x + (x *(r+y)+y))
+ *
+ *      4. For x in [0.67434,pi/4],  let y = pi/4 - x, then
+ *              tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
+ *                     = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
+ */
+double __kernel_tan(double x, double y, int iy) {
+  static const double xxx[] = {
+      3.33333333333334091986e-01,             /* 3FD55555, 55555563 */
+      1.33333333333201242699e-01,             /* 3FC11111, 1110FE7A */
+      5.39682539762260521377e-02,             /* 3FABA1BA, 1BB341FE */
+      2.18694882948595424599e-02,             /* 3F9664F4, 8406D637 */
+      8.86323982359930005737e-03,             /* 3F8226E3, E96E8493 */
+      3.59207910759131235356e-03,             /* 3F6D6D22, C9560328 */
+      1.45620945432529025516e-03,             /* 3F57DBC8, FEE08315 */
+      5.88041240820264096874e-04,             /* 3F4344D8, F2F26501 */
+      2.46463134818469906812e-04,             /* 3F3026F7, 1A8D1068 */
+      7.81794442939557092300e-05,             /* 3F147E88, A03792A6 */
+      7.14072491382608190305e-05,             /* 3F12B80F, 32F0A7E9 */
+      -1.85586374855275456654e-05,            /* BEF375CB, DB605373 */
+      2.59073051863633712884e-05,             /* 3EFB2A70, 74BF7AD4 */
+      /* one */ 1.00000000000000000000e+00,   /* 3FF00000, 00000000 */
+      /* pio4 */ 7.85398163397448278999e-01,  /* 3FE921FB, 54442D18 */
+      /* pio4lo */ 3.06161699786838301793e-17 /* 3C81A626, 33145C07 */
+  };
+#define one xxx[13]
+#define pio4 xxx[14]
+#define pio4lo xxx[15]
+#define T xxx
+
+  double z, r, v, w, s;
+  int32_t ix, hx;
+
+  GET_HIGH_WORD(hx, x);             /* high word of x */
+  ix = hx & 0x7fffffff;             /* high word of |x| */
+  if (ix < 0x3e300000) {            /* x < 2**-28 */
+    if (static_cast<int>(x) == 0) { /* generate inexact */
+      uint32_t low;
+      GET_LOW_WORD(low, x);
+      if (((ix | low) | (iy + 1)) == 0) {
+        return one / fabs(x);
+      } else {
+        if (iy == 1) {
+          return x;
+        } else { /* compute -1 / (x+y) carefully */
+          double a, t;
+
+          z = w = x + y;
+          SET_LOW_WORD(z, 0);
+          v = y - (z - x);
+          t = a = -one / w;
+          SET_LOW_WORD(t, 0);
+          s = one + t * z;
+          return t + a * (s + t * v);
+        }
+      }
+    }
+  }
+  if (ix >= 0x3FE59428) { /* |x| >= 0.6744 */
+    if (hx < 0) {
+      x = -x;
+      y = -y;
+    }
+    z = pio4 - x;
+    w = pio4lo - y;
+    x = z + w;
+    y = 0.0;
+  }
+  z = x * x;
+  w = z * z;
+  /*
+   * Break x^5*(T[1]+x^2*T[2]+...) into
+   * x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
+   * x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
+   */
+  r = T[1] + w * (T[3] + w * (T[5] + w * (T[7] + w * (T[9] + w * T[11]))));
+  v = z *
+      (T[2] + w * (T[4] + w * (T[6] + w * (T[8] + w * (T[10] + w * T[12])))));
+  s = z * x;
+  r = y + z * (s * (r + v) + y);
+  r += T[0] * s;
+  w = x + r;
+  if (ix >= 0x3FE59428) {
+    v = iy;
+    return (1 - ((hx >> 30) & 2)) * (v - 2.0 * (x - (w * w / (w + v) - r)));
+  }
+  if (iy == 1) {
+    return w;
+  } else {
+    /*
+     * if allow error up to 2 ulp, simply return
+     * -1.0 / (x+r) here
+     */
+    /* compute -1.0 / (x+r) accurately */
+    double a, t;
+    z = w;
+    SET_LOW_WORD(z, 0);
+    v = r - (z - x);  /* z+v = r+x */
+    t = a = -1.0 / w; /* a = -1.0/w */
+    SET_LOW_WORD(t, 0);
+    s = 1.0 + t * z;
+    return t + a * (s + t * v);
+  }
+
+#undef one
+#undef pio4
+#undef pio4lo
+#undef T
+}
+
+}  // namespace
+
+/* atan(x)
+ * Method
+ *   1. Reduce x to positive by atan(x) = -atan(-x).
+ *   2. According to the integer k=4t+0.25 chopped, t=x, the argument
+ *      is further reduced to one of the following intervals and the
+ *      arctangent of t is evaluated by the corresponding formula:
+ *
+ *      [0,7/16]      atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
+ *      [7/16,11/16]  atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) )
+ *      [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) )
+ *      [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) )
+ *      [39/16,INF]   atan(x) = atan(INF) + atan( -1/t )
+ *
+ * Constants:
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
+ * to produce the hexadecimal values shown.
+ */
+double atan(double x) {
+  static const double atanhi[] = {
+      4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */
+      7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */
+      9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */
+      1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */
+  };
+
+  static const double atanlo[] = {
+      2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */
+      3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */
+      1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */
+      6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */
+  };
+
+  static const double aT[] = {
+      3.33333333333329318027e-01,  /* 0x3FD55555, 0x5555550D */
+      -1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */
+      1.42857142725034663711e-01,  /* 0x3FC24924, 0x920083FF */
+      -1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */
+      9.09088713343650656196e-02,  /* 0x3FB745CD, 0xC54C206E */
+      -7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */
+      6.66107313738753120669e-02,  /* 0x3FB10D66, 0xA0D03D51 */
+      -5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */
+      4.97687799461593236017e-02,  /* 0x3FA97B4B, 0x24760DEB */
+      -3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */
+      1.62858201153657823623e-02,  /* 0x3F90AD3A, 0xE322DA11 */
+  };
+
+  static const double one = 1.0, huge = 1.0e300;
+
+  double w, s1, s2, z;
+  int32_t ix, hx, id;
+
+  GET_HIGH_WORD(hx, x);
+  ix = hx & 0x7fffffff;
+  if (ix >= 0x44100000) { /* if |x| >= 2^66 */
+    uint32_t low;
+    GET_LOW_WORD(low, x);
+    if (ix > 0x7ff00000 || (ix == 0x7ff00000 && (low != 0)))
+      return x + x; /* NaN */
+    if (hx > 0)
+      return atanhi[3] + *(volatile double *)&atanlo[3];
+    else
+      return -atanhi[3] - *(volatile double *)&atanlo[3];
+  }
+  if (ix < 0x3fdc0000) {            /* |x| < 0.4375 */
+    if (ix < 0x3e400000) {          /* |x| < 2^-27 */
+      if (huge + x > one) return x; /* raise inexact */
+    }
+    id = -1;
+  } else {
+    x = fabs(x);
+    if (ix < 0x3ff30000) {   /* |x| < 1.1875 */
+      if (ix < 0x3fe60000) { /* 7/16 <=|x|<11/16 */
+        id = 0;
+        x = (2.0 * x - one) / (2.0 + x);
+      } else { /* 11/16<=|x|< 19/16 */
+        id = 1;
+        x = (x - one) / (x + one);
+      }
+    } else {
+      if (ix < 0x40038000) { /* |x| < 2.4375 */
+        id = 2;
+        x = (x - 1.5) / (one + 1.5 * x);
+      } else { /* 2.4375 <= |x| < 2^66 */
+        id = 3;
+        x = -1.0 / x;
+      }
+    }
+  }
+  /* end of argument reduction */
+  z = x * x;
+  w = z * z;
+  /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
+  s1 = z * (aT[0] +
+            w * (aT[2] + w * (aT[4] + w * (aT[6] + w * (aT[8] + w * aT[10])))));
+  s2 = w * (aT[1] + w * (aT[3] + w * (aT[5] + w * (aT[7] + w * aT[9]))));
+  if (id < 0) {
+    return x - x * (s1 + s2);
+  } else {
+    z = atanhi[id] - ((x * (s1 + s2) - atanlo[id]) - x);
+    return (hx < 0) ? -z : z;
+  }
+}
+
+/* atan2(y,x)
+ * Method :
+ *  1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
+ *  2. Reduce x to positive by (if x and y are unexceptional):
+ *    ARG (x+iy) = arctan(y/x)       ... if x > 0,
+ *    ARG (x+iy) = pi - arctan[y/(-x)]   ... if x < 0,
+ *
+ * Special cases:
+ *
+ *  ATAN2((anything), NaN ) is NaN;
+ *  ATAN2(NAN , (anything) ) is NaN;
+ *  ATAN2(+-0, +(anything but NaN)) is +-0  ;
+ *  ATAN2(+-0, -(anything but NaN)) is +-pi ;
+ *  ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
+ *  ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
+ *  ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
+ *  ATAN2(+-INF,+INF ) is +-pi/4 ;
+ *  ATAN2(+-INF,-INF ) is +-3pi/4;
+ *  ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
+ *
+ * Constants:
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
+ * to produce the hexadecimal values shown.
+ */
+double atan2(double y, double x) {
+  static volatile double tiny = 1.0e-300;
+  static const double
+      zero = 0.0,
+      pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */
+      pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */
+      pi = 3.1415926535897931160E+00;     /* 0x400921FB, 0x54442D18 */
+  static volatile double pi_lo =
+      1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
+
+  double z;
+  int32_t k, m, hx, hy, ix, iy;
+  uint32_t lx, ly;
+
+  EXTRACT_WORDS(hx, lx, x);
+  ix = hx & 0x7fffffff;
+  EXTRACT_WORDS(hy, ly, y);
+  iy = hy & 0x7fffffff;
+  if (((ix | ((lx | -static_cast<int32_t>(lx)) >> 31)) > 0x7ff00000) ||
+      ((iy | ((ly | -static_cast<int32_t>(ly)) >> 31)) > 0x7ff00000)) {
+    return x + y; /* x or y is NaN */
+  }
+  if (((hx - 0x3ff00000) | lx) == 0) return atan(y); /* x=1.0 */
+  m = ((hy >> 31) & 1) | ((hx >> 30) & 2);           /* 2*sign(x)+sign(y) */
+
+  /* when y = 0 */
+  if ((iy | ly) == 0) {
+    switch (m) {
+      case 0:
+      case 1:
+        return y; /* atan(+-0,+anything)=+-0 */
+      case 2:
+        return pi + tiny; /* atan(+0,-anything) = pi */
+      case 3:
+        return -pi - tiny; /* atan(-0,-anything) =-pi */
+    }
+  }
+  /* when x = 0 */
+  if ((ix | lx) == 0) return (hy < 0) ? -pi_o_2 - tiny : pi_o_2 + tiny;
+
+  /* when x is INF */
+  if (ix == 0x7ff00000) {
+    if (iy == 0x7ff00000) {
+      switch (m) {
+        case 0:
+          return pi_o_4 + tiny; /* atan(+INF,+INF) */
+        case 1:
+          return -pi_o_4 - tiny; /* atan(-INF,+INF) */
+        case 2:
+          return 3.0 * pi_o_4 + tiny; /*atan(+INF,-INF)*/
+        case 3:
+          return -3.0 * pi_o_4 - tiny; /*atan(-INF,-INF)*/
+      }
+    } else {
+      switch (m) {
+        case 0:
+          return zero; /* atan(+...,+INF) */
+        case 1:
+          return -zero; /* atan(-...,+INF) */
+        case 2:
+          return pi + tiny; /* atan(+...,-INF) */
+        case 3:
+          return -pi - tiny; /* atan(-...,-INF) */
+      }
+    }
+  }
+  /* when y is INF */
+  if (iy == 0x7ff00000) return (hy < 0) ? -pi_o_2 - tiny : pi_o_2 + tiny;
+
+  /* compute y/x */
+  k = (iy - ix) >> 20;
+  if (k > 60) { /* |y/x| >  2**60 */
+    z = pi_o_2 + 0.5 * pi_lo;
+    m &= 1;
+  } else if (hx < 0 && k < -60) {
+    z = 0.0; /* 0 > |y|/x > -2**-60 */
+  } else {
+    z = atan(fabs(y / x)); /* safe to do y/x */
+  }
+  switch (m) {
+    case 0:
+      return z; /* atan(+,+) */
+    case 1:
+      return -z; /* atan(-,+) */
+    case 2:
+      return pi - (z - pi_lo); /* atan(+,-) */
+    default:                   /* case 3 */
+      return (z - pi_lo) - pi; /* atan(-,-) */
+  }
+}
+
+/* cos(x)
+ * Return cosine function of x.
+ *
+ * kernel function:
+ *      __kernel_sin            ... sine function on [-pi/4,pi/4]
+ *      __kernel_cos            ... cosine function on [-pi/4,pi/4]
+ *      __ieee754_rem_pio2      ... argument reduction routine
+ *
+ * Method.
+ *      Let S,C and T denote the sin, cos and tan respectively on
+ *      [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
+ *      in [-pi/4 , +pi/4], and let n = k mod 4.
+ *      We have
+ *
+ *          n        sin(x)      cos(x)        tan(x)
+ *     ----------------------------------------------------------
+ *          0          S           C             T
+ *          1          C          -S            -1/T
+ *          2         -S          -C             T
+ *          3         -C           S            -1/T
+ *     ----------------------------------------------------------
+ *
+ * Special cases:
+ *      Let trig be any of sin, cos, or tan.
+ *      trig(+-INF)  is NaN, with signals;
+ *      trig(NaN)    is that NaN;
+ *
+ * Accuracy:
+ *      TRIG(x) returns trig(x) nearly rounded
+ */
+double cos(double x) {
+  double y[2], z = 0.0;
+  int32_t n, ix;
+
+  /* High word of x. */
+  GET_HIGH_WORD(ix, x);
+
+  /* |x| ~< pi/4 */
+  ix &= 0x7fffffff;
+  if (ix <= 0x3fe921fb) {
+    return __kernel_cos(x, z);
+  } else if (ix >= 0x7ff00000) {
+    /* cos(Inf or NaN) is NaN */
+    return x - x;
+  } else {
+    /* argument reduction needed */
+    n = __ieee754_rem_pio2(x, y);
+    switch (n & 3) {
+      case 0:
+        return __kernel_cos(y[0], y[1]);
+      case 1:
+        return -__kernel_sin(y[0], y[1], 1);
+      case 2:
+        return -__kernel_cos(y[0], y[1]);
+      default:
+        return __kernel_sin(y[0], y[1], 1);
+    }
+  }
+}
+
+/* exp(x)
+ * Returns the exponential of x.
+ *
+ * Method
+ *   1. Argument reduction:
+ *      Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
+ *      Given x, find r and integer k such that
+ *
+ *               x = k*ln2 + r,  |r| <= 0.5*ln2.
+ *
+ *      Here r will be represented as r = hi-lo for better
+ *      accuracy.
+ *
+ *   2. Approximation of exp(r) by a special rational function on
+ *      the interval [0,0.34658]:
+ *      Write
+ *          R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
+ *      We use a special Remes algorithm on [0,0.34658] to generate
+ *      a polynomial of degree 5 to approximate R. The maximum error
+ *      of this polynomial approximation is bounded by 2**-59. In
+ *      other words,
+ *          R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
+ *      (where z=r*r, and the values of P1 to P5 are listed below)
+ *      and
+ *          |                  5          |     -59
+ *          | 2.0+P1*z+...+P5*z   -  R(z) | <= 2
+ *          |                             |
+ *      The computation of exp(r) thus becomes
+ *                             2*r
+ *              exp(r) = 1 + -------
+ *                            R - r
+ *                                 r*R1(r)
+ *                     = 1 + r + ----------- (for better accuracy)
+ *                                2 - R1(r)
+ *      where
+ *                               2       4             10
+ *              R1(r) = r - (P1*r  + P2*r  + ... + P5*r   ).
+ *
+ *   3. Scale back to obtain exp(x):
+ *      From step 1, we have
+ *         exp(x) = 2^k * exp(r)
+ *
+ * Special cases:
+ *      exp(INF) is INF, exp(NaN) is NaN;
+ *      exp(-INF) is 0, and
+ *      for finite argument, only exp(0)=1 is exact.
+ *
+ * Accuracy:
+ *      according to an error analysis, the error is always less than
+ *      1 ulp (unit in the last place).
+ *
+ * Misc. info.
+ *      For IEEE double
+ *          if x >  7.09782712893383973096e+02 then exp(x) overflow
+ *          if x < -7.45133219101941108420e+02 then exp(x) underflow
+ *
+ * Constants:
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
+ * to produce the hexadecimal values shown.
+ */
+double exp(double x) {
+  static const double
+      one = 1.0,
+      halF[2] = {0.5, -0.5},
+      o_threshold = 7.09782712893383973096e+02,  /* 0x40862E42, 0xFEFA39EF */
+      u_threshold = -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */
+      ln2HI[2] = {6.93147180369123816490e-01,    /* 0x3fe62e42, 0xfee00000 */
+                  -6.93147180369123816490e-01},  /* 0xbfe62e42, 0xfee00000 */
+      ln2LO[2] = {1.90821492927058770002e-10,    /* 0x3dea39ef, 0x35793c76 */
+                  -1.90821492927058770002e-10},  /* 0xbdea39ef, 0x35793c76 */
+      invln2 = 1.44269504088896338700e+00,       /* 0x3ff71547, 0x652b82fe */
+      P1 = 1.66666666666666019037e-01,           /* 0x3FC55555, 0x5555553E */
+      P2 = -2.77777777770155933842e-03,          /* 0xBF66C16C, 0x16BEBD93 */
+      P3 = 6.61375632143793436117e-05,           /* 0x3F11566A, 0xAF25DE2C */
+      P4 = -1.65339022054652515390e-06,          /* 0xBEBBBD41, 0xC5D26BF1 */
+      P5 = 4.13813679705723846039e-08,           /* 0x3E663769, 0x72BEA4D0 */
+      E = 2.718281828459045;                     /* 0x4005bf0a, 0x8b145769 */
+
+  static volatile double
+      huge = 1.0e+300,
+      twom1000 = 9.33263618503218878990e-302, /* 2**-1000=0x01700000,0*/
+      two1023 = 8.988465674311579539e307;     /* 0x1p1023 */
+
+  double y, hi = 0.0, lo = 0.0, c, t, twopk;
+  int32_t k = 0, xsb;
+  uint32_t hx;
+
+  GET_HIGH_WORD(hx, x);
+  xsb = (hx >> 31) & 1; /* sign bit of x */
+  hx &= 0x7fffffff;     /* high word of |x| */
+
+  /* filter out non-finite argument */
+  if (hx >= 0x40862E42) { /* if |x|>=709.78... */
+    if (hx >= 0x7ff00000) {
+      uint32_t lx;
+      GET_LOW_WORD(lx, x);
+      if (((hx & 0xfffff) | lx) != 0)
+        return x + x; /* NaN */
+      else
+        return (xsb == 0) ? x : 0.0; /* exp(+-inf)={inf,0} */
+    }
+    if (x > o_threshold) return huge * huge;         /* overflow */
+    if (x < u_threshold) return twom1000 * twom1000; /* underflow */
+  }
+
+  /* argument reduction */
+  if (hx > 0x3fd62e42) {   /* if  |x| > 0.5 ln2 */
+    if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
+      /* TODO(rtoy): We special case exp(1) here to return the correct
+       * value of E, as the computation below would get the last bit
+       * wrong. We should probably fix the algorithm instead.
+       */
+      if (x == 1.0) return E;
+      hi = x - ln2HI[xsb];
+      lo = ln2LO[xsb];
+      k = 1 - xsb - xsb;
+    } else {
+      k = static_cast<int>(invln2 * x + halF[xsb]);
+      t = k;
+      hi = x - t * ln2HI[0]; /* t*ln2HI is exact here */
+      lo = t * ln2LO[0];
+    }
+    STRICT_ASSIGN(double, x, hi - lo);
+  } else if (hx < 0x3e300000) {         /* when |x|<2**-28 */
+    if (huge + x > one) return one + x; /* trigger inexact */
+  } else {
+    k = 0;
+  }
+
+  /* x is now in primary range */
+  t = x * x;
+  if (k >= -1021) {
+    INSERT_WORDS(twopk, 0x3ff00000 + (k << 20), 0);
+  } else {
+    INSERT_WORDS(twopk, 0x3ff00000 + ((k + 1000) << 20), 0);
+  }
+  c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
+  if (k == 0) {
+    return one - ((x * c) / (c - 2.0) - x);
+  } else {
+    y = one - ((lo - (x * c) / (2.0 - c)) - hi);
+  }
+  if (k >= -1021) {
+    if (k == 1024) return y * 2.0 * two1023;
+    return y * twopk;
+  } else {
+    return y * twopk * twom1000;
+  }
+}
+
+/*
+ * Method :
+ *    1.Reduced x to positive by atanh(-x) = -atanh(x)
+ *    2.For x>=0.5
+ *              1              2x                          x
+ *  atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
+ *              2             1 - x                      1 - x
+ *
+ *   For x<0.5
+ *  atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
+ *
+ * Special cases:
+ *  atanh(x) is NaN if |x| > 1 with signal;
+ *  atanh(NaN) is that NaN with no signal;
+ *  atanh(+-1) is +-INF with signal.
+ *
+ */
+double atanh(double x) {
+  static const double one = 1.0, huge = 1e300;
+  static const double zero = 0.0;
+
+  double t;
+  int32_t hx, ix;
+  uint32_t lx;
+  EXTRACT_WORDS(hx, lx, x);
+  ix = hx & 0x7fffffff;
+  if ((ix | ((lx | -static_cast<int32_t>(lx)) >> 31)) > 0x3ff00000) /* |x|>1 */
+    return (x - x) / (x - x);
+  if (ix == 0x3ff00000) return x / zero;
+  if (ix < 0x3e300000 && (huge + x) > zero) return x; /* x<2**-28 */
+  SET_HIGH_WORD(x, ix);
+  if (ix < 0x3fe00000) { /* x < 0.5 */
+    t = x + x;
+    t = 0.5 * log1p(t + t * x / (one - x));
+  } else {
+    t = 0.5 * log1p((x + x) / (one - x));
+  }
+  if (hx >= 0)
+    return t;
+  else
+    return -t;
+}
+
+/* log(x)
+ * Return the logrithm of x
+ *
+ * Method :
+ *   1. Argument Reduction: find k and f such that
+ *     x = 2^k * (1+f),
+ *     where  sqrt(2)/2 < 1+f < sqrt(2) .
+ *
+ *   2. Approximation of log(1+f).
+ *  Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
+ *     = 2s + 2/3 s**3 + 2/5 s**5 + .....,
+ *         = 2s + s*R
+ *      We use a special Reme algorithm on [0,0.1716] to generate
+ *  a polynomial of degree 14 to approximate R The maximum error
+ *  of this polynomial approximation is bounded by 2**-58.45. In
+ *  other words,
+ *            2      4      6      8      10      12      14
+ *      R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s  +Lg6*s  +Lg7*s
+ *    (the values of Lg1 to Lg7 are listed in the program)
+ *  and
+ *      |      2          14          |     -58.45
+ *      | Lg1*s +...+Lg7*s    -  R(z) | <= 2
+ *      |                             |
+ *  Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
+ *  In order to guarantee error in log below 1ulp, we compute log
+ *  by
+ *    log(1+f) = f - s*(f - R)  (if f is not too large)
+ *    log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy)
+ *
+ *  3. Finally,  log(x) = k*ln2 + log(1+f).
+ *          = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
+ *     Here ln2 is split into two floating point number:
+ *      ln2_hi + ln2_lo,
+ *     where n*ln2_hi is always exact for |n| < 2000.
+ *
+ * Special cases:
+ *  log(x) is NaN with signal if x < 0 (including -INF) ;
+ *  log(+INF) is +INF; log(0) is -INF with signal;
+ *  log(NaN) is that NaN with no signal.
+ *
+ * Accuracy:
+ *  according to an error analysis, the error is always less than
+ *  1 ulp (unit in the last place).
+ *
+ * Constants:
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
+ * to produce the hexadecimal values shown.
+ */
+double log(double x) {
+  static const double                      /* -- */
+      ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */
+      ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */
+      two54 = 1.80143985094819840000e+16,  /* 43500000 00000000 */
+      Lg1 = 6.666666666666735130e-01,      /* 3FE55555 55555593 */
+      Lg2 = 3.999999999940941908e-01,      /* 3FD99999 9997FA04 */
+      Lg3 = 2.857142874366239149e-01,      /* 3FD24924 94229359 */
+      Lg4 = 2.222219843214978396e-01,      /* 3FCC71C5 1D8E78AF */
+      Lg5 = 1.818357216161805012e-01,      /* 3FC74664 96CB03DE */
+      Lg6 = 1.531383769920937332e-01,      /* 3FC39A09 D078C69F */
+      Lg7 = 1.479819860511658591e-01;      /* 3FC2F112 DF3E5244 */
+
+  static const double zero = 0.0;
+  static volatile double vzero = 0.0;
+
+  double hfsq, f, s, z, R, w, t1, t2, dk;
+  int32_t k, hx, i, j;
+  uint32_t lx;
+
+  EXTRACT_WORDS(hx, lx, x);
+
+  k = 0;
+  if (hx < 0x00100000) { /* x < 2**-1022  */
+    if (((hx & 0x7fffffff) | lx) == 0)
+      return -two54 / vzero;           /* log(+-0)=-inf */
+    if (hx < 0) return (x - x) / zero; /* log(-#) = NaN */
+    k -= 54;
+    x *= two54; /* subnormal number, scale up x */
+    GET_HIGH_WORD(hx, x);
+  }
+  if (hx >= 0x7ff00000) return x + x;
+  k += (hx >> 20) - 1023;
+  hx &= 0x000fffff;
+  i = (hx + 0x95f64) & 0x100000;
+  SET_HIGH_WORD(x, hx | (i ^ 0x3ff00000)); /* normalize x or x/2 */
+  k += (i >> 20);
+  f = x - 1.0;
+  if ((0x000fffff & (2 + hx)) < 3) { /* -2**-20 <= f < 2**-20 */
+    if (f == zero) {
+      if (k == 0) {
+        return zero;
+      } else {
+        dk = static_cast<double>(k);
+        return dk * ln2_hi + dk * ln2_lo;
+      }
+    }
+    R = f * f * (0.5 - 0.33333333333333333 * f);
+    if (k == 0) {
+      return f - R;
+    } else {
+      dk = static_cast<double>(k);
+      return dk * ln2_hi - ((R - dk * ln2_lo) - f);
+    }
+  }
+  s = f / (2.0 + f);
+  dk = static_cast<double>(k);
+  z = s * s;
+  i = hx - 0x6147a;
+  w = z * z;
+  j = 0x6b851 - hx;
+  t1 = w * (Lg2 + w * (Lg4 + w * Lg6));
+  t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7)));
+  i |= j;
+  R = t2 + t1;
+  if (i > 0) {
+    hfsq = 0.5 * f * f;
+    if (k == 0)
+      return f - (hfsq - s * (hfsq + R));
+    else
+      return dk * ln2_hi - ((hfsq - (s * (hfsq + R) + dk * ln2_lo)) - f);
+  } else {
+    if (k == 0)
+      return f - s * (f - R);
+    else
+      return dk * ln2_hi - ((s * (f - R) - dk * ln2_lo) - f);
+  }
+}
+
+/* double log1p(double x)
+ *
+ * Method :
+ *   1. Argument Reduction: find k and f such that
+ *      1+x = 2^k * (1+f),
+ *     where  sqrt(2)/2 < 1+f < sqrt(2) .
+ *
+ *      Note. If k=0, then f=x is exact. However, if k!=0, then f
+ *  may not be representable exactly. In that case, a correction
+ *  term is need. Let u=1+x rounded. Let c = (1+x)-u, then
+ *  log(1+x) - log(u) ~ c/u. Thus, we proceed to compute log(u),
+ *  and add back the correction term c/u.
+ *  (Note: when x > 2**53, one can simply return log(x))
+ *
+ *   2. Approximation of log1p(f).
+ *  Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
+ *     = 2s + 2/3 s**3 + 2/5 s**5 + .....,
+ *         = 2s + s*R
+ *      We use a special Reme algorithm on [0,0.1716] to generate
+ *  a polynomial of degree 14 to approximate R The maximum error
+ *  of this polynomial approximation is bounded by 2**-58.45. In
+ *  other words,
+ *            2      4      6      8      10      12      14
+ *      R(z) ~ Lp1*s +Lp2*s +Lp3*s +Lp4*s +Lp5*s  +Lp6*s  +Lp7*s
+ *    (the values of Lp1 to Lp7 are listed in the program)
+ *  and
+ *      |      2          14          |     -58.45
+ *      | Lp1*s +...+Lp7*s    -  R(z) | <= 2
+ *      |                             |
+ *  Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
+ *  In order to guarantee error in log below 1ulp, we compute log
+ *  by
+ *    log1p(f) = f - (hfsq - s*(hfsq+R)).
+ *
+ *  3. Finally, log1p(x) = k*ln2 + log1p(f).
+ *           = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
+ *     Here ln2 is split into two floating point number:
+ *      ln2_hi + ln2_lo,
+ *     where n*ln2_hi is always exact for |n| < 2000.
+ *
+ * Special cases:
+ *  log1p(x) is NaN with signal if x < -1 (including -INF) ;
+ *  log1p(+INF) is +INF; log1p(-1) is -INF with signal;
+ *  log1p(NaN) is that NaN with no signal.
+ *
+ * Accuracy:
+ *  according to an error analysis, the error is always less than
+ *  1 ulp (unit in the last place).
+ *
+ * Constants:
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
+ * to produce the hexadecimal values shown.
+ *
+ * Note: Assuming log() return accurate answer, the following
+ *   algorithm can be used to compute log1p(x) to within a few ULP:
+ *
+ *    u = 1+x;
+ *    if(u==1.0) return x ; else
+ *         return log(u)*(x/(u-1.0));
+ *
+ *   See HP-15C Advanced Functions Handbook, p.193.
+ */
+double log1p(double x) {
+  static const double                      /* -- */
+      ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */
+      ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */
+      two54 = 1.80143985094819840000e+16,  /* 43500000 00000000 */
+      Lp1 = 6.666666666666735130e-01,      /* 3FE55555 55555593 */
+      Lp2 = 3.999999999940941908e-01,      /* 3FD99999 9997FA04 */
+      Lp3 = 2.857142874366239149e-01,      /* 3FD24924 94229359 */
+      Lp4 = 2.222219843214978396e-01,      /* 3FCC71C5 1D8E78AF */
+      Lp5 = 1.818357216161805012e-01,      /* 3FC74664 96CB03DE */
+      Lp6 = 1.531383769920937332e-01,      /* 3FC39A09 D078C69F */
+      Lp7 = 1.479819860511658591e-01;      /* 3FC2F112 DF3E5244 */
+
+  static const double zero = 0.0;
+  static volatile double vzero = 0.0;
+
+  double hfsq, f, c, s, z, R, u;
+  int32_t k, hx, hu, ax;
+
+  GET_HIGH_WORD(hx, x);
+  ax = hx & 0x7fffffff;
+
+  k = 1;
+  if (hx < 0x3FDA827A) {    /* 1+x < sqrt(2)+ */
+    if (ax >= 0x3ff00000) { /* x <= -1.0 */
+      if (x == -1.0)
+        return -two54 / vzero; /* log1p(-1)=+inf */
+      else
+        return (x - x) / (x - x); /* log1p(x<-1)=NaN */
+    }
+    if (ax < 0x3e200000) {    /* |x| < 2**-29 */
+      if (two54 + x > zero    /* raise inexact */
+          && ax < 0x3c900000) /* |x| < 2**-54 */
+        return x;
+      else
+        return x - x * x * 0.5;
+    }
+    if (hx > 0 || hx <= static_cast<int32_t>(0xbfd2bec4)) {
+      k = 0;
+      f = x;
+      hu = 1;
+    } /* sqrt(2)/2- <= 1+x < sqrt(2)+ */
+  }
+  if (hx >= 0x7ff00000) return x + x;
+  if (k != 0) {
+    if (hx < 0x43400000) {
+      STRICT_ASSIGN(double, u, 1.0 + x);
+      GET_HIGH_WORD(hu, u);
+      k = (hu >> 20) - 1023;
+      c = (k > 0) ? 1.0 - (u - x) : x - (u - 1.0); /* correction term */
+      c /= u;
+    } else {
+      u = x;
+      GET_HIGH_WORD(hu, u);
+      k = (hu >> 20) - 1023;
+      c = 0;
+    }
+    hu &= 0x000fffff;
+    /*
+     * The approximation to sqrt(2) used in thresholds is not
+     * critical.  However, the ones used above must give less
+     * strict bounds than the one here so that the k==0 case is
+     * never reached from here, since here we have committed to
+     * using the correction term but don't use it if k==0.
+     */
+    if (hu < 0x6a09e) {                  /* u ~< sqrt(2) */
+      SET_HIGH_WORD(u, hu | 0x3ff00000); /* normalize u */
+    } else {
+      k += 1;
+      SET_HIGH_WORD(u, hu | 0x3fe00000); /* normalize u/2 */
+      hu = (0x00100000 - hu) >> 2;
+    }
+    f = u - 1.0;
+  }
+  hfsq = 0.5 * f * f;
+  if (hu == 0) { /* |f| < 2**-20 */
+    if (f == zero) {
+      if (k == 0) {
+        return zero;
+      } else {
+        c += k * ln2_lo;
+        return k * ln2_hi + c;
+      }
+    }
+    R = hfsq * (1.0 - 0.66666666666666666 * f);
+    if (k == 0)
+      return f - R;
+    else
+      return k * ln2_hi - ((R - (k * ln2_lo + c)) - f);
+  }
+  s = f / (2.0 + f);
+  z = s * s;
+  R = z * (Lp1 +
+           z * (Lp2 + z * (Lp3 + z * (Lp4 + z * (Lp5 + z * (Lp6 + z * Lp7))))));
+  if (k == 0)
+    return f - (hfsq - s * (hfsq + R));
+  else
+    return k * ln2_hi - ((hfsq - (s * (hfsq + R) + (k * ln2_lo + c))) - f);
+}
+
+/*
+ * k_log1p(f):
+ * Return log(1+f) - f for 1+f in ~[sqrt(2)/2, sqrt(2)].
+ *
+ * The following describes the overall strategy for computing
+ * logarithms in base e.  The argument reduction and adding the final
+ * term of the polynomial are done by the caller for increased accuracy
+ * when different bases are used.
+ *
+ * Method :
+ *   1. Argument Reduction: find k and f such that
+ *         x = 2^k * (1+f),
+ *         where  sqrt(2)/2 < 1+f < sqrt(2) .
+ *
+ *   2. Approximation of log(1+f).
+ *      Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
+ *            = 2s + 2/3 s**3 + 2/5 s**5 + .....,
+ *            = 2s + s*R
+ *      We use a special Reme algorithm on [0,0.1716] to generate
+ *      a polynomial of degree 14 to approximate R The maximum error
+ *      of this polynomial approximation is bounded by 2**-58.45. In
+ *      other words,
+ *          2      4      6      8      10      12      14
+ *          R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s  +Lg6*s  +Lg7*s
+ *      (the values of Lg1 to Lg7 are listed in the program)
+ *      and
+ *          |      2          14          |     -58.45
+ *          | Lg1*s +...+Lg7*s    -  R(z) | <= 2
+ *          |                             |
+ *      Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
+ *      In order to guarantee error in log below 1ulp, we compute log
+ *      by
+ *          log(1+f) = f - s*(f - R)            (if f is not too large)
+ *          log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy)
+ *
+ *   3. Finally,  log(x) = k*ln2 + log(1+f).
+ *          = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
+ *      Here ln2 is split into two floating point number:
+ *          ln2_hi + ln2_lo,
+ *      where n*ln2_hi is always exact for |n| < 2000.
+ *
+ * Special cases:
+ *      log(x) is NaN with signal if x < 0 (including -INF) ;
+ *      log(+INF) is +INF; log(0) is -INF with signal;
+ *      log(NaN) is that NaN with no signal.
+ *
+ * Accuracy:
+ *      according to an error analysis, the error is always less than
+ *      1 ulp (unit in the last place).
+ *
+ * Constants:
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
+ * to produce the hexadecimal values shown.
+ */
+
+static const double Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */
+    Lg2 = 3.999999999940941908e-01,                 /* 3FD99999 9997FA04 */
+    Lg3 = 2.857142874366239149e-01,                 /* 3FD24924 94229359 */
+    Lg4 = 2.222219843214978396e-01,                 /* 3FCC71C5 1D8E78AF */
+    Lg5 = 1.818357216161805012e-01,                 /* 3FC74664 96CB03DE */
+    Lg6 = 1.531383769920937332e-01,                 /* 3FC39A09 D078C69F */
+    Lg7 = 1.479819860511658591e-01;                 /* 3FC2F112 DF3E5244 */
+
+/*
+ * We always inline k_log1p(), since doing so produces a
+ * substantial performance improvement (~40% on amd64).
+ */
+static inline double k_log1p(double f) {
+  double hfsq, s, z, R, w, t1, t2;
+
+  s = f / (2.0 + f);
+  z = s * s;
+  w = z * z;
+  t1 = w * (Lg2 + w * (Lg4 + w * Lg6));
+  t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7)));
+  R = t2 + t1;
+  hfsq = 0.5 * f * f;
+  return s * (hfsq + R);
+}
+
+/*
+ * Return the base 2 logarithm of x.  See e_log.c and k_log.h for most
+ * comments.
+ *
+ * This reduces x to {k, 1+f} exactly as in e_log.c, then calls the kernel,
+ * then does the combining and scaling steps
+ *    log2(x) = (f - 0.5*f*f + k_log1p(f)) / ln2 + k
+ * in not-quite-routine extra precision.
+ */
+double log2(double x) {
+  static const double
+      two54 = 1.80143985094819840000e+16,   /* 0x43500000, 0x00000000 */
+      ivln2hi = 1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */
+      ivln2lo = 1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */
+
+  static const double zero = 0.0;
+  static volatile double vzero = 0.0;
+
+  double f, hfsq, hi, lo, r, val_hi, val_lo, w, y;
+  int32_t i, k, hx;
+  uint32_t lx;
+
+  EXTRACT_WORDS(hx, lx, x);
+
+  k = 0;
+  if (hx < 0x00100000) { /* x < 2**-1022  */
+    if (((hx & 0x7fffffff) | lx) == 0)
+      return -two54 / vzero;           /* log(+-0)=-inf */
+    if (hx < 0) return (x - x) / zero; /* log(-#) = NaN */
+    k -= 54;
+    x *= two54; /* subnormal number, scale up x */
+    GET_HIGH_WORD(hx, x);
+  }
+  if (hx >= 0x7ff00000) return x + x;
+  if (hx == 0x3ff00000 && lx == 0) return zero; /* log(1) = +0 */
+  k += (hx >> 20) - 1023;
+  hx &= 0x000fffff;
+  i = (hx + 0x95f64) & 0x100000;
+  SET_HIGH_WORD(x, hx | (i ^ 0x3ff00000)); /* normalize x or x/2 */
+  k += (i >> 20);
+  y = static_cast<double>(k);
+  f = x - 1.0;
+  hfsq = 0.5 * f * f;
+  r = k_log1p(f);
+
+  /*
+   * f-hfsq must (for args near 1) be evaluated in extra precision
+   * to avoid a large cancellation when x is near sqrt(2) or 1/sqrt(2).
+   * This is fairly efficient since f-hfsq only depends on f, so can
+   * be evaluated in parallel with R.  Not combining hfsq with R also
+   * keeps R small (though not as small as a true `lo' term would be),
+   * so that extra precision is not needed for terms involving R.
+   *
+   * Compiler bugs involving extra precision used to break Dekker's
+   * theorem for spitting f-hfsq as hi+lo, unless double_t was used
+   * or the multi-precision calculations were avoided when double_t
+   * has extra precision.  These problems are now automatically
+   * avoided as a side effect of the optimization of combining the
+   * Dekker splitting step with the clear-low-bits step.
+   *
+   * y must (for args near sqrt(2) and 1/sqrt(2)) be added in extra
+   * precision to avoid a very large cancellation when x is very near
+   * these values.  Unlike the above cancellations, this problem is
+   * specific to base 2.  It is strange that adding +-1 is so much
+   * harder than adding +-ln2 or +-log10_2.
+   *
+   * This uses Dekker's theorem to normalize y+val_hi, so the
+   * compiler bugs are back in some configurations, sigh.  And I
+   * don't want to used double_t to avoid them, since that gives a
+   * pessimization and the support for avoiding the pessimization
+   * is not yet available.
+   *
+   * The multi-precision calculations for the multiplications are
+   * routine.
+   */
+  hi = f - hfsq;
+  SET_LOW_WORD(hi, 0);
+  lo = (f - hi) - hfsq + r;
+  val_hi = hi * ivln2hi;
+  val_lo = (lo + hi) * ivln2lo + lo * ivln2hi;
+
+  /* spadd(val_hi, val_lo, y), except for not using double_t: */
+  w = y + val_hi;
+  val_lo += (y - w) + val_hi;
+  val_hi = w;
+
+  return val_lo + val_hi;
+}
+
+/*
+ * Return the base 10 logarithm of x
+ *
+ * Method :
+ *      Let log10_2hi = leading 40 bits of log10(2) and
+ *          log10_2lo = log10(2) - log10_2hi,
+ *          ivln10   = 1/log(10) rounded.
+ *      Then
+ *              n = ilogb(x),
+ *              if(n<0)  n = n+1;
+ *              x = scalbn(x,-n);
+ *              log10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x))
+ *
+ *  Note 1:
+ *     To guarantee log10(10**n)=n, where 10**n is normal, the rounding
+ *     mode must set to Round-to-Nearest.
+ *  Note 2:
+ *      [1/log(10)] rounded to 53 bits has error .198 ulps;
+ *      log10 is monotonic at all binary break points.
+ *
+ *  Special cases:
+ *      log10(x) is NaN if x < 0;
+ *      log10(+INF) is +INF; log10(0) is -INF;
+ *      log10(NaN) is that NaN;
+ *      log10(10**N) = N  for N=0,1,...,22.
+ */
+double log10(double x) {
+  static const double
+      two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
+      ivln10 = 4.34294481903251816668e-01,
+      log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */
+      log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
+
+  static const double zero = 0.0;
+  static volatile double vzero = 0.0;
+
+  double y;
+  int32_t i, k, hx;
+  uint32_t lx;
+
+  EXTRACT_WORDS(hx, lx, x);
+
+  k = 0;
+  if (hx < 0x00100000) { /* x < 2**-1022  */
+    if (((hx & 0x7fffffff) | lx) == 0)
+      return -two54 / vzero;           /* log(+-0)=-inf */
+    if (hx < 0) return (x - x) / zero; /* log(-#) = NaN */
+    k -= 54;
+    x *= two54; /* subnormal number, scale up x */
+    GET_HIGH_WORD(hx, x);
+    GET_LOW_WORD(lx, x);
+  }
+  if (hx >= 0x7ff00000) return x + x;
+  if (hx == 0x3ff00000 && lx == 0) return zero; /* log(1) = +0 */
+  k += (hx >> 20) - 1023;
+
+  i = (k & 0x80000000) >> 31;
+  hx = (hx & 0x000fffff) | ((0x3ff - i) << 20);
+  y = k + i;
+  SET_HIGH_WORD(x, hx);
+  SET_LOW_WORD(x, lx);
+
+  double z = y * log10_2lo + ivln10 * log(x);
+  return z + y * log10_2hi;
+}
+
+/* expm1(x)
+ * Returns exp(x)-1, the exponential of x minus 1.
+ *
+ * Method
+ *   1. Argument reduction:
+ *  Given x, find r and integer k such that
+ *
+ *               x = k*ln2 + r,  |r| <= 0.5*ln2 ~ 0.34658
+ *
+ *      Here a correction term c will be computed to compensate
+ *  the error in r when rounded to a floating-point number.
+ *
+ *   2. Approximating expm1(r) by a special rational function on
+ *  the interval [0,0.34658]:
+ *  Since
+ *      r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 - r^4/360 + ...
+ *  we define R1(r*r) by
+ *      r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 * R1(r*r)
+ *  That is,
+ *      R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r)
+ *         = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r))
+ *         = 1 - r^2/60 + r^4/2520 - r^6/100800 + ...
+ *      We use a special Reme algorithm on [0,0.347] to generate
+ *   a polynomial of degree 5 in r*r to approximate R1. The
+ *  maximum error of this polynomial approximation is bounded
+ *  by 2**-61. In other words,
+ *      R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5
+ *  where   Q1  =  -1.6666666666666567384E-2,
+ *     Q2  =   3.9682539681370365873E-4,
+ *     Q3  =  -9.9206344733435987357E-6,
+ *     Q4  =   2.5051361420808517002E-7,
+ *     Q5  =  -6.2843505682382617102E-9;
+ *    z   =  r*r,
+ *  with error bounded by
+ *      |                  5           |     -61
+ *      | 1.0+Q1*z+...+Q5*z   -  R1(z) | <= 2
+ *      |                              |
+ *
+ *  expm1(r) = exp(r)-1 is then computed by the following
+ *   specific way which minimize the accumulation rounding error:
+ *             2     3
+ *            r     r    [ 3 - (R1 + R1*r/2)  ]
+ *        expm1(r) = r + --- + --- * [--------------------]
+ *                  2     2    [ 6 - r*(3 - R1*r/2) ]
+ *
+ *  To compensate the error in the argument reduction, we use
+ *    expm1(r+c) = expm1(r) + c + expm1(r)*c
+ *         ~ expm1(r) + c + r*c
+ *  Thus c+r*c will be added in as the correction terms for
+ *  expm1(r+c). Now rearrange the term to avoid optimization
+ *   screw up:
+ *            (      2                                    2 )
+ *            ({  ( r    [ R1 -  (3 - R1*r/2) ]  )  }    r  )
+ *   expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- )
+ *                  ({  ( 2    [ 6 - r*(3 - R1*r/2) ]  )  }    2  )
+ *                      (                                             )
+ *
+ *       = r - E
+ *   3. Scale back to obtain expm1(x):
+ *  From step 1, we have
+ *     expm1(x) = either 2^k*[expm1(r)+1] - 1
+ *        = or     2^k*[expm1(r) + (1-2^-k)]
+ *   4. Implementation notes:
+ *  (A). To save one multiplication, we scale the coefficient Qi
+ *       to Qi*2^i, and replace z by (x^2)/2.
+ *  (B). To achieve maximum accuracy, we compute expm1(x) by
+ *    (i)   if x < -56*ln2, return -1.0, (raise inexact if x!=inf)
+ *    (ii)  if k=0, return r-E
+ *    (iii) if k=-1, return 0.5*(r-E)-0.5
+ *        (iv)  if k=1 if r < -0.25, return 2*((r+0.5)- E)
+ *                  else       return  1.0+2.0*(r-E);
+ *    (v)   if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1)
+ *    (vi)  if k <= 20, return 2^k((1-2^-k)-(E-r)), else
+ *    (vii) return 2^k(1-((E+2^-k)-r))
+ *
+ * Special cases:
+ *  expm1(INF) is INF, expm1(NaN) is NaN;
+ *  expm1(-INF) is -1, and
+ *  for finite argument, only expm1(0)=0 is exact.
+ *
+ * Accuracy:
+ *  according to an error analysis, the error is always less than
+ *  1 ulp (unit in the last place).
+ *
+ * Misc. info.
+ *  For IEEE double
+ *      if x >  7.09782712893383973096e+02 then expm1(x) overflow
+ *
+ * Constants:
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
+ * to produce the hexadecimal values shown.
+ */
+double expm1(double x) {
+  static const double
+      one = 1.0,
+      tiny = 1.0e-300,
+      o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
+      ln2_hi = 6.93147180369123816490e-01,      /* 0x3fe62e42, 0xfee00000 */
+      ln2_lo = 1.90821492927058770002e-10,      /* 0x3dea39ef, 0x35793c76 */
+      invln2 = 1.44269504088896338700e+00,      /* 0x3ff71547, 0x652b82fe */
+      /* Scaled Q's: Qn_here = 2**n * Qn_above, for R(2*z) where z = hxs =
+         x*x/2: */
+      Q1 = -3.33333333333331316428e-02, /* BFA11111 111110F4 */
+      Q2 = 1.58730158725481460165e-03,  /* 3F5A01A0 19FE5585 */
+      Q3 = -7.93650757867487942473e-05, /* BF14CE19 9EAADBB7 */
+      Q4 = 4.00821782732936239552e-06,  /* 3ED0CFCA 86E65239 */
+      Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
+
+  static volatile double huge = 1.0e+300;
+
+  double y, hi, lo, c, t, e, hxs, hfx, r1, twopk;
+  int32_t k, xsb;
+  uint32_t hx;
+
+  GET_HIGH_WORD(hx, x);
+  xsb = hx & 0x80000000; /* sign bit of x */
+  hx &= 0x7fffffff;      /* high word of |x| */
+
+  /* filter out huge and non-finite argument */
+  if (hx >= 0x4043687A) {   /* if |x|>=56*ln2 */
+    if (hx >= 0x40862E42) { /* if |x|>=709.78... */
+      if (hx >= 0x7ff00000) {
+        uint32_t low;
+        GET_LOW_WORD(low, x);
+        if (((hx & 0xfffff) | low) != 0)
+          return x + x; /* NaN */
+        else
+          return (xsb == 0) ? x : -1.0; /* exp(+-inf)={inf,-1} */
+      }
+      if (x > o_threshold) return huge * huge; /* overflow */
+    }
+    if (xsb != 0) {        /* x < -56*ln2, return -1.0 with inexact */
+      if (x + tiny < 0.0)  /* raise inexact */
+        return tiny - one; /* return -1 */
+    }
+  }
+
+  /* argument reduction */
+  if (hx > 0x3fd62e42) {   /* if  |x| > 0.5 ln2 */
+    if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
+      if (xsb == 0) {
+        hi = x - ln2_hi;
+        lo = ln2_lo;
+        k = 1;
+      } else {
+        hi = x + ln2_hi;
+        lo = -ln2_lo;
+        k = -1;
+      }
+    } else {
+      k = invln2 * x + ((xsb == 0) ? 0.5 : -0.5);
+      t = k;
+      hi = x - t * ln2_hi; /* t*ln2_hi is exact here */
+      lo = t * ln2_lo;
+    }
+    STRICT_ASSIGN(double, x, hi - lo);
+    c = (hi - x) - lo;
+  } else if (hx < 0x3c900000) { /* when |x|<2**-54, return x */
+    t = huge + x;               /* return x with inexact flags when x!=0 */
+    return x - (t - (huge + x));
+  } else {
+    k = 0;
+  }
+
+  /* x is now in primary range */
+  hfx = 0.5 * x;
+  hxs = x * hfx;
+  r1 = one + hxs * (Q1 + hxs * (Q2 + hxs * (Q3 + hxs * (Q4 + hxs * Q5))));
+  t = 3.0 - r1 * hfx;
+  e = hxs * ((r1 - t) / (6.0 - x * t));
+  if (k == 0) {
+    return x - (x * e - hxs); /* c is 0 */
+  } else {
+    INSERT_WORDS(twopk, 0x3ff00000 + (k << 20), 0); /* 2^k */
+    e = (x * (e - c) - c);
+    e -= hxs;
+    if (k == -1) return 0.5 * (x - e) - 0.5;
+    if (k == 1) {
+      if (x < -0.25)
+        return -2.0 * (e - (x + 0.5));
+      else
+        return one + 2.0 * (x - e);
+    }
+    if (k <= -2 || k > 56) { /* suffice to return exp(x)-1 */
+      y = one - (e - x);
+      // TODO(mvstanton): is this replacement for the hex float
+      // sufficient?
+      // if (k == 1024) y = y*2.0*0x1p1023;
+      if (k == 1024)
+        y = y * 2.0 * 8.98846567431158e+307;
+      else
+        y = y * twopk;
+      return y - one;
+    }
+    t = one;
+    if (k < 20) {
+      SET_HIGH_WORD(t, 0x3ff00000 - (0x200000 >> k)); /* t=1-2^-k */
+      y = t - (e - x);
+      y = y * twopk;
+    } else {
+      SET_HIGH_WORD(t, ((0x3ff - k) << 20)); /* 2^-k */
+      y = x - (e + t);
+      y += one;
+      y = y * twopk;
+    }
+  }
+  return y;
+}
+
+double cbrt(double x) {
+  static const uint32_t
+      B1 = 715094163, /* B1 = (1023-1023/3-0.03306235651)*2**20 */
+      B2 = 696219795; /* B2 = (1023-1023/3-54/3-0.03306235651)*2**20 */
+
+  /* |1/cbrt(x) - p(x)| < 2**-23.5 (~[-7.93e-8, 7.929e-8]). */
+  static const double P0 = 1.87595182427177009643, /* 0x3ffe03e6, 0x0f61e692 */
+      P1 = -1.88497979543377169875,                /* 0xbffe28e0, 0x92f02420 */
+      P2 = 1.621429720105354466140,                /* 0x3ff9f160, 0x4a49d6c2 */
+      P3 = -0.758397934778766047437,               /* 0xbfe844cb, 0xbee751d9 */
+      P4 = 0.145996192886612446982;                /* 0x3fc2b000, 0xd4e4edd7 */
+
+  int32_t hx;
+  union {
+    double value;
+    uint64_t bits;
+  } u;
+  double r, s, t = 0.0, w;
+  uint32_t sign;
+  uint32_t high, low;
+
+  EXTRACT_WORDS(hx, low, x);
+  sign = hx & 0x80000000; /* sign= sign(x) */
+  hx ^= sign;
+  if (hx >= 0x7ff00000) return (x + x); /* cbrt(NaN,INF) is itself */
+
+  /*
+   * Rough cbrt to 5 bits:
+   *    cbrt(2**e*(1+m) ~= 2**(e/3)*(1+(e%3+m)/3)
+   * where e is integral and >= 0, m is real and in [0, 1), and "/" and
+   * "%" are integer division and modulus with rounding towards minus
+   * infinity.  The RHS is always >= the LHS and has a maximum relative
+   * error of about 1 in 16.  Adding a bias of -0.03306235651 to the
+   * (e%3+m)/3 term reduces the error to about 1 in 32. With the IEEE
+   * floating point representation, for finite positive normal values,
+   * ordinary integer divison of the value in bits magically gives
+   * almost exactly the RHS of the above provided we first subtract the
+   * exponent bias (1023 for doubles) and later add it back.  We do the
+   * subtraction virtually to keep e >= 0 so that ordinary integer
+   * division rounds towards minus infinity; this is also efficient.
+   */
+  if (hx < 0x00100000) {             /* zero or subnormal? */
+    if ((hx | low) == 0) return (x); /* cbrt(0) is itself */
+    SET_HIGH_WORD(t, 0x43500000);    /* set t= 2**54 */
+    t *= x;
+    GET_HIGH_WORD(high, t);
+    INSERT_WORDS(t, sign | ((high & 0x7fffffff) / 3 + B2), 0);
+  } else {
+    INSERT_WORDS(t, sign | (hx / 3 + B1), 0);
+  }
+
+  /*
+   * New cbrt to 23 bits:
+   *    cbrt(x) = t*cbrt(x/t**3) ~= t*P(t**3/x)
+   * where P(r) is a polynomial of degree 4 that approximates 1/cbrt(r)
+   * to within 2**-23.5 when |r - 1| < 1/10.  The rough approximation
+   * has produced t such than |t/cbrt(x) - 1| ~< 1/32, and cubing this
+   * gives us bounds for r = t**3/x.
+   *
+   * Try to optimize for parallel evaluation as in k_tanf.c.
+   */
+  r = (t * t) * (t / x);
+  t = t * ((P0 + r * (P1 + r * P2)) + ((r * r) * r) * (P3 + r * P4));
+
+  /*
+   * Round t away from zero to 23 bits (sloppily except for ensuring that
+   * the result is larger in magnitude than cbrt(x) but not much more than
+   * 2 23-bit ulps larger).  With rounding towards zero, the error bound
+   * would be ~5/6 instead of ~4/6.  With a maximum error of 2 23-bit ulps
+   * in the rounded t, the infinite-precision error in the Newton
+   * approximation barely affects third digit in the final error
+   * 0.667; the error in the rounded t can be up to about 3 23-bit ulps
+   * before the final error is larger than 0.667 ulps.
+   */
+  u.value = t;
+  u.bits = (u.bits + 0x80000000) & 0xffffffffc0000000ULL;
+  t = u.value;
+
+  /* one step Newton iteration to 53 bits with error < 0.667 ulps */
+  s = t * t;             /* t*t is exact */
+  r = x / s;             /* error <= 0.5 ulps; |r| < |t| */
+  w = t + t;             /* t+t is exact */
+  r = (r - t) / (w + r); /* r-t is exact; w+r ~= 3*t */
+  t = t + t * r;         /* error <= 0.5 + 0.5/3 + epsilon */
+
+  return (t);
+}
+
+/* sin(x)
+ * Return sine function of x.
+ *
+ * kernel function:
+ *      __kernel_sin            ... sine function on [-pi/4,pi/4]
+ *      __kernel_cos            ... cose function on [-pi/4,pi/4]
+ *      __ieee754_rem_pio2      ... argument reduction routine
+ *
+ * Method.
+ *      Let S,C and T denote the sin, cos and tan respectively on
+ *      [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
+ *      in [-pi/4 , +pi/4], and let n = k mod 4.
+ *      We have
+ *
+ *          n        sin(x)      cos(x)        tan(x)
+ *     ----------------------------------------------------------
+ *          0          S           C             T
+ *          1          C          -S            -1/T
+ *          2         -S          -C             T
+ *          3         -C           S            -1/T
+ *     ----------------------------------------------------------
+ *
+ * Special cases:
+ *      Let trig be any of sin, cos, or tan.
+ *      trig(+-INF)  is NaN, with signals;
+ *      trig(NaN)    is that NaN;
+ *
+ * Accuracy:
+ *      TRIG(x) returns trig(x) nearly rounded
+ */
+double sin(double x) {
+  double y[2], z = 0.0;
+  int32_t n, ix;
+
+  /* High word of x. */
+  GET_HIGH_WORD(ix, x);
+
+  /* |x| ~< pi/4 */
+  ix &= 0x7fffffff;
+  if (ix <= 0x3fe921fb) {
+    return __kernel_sin(x, z, 0);
+  } else if (ix >= 0x7ff00000) {
+    /* sin(Inf or NaN) is NaN */
+    return x - x;
+  } else {
+    /* argument reduction needed */
+    n = __ieee754_rem_pio2(x, y);
+    switch (n & 3) {
+      case 0:
+        return __kernel_sin(y[0], y[1], 1);
+      case 1:
+        return __kernel_cos(y[0], y[1]);
+      case 2:
+        return -__kernel_sin(y[0], y[1], 1);
+      default:
+        return -__kernel_cos(y[0], y[1]);
+    }
+  }
+}
+
+/* tan(x)
+ * Return tangent function of x.
+ *
+ * kernel function:
+ *      __kernel_tan            ... tangent function on [-pi/4,pi/4]
+ *      __ieee754_rem_pio2      ... argument reduction routine
+ *
+ * Method.
+ *      Let S,C and T denote the sin, cos and tan respectively on
+ *      [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
+ *      in [-pi/4 , +pi/4], and let n = k mod 4.
+ *      We have
+ *
+ *          n        sin(x)      cos(x)        tan(x)
+ *     ----------------------------------------------------------
+ *          0          S           C             T
+ *          1          C          -S            -1/T
+ *          2         -S          -C             T
+ *          3         -C           S            -1/T
+ *     ----------------------------------------------------------
+ *
+ * Special cases:
+ *      Let trig be any of sin, cos, or tan.
+ *      trig(+-INF)  is NaN, with signals;
+ *      trig(NaN)    is that NaN;
+ *
+ * Accuracy:
+ *      TRIG(x) returns trig(x) nearly rounded
+ */
+double tan(double x) {
+  double y[2], z = 0.0;
+  int32_t n, ix;
+
+  /* High word of x. */
+  GET_HIGH_WORD(ix, x);
+
+  /* |x| ~< pi/4 */
+  ix &= 0x7fffffff;
+  if (ix <= 0x3fe921fb) {
+    return __kernel_tan(x, z, 1);
+  } else if (ix >= 0x7ff00000) {
+    /* tan(Inf or NaN) is NaN */
+    return x - x; /* NaN */
+  } else {
+    /* argument reduction needed */
+    n = __ieee754_rem_pio2(x, y);
+    /* 1 -> n even, -1 -> n odd */
+    return __kernel_tan(y[0], y[1], 1 - ((n & 1) << 1));
+  }
+}
+
+}  // namespace ieee754
+}  // namespace base
+}  // namespace v8
diff --git a/src/base/ieee754.h b/src/base/ieee754.h
new file mode 100644
index 0000000..cf33580
--- /dev/null
+++ b/src/base/ieee754.h
@@ -0,0 +1,57 @@
+// Copyright 2016 the V8 project authors. All rights reserved.
+// Use of this source code is governed by a BSD-style license that can be
+// found in the LICENSE file.
+
+#ifndef V8_BASE_IEEE754_H_
+#define V8_BASE_IEEE754_H_
+
+namespace v8 {
+namespace base {
+namespace ieee754 {
+
+// Returns the principal value of the arc tangent of |x|; that is the value
+// whose tangent is |x|.
+double atan(double x);
+
+// Returns the principal value of the arc tangent of |y/x|, using the signs of
+// the two arguments to determine the quadrant of the result.
+double atan2(double y, double x);
+
+// Returns the cosine of |x|, where |x| is given in radians.
+double cos(double x);
+
+// Returns the base-e exponential of |x|.
+double exp(double x);
+
+double atanh(double x);
+
+// Returns the natural logarithm of |x|.
+double log(double x);
+
+// Returns a value equivalent to |log(1+x)|, but computed in a way that is
+// accurate even if the value of |x| is near zero.
+double log1p(double x);
+
+// Returns the base 2 logarithm of |x|.
+double log2(double x);
+
+// Returns the base 10 logarithm of |x|.
+double log10(double x);
+
+// Returns the cube root of |x|.
+double cbrt(double x);
+
+// Returns exp(x)-1, the exponential of |x| minus 1.
+double expm1(double x);
+
+// Returns the sine of |x|, where |x| is given in radians.
+double sin(double x);
+
+// Returns the tangent of |x|, where |x| is given in radians.
+double tan(double x);
+
+}  // namespace ieee754
+}  // namespace base
+}  // namespace v8
+
+#endif  // V8_BASE_IEEE754_H_
diff --git a/src/base/platform/semaphore.cc b/src/base/platform/semaphore.cc
index a7e522a..7bf5986 100644
--- a/src/base/platform/semaphore.cc
+++ b/src/base/platform/semaphore.cc
@@ -87,10 +87,6 @@
       0, reinterpret_cast<uintptr_t>(&native_handle_) &
       kSemaphoreAlignmentMask);
   DCHECK(count >= 0);
-#if V8_LIBC_GLIBC
-  // sem_init in glibc prior to 2.1 does not zero out semaphores.
-  memset(&native_handle_, 0, sizeof(native_handle_));
-#endif
   int result = sem_init(&native_handle_, 0, count);
   DCHECK_EQ(0, result);
   USE(result);
@@ -105,6 +101,9 @@
 
 void Semaphore::Signal() {
   int result = sem_post(&native_handle_);
+  // This check may fail with <libc-2.21, which we use on the try bots, if the
+  // semaphore is destroyed while sem_post is still executed. A work around is
+  // to extend the lifetime of the semaphore.
   CHECK_EQ(0, result);
 }
 
diff --git a/src/base/platform/time.cc b/src/base/platform/time.cc
index b2355a3..786ef2e 100644
--- a/src/base/platform/time.cc
+++ b/src/base/platform/time.cc
@@ -41,9 +41,11 @@
   CHECK(kr == KERN_SUCCESS);
 
   v8::base::CheckedNumeric<int64_t> absolute_micros(
-      thread_info_data.user_time.seconds);
+      thread_info_data.user_time.seconds +
+      thread_info_data.system_time.seconds);
   absolute_micros *= v8::base::Time::kMicrosecondsPerSecond;
-  absolute_micros += thread_info_data.user_time.microseconds;
+  absolute_micros += (thread_info_data.user_time.microseconds +
+                      thread_info_data.system_time.microseconds);
   return absolute_micros.ValueOrDie();
 }
 #elif V8_OS_POSIX
@@ -51,9 +53,20 @@
 // microsecond timebase. Minimum requirement is MONOTONIC_CLOCK to be supported
 // on the system. FreeBSD 6 has CLOCK_MONOTONIC but defines
 // _POSIX_MONOTONIC_CLOCK to -1.
-inline int64_t ClockNow(clockid_t clk_id) {
+V8_INLINE int64_t ClockNow(clockid_t clk_id) {
 #if (defined(_POSIX_MONOTONIC_CLOCK) && _POSIX_MONOTONIC_CLOCK >= 0) || \
   defined(V8_OS_BSD) || defined(V8_OS_ANDROID)
+// On AIX clock_gettime for CLOCK_THREAD_CPUTIME_ID outputs time with
+// resolution of 10ms. thread_cputime API provides the time in ns
+#if defined(V8_OS_AIX)
+  thread_cputime_t tc;
+  if (clk_id == CLOCK_THREAD_CPUTIME_ID) {
+    if (thread_cputime(-1, &tc) != 0) {
+      UNREACHABLE();
+      return 0;
+    }
+  }
+#endif
   struct timespec ts;
   if (clock_gettime(clk_id, &ts) != 0) {
     UNREACHABLE();
@@ -61,12 +74,38 @@
   }
   v8::base::internal::CheckedNumeric<int64_t> result(ts.tv_sec);
   result *= v8::base::Time::kMicrosecondsPerSecond;
+#if defined(V8_OS_AIX)
+  if (clk_id == CLOCK_THREAD_CPUTIME_ID) {
+    result += (tc.stime / v8::base::Time::kNanosecondsPerMicrosecond);
+  } else {
+    result += (ts.tv_nsec / v8::base::Time::kNanosecondsPerMicrosecond);
+  }
+#else
   result += (ts.tv_nsec / v8::base::Time::kNanosecondsPerMicrosecond);
+#endif
   return result.ValueOrDie();
 #else  // Monotonic clock not supported.
   return 0;
 #endif
 }
+#elif V8_OS_WIN
+V8_INLINE bool IsQPCReliable() {
+  v8::base::CPU cpu;
+  // On Athlon X2 CPUs (e.g. model 15) QueryPerformanceCounter is unreliable.
+  return strcmp(cpu.vendor(), "AuthenticAMD") == 0 && cpu.family() == 15;
+}
+
+// Returns the current value of the performance counter.
+V8_INLINE uint64_t QPCNowRaw() {
+  LARGE_INTEGER perf_counter_now = {};
+  // According to the MSDN documentation for QueryPerformanceCounter(), this
+  // will never fail on systems that run XP or later.
+  // https://msdn.microsoft.com/library/windows/desktop/ms644904.aspx
+  BOOL result = ::QueryPerformanceCounter(&perf_counter_now);
+  DCHECK(result);
+  USE(result);
+  return perf_counter_now.QuadPart;
+}
 #endif  // V8_OS_MACOSX
 
 
@@ -456,15 +495,12 @@
   virtual ~HighResolutionTickClock() {}
 
   int64_t Now() override {
-    LARGE_INTEGER now;
-    BOOL result = QueryPerformanceCounter(&now);
-    DCHECK(result);
-    USE(result);
+    uint64_t now = QPCNowRaw();
 
     // Intentionally calculate microseconds in a round about manner to avoid
     // overflow and precision issues. Think twice before simplifying!
-    int64_t whole_seconds = now.QuadPart / ticks_per_second_;
-    int64_t leftover_ticks = now.QuadPart % ticks_per_second_;
+    int64_t whole_seconds = now / ticks_per_second_;
+    int64_t leftover_ticks = now % ticks_per_second_;
     int64_t ticks = (whole_seconds * Time::kMicrosecondsPerSecond) +
         ((leftover_ticks * Time::kMicrosecondsPerSecond) / ticks_per_second_);
 
@@ -529,10 +565,8 @@
       return tick_clock.Pointer();
     }
 
-    // On Athlon X2 CPUs (e.g. model 15) the QueryPerformanceCounter
-    // is unreliable, fallback to the low-resolution tick clock.
-    CPU cpu;
-    if (strcmp(cpu.vendor(), "AuthenticAMD") == 0 && cpu.family() == 15) {
+    // If QPC not reliable, fallback to low-resolution tick clock.
+    if (IsQPCReliable()) {
       return tick_clock.Pointer();
     }
 
@@ -621,11 +655,106 @@
 #elif(defined(_POSIX_THREAD_CPUTIME) && (_POSIX_THREAD_CPUTIME >= 0)) || \
   defined(V8_OS_ANDROID)
   return ThreadTicks(ClockNow(CLOCK_THREAD_CPUTIME_ID));
+#elif V8_OS_WIN
+  return ThreadTicks::GetForThread(::GetCurrentThread());
 #else
   UNREACHABLE();
   return ThreadTicks();
 #endif
 }
 
+
+#if V8_OS_WIN
+ThreadTicks ThreadTicks::GetForThread(const HANDLE& thread_handle) {
+  DCHECK(IsSupported());
+
+  // Get the number of TSC ticks used by the current thread.
+  ULONG64 thread_cycle_time = 0;
+  ::QueryThreadCycleTime(thread_handle, &thread_cycle_time);
+
+  // Get the frequency of the TSC.
+  double tsc_ticks_per_second = TSCTicksPerSecond();
+  if (tsc_ticks_per_second == 0)
+    return ThreadTicks();
+
+  // Return the CPU time of the current thread.
+  double thread_time_seconds = thread_cycle_time / tsc_ticks_per_second;
+  return ThreadTicks(
+      static_cast<int64_t>(thread_time_seconds * Time::kMicrosecondsPerSecond));
+}
+
+// static
+bool ThreadTicks::IsSupportedWin() {
+  static bool is_supported = base::CPU().has_non_stop_time_stamp_counter() &&
+                             !IsQPCReliable();
+  return is_supported;
+}
+
+// static
+void ThreadTicks::WaitUntilInitializedWin() {
+  while (TSCTicksPerSecond() == 0)
+    ::Sleep(10);
+}
+
+double ThreadTicks::TSCTicksPerSecond() {
+  DCHECK(IsSupported());
+
+  // The value returned by QueryPerformanceFrequency() cannot be used as the TSC
+  // frequency, because there is no guarantee that the TSC frequency is equal to
+  // the performance counter frequency.
+
+  // The TSC frequency is cached in a static variable because it takes some time
+  // to compute it.
+  static double tsc_ticks_per_second = 0;
+  if (tsc_ticks_per_second != 0)
+    return tsc_ticks_per_second;
+
+  // Increase the thread priority to reduces the chances of having a context
+  // switch during a reading of the TSC and the performance counter.
+  int previous_priority = ::GetThreadPriority(::GetCurrentThread());
+  ::SetThreadPriority(::GetCurrentThread(), THREAD_PRIORITY_HIGHEST);
+
+  // The first time that this function is called, make an initial reading of the
+  // TSC and the performance counter.
+  static const uint64_t tsc_initial = __rdtsc();
+  static const uint64_t perf_counter_initial = QPCNowRaw();
+
+  // Make a another reading of the TSC and the performance counter every time
+  // that this function is called.
+  uint64_t tsc_now = __rdtsc();
+  uint64_t perf_counter_now = QPCNowRaw();
+
+  // Reset the thread priority.
+  ::SetThreadPriority(::GetCurrentThread(), previous_priority);
+
+  // Make sure that at least 50 ms elapsed between the 2 readings. The first
+  // time that this function is called, we don't expect this to be the case.
+  // Note: The longer the elapsed time between the 2 readings is, the more
+  //   accurate the computed TSC frequency will be. The 50 ms value was
+  //   chosen because local benchmarks show that it allows us to get a
+  //   stddev of less than 1 tick/us between multiple runs.
+  // Note: According to the MSDN documentation for QueryPerformanceFrequency(),
+  //   this will never fail on systems that run XP or later.
+  //   https://msdn.microsoft.com/library/windows/desktop/ms644905.aspx
+  LARGE_INTEGER perf_counter_frequency = {};
+  ::QueryPerformanceFrequency(&perf_counter_frequency);
+  DCHECK_GE(perf_counter_now, perf_counter_initial);
+  uint64_t perf_counter_ticks = perf_counter_now - perf_counter_initial;
+  double elapsed_time_seconds =
+      perf_counter_ticks / static_cast<double>(perf_counter_frequency.QuadPart);
+
+  const double kMinimumEvaluationPeriodSeconds = 0.05;
+  if (elapsed_time_seconds < kMinimumEvaluationPeriodSeconds)
+    return 0;
+
+  // Compute the frequency of the TSC.
+  DCHECK_GE(tsc_now, tsc_initial);
+  uint64_t tsc_ticks = tsc_now - tsc_initial;
+  tsc_ticks_per_second = tsc_ticks / elapsed_time_seconds;
+
+  return tsc_ticks_per_second;
+}
+#endif  // V8_OS_WIN
+
 }  // namespace base
 }  // namespace v8
diff --git a/src/base/platform/time.h b/src/base/platform/time.h
index e17fc1d..be62014 100644
--- a/src/base/platform/time.h
+++ b/src/base/platform/time.h
@@ -12,6 +12,9 @@
 #include "src/base/bits.h"
 #include "src/base/macros.h"
 #include "src/base/safe_math.h"
+#if V8_OS_WIN
+#include "src/base/win32-headers.h"
+#endif
 
 // Forward declarations.
 extern "C" {
@@ -380,17 +383,42 @@
   // Returns true if ThreadTicks::Now() is supported on this system.
   static bool IsSupported();
 
+  // Waits until the initialization is completed. Needs to be guarded with a
+  // call to IsSupported().
+  static void WaitUntilInitialized() {
+#if V8_OS_WIN
+    WaitUntilInitializedWin();
+#endif
+  }
+
   // Returns thread-specific CPU-time on systems that support this feature.
   // Needs to be guarded with a call to IsSupported(). Use this timer
   // to (approximately) measure how much time the calling thread spent doing
   // actual work vs. being de-scheduled. May return bogus results if the thread
   // migrates to another CPU between two calls. Returns an empty ThreadTicks
-  // object until the initialization is completed.
+  // object until the initialization is completed. If a clock reading is
+  // absolutely needed, call WaitUntilInitialized() before this method.
   static ThreadTicks Now();
 
+#if V8_OS_WIN
+  // Similar to Now() above except this returns thread-specific CPU time for an
+  // arbitrary thread. All comments for Now() method above apply apply to this
+  // method as well.
+  static ThreadTicks GetForThread(const HANDLE& thread_handle);
+#endif
+
  private:
-  // This is for internal use and testing. Ticks are in microseconds.
+  // Please use Now() or GetForThread() to create a new object. This is for
+  // internal use and testing. Ticks are in microseconds.
   explicit ThreadTicks(int64_t ticks) : TimeBase(ticks) {}
+
+#if V8_OS_WIN
+  // Returns the frequency of the TSC in ticks per second, or 0 if it hasn't
+  // been measured yet. Needs to be guarded with a call to IsSupported().
+  static double TSCTicksPerSecond();
+  static bool IsSupportedWin();
+  static void WaitUntilInitializedWin();
+#endif
 };
 
 }  // namespace base
diff --git a/src/base/utils/random-number-generator.cc b/src/base/utils/random-number-generator.cc
index ff42840..3a6f2c6 100644
--- a/src/base/utils/random-number-generator.cc
+++ b/src/base/utils/random-number-generator.cc
@@ -124,10 +124,10 @@
 
 
 void RandomNumberGenerator::SetSeed(int64_t seed) {
-  if (seed == 0) seed = 1;
   initial_seed_ = seed;
   state0_ = MurmurHash3(bit_cast<uint64_t>(seed));
-  state1_ = MurmurHash3(state0_);
+  state1_ = MurmurHash3(~state0_);
+  CHECK(state0_ != 0 || state1_ != 0);
 }
 
 
diff --git a/src/base/win32-headers.h b/src/base/win32-headers.h
index 20ec8e0..b61ce71 100644
--- a/src/base/win32-headers.h
+++ b/src/base/win32-headers.h
@@ -27,10 +27,10 @@
 #ifndef NOMCX
 #define NOMCX
 #endif
-// Require Windows XP or higher (this is required for the RtlCaptureContext
-// function to be present).
+// Require Windows Vista or higher (this is required for the
+// QueryThreadCycleTime function to be present).
 #ifndef _WIN32_WINNT
-#define _WIN32_WINNT 0x501
+#define _WIN32_WINNT 0x0600
 #endif
 
 #include <windows.h>