Merge V8 5.3.332.45. DO NOT MERGE
Test: Manual
FPIIM-449
Change-Id: Id3254828b068abdea3cb10442e0172a8c9a98e03
(cherry picked from commit 13e2dadd00298019ed862f2b2fc5068bba730bcf)
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>