- migrated to GZDoom's utility code.

This commit is contained in:
Christoph Oelckers 2020-04-11 23:50:43 +02:00
parent 2a9813eb5c
commit f671eb622f
168 changed files with 10543 additions and 737 deletions

View file

@ -603,6 +603,11 @@ file( GLOB HEADER_FILES
core/rendering/hwrenderer/postprocessing/*.h
core/rendering/hwrenderer/utility/*.h
common/utility/*.h
common/thirdparty/*.h
common/thirdparty/rapidjson/*.h
common/thirdparty/math./*h
build/src/*.h
platform/win32/*.h
platform/posix/*.h
@ -634,8 +639,9 @@ set( FASTMATH_SOURCES
libsmackerdec/src/LogError.cpp
libsmackerdec/src/SmackerDecoder.cpp
common/utility/matrix.cpp
# The rest is only here because it is C, not C++
core/utility/strnatcmp.c
core/rendering/gl_load/gl_load.c
gitinfo.cpp
@ -714,26 +720,29 @@ set (PCH_SOURCES
core/console/d_event.cpp
core/console/c_con.cpp
core/utility/i_module.cpp
core/utility/i_time.cpp
common/thirdparty/sfmt/SFMT.cpp
common/utility/engineerrors.cpp
common/utility/i_module.cpp
common/utility/m_alloc.cpp
common/utility/utf8.cpp
common/utility/palette.cpp
common/utility/files.cpp
common/utility/files_decompress.cpp
common/utility/memarena.cpp
common/utility/cmdlib.cpp
common/utility/configfile.cpp
common/utility/i_time.cpp
common/utility/m_argv.cpp
common/utility/s_playlist.cpp
common/utility/zstrformat.cpp
common/thirdparty/md5.cpp
common/thirdparty/superfasthash.cpp
core/utility/name.cpp
core/utility/m_alloc.cpp
core/utility/cmdlib.cpp
core/utility/m_argv.cpp
core/utility/files.cpp
core/utility/files_decompress.cpp
core/utility/zstring.cpp
core/utility/zstrformat.cpp
core/utility/utf8.cpp
core/utility/superfasthash.cpp
core/utility/configfile.cpp
core/utility/matrix.cpp
core/utility/m_png.cpp
core/utility/memarena.cpp
core/utility/sc_man.cpp
core/utility/stringtable.cpp
core/utility/stats.cpp
core/utility/palette.cpp
core/filesystem/filesystem.cpp
core/filesystem/ancientzip.cpp
@ -826,9 +835,28 @@ add_executable( ${PROJECT_NAME} WIN32 MACOSX_BUNDLE
${SYSTEM_SOURCES}
${FASTMATH_SOURCES}
${PCH_SOURCES}
#utility/strnatcmp.c
#utility/zstring.cpp
#zzautozend.cpp
common/utility/x86.cpp
common/thirdparty/strnatcmp.c
common/utility/zstring.cpp
common/utility/findfile.cpp
common/thirdparty/math/asin.c
common/thirdparty/math/atan.c
common/thirdparty/math/const.c
common/thirdparty/math/cosh.c
common/thirdparty/math/exp.c
common/thirdparty/math/isnan.c
common/thirdparty/math/log.c
common/thirdparty/math/log10.c
common/thirdparty/math/mtherr.c
common/thirdparty/math/polevl.c
common/thirdparty/math/pow.c
common/thirdparty/math/powi.c
common/thirdparty/math/sin.c
common/thirdparty/math/sinh.c
common/thirdparty/math/sqrt.c
common/thirdparty/math/tan.c
common/thirdparty/math/tanh.c
common/thirdparty/math/fastsin.cpp
)
#set_source_files_properties( ${FASTMATH_SOURCES} PROPERTIES COMPILE_FLAGS ${DEM_FASTMATH_FLAG} )
@ -881,6 +909,8 @@ include_directories(
core/rendering/hwrenderer/utility
core/rendering
platform
common/thirdparty
common/utility
${CMAKE_BINARY_DIR}/libraries/gdtoa
@ -990,6 +1020,12 @@ source_group("Core\\Rendering\\GL\\System" REGULAR_EXPRESSION "^${CMAKE_CURRENT_
source_group("Platform" REGULAR_EXPRESSION "^${CMAKE_CURRENT_SOURCE_DIR}/platform/.+")
source_group("Platform\\Win32" REGULAR_EXPRESSION "^${CMAKE_CURRENT_SOURCE_DIR}/platform/win32/.+")
source_group("Platform\\POSIX" REGULAR_EXPRESSION "^${CMAKE_CURRENT_SOURCE_DIR}/platform/posix/.+")
source_group("Common" REGULAR_EXPRESSION "^${CMAKE_CURRENT_SOURCE_DIR}/common/.+")
source_group("Common\\Utility" REGULAR_EXPRESSION "^${CMAKE_CURRENT_SOURCE_DIR}/utility/.+")
source_group("Common\\Third Party" REGULAR_EXPRESSION "^${CMAKE_CURRENT_SOURCE_DIR}/common/thirdparty/.+")
source_group("Common\\Third Party\\Math" REGULAR_EXPRESSION "^${CMAKE_CURRENT_SOURCE_DIR}/common/thirdparty/math/.+")
source_group("Common\\Third Party\\RapidJSON" REGULAR_EXPRESSION "^${CMAKE_CURRENT_SOURCE_DIR}/common/thirdparty/rapidjson/.+")
source_group("Common\\Third Party\\SFMT" REGULAR_EXPRESSION "^${CMAKE_CURRENT_SOURCE_DIR}/common/thirdparty/sfmt/.+")
source_group("Utility\\Smackerdec" REGULAR_EXPRESSION "^${CMAKE_CURRENT_SOURCE_DIR}/smackerdec/.+")
source_group("Utility\\Smackerdec\\Headers" REGULAR_EXPRESSION "^${CMAKE_CURRENT_SOURCE_DIR}/libsmackerdec/include/.+")
source_group("Utility\\Smackerdec\\Sources" REGULAR_EXPRESSION "^${CMAKE_CURRENT_SOURCE_DIR}/libsmackerdec/src/.+")

View file

@ -614,7 +614,7 @@ bool GameInterface::SaveGame(FSaveGameNode* node)
dword_27AA38 = 0;
}
}
catch (std::runtime_error & err)
catch (CRecoverableError & err)
{
// Let's not abort for write errors.
Printf(TEXTCOLOR_RED "%s\n", err.what());

View file

@ -432,6 +432,7 @@ defined __x86_64__ || defined __amd64__ || defined _M_X64 || defined _M_IA64 ||
#endif
#endif
#include "engineerrors.h"
////////// DEPRECATED: Standard library prefixing //////////
@ -508,17 +509,9 @@ static inline int Blrintf(const double x)
# define Bsqrtf sqrtf
#endif
class ExitEvent : public std::exception
{
int reason;
public:
ExitEvent(int a) { reason = a; }
int Reason() const { return reason; }
};
inline void Bexit(int a)
{
throw ExitEvent(a);
throw CExitEvent(a);
}

251
source/common/thirdparty/ctpl.h vendored Normal file
View file

@ -0,0 +1,251 @@
/*********************************************************
*
* Copyright (C) 2014 by Vitaliy Vitsentiy
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*********************************************************/
#ifndef __ctpl_stl_thread_pool_H__
#define __ctpl_stl_thread_pool_H__
#include <functional>
#include <thread>
#include <atomic>
#include <vector>
#include <memory>
#include <exception>
#include <future>
#include <mutex>
#include <queue>
// thread pool to run user's functors with signature
// ret func(int id, other_params)
// where id is the index of the thread that runs the functor
// ret is some return type
namespace ctpl {
namespace detail {
template <typename T>
class Queue {
public:
bool push(T const & value) {
std::unique_lock<std::mutex> lock(this->mutex);
this->q.push(value);
return true;
}
// deletes the retrieved element, do not use for non integral types
bool pop(T & v) {
std::unique_lock<std::mutex> lock(this->mutex);
if (this->q.empty())
return false;
v = this->q.front();
this->q.pop();
return true;
}
bool empty() {
std::unique_lock<std::mutex> lock(this->mutex);
return this->q.empty();
}
private:
std::queue<T> q;
std::mutex mutex;
};
}
class thread_pool {
public:
thread_pool() { this->init(); }
thread_pool(int nThreads) { this->init(); this->resize(nThreads); }
// the destructor waits for all the functions in the queue to be finished
~thread_pool() {
this->stop(true);
}
// get the number of running threads in the pool
int size() { return static_cast<int>(this->threads.size()); }
// number of idle threads
int n_idle() { return this->nWaiting; }
std::thread & get_thread(int i) { return *this->threads[i]; }
// change the number of threads in the pool
// should be called from one thread, otherwise be careful to not interleave, also with this->stop()
// nThreads must be >= 0
void resize(int nThreads) {
if (!this->isStop && !this->isDone) {
int oldNThreads = static_cast<int>(this->threads.size());
if (oldNThreads <= nThreads) { // if the number of threads is increased
this->threads.resize(nThreads);
this->flags.resize(nThreads);
for (int i = oldNThreads; i < nThreads; ++i) {
this->flags[i] = std::make_shared<std::atomic<bool>>(false);
this->set_thread(i);
}
}
else { // the number of threads is decreased
for (int i = oldNThreads - 1; i >= nThreads; --i) {
*this->flags[i] = true; // this thread will finish
this->threads[i]->detach();
}
{
// stop the detached threads that were waiting
std::unique_lock<std::mutex> lock(this->mutex);
this->cv.notify_all();
}
this->threads.resize(nThreads); // safe to delete because the threads are detached
this->flags.resize(nThreads); // safe to delete because the threads have copies of shared_ptr of the flags, not originals
}
}
}
// empty the queue
void clear_queue() {
std::function<void(int id)> * _f;
while (this->q.pop(_f))
delete _f; // empty the queue
}
// pops a functional wrapper to the original function
std::function<void(int)> pop() {
std::function<void(int id)> * _f = nullptr;
this->q.pop(_f);
std::unique_ptr<std::function<void(int id)>> func(_f); // at return, delete the function even if an exception occurred
std::function<void(int)> f;
if (_f)
f = *_f;
return f;
}
// wait for all computing threads to finish and stop all threads
// may be called asynchronously to not pause the calling thread while waiting
// if isWait == true, all the functions in the queue are run, otherwise the queue is cleared without running the functions
void stop(bool isWait = false) {
if (!isWait) {
if (this->isStop)
return;
this->isStop = true;
for (int i = 0, n = this->size(); i < n; ++i) {
*this->flags[i] = true; // command the threads to stop
}
this->clear_queue(); // empty the queue
}
else {
if (this->isDone || this->isStop)
return;
this->isDone = true; // give the waiting threads a command to finish
}
{
std::unique_lock<std::mutex> lock(this->mutex);
this->cv.notify_all(); // stop all waiting threads
}
for (int i = 0; i < static_cast<int>(this->threads.size()); ++i) { // wait for the computing threads to finish
if (this->threads[i]->joinable())
this->threads[i]->join();
}
// if there were no threads in the pool but some functors in the queue, the functors are not deleted by the threads
// therefore delete them here
this->clear_queue();
this->threads.clear();
this->flags.clear();
}
template<typename F, typename... Rest>
auto push(F && f, Rest&&... rest) ->std::future<decltype(f(0, rest...))> {
auto pck = std::make_shared<std::packaged_task<decltype(f(0, rest...))(int)>>(
std::bind(std::forward<F>(f), std::placeholders::_1, std::forward<Rest>(rest)...)
);
auto _f = new std::function<void(int id)>([pck](int id) {
(*pck)(id);
});
this->q.push(_f);
std::unique_lock<std::mutex> lock(this->mutex);
this->cv.notify_one();
return pck->get_future();
}
// run the user's function that excepts argument int - id of the running thread. returned value is templatized
// operator returns std::future, where the user can get the result and rethrow the catched exceptins
template<typename F>
auto push(F && f) ->std::future<decltype(f(0))> {
auto pck = std::make_shared<std::packaged_task<decltype(f(0))(int)>>(std::forward<F>(f));
auto _f = new std::function<void(int id)>([pck](int id) {
(*pck)(id);
});
this->q.push(_f);
std::unique_lock<std::mutex> lock(this->mutex);
this->cv.notify_one();
return pck->get_future();
}
private:
// deleted
thread_pool(const thread_pool &);// = delete;
thread_pool(thread_pool &&);// = delete;
thread_pool & operator=(const thread_pool &);// = delete;
thread_pool & operator=(thread_pool &&);// = delete;
void set_thread(int i) {
std::shared_ptr<std::atomic<bool>> flag(this->flags[i]); // a copy of the shared ptr to the flag
auto f = [this, i, flag/* a copy of the shared ptr to the flag */]() {
std::atomic<bool> & _flag = *flag;
std::function<void(int id)> * _f;
bool isPop = this->q.pop(_f);
while (true) {
while (isPop) { // if there is anything in the queue
std::unique_ptr<std::function<void(int id)>> func(_f); // at return, delete the function even if an exception occurred
(*_f)(i);
if (_flag)
return; // the thread is wanted to stop, return even if the queue is not empty yet
else
isPop = this->q.pop(_f);
}
// the queue is empty here, wait for the next command
std::unique_lock<std::mutex> lock(this->mutex);
++this->nWaiting;
this->cv.wait(lock, [this, &_f, &isPop, &_flag](){ isPop = this->q.pop(_f); return isPop || this->isDone || _flag; });
--this->nWaiting;
if (!isPop)
return; // if the queue is empty and this->isDone == true or *flag then return
}
};
this->threads[i].reset(new std::thread(f)); // compiler may not support std::make_unique()
}
void init() { this->nWaiting = 0; this->isStop = false; this->isDone = false; }
std::vector<std::unique_ptr<std::thread>> threads;
std::vector<std::shared_ptr<std::atomic<bool>>> flags;
detail::Queue<std::function<void(int id)> *> q;
std::atomic<bool> isDone;
std::atomic<bool> isStop;
std::atomic<int> nWaiting; // how many threads are waiting
std::mutex mutex;
std::condition_variable cv;
};
}
#endif // __ctpl_stl_thread_pool_H__

790
source/common/thirdparty/earcut.hpp vendored Normal file
View file

@ -0,0 +1,790 @@
/*
ISC License
Copyright (c) 2015, Mapbox
Permission to use, copy, modify, and/or distribute this software for any purpose
with or without fee is hereby granted, provided that the above copyright notice
and this permission notice appear in all copies.
THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH
REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,
INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.
*/
#pragma once
#include <algorithm>
#include <cassert>
#include <cmath>
#include <memory>
#include <vector>
namespace mapbox {
namespace util {
template <std::size_t I, typename T> struct nth {
inline static typename std::tuple_element<I, T>::type
get(const T& t) { return std::get<I>(t); };
};
}
namespace detail {
template <typename N = uint32_t>
class Earcut {
public:
std::vector<N> indices;
std::size_t vertices = 0;
template <typename Polygon>
void operator()(const Polygon& points);
private:
struct Node {
Node(N index, double x_, double y_) : i(index), x(x_), y(y_) {}
Node(const Node&) = delete;
Node& operator=(const Node&) = delete;
Node(Node&&) = delete;
Node& operator=(Node&&) = delete;
const N i;
const double x;
const double y;
// previous and next vertice nodes in a polygon ring
Node* prev = nullptr;
Node* next = nullptr;
// z-order curve value
int32_t z = 0;
// previous and next nodes in z-order
Node* prevZ = nullptr;
Node* nextZ = nullptr;
// indicates whether this is a steiner point
bool steiner = false;
};
template <typename Ring> Node* linkedList(const Ring& points, const bool clockwise);
Node* filterPoints(Node* start, Node* end = nullptr);
void earcutLinked(Node* ear, int pass = 0);
bool isEar(Node* ear);
bool isEarHashed(Node* ear);
Node* cureLocalIntersections(Node* start);
void splitEarcut(Node* start);
template <typename Polygon> Node* eliminateHoles(const Polygon& points, Node* outerNode);
void eliminateHole(Node* hole, Node* outerNode);
Node* findHoleBridge(Node* hole, Node* outerNode);
void indexCurve(Node* start);
Node* sortLinked(Node* list);
int32_t zOrder(const double x_, const double y_);
Node* getLeftmost(Node* start);
bool pointInTriangle(double ax, double ay, double bx, double by, double cx, double cy, double px, double py) const;
bool isValidDiagonal(Node* a, Node* b);
double area(const Node* p, const Node* q, const Node* r) const;
bool equals(const Node* p1, const Node* p2);
bool intersects(const Node* p1, const Node* q1, const Node* p2, const Node* q2);
bool intersectsPolygon(const Node* a, const Node* b);
bool locallyInside(const Node* a, const Node* b);
bool middleInside(const Node* a, const Node* b);
Node* splitPolygon(Node* a, Node* b);
template <typename Point> Node* insertNode(std::size_t i, const Point& p, Node* last);
void removeNode(Node* p);
bool hashing;
double minX, maxX;
double minY, maxY;
double inv_size = 0;
template <typename T, typename Alloc = std::allocator<T>>
class ObjectPool {
public:
ObjectPool() { }
ObjectPool(std::size_t blockSize_) {
reset(blockSize_);
}
~ObjectPool() {
clear();
}
template <typename... Args>
T* construct(Args&&... args) {
if (currentIndex >= blockSize) {
currentBlock = alloc.allocate(blockSize);
allocations.emplace_back(currentBlock);
currentIndex = 0;
}
T* object = &currentBlock[currentIndex++];
alloc.construct(object, std::forward<Args>(args)...);
return object;
}
void reset(std::size_t newBlockSize) {
for (auto allocation : allocations) alloc.deallocate(allocation, blockSize);
allocations.clear();
blockSize = std::max<std::size_t>(1, newBlockSize);
currentBlock = nullptr;
currentIndex = blockSize;
}
void clear() { reset(blockSize); }
private:
T* currentBlock = nullptr;
std::size_t currentIndex = 1;
std::size_t blockSize = 1;
std::vector<T*> allocations;
Alloc alloc;
};
ObjectPool<Node> nodes;
};
template <typename N> template <typename Polygon>
void Earcut<N>::operator()(const Polygon& points) {
// reset
indices.clear();
vertices = 0;
if (points.empty()) return;
double x;
double y;
int threshold = 80;
std::size_t len = 0;
for (size_t i = 0; threshold >= 0 && i < points.size(); i++) {
threshold -= static_cast<int>(points[i].size());
len += points[i].size();
}
//estimate size of nodes and indices
nodes.reset(len * 3 / 2);
indices.reserve(len + points[0].size());
Node* outerNode = linkedList(points[0], true);
if (!outerNode) return;
if (points.size() > 1) outerNode = eliminateHoles(points, outerNode);
// if the shape is not too simple, we'll use z-order curve hash later; calculate polygon bbox
hashing = threshold < 0;
if (hashing) {
Node* p = outerNode->next;
minX = maxX = outerNode->x;
minY = maxY = outerNode->y;
do {
x = p->x;
y = p->y;
minX = std::min<double>(minX, x);
minY = std::min<double>(minY, y);
maxX = std::max<double>(maxX, x);
maxY = std::max<double>(maxY, y);
p = p->next;
} while (p != outerNode);
// minX, minY and size are later used to transform coords into integers for z-order calculation
inv_size = std::max<double>(maxX - minX, maxY - minY);
inv_size = inv_size != .0 ? (1. / inv_size) : .0;
}
earcutLinked(outerNode);
nodes.clear();
}
// create a circular doubly linked list from polygon points in the specified winding order
template <typename N> template <typename Ring>
typename Earcut<N>::Node*
Earcut<N>::linkedList(const Ring& points, const bool clockwise) {
using Point = typename Ring::value_type;
double sum = 0;
const std::size_t len = points.size();
std::size_t i, j;
Node* last = nullptr;
// calculate original winding order of a polygon ring
for (i = 0, j = len > 0 ? len - 1 : 0; i < len; j = i++) {
const auto& p1 = points[i];
const auto& p2 = points[j];
const double p20 = util::nth<0, Point>::get(p2);
const double p10 = util::nth<0, Point>::get(p1);
const double p11 = util::nth<1, Point>::get(p1);
const double p21 = util::nth<1, Point>::get(p2);
sum += (p20 - p10) * (p11 + p21);
}
// link points into circular doubly-linked list in the specified winding order
if (clockwise == (sum > 0)) {
for (i = 0; i < len; i++) last = insertNode(vertices + i, points[i], last);
} else {
for (i = len; i-- > 0;) last = insertNode(vertices + i, points[i], last);
}
if (last && equals(last, last->next)) {
removeNode(last);
last = last->next;
}
vertices += len;
return last;
}
// eliminate colinear or duplicate points
template <typename N>
typename Earcut<N>::Node*
Earcut<N>::filterPoints(Node* start, Node* end) {
if (!end) end = start;
Node* p = start;
bool again;
do {
again = false;
if (!p->steiner && (equals(p, p->next) /*|| area(p->prev, p, p->next) == 0*/))
{
removeNode(p);
p = end = p->prev;
if (p == p->next) break;
again = true;
} else {
p = p->next;
}
} while (again || p != end);
return end;
}
// main ear slicing loop which triangulates a polygon (given as a linked list)
template <typename N>
void Earcut<N>::earcutLinked(Node* ear, int pass) {
if (!ear) return;
// interlink polygon nodes in z-order
if (!pass && hashing) indexCurve(ear);
Node* stop = ear;
Node* prev;
Node* next;
int iterations = 0;
// iterate through ears, slicing them one by one
while (ear->prev != ear->next) {
iterations++;
prev = ear->prev;
next = ear->next;
if (hashing ? isEarHashed(ear) : isEar(ear)) {
// cut off the triangle
indices.emplace_back(prev->i);
indices.emplace_back(ear->i);
indices.emplace_back(next->i);
removeNode(ear);
// skipping the next vertice leads to less sliver triangles
ear = next->next;
stop = next->next;
continue;
}
ear = next;
// if we looped through the whole remaining polygon and can't find any more ears
if (ear == stop) {
// try filtering points and slicing again
if (!pass) earcutLinked(filterPoints(ear), 1);
// if this didn't work, try curing all small self-intersections locally
else if (pass == 1) {
ear = cureLocalIntersections(ear);
earcutLinked(ear, 2);
// as a last resort, try splitting the remaining polygon into two
} else if (pass == 2) splitEarcut(ear);
break;
}
}
}
// check whether a polygon node forms a valid ear with adjacent nodes
template <typename N>
bool Earcut<N>::isEar(Node* ear) {
const Node* a = ear->prev;
const Node* b = ear;
const Node* c = ear->next;
if (area(a, b, c) >= 0) return false; // reflex, can't be an ear
// now make sure we don't have other points inside the potential ear
Node* p = ear->next->next;
while (p != ear->prev) {
if (pointInTriangle(a->x, a->y, b->x, b->y, c->x, c->y, p->x, p->y) &&
area(p->prev, p, p->next) >= 0) return false;
p = p->next;
}
return true;
}
template <typename N>
bool Earcut<N>::isEarHashed(Node* ear) {
const Node* a = ear->prev;
const Node* b = ear;
const Node* c = ear->next;
if (area(a, b, c) >= 0) return false; // reflex, can't be an ear
// triangle bbox; min & max are calculated like this for speed
const double minTX = std::min<double>(a->x, std::min<double>(b->x, c->x));
const double minTY = std::min<double>(a->y, std::min<double>(b->y, c->y));
const double maxTX = std::max<double>(a->x, std::max<double>(b->x, c->x));
const double maxTY = std::max<double>(a->y, std::max<double>(b->y, c->y));
// z-order range for the current triangle bbox;
const int32_t minZ = zOrder(minTX, minTY);
const int32_t maxZ = zOrder(maxTX, maxTY);
// first look for points inside the triangle in increasing z-order
Node* p = ear->nextZ;
while (p && p->z <= maxZ) {
if (p != ear->prev && p != ear->next &&
pointInTriangle(a->x, a->y, b->x, b->y, c->x, c->y, p->x, p->y) &&
area(p->prev, p, p->next) >= 0) return false;
p = p->nextZ;
}
// then look for points in decreasing z-order
p = ear->prevZ;
while (p && p->z >= minZ) {
if (p != ear->prev && p != ear->next &&
pointInTriangle(a->x, a->y, b->x, b->y, c->x, c->y, p->x, p->y) &&
area(p->prev, p, p->next) >= 0) return false;
p = p->prevZ;
}
return true;
}
// go through all polygon nodes and cure small local self-intersections
template <typename N>
typename Earcut<N>::Node*
Earcut<N>::cureLocalIntersections(Node* start) {
Node* p = start;
do {
Node* a = p->prev;
Node* b = p->next->next;
// a self-intersection where edge (v[i-1],v[i]) intersects (v[i+1],v[i+2])
if (!equals(a, b) && intersects(a, p, p->next, b) && locallyInside(a, b) && locallyInside(b, a)) {
indices.emplace_back(a->i);
indices.emplace_back(p->i);
indices.emplace_back(b->i);
// remove two nodes involved
removeNode(p);
removeNode(p->next);
p = start = b;
}
p = p->next;
} while (p != start);
return p;
}
// try splitting polygon into two and triangulate them independently
template <typename N>
void Earcut<N>::splitEarcut(Node* start) {
// look for a valid diagonal that divides the polygon into two
Node* a = start;
do {
Node* b = a->next->next;
while (b != a->prev) {
if (a->i != b->i && isValidDiagonal(a, b)) {
// split the polygon in two by the diagonal
Node* c = splitPolygon(a, b);
// filter colinear points around the cuts
a = filterPoints(a, a->next);
c = filterPoints(c, c->next);
// run earcut on each half
earcutLinked(a);
earcutLinked(c);
return;
}
b = b->next;
}
a = a->next;
} while (a != start);
}
// link every hole into the outer loop, producing a single-ring polygon without holes
template <typename N> template <typename Polygon>
typename Earcut<N>::Node*
Earcut<N>::eliminateHoles(const Polygon& points, Node* outerNode) {
const size_t len = points.size();
std::vector<Node*> queue;
for (size_t i = 1; i < len; i++) {
Node* list = linkedList(points[i], false);
if (list) {
if (list == list->next) list->steiner = true;
queue.push_back(getLeftmost(list));
}
}
std::sort(queue.begin(), queue.end(), [](const Node* a, const Node* b) {
return a->x < b->x;
});
// process holes from left to right
for (size_t i = 0; i < queue.size(); i++) {
eliminateHole(queue[i], outerNode);
outerNode = filterPoints(outerNode, outerNode->next);
}
return outerNode;
}
// find a bridge between vertices that connects hole with an outer ring and and link it
template <typename N>
void Earcut<N>::eliminateHole(Node* hole, Node* outerNode) {
outerNode = findHoleBridge(hole, outerNode);
if (outerNode) {
Node* b = splitPolygon(outerNode, hole);
filterPoints(b, b->next);
}
}
// David Eberly's algorithm for finding a bridge between hole and outer polygon
template <typename N>
typename Earcut<N>::Node*
Earcut<N>::findHoleBridge(Node* hole, Node* outerNode) {
Node* p = outerNode;
double hx = hole->x;
double hy = hole->y;
double qx = -std::numeric_limits<double>::infinity();
Node* m = nullptr;
// find a segment intersected by a ray from the hole's leftmost Vertex to the left;
// segment's endpoint with lesser x will be potential connection Vertex
do {
if (hy <= p->y && hy >= p->next->y && p->next->y != p->y) {
double x = p->x + (hy - p->y) * (p->next->x - p->x) / (p->next->y - p->y);
if (x <= hx && x > qx) {
qx = x;
if (x == hx) {
if (hy == p->y) return p;
if (hy == p->next->y) return p->next;
}
m = p->x < p->next->x ? p : p->next;
}
}
p = p->next;
} while (p != outerNode);
if (!m) return 0;
if (hx == qx) return m->prev;
// look for points inside the triangle of hole Vertex, segment intersection and endpoint;
// if there are no points found, we have a valid connection;
// otherwise choose the Vertex of the minimum angle with the ray as connection Vertex
const Node* stop = m;
double tanMin = std::numeric_limits<double>::infinity();
double tanCur = 0;
p = m->next;
double mx = m->x;
double my = m->y;
while (p != stop) {
if (hx >= p->x && p->x >= mx && hx != p->x &&
pointInTriangle(hy < my ? hx : qx, hy, mx, my, hy < my ? qx : hx, hy, p->x, p->y)) {
tanCur = std::abs(hy - p->y) / (hx - p->x); // tangential
if ((tanCur < tanMin || (tanCur == tanMin && p->x > m->x)) && locallyInside(p, hole)) {
m = p;
tanMin = tanCur;
}
}
p = p->next;
}
return m;
}
// interlink polygon nodes in z-order
template <typename N>
void Earcut<N>::indexCurve(Node* start) {
assert(start);
Node* p = start;
do {
p->z = p->z ? p->z : zOrder(p->x, p->y);
p->prevZ = p->prev;
p->nextZ = p->next;
p = p->next;
} while (p != start);
p->prevZ->nextZ = nullptr;
p->prevZ = nullptr;
sortLinked(p);
}
// Simon Tatham's linked list merge sort algorithm
// http://www.chiark.greenend.org.uk/~sgtatham/algorithms/listsort.html
template <typename N>
typename Earcut<N>::Node*
Earcut<N>::sortLinked(Node* list) {
assert(list);
Node* p;
Node* q;
Node* e;
Node* tail;
int i, numMerges, pSize, qSize;
int inSize = 1;
for (;;) {
p = list;
list = nullptr;
tail = nullptr;
numMerges = 0;
while (p) {
numMerges++;
q = p;
pSize = 0;
for (i = 0; i < inSize; i++) {
pSize++;
q = q->nextZ;
if (!q) break;
}
qSize = inSize;
while (pSize > 0 || (qSize > 0 && q)) {
if (pSize == 0) {
e = q;
q = q->nextZ;
qSize--;
} else if (qSize == 0 || !q) {
e = p;
p = p->nextZ;
pSize--;
} else if (p->z <= q->z) {
e = p;
p = p->nextZ;
pSize--;
} else {
e = q;
q = q->nextZ;
qSize--;
}
if (tail) tail->nextZ = e;
else list = e;
e->prevZ = tail;
tail = e;
}
p = q;
}
tail->nextZ = nullptr;
if (numMerges <= 1) return list;
inSize *= 2;
}
}
// z-order of a Vertex given coords and size of the data bounding box
template <typename N>
int32_t Earcut<N>::zOrder(const double x_, const double y_) {
// coords are transformed into non-negative 15-bit integer range
int32_t x = static_cast<int32_t>(32767.0 * (x_ - minX) * inv_size);
int32_t y = static_cast<int32_t>(32767.0 * (y_ - minY) * inv_size);
x = (x | (x << 8)) & 0x00FF00FF;
x = (x | (x << 4)) & 0x0F0F0F0F;
x = (x | (x << 2)) & 0x33333333;
x = (x | (x << 1)) & 0x55555555;
y = (y | (y << 8)) & 0x00FF00FF;
y = (y | (y << 4)) & 0x0F0F0F0F;
y = (y | (y << 2)) & 0x33333333;
y = (y | (y << 1)) & 0x55555555;
return x | (y << 1);
}
// find the leftmost node of a polygon ring
template <typename N>
typename Earcut<N>::Node*
Earcut<N>::getLeftmost(Node* start) {
Node* p = start;
Node* leftmost = start;
do {
if (p->x < leftmost->x) leftmost = p;
p = p->next;
} while (p != start);
return leftmost;
}
// check if a point lies within a convex triangle
template <typename N>
bool Earcut<N>::pointInTriangle(double ax, double ay, double bx, double by, double cx, double cy, double px, double py) const {
return (cx - px) * (ay - py) - (ax - px) * (cy - py) >= 0 &&
(ax - px) * (by - py) - (bx - px) * (ay - py) >= 0 &&
(bx - px) * (cy - py) - (cx - px) * (by - py) >= 0;
}
// check if a diagonal between two polygon nodes is valid (lies in polygon interior)
template <typename N>
bool Earcut<N>::isValidDiagonal(Node* a, Node* b) {
return a->next->i != b->i && a->prev->i != b->i && !intersectsPolygon(a, b) &&
locallyInside(a, b) && locallyInside(b, a) && middleInside(a, b);
}
// signed area of a triangle
template <typename N>
double Earcut<N>::area(const Node* p, const Node* q, const Node* r) const {
return (q->y - p->y) * (r->x - q->x) - (q->x - p->x) * (r->y - q->y);
}
// check if two points are equal
template <typename N>
bool Earcut<N>::equals(const Node* p1, const Node* p2) {
return p1->x == p2->x && p1->y == p2->y;
}
// check if two segments intersect
template <typename N>
bool Earcut<N>::intersects(const Node* p1, const Node* q1, const Node* p2, const Node* q2) {
if ((equals(p1, q1) && equals(p2, q2)) ||
(equals(p1, q2) && equals(p2, q1))) return true;
return (area(p1, q1, p2) > 0) != (area(p1, q1, q2) > 0) &&
(area(p2, q2, p1) > 0) != (area(p2, q2, q1) > 0);
}
// check if a polygon diagonal intersects any polygon segments
template <typename N>
bool Earcut<N>::intersectsPolygon(const Node* a, const Node* b) {
const Node* p = a;
do {
if (p->i != a->i && p->next->i != a->i && p->i != b->i && p->next->i != b->i &&
intersects(p, p->next, a, b)) return true;
p = p->next;
} while (p != a);
return false;
}
// check if a polygon diagonal is locally inside the polygon
template <typename N>
bool Earcut<N>::locallyInside(const Node* a, const Node* b) {
return area(a->prev, a, a->next) < 0 ?
area(a, b, a->next) >= 0 && area(a, a->prev, b) >= 0 :
area(a, b, a->prev) < 0 || area(a, a->next, b) < 0;
}
// check if the middle Vertex of a polygon diagonal is inside the polygon
template <typename N>
bool Earcut<N>::middleInside(const Node* a, const Node* b) {
const Node* p = a;
bool inside = false;
double px = (a->x + b->x) / 2;
double py = (a->y + b->y) / 2;
do {
if (((p->y > py) != (p->next->y > py)) && p->next->y != p->y &&
(px < (p->next->x - p->x) * (py - p->y) / (p->next->y - p->y) + p->x))
inside = !inside;
p = p->next;
} while (p != a);
return inside;
}
// link two polygon vertices with a bridge; if the vertices belong to the same ring, it splits
// polygon into two; if one belongs to the outer ring and another to a hole, it merges it into a
// single ring
template <typename N>
typename Earcut<N>::Node*
Earcut<N>::splitPolygon(Node* a, Node* b) {
Node* a2 = nodes.construct(a->i, a->x, a->y);
Node* b2 = nodes.construct(b->i, b->x, b->y);
Node* an = a->next;
Node* bp = b->prev;
a->next = b;
b->prev = a;
a2->next = an;
an->prev = a2;
b2->next = a2;
a2->prev = b2;
bp->next = b2;
b2->prev = bp;
return b2;
}
// create a node and util::optionally link it with previous one (in a circular doubly linked list)
template <typename N> template <typename Point>
typename Earcut<N>::Node*
Earcut<N>::insertNode(std::size_t i, const Point& pt, Node* last) {
Node* p = nodes.construct(static_cast<N>(i), util::nth<0, Point>::get(pt), util::nth<1, Point>::get(pt));
if (!last) {
p->prev = p;
p->next = p;
} else {
assert(last);
p->next = last->next;
p->prev = last;
last->next->prev = p;
last->next = p;
}
return p;
}
template <typename N>
void Earcut<N>::removeNode(Node* p) {
p->next->prev = p->prev;
p->prev->next = p->next;
if (p->prevZ) p->prevZ->nextZ = p->nextZ;
if (p->nextZ) p->nextZ->prevZ = p->prevZ;
}
}
template <typename N = uint32_t, typename Polygon>
std::vector<N> earcut(const Polygon& poly) {
mapbox::detail::Earcut<N> earcut;
earcut(poly);
return std::move(earcut.indices);
}
}

339
source/common/thirdparty/math/asin.c vendored Normal file
View file

@ -0,0 +1,339 @@
/* asin.c
*
* Inverse circular sine
*
*
*
* SYNOPSIS:
*
* double x, y, asin();
*
* y = asin( x );
*
*
*
* DESCRIPTION:
*
* Returns radian angle between -pi/2 and +pi/2 whose sine is x.
*
* A rational function of the form x + x**3 P(x**2)/Q(x**2)
* is used for |x| in the interval [0, 0.5]. If |x| > 0.5 it is
* transformed by the identity
*
* asin(x) = pi/2 - 2 asin( sqrt( (1-x)/2 ) ).
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC -1, 1 40000 2.6e-17 7.1e-18
* IEEE -1, 1 10^6 1.9e-16 5.4e-17
*
*
* ERROR MESSAGES:
*
* message condition value returned
* asin domain |x| > 1 NAN
*
*/
/* acos()
*
* Inverse circular cosine
*
*
*
* SYNOPSIS:
*
* double x, y, acos();
*
* y = acos( x );
*
*
*
* DESCRIPTION:
*
* Returns radian angle between 0 and pi whose cosine
* is x.
*
* Analytically, acos(x) = pi/2 - asin(x). However if |x| is
* near 1, there is cancellation error in subtracting asin(x)
* from pi/2. Hence if x < -0.5,
*
* acos(x) = pi - 2.0 * asin( sqrt((1+x)/2) );
*
* or if x > +0.5,
*
* acos(x) = 2.0 * asin( sqrt((1-x)/2) ).
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC -1, 1 50000 3.3e-17 8.2e-18
* IEEE -1, 1 10^6 2.2e-16 6.5e-17
*
*
* ERROR MESSAGES:
*
* message condition value returned
* asin domain |x| > 1 NAN
*/
/* asin.c */
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. Neither the name of the <ORGANIZATION> nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
*/
#include "mconf.h"
/* arcsin(x) = x + x^3 P(x^2)/Q(x^2)
0 <= x <= 0.625
Peak relative error = 1.2e-18 */
#if UNK
static double P[6] = {
4.253011369004428248960E-3,
-6.019598008014123785661E-1,
5.444622390564711410273E0,
-1.626247967210700244449E1,
1.956261983317594739197E1,
-8.198089802484824371615E0,
};
static double Q[5] = {
/* 1.000000000000000000000E0, */
-1.474091372988853791896E1,
7.049610280856842141659E1,
-1.471791292232726029859E2,
1.395105614657485689735E2,
-4.918853881490881290097E1,
};
#endif
#if DEC
static short P[24] = {
0036213,0056330,0057244,0053234,
0140032,0015011,0114762,0160255,
0040656,0035130,0136121,0067313,
0141202,0014616,0170474,0101731,
0041234,0100076,0151674,0111310,
0141003,0025540,0033165,0077246,
};
static short Q[20] = {
/* 0040200,0000000,0000000,0000000, */
0141153,0155310,0055360,0072530,
0041614,0177001,0027764,0101237,
0142023,0026733,0064653,0133266,
0042013,0101264,0023775,0176351,
0141504,0140420,0050660,0036543,
};
#endif
#if IBMPC
static short P[24] = {
0x8ad3,0x0bd4,0x6b9b,0x3f71,
0x5c16,0x333e,0x4341,0xbfe3,
0x2dd9,0x178a,0xc74b,0x4015,
0x907b,0xde27,0x4331,0xc030,
0x9259,0xda77,0x9007,0x4033,
0xafd5,0x06ce,0x656c,0xc020,
};
static short Q[20] = {
/* 0x0000,0x0000,0x0000,0x3ff0, */
0x0eab,0x0b5e,0x7b59,0xc02d,
0x9054,0x25fe,0x9fc0,0x4051,
0x76d7,0x6d35,0x65bb,0xc062,
0xbf9d,0x84ff,0x7056,0x4061,
0x07ac,0x0a36,0x9822,0xc048,
};
#endif
#if MIEEE
static short P[24] = {
0x3f71,0x6b9b,0x0bd4,0x8ad3,
0xbfe3,0x4341,0x333e,0x5c16,
0x4015,0xc74b,0x178a,0x2dd9,
0xc030,0x4331,0xde27,0x907b,
0x4033,0x9007,0xda77,0x9259,
0xc020,0x656c,0x06ce,0xafd5,
};
static short Q[20] = {
/* 0x3ff0,0x0000,0x0000,0x0000, */
0xc02d,0x7b59,0x0b5e,0x0eab,
0x4051,0x9fc0,0x25fe,0x9054,
0xc062,0x65bb,0x6d35,0x76d7,
0x4061,0x7056,0x84ff,0xbf9d,
0xc048,0x9822,0x0a36,0x07ac,
};
#endif
/* arcsin(1-x) = pi/2 - sqrt(2x)(1+R(x))
0 <= x <= 0.5
Peak relative error = 4.2e-18 */
#if UNK
static double R[5] = {
2.967721961301243206100E-3,
-5.634242780008963776856E-1,
6.968710824104713396794E0,
-2.556901049652824852289E1,
2.853665548261061424989E1,
};
static double S[4] = {
/* 1.000000000000000000000E0, */
-2.194779531642920639778E1,
1.470656354026814941758E2,
-3.838770957603691357202E2,
3.424398657913078477438E2,
};
#endif
#if DEC
static short R[20] = {
0036102,0077034,0142164,0174103,
0140020,0036222,0147711,0044173,
0040736,0177655,0153631,0171523,
0141314,0106525,0060015,0055474,
0041344,0045422,0003630,0040344,
};
static short S[16] = {
/* 0040200,0000000,0000000,0000000, */
0141257,0112425,0132772,0166136,
0042023,0010315,0075523,0175020,
0142277,0170104,0126203,0017563,
0042253,0034115,0102662,0022757,
};
#endif
#if IBMPC
static short R[20] = {
0x9f08,0x988e,0x4fc3,0x3f68,
0x290f,0x59f9,0x0792,0xbfe2,
0x3e6a,0xbaf3,0xdff5,0x401b,
0xab68,0xac01,0x91aa,0xc039,
0x081d,0x40f3,0x8962,0x403c,
};
static short S[16] = {
/* 0x0000,0x0000,0x0000,0x3ff0, */
0x5d8c,0xb6bf,0xf2a2,0xc035,
0x7f42,0xaf6a,0x6219,0x4062,
0x63ee,0x9590,0xfe08,0xc077,
0x44be,0xb0b6,0x6709,0x4075,
};
#endif
#if MIEEE
static short R[20] = {
0x3f68,0x4fc3,0x988e,0x9f08,
0xbfe2,0x0792,0x59f9,0x290f,
0x401b,0xdff5,0xbaf3,0x3e6a,
0xc039,0x91aa,0xac01,0xab68,
0x403c,0x8962,0x40f3,0x081d,
};
static short S[16] = {
/* 0x3ff0,0x0000,0x0000,0x0000, */
0xc035,0xf2a2,0xb6bf,0x5d8c,
0x4062,0x6219,0xaf6a,0x7f42,
0xc077,0xfe08,0x9590,0x63ee,
0x4075,0x6709,0xb0b6,0x44be,
};
#endif
/* pi/2 = PIO2 + MOREBITS. */
#ifdef DEC
#define MOREBITS 5.721188726109831840122E-18
#else
#define MOREBITS 6.123233995736765886130E-17
#endif
#ifdef ANSIPROT
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern double c_sqrt ( double );
double c_asin ( double );
#else
double c_sqrt(), polevl(), p1evl();
double c_asin();
#endif
extern double PIO2, PIO4, NAN;
double c_asin(x)
double x;
{
double a, p, z, zz;
short sign;
if( x > 0 )
{
sign = 1;
a = x;
}
else
{
sign = -1;
a = -x;
}
if( a > 1.0 )
{
mtherr( "asin", DOMAIN );
return( NAN );
}
if( a > 0.625 )
{
/* arcsin(1-x) = pi/2 - sqrt(2x)(1+R(x)) */
zz = 1.0 - a;
p = zz * polevl( zz, R, 4)/p1evl( zz, S, 4);
zz = c_sqrt(zz+zz);
z = PIO4 - zz;
zz = zz * p - MOREBITS;
z = z - zz;
z = z + PIO4;
}
else
{
if( a < 1.0e-8 )
{
return(x);
}
zz = a * a;
z = zz * polevl( zz, P, 5)/p1evl( zz, Q, 5);
z = a * z + a;
}
if( sign < 0 )
z = -z;
return(z);
}
double c_acos(x)
double x;
{
if( (x < -1.0) || (x > 1.0) )
{
mtherr( "acos", DOMAIN );
return( NAN );
}
return PIO2 - c_asin(x) + MOREBITS;
}

417
source/common/thirdparty/math/atan.c vendored Normal file
View file

@ -0,0 +1,417 @@
/* atan.c
*
* Inverse circular tangent
* (arctangent)
*
*
*
* SYNOPSIS:
*
* double x, y, atan();
*
* y = atan( x );
*
*
*
* DESCRIPTION:
*
* Returns radian angle between -pi/2 and +pi/2 whose tangent
* is x.
*
* Range reduction is from three intervals into the interval
* from zero to 0.66. The approximant uses a rational
* function of degree 4/5 of the form x + x**3 P(x)/Q(x).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC -10, 10 50000 2.4e-17 8.3e-18
* IEEE -10, 10 10^6 1.8e-16 5.0e-17
*
*/
/* atan2()
*
* Quadrant correct inverse circular tangent
*
*
*
* SYNOPSIS:
*
* double x, y, z, atan2();
*
* z = atan2( y, x );
*
*
*
* DESCRIPTION:
*
* Returns radian angle whose tangent is y/x.
* Define compile time symbol ANSIC = 1 for ANSI standard,
* range -PI < z <= +PI, args (y,x); else ANSIC = 0 for range
* 0 to 2PI, args (x,y).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -10, 10 10^6 2.5e-16 6.9e-17
* See atan.c.
*
*/
/* atan.c */
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. Neither the name of the <ORGANIZATION> nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
*/
#include "mconf.h"
/* arctan(x) = x + x^3 P(x^2)/Q(x^2)
0 <= x <= 0.66
Peak relative error = 2.6e-18 */
#ifdef UNK
static double P[5] = {
-8.750608600031904122785E-1,
-1.615753718733365076637E1,
-7.500855792314704667340E1,
-1.228866684490136173410E2,
-6.485021904942025371773E1,
};
static double Q[5] = {
/* 1.000000000000000000000E0, */
2.485846490142306297962E1,
1.650270098316988542046E2,
4.328810604912902668951E2,
4.853903996359136964868E2,
1.945506571482613964425E2,
};
/* tan( 3*pi/8 ) */
static double T3P8 = 2.41421356237309504880;
#endif
#ifdef DEC
static short P[20] = {
0140140,0001775,0007671,0026242,
0141201,0041242,0155534,0001715,
0141626,0002141,0132100,0011625,
0141765,0142771,0064055,0150453,
0141601,0131517,0164507,0062164,
};
static short Q[20] = {
/* 0040200,0000000,0000000,0000000, */
0041306,0157042,0154243,0000742,
0042045,0003352,0016707,0150452,
0042330,0070306,0113425,0170730,
0042362,0130770,0116602,0047520,
0042102,0106367,0156753,0013541,
};
/* tan( 3*pi/8 ) = 2.41421356237309504880 */
static unsigned short T3P8A[] = {040432,0101171,0114774,0167462,};
#define T3P8 *(double *)T3P8A
#endif
#ifdef IBMPC
static short P[20] = {
0x2594,0xa1f7,0x007f,0xbfec,
0x807a,0x5b6b,0x2854,0xc030,
0x0273,0x3688,0xc08c,0xc052,
0xba25,0x2d05,0xb8bf,0xc05e,
0xec8e,0xfd28,0x3669,0xc050,
};
static short Q[20] = {
/* 0x0000,0x0000,0x0000,0x3ff0, */
0x603c,0x5b14,0xdbc4,0x4038,
0xfa25,0x43b8,0xa0dd,0x4064,
0xbe3b,0xd2e2,0x0e18,0x407b,
0x49ea,0x13b0,0x563f,0x407e,
0x62ec,0xfbbd,0x519e,0x4068,
};
/* tan( 3*pi/8 ) = 2.41421356237309504880 */
static unsigned short T3P8A[] = {0x9de6,0x333f,0x504f,0x4003};
#define T3P8 *(double *)T3P8A
#endif
#ifdef MIEEE
static short P[20] = {
0xbfec,0x007f,0xa1f7,0x2594,
0xc030,0x2854,0x5b6b,0x807a,
0xc052,0xc08c,0x3688,0x0273,
0xc05e,0xb8bf,0x2d05,0xba25,
0xc050,0x3669,0xfd28,0xec8e,
};
static short Q[20] = {
/* 0x3ff0,0x0000,0x0000,0x0000, */
0x4038,0xdbc4,0x5b14,0x603c,
0x4064,0xa0dd,0x43b8,0xfa25,
0x407b,0x0e18,0xd2e2,0xbe3b,
0x407e,0x563f,0x13b0,0x49ea,
0x4068,0x519e,0xfbbd,0x62ec,
};
/* tan( 3*pi/8 ) = 2.41421356237309504880 */
static unsigned short T3P8A[] = {
0x4003,0x504f,0x333f,0x9de6
};
#define T3P8 *(double *)T3P8A
#endif
#ifdef ANSIPROT
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern double atan ( double );
extern double fabs ( double );
extern int signbit ( double );
extern int isnan ( double );
#else
double polevl(), p1evl(), atan(), fabs();
int signbit(), isnan();
#endif
extern double PI, PIO2, PIO4, INFINITY, NEGZERO, MAXNUM;
/* pi/2 = PIO2 + MOREBITS. */
#ifdef DEC
#define MOREBITS 5.721188726109831840122E-18
#else
#define MOREBITS 6.123233995736765886130E-17
#endif
double c_atan(x)
double x;
{
double y, z;
short sign, flag;
#ifdef MINUSZERO
if( x == 0.0 )
return(x);
#endif
#ifdef INFINITIES
if(x == INFINITY)
return(PIO2);
if(x == -INFINITY)
return(-PIO2);
#endif
/* make argument positive and save the sign */
sign = 1;
if( x < 0.0 )
{
sign = -1;
x = -x;
}
/* range reduction */
flag = 0;
if( x > T3P8 )
{
y = PIO2;
flag = 1;
x = -( 1.0/x );
}
else if( x <= 0.66 )
{
y = 0.0;
}
else
{
y = PIO4;
flag = 2;
x = (x-1.0)/(x+1.0);
}
z = x * x;
z = z * polevl( z, P, 4 ) / p1evl( z, Q, 5 );
z = x * z + x;
if( flag == 2 )
z += 0.5 * MOREBITS;
else if( flag == 1 )
z += MOREBITS;
y = y + z;
if( sign < 0 )
y = -y;
return(y);
}
/* atan2 */
#ifdef ANSIC
double c_atan2( y, x )
#else
double c_atan2( x, y )
#endif
double x, y;
{
double z, w;
short code;
code = 0;
#ifdef NANS
if( isnan(x) )
return(x);
if( isnan(y) )
return(y);
#endif
#ifdef MINUSZERO
if( y == 0.0 )
{
if( signbit(y) )
{
if( x > 0.0 )
z = y;
else if( x < 0.0 )
z = -PI;
else
{
if( signbit(x) )
z = -PI;
else
z = y;
}
}
else /* y is +0 */
{
if( x == 0.0 )
{
if( signbit(x) )
z = PI;
else
z = 0.0;
}
else if( x > 0.0 )
z = 0.0;
else
z = PI;
}
return z;
}
if( x == 0.0 )
{
if( y > 0.0 )
z = PIO2;
else
z = -PIO2;
return z;
}
#endif /* MINUSZERO */
#ifdef INFINITIES
if( x == INFINITY )
{
if( y == INFINITY )
z = 0.25 * PI;
else if( y == -INFINITY )
z = -0.25 * PI;
else if( y < 0.0 )
z = NEGZERO;
else
z = 0.0;
return z;
}
if( x == -INFINITY )
{
if( y == INFINITY )
z = 0.75 * PI;
else if( y <= -INFINITY )
z = -0.75 * PI;
else if( y >= 0.0 )
z = PI;
else
z = -PI;
return z;
}
if( y == INFINITY )
return( PIO2 );
if( y == -INFINITY )
return( -PIO2 );
#endif
if( x < 0.0 )
code = 2;
if( y < 0.0 )
code |= 1;
#ifdef INFINITIES
if( x == 0.0 )
#else
if( fabs(x) <= (fabs(y) / MAXNUM) )
#endif
{
if( code & 1 )
{
#if ANSIC
return( -PIO2 );
#else
return( 3.0*PIO2 );
#endif
}
if( y == 0.0 )
return( 0.0 );
return( PIO2 );
}
if( y == 0.0 )
{
if( code & 2 )
return( PI );
return( 0.0 );
}
switch( code )
{
#if ANSIC
default:
case 0:
case 1: w = 0.0; break;
case 2: w = PI; break;
case 3: w = -PI; break;
#else
default:
case 0: w = 0.0; break;
case 1: w = 2.0 * PI; break;
case 2:
case 3: w = PI; break;
#endif
}
z = w + c_atan( y/x );
#ifdef MINUSZERO
if( z == 0.0 && y < 0 )
z = NEGZERO;
#endif
return( z );
}

149
source/common/thirdparty/math/cmath.h vendored Normal file
View file

@ -0,0 +1,149 @@
#ifndef __CMATH_H
#define __CMATH_H
#include "xs_Float.h"
#define USE_CUSTOM_MATH // we want repreducably reliable results, even at the cost of performance
#define USE_FAST_MATH // use faster table-based sin and cos variants with limited precision (sufficient for Doom gameplay)
extern"C"
{
double c_asin(double);
double c_acos(double);
double c_atan(double);
double c_atan2(double, double);
double c_sin(double);
double c_cos(double);
double c_tan(double);
double c_cot(double);
double c_sqrt(double);
double c_sinh(double);
double c_cosh(double);
double c_tanh(double);
double c_exp(double);
double c_log(double);
double c_log10(double);
double c_pow(double, double);
}
// This uses a sine table with linear interpolation
// For in-game calculations this is precise enough
// and this code is more than 10x faster than the
// Cephes sin and cos function.
struct FFastTrig
{
static const int TBLPERIOD = 8192;
static const int BITSHIFT = 19;
static const int REMAINDER = (1 << BITSHIFT) - 1;
float sinetable[2049];
double sinq1(unsigned);
public:
FFastTrig();
double sin(unsigned);
double cos(unsigned);
};
extern FFastTrig fasttrig;
// This must use xs_Float to guarantee proper integer wraparound.
#define DEG2BAM(f) ((unsigned)xs_CRoundToInt((f) * (0x40000000/90.)))
#define RAD2BAM(f) ((unsigned)xs_CRoundToInt((f) * (0x80000000/3.14159265358979323846)))
inline double fastcosdeg(double v)
{
return fasttrig.cos(DEG2BAM(v));
}
inline double fastsindeg(double v)
{
return fasttrig.sin(DEG2BAM(v));
}
inline double fastcos(double v)
{
return fasttrig.cos(RAD2BAM(v));
}
inline double fastsin(double v)
{
return fasttrig.sin(RAD2BAM(v));
}
// these are supposed to be local to this file.
#undef DEG2BAM
#undef RAD2BAM
inline double sindeg(double v)
{
#ifdef USE_CUSTOM_MATH
return c_sin(v * (3.14159265358979323846 / 180.));
#else
return sin(v * (3.14159265358979323846 / 180.));
#endif
}
inline double cosdeg(double v)
{
#ifdef USE_CUSTOM_MATH
return c_cos(v * (3.14159265358979323846 / 180.));
#else
return cos(v * (3.14159265358979323846 / 180.));
#endif
}
#ifndef USE_CUSTOM_MATH
#define g_asin asin
#define g_acos acos
#define g_atan atan
#define g_atan2 atan2
#define g_sin sin
#define g_cos cos
#define g_sindeg sindeg
#define g_cosdeg cosdeg
#define g_tan tan
#define g_cot cot
#define g_sqrt sqrt
#define g_sinh sinh
#define g_cosh cosh
#define g_tanh tanh
#define g_exp exp
#define g_log log
#define g_log10 log10
#define g_pow pow
#else
#define g_asin c_asin
#define g_acos c_acos
#define g_atan c_atan
#define g_atan2 c_atan2
#ifndef USE_FAST_MATH
#define g_sindeg sindeg
#define g_cosdeg cosdeg
#define g_sin c_sin
#define g_cos c_cos
#else
#define g_sindeg fastsindeg
#define g_cosdeg fastcosdeg
#define g_sin fastsin
#define g_cos fastcos
#endif
#define g_tan c_tan
#define g_cot c_cot
#define g_sqrt c_sqrt
#define g_sinh c_sinh
#define g_cosh c_cosh
#define g_tanh c_tanh
#define g_exp c_exp
#define g_log c_log
#define g_log10 c_log10
#define g_pow c_pow
#endif
#endif

276
source/common/thirdparty/math/const.c vendored Normal file
View file

@ -0,0 +1,276 @@
/* const.c
*
* Globally declared constants
*
*
*
* SYNOPSIS:
*
* extern double nameofconstant;
*
*
*
*
* DESCRIPTION:
*
* This file contains a number of mathematical constants and
* also some needed size parameters of the computer arithmetic.
* The values are supplied as arrays of hexadecimal integers
* for IEEE arithmetic; arrays of octal constants for DEC
* arithmetic; and in a normal decimal scientific notation for
* other machines. The particular notation used is determined
* by a symbol (DEC, IBMPC, or UNK) defined in the include file
* mconf.h.
*
* The default size parameters are as follows.
*
* For DEC and UNK modes:
* MACHEP = 1.38777878078144567553E-17 2**-56
* MAXLOG = 8.8029691931113054295988E1 log(2**127)
* MINLOG = -8.872283911167299960540E1 log(2**-128)
* MAXNUM = 1.701411834604692317316873e38 2**127
*
* For IEEE arithmetic (IBMPC):
* MACHEP = 1.11022302462515654042E-16 2**-53
* MAXLOG = 7.09782712893383996843E2 log(2**1024)
* MINLOG = -7.08396418532264106224E2 log(2**-1022)
* MAXNUM = 1.7976931348623158E308 2**1024
*
* The global symbols for mathematical constants are
* PI = 3.14159265358979323846 pi
* PIO2 = 1.57079632679489661923 pi/2
* PIO4 = 7.85398163397448309616E-1 pi/4
* SQRT2 = 1.41421356237309504880 sqrt(2)
* SQRTH = 7.07106781186547524401E-1 sqrt(2)/2
* LOG2E = 1.4426950408889634073599 1/log(2)
* SQ2OPI = 7.9788456080286535587989E-1 sqrt( 2/pi )
* LOGE2 = 6.93147180559945309417E-1 log(2)
* LOGSQ2 = 3.46573590279972654709E-1 log(2)/2
* THPIO4 = 2.35619449019234492885 3*pi/4
* TWOOPI = 6.36619772367581343075535E-1 2/pi
*
* These lists are subject to change.
*/
/* const.c */
/*
Cephes Math Library Release 2.3: March, 1995
Copyright 1984, 1995 by Stephen L. Moshier
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. Neither the name of the <ORGANIZATION> nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
*/
#include "mconf.h"
#ifdef UNK
#if 1
double MACHEP = 1.11022302462515654042E-16; /* 2**-53 */
#else
double MACHEP = 1.38777878078144567553E-17; /* 2**-56 */
#endif
double UFLOWTHRESH = 2.22507385850720138309E-308; /* 2**-1022 */
#ifdef DENORMAL
double MAXLOG = 7.09782712893383996732E2; /* log(MAXNUM) */
/* double MINLOG = -7.44440071921381262314E2; */ /* log(2**-1074) */
double MINLOG = -7.451332191019412076235E2; /* log(2**-1075) */
#else
double MAXLOG = 7.08396418532264106224E2; /* log 2**1022 */
double MINLOG = -7.08396418532264106224E2; /* log 2**-1022 */
#endif
double MAXNUM = 1.79769313486231570815E308; /* 2**1024*(1-MACHEP) */
double PI = 3.14159265358979323846; /* pi */
double PIO2 = 1.57079632679489661923; /* pi/2 */
double PIO4 = 7.85398163397448309616E-1; /* pi/4 */
double SQRT2 = 1.41421356237309504880; /* sqrt(2) */
double SQRTH = 7.07106781186547524401E-1; /* sqrt(2)/2 */
double LOG2E = 1.4426950408889634073599; /* 1/log(2) */
double SQ2OPI = 7.9788456080286535587989E-1; /* sqrt( 2/pi ) */
double LOGE2 = 6.93147180559945309417E-1; /* log(2) */
double LOGSQ2 = 3.46573590279972654709E-1; /* log(2)/2 */
double THPIO4 = 2.35619449019234492885; /* 3*pi/4 */
double TWOOPI = 6.36619772367581343075535E-1; /* 2/pi */
#ifdef INFINITIES
double INFINITY = 1.0/0.0; /* 99e999; */
#else
double INFINITY = 1.79769313486231570815E308; /* 2**1024*(1-MACHEP) */
#endif
#ifdef NANS
double NAN = 1.0/0.0 - 1.0/0.0;
#else
double NAN = 0.0;
#endif
#ifdef MINUSZERO
double NEGZERO = -0.0;
#else
double NEGZERO = 0.0;
#endif
#endif
#ifdef IBMPC
/* 2**-53 = 1.11022302462515654042E-16 */
unsigned short MACHEP[4] = {0x0000,0x0000,0x0000,0x3ca0};
unsigned short UFLOWTHRESH[4] = {0x0000,0x0000,0x0000,0x0010};
#ifdef DENORMAL
/* log(MAXNUM) = 7.09782712893383996732224E2 */
unsigned short MAXLOG[4] = {0x39ef,0xfefa,0x2e42,0x4086};
/* log(2**-1074) = - -7.44440071921381262314E2 */
/*unsigned short MINLOG[4] = {0x71c3,0x446d,0x4385,0xc087};*/
unsigned short MINLOG[4] = {0x3052,0xd52d,0x4910,0xc087};
#else
/* log(2**1022) = 7.08396418532264106224E2 */
unsigned short MAXLOG[4] = {0xbcd2,0xdd7a,0x232b,0x4086};
/* log(2**-1022) = - 7.08396418532264106224E2 */
unsigned short MINLOG[4] = {0xbcd2,0xdd7a,0x232b,0xc086};
#endif
/* 2**1024*(1-MACHEP) = 1.7976931348623158E308 */
unsigned short MAXNUM[4] = {0xffff,0xffff,0xffff,0x7fef};
unsigned short PI[4] = {0x2d18,0x5444,0x21fb,0x4009};
unsigned short PIO2[4] = {0x2d18,0x5444,0x21fb,0x3ff9};
unsigned short PIO4[4] = {0x2d18,0x5444,0x21fb,0x3fe9};
unsigned short SQRT2[4] = {0x3bcd,0x667f,0xa09e,0x3ff6};
unsigned short SQRTH[4] = {0x3bcd,0x667f,0xa09e,0x3fe6};
unsigned short LOG2E[4] = {0x82fe,0x652b,0x1547,0x3ff7};
unsigned short SQ2OPI[4] = {0x3651,0x33d4,0x8845,0x3fe9};
unsigned short LOGE2[4] = {0x39ef,0xfefa,0x2e42,0x3fe6};
unsigned short LOGSQ2[4] = {0x39ef,0xfefa,0x2e42,0x3fd6};
unsigned short THPIO4[4] = {0x21d2,0x7f33,0xd97c,0x4002};
unsigned short TWOOPI[4] = {0xc883,0x6dc9,0x5f30,0x3fe4};
#ifdef INFINITIES
unsigned short INFINITY[4] = {0x0000,0x0000,0x0000,0x7ff0};
#else
unsigned short INFINITY[4] = {0xffff,0xffff,0xffff,0x7fef};
#endif
#ifdef NANS
unsigned short NAN[4] = {0x0000,0x0000,0x0000,0x7ffc};
#else
unsigned short NAN[4] = {0x0000,0x0000,0x0000,0x0000};
#endif
#ifdef MINUSZERO
unsigned short NEGZERO[4] = {0x0000,0x0000,0x0000,0x8000};
#else
unsigned short NEGZERO[4] = {0x0000,0x0000,0x0000,0x0000};
#endif
#endif
#ifdef MIEEE
/* 2**-53 = 1.11022302462515654042E-16 */
unsigned short MACHEP[4] = {0x3ca0,0x0000,0x0000,0x0000};
unsigned short UFLOWTHRESH[4] = {0x0010,0x0000,0x0000,0x0000};
#ifdef DENORMAL
/* log(2**1024) = 7.09782712893383996843E2 */
unsigned short MAXLOG[4] = {0x4086,0x2e42,0xfefa,0x39ef};
/* log(2**-1074) = - -7.44440071921381262314E2 */
/* unsigned short MINLOG[4] = {0xc087,0x4385,0x446d,0x71c3}; */
unsigned short MINLOG[4] = {0xc087,0x4910,0xd52d,0x3052};
#else
/* log(2**1022) = 7.08396418532264106224E2 */
unsigned short MAXLOG[4] = {0x4086,0x232b,0xdd7a,0xbcd2};
/* log(2**-1022) = - 7.08396418532264106224E2 */
unsigned short MINLOG[4] = {0xc086,0x232b,0xdd7a,0xbcd2};
#endif
/* 2**1024*(1-MACHEP) = 1.7976931348623158E308 */
unsigned short MAXNUM[4] = {0x7fef,0xffff,0xffff,0xffff};
unsigned short PI[4] = {0x4009,0x21fb,0x5444,0x2d18};
unsigned short PIO2[4] = {0x3ff9,0x21fb,0x5444,0x2d18};
unsigned short PIO4[4] = {0x3fe9,0x21fb,0x5444,0x2d18};
unsigned short SQRT2[4] = {0x3ff6,0xa09e,0x667f,0x3bcd};
unsigned short SQRTH[4] = {0x3fe6,0xa09e,0x667f,0x3bcd};
unsigned short LOG2E[4] = {0x3ff7,0x1547,0x652b,0x82fe};
unsigned short SQ2OPI[4] = {0x3fe9,0x8845,0x33d4,0x3651};
unsigned short LOGE2[4] = {0x3fe6,0x2e42,0xfefa,0x39ef};
unsigned short LOGSQ2[4] = {0x3fd6,0x2e42,0xfefa,0x39ef};
unsigned short THPIO4[4] = {0x4002,0xd97c,0x7f33,0x21d2};
unsigned short TWOOPI[4] = {0x3fe4,0x5f30,0x6dc9,0xc883};
#ifdef INFINITIES
unsigned short INFINITY[4] = {0x7ff0,0x0000,0x0000,0x0000};
#else
unsigned short INFINITY[4] = {0x7fef,0xffff,0xffff,0xffff};
#endif
#ifdef NANS
unsigned short NAN[4] = {0x7ff8,0x0000,0x0000,0x0000};
#else
unsigned short NAN[4] = {0x0000,0x0000,0x0000,0x0000};
#endif
#ifdef MINUSZERO
unsigned short NEGZERO[4] = {0x8000,0x0000,0x0000,0x0000};
#else
unsigned short NEGZERO[4] = {0x0000,0x0000,0x0000,0x0000};
#endif
#endif
#ifdef DEC
/* 2**-56 = 1.38777878078144567553E-17 */
unsigned short MACHEP[4] = {0022200,0000000,0000000,0000000};
unsigned short UFLOWTHRESH[4] = {0x0080,0x0000,0x0000,0x0000};
/* log 2**127 = 88.029691931113054295988 */
unsigned short MAXLOG[4] = {041660,007463,0143742,025733,};
/* log 2**-128 = -88.72283911167299960540 */
unsigned short MINLOG[4] = {0141661,071027,0173721,0147572,};
/* 2**127 = 1.701411834604692317316873e38 */
unsigned short MAXNUM[4] = {077777,0177777,0177777,0177777,};
unsigned short PI[4] = {040511,007732,0121041,064302,};
unsigned short PIO2[4] = {040311,007732,0121041,064302,};
unsigned short PIO4[4] = {040111,007732,0121041,064302,};
unsigned short SQRT2[4] = {040265,002363,031771,0157145,};
unsigned short SQRTH[4] = {040065,002363,031771,0157144,};
unsigned short LOG2E[4] = {040270,0125073,024534,013761,};
unsigned short SQ2OPI[4] = {040114,041051,0117241,0131204,};
unsigned short LOGE2[4] = {040061,071027,0173721,0147572,};
unsigned short LOGSQ2[4] = {037661,071027,0173721,0147572,};
unsigned short THPIO4[4] = {040426,0145743,0174631,007222,};
unsigned short TWOOPI[4] = {040042,0174603,067116,042025,};
/* Approximate infinity by MAXNUM. */
unsigned short INFINITY[4] = {077777,0177777,0177777,0177777,};
unsigned short NAN[4] = {0000000,0000000,0000000,0000000};
#ifdef MINUSZERO
unsigned short NEGZERO[4] = {0000000,0000000,0000000,0100000};
#else
unsigned short NEGZERO[4] = {0000000,0000000,0000000,0000000};
#endif
#endif
#ifndef UNK
extern unsigned short MACHEP[];
extern unsigned short UFLOWTHRESH[];
extern unsigned short MAXLOG[];
extern unsigned short UNDLOG[];
extern unsigned short MINLOG[];
extern unsigned short MAXNUM[];
extern unsigned short PI[];
extern unsigned short PIO2[];
extern unsigned short PIO4[];
extern unsigned short SQRT2[];
extern unsigned short SQRTH[];
extern unsigned short LOG2E[];
extern unsigned short SQ2OPI[];
extern unsigned short LOGE2[];
extern unsigned short LOGSQ2[];
extern unsigned short THPIO4[];
extern unsigned short TWOOPI[];
extern unsigned short INFINITY[];
extern unsigned short NAN[];
extern unsigned short NEGZERO[];
#endif

107
source/common/thirdparty/math/cosh.c vendored Normal file
View file

@ -0,0 +1,107 @@
/* cosh.c
*
* Hyperbolic cosine
*
*
*
* SYNOPSIS:
*
* double x, y, cosh();
*
* y = cosh( x );
*
*
*
* DESCRIPTION:
*
* Returns hyperbolic cosine of argument in the range MINLOG to
* MAXLOG.
*
* cosh(x) = ( exp(x) + exp(-x) )/2.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC +- 88 50000 4.0e-17 7.7e-18
* IEEE +-MAXLOG 30000 2.6e-16 5.7e-17
*
*
* ERROR MESSAGES:
*
* message condition value returned
* cosh overflow |x| > MAXLOG MAXNUM
*
*
*/
/* cosh.c */
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1985, 1995, 2000 by Stephen L. Moshier
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. Neither the name of the <ORGANIZATION> nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
*/
#include "mconf.h"
#ifdef ANSIPROT
extern double c_exp ( double );
extern int isnan ( double );
extern int isfinite ( double );
#else
double c_exp();
int isnan(), isfinite();
#endif
extern double MAXLOG, INFINITY, LOGE2;
double c_cosh(x)
double x;
{
double y;
#ifdef NANS
if( isnan(x) )
return(x);
#endif
if( x < 0 )
x = -x;
if( x > (MAXLOG + LOGE2) )
{
mtherr( "cosh", OVERFLOW );
return( INFINITY );
}
if( x >= (MAXLOG - LOGE2) )
{
y = c_exp(0.5 * x);
y = (0.5 * y) * y;
return(y);
}
y = c_exp(x);
y = 0.5 * (y + 1.0 / y);
return( y );
}

207
source/common/thirdparty/math/exp.c vendored Normal file
View file

@ -0,0 +1,207 @@
/* exp.c
*
* Exponential function
*
*
*
* SYNOPSIS:
*
* double x, y, exp();
*
* y = exp( x );
*
*
*
* DESCRIPTION:
*
* Returns e (2.71828...) raised to the x power.
*
* Range reduction is accomplished by separating the argument
* into an integer k and fraction f such that
*
* x k f
* e = 2 e.
*
* A Pade' form 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
* of degree 2/3 is used to approximate exp(f) in the basic
* interval [-0.5, 0.5].
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC 0, MAXLOG 38000 3.0e-17 6.2e-18
* IEEE +- 708 40000 2.0e-16 5.6e-17
*
*
* Error amplification in the exponential function can be
* a serious matter. The error propagation involves
* exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ),
* which shows that a 1 lsb error in representing X produces
* a relative error of X times 1 lsb in the function.
* While the routine gives an accurate result for arguments
* that are exactly represented by a double precision
* computer number, the result contains amplified roundoff
* error for large arguments not exactly represented.
*
*
* ERROR MESSAGES:
*
* message condition value returned
* exp underflow x < MINLOG 0.0
* exp overflow x > MAXLOG MAXNUM
*
*/
/*
Cephes Math Library Release 2.2: January, 1991
Copyright 1984, 1991 by Stephen L. Moshier
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. Neither the name of the <ORGANIZATION> nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
/* Exponential function */
#include "mconf.h"
static char fname[] = {"exp"};
#ifdef UNK
static double P[] = {
1.26177193074810590878E-4,
3.02994407707441961300E-2,
9.99999999999999999910E-1,
};
static double Q[] = {
3.00198505138664455042E-6,
2.52448340349684104192E-3,
2.27265548208155028766E-1,
2.00000000000000000009E0,
};
static double C1 = 6.93145751953125E-1;
static double C2 = 1.42860682030941723212E-6;
#endif
#ifdef DEC
static short P[] = {
0035004,0047156,0127442,0057502,
0036770,0033210,0063121,0061764,
0040200,0000000,0000000,0000000,
};
static short Q[] = {
0033511,0072665,0160662,0176377,
0036045,0070715,0124105,0132777,
0037550,0134114,0142077,0001637,
0040400,0000000,0000000,0000000,
};
static short sc1[] = {0040061,0071000,0000000,0000000};
#define C1 (*(double *)sc1)
static short sc2[] = {0033277,0137216,0075715,0057117};
#define C2 (*(double *)sc2)
#endif
#ifdef IBMPC
static short P[] = {
0x4be8,0xd5e4,0x89cd,0x3f20,
0x2c7e,0x0cca,0x06d1,0x3f9f,
0x0000,0x0000,0x0000,0x3ff0,
};
static short Q[] = {
0x5fa0,0xbc36,0x2eb6,0x3ec9,
0xb6c0,0xb508,0xae39,0x3f64,
0xe074,0x9887,0x1709,0x3fcd,
0x0000,0x0000,0x0000,0x4000,
};
static short sc1[] = {0x0000,0x0000,0x2e40,0x3fe6};
#define C1 (*(double *)sc1)
static short sc2[] = {0xabca,0xcf79,0xf7d1,0x3eb7};
#define C2 (*(double *)sc2)
#endif
#ifdef MIEEE
static short P[] = {
0x3f20,0x89cd,0xd5e4,0x4be8,
0x3f9f,0x06d1,0x0cca,0x2c7e,
0x3ff0,0x0000,0x0000,0x0000,
};
static short Q[] = {
0x3ec9,0x2eb6,0xbc36,0x5fa0,
0x3f64,0xae39,0xb508,0xb6c0,
0x3fcd,0x1709,0x9887,0xe074,
0x4000,0x0000,0x0000,0x0000,
};
static short sc1[] = {0x3fe6,0x2e40,0x0000,0x0000};
#define C1 (*(double *)sc1)
static short sc2[] = {0x3eb7,0xf7d1,0xcf79,0xabca};
#define C2 (*(double *)sc2)
#endif
extern double LOGE2, LOG2E, MAXLOG, MINLOG, MAXNUM;
double c_exp(x)
double x;
{
double px, xx;
int n;
double polevl(), floor(), ldexp();
if( x > MAXLOG)
{
mtherr( fname, OVERFLOW );
return( MAXNUM );
}
if( x < MINLOG )
{
mtherr( fname, UNDERFLOW );
return(0.0);
}
/* Express e**x = e**g 2**n
* = e**g e**( n loge(2) )
* = e**( g + n loge(2) )
*/
px = floor( LOG2E * x + 0.5 ); /* floor() truncates toward -infinity. */
n = (int)px;
x -= px * C1;
x -= px * C2;
/* rational approximation for exponential
* of the fractional part:
* e**x = 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
*/
xx = x * x;
px = x * polevl( xx, P, 2 );
x = px/( polevl( xx, Q, 3 ) - px );
x = 1.0 + ldexp( x, 1 );
/* multiply by power of 2 */
x = ldexp( x, n );
return(x);
}

View file

@ -0,0 +1,106 @@
/*
** fastsin.cpp
** a table/linear interpolation-based sine function that is both
** precise and fast enough for most purposes.
**
**---------------------------------------------------------------------------
** Copyright 2015 Christoph Oelckers
** All rights reserved.
**
** Redistribution and use in source and binary forms, with or without
** modification, are permitted provided that the following conditions
** are met:
**
** 1. Redistributions of source code must retain the above copyright
** notice, this list of conditions and the following disclaimer.
** 2. Redistributions in binary form must reproduce the above copyright
** notice, this list of conditions and the following disclaimer in the
** documentation and/or other materials provided with the distribution.
** 3. The name of the author may not be used to endorse or promote products
** derived from this software without specific prior written permission.
**
** THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
** IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
** OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
** IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
** INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
** NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
** DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
** THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
** (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
** THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
**---------------------------------------------------------------------------
**
*/
#include <math.h>
#include "cmath.h"
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
FFastTrig fasttrig;
FFastTrig::FFastTrig()
{
const double pimul = M_PI * 2 / TBLPERIOD;
for (int i = 0; i < 2049; i++)
{
sinetable[i] = (float)c_sin(i*pimul);
}
}
__forceinline double FFastTrig::sinq1(unsigned bangle)
{
unsigned int index = bangle >> BITSHIFT;
if ((bangle &= (REMAINDER)) == 0) // This is to avoid precision problems at 180°
{
return double(sinetable[index]);
}
else
{
return (double(sinetable[index]) * (REMAINDER - bangle) + double(sinetable[index + 1]) * bangle) * (1. / REMAINDER);
}
}
double FFastTrig::sin(unsigned bangle)
{
switch (bangle & 0xc0000000)
{
default:
return sinq1(bangle);
case 0x40000000:
return sinq1(0x80000000 - bangle);
case 0x80000000:
return -sinq1(bangle - 0x80000000);
case 0xc0000000:
return -sinq1(0 - bangle);
}
}
double FFastTrig::cos(unsigned bangle)
{
switch (bangle & 0xc0000000)
{
default:
return sinq1(0x40000000 - bangle);
case 0x40000000:
return -sinq1(bangle - 0x40000000);
case 0x80000000:
return -sinq1(0xc0000000 - bangle);
case 0xc0000000:
return sinq1(bangle - 0xc0000000);
}
}

261
source/common/thirdparty/math/isnan.c vendored Normal file
View file

@ -0,0 +1,261 @@
/* isnan()
* signbit()
* isfinite()
*
* Floating point numeric utilities
*
*
*
* SYNOPSIS:
*
* double ceil(), floor(), frexp(), ldexp();
* int signbit(), isnan(), isfinite();
* double x, y;
* int expnt, n;
*
* y = floor(x);
* y = ceil(x);
* y = frexp( x, &expnt );
* y = ldexp( x, n );
* n = signbit(x);
* n = isnan(x);
* n = isfinite(x);
*
*
*
* DESCRIPTION:
*
* All four routines return a double precision floating point
* result.
*
* floor() returns the largest integer less than or equal to x.
* It truncates toward minus infinity.
*
* ceil() returns the smallest integer greater than or equal
* to x. It truncates toward plus infinity.
*
* frexp() extracts the exponent from x. It returns an integer
* power of two to expnt and the significand between 0.5 and 1
* to y. Thus x = y * 2**expn.
*
* ldexp() multiplies x by 2**n.
*
* signbit(x) returns 1 if the sign bit of x is 1, else 0.
*
* These functions are part of the standard C run time library
* for many but not all C compilers. The ones supplied are
* written in C for either DEC or IEEE arithmetic. They should
* be used only if your compiler library does not already have
* them.
*
* The IEEE versions assume that denormal numbers are implemented
* in the arithmetic. Some modifications will be required if
* the arithmetic has abrupt rather than gradual underflow.
*/
/*
Cephes Math Library Release 2.3: March, 1995
Copyright 1984, 1995 by Stephen L. Moshier
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. Neither the name of the <ORGANIZATION> nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
*/
#include "mconf.h"
#ifdef UNK
/* ceil(), floor(), frexp(), ldexp() may need to be rewritten. */
#undef UNK
#if BIGENDIAN
#define MIEEE 1
#else
#define IBMPC 1
#endif
#endif
/* Return 1 if the sign bit of x is 1, else 0. */
int signbit(x)
double x;
{
union
{
double d;
short s[4];
int i[2];
} u;
u.d = x;
if( sizeof(int) == 4 )
{
#ifdef IBMPC
return( u.i[1] < 0 );
#endif
#ifdef DEC
return( u.s[3] < 0 );
#endif
#ifdef MIEEE
return( u.i[0] < 0 );
#endif
}
else
{
#ifdef IBMPC
return( u.s[3] < 0 );
#endif
#ifdef DEC
return( u.s[3] < 0 );
#endif
#ifdef MIEEE
return( u.s[0] < 0 );
#endif
}
}
/* Return 1 if x is a number that is Not a Number, else return 0. */
int isnan(x)
double x;
{
#ifdef NANS
union
{
double d;
unsigned short s[4];
unsigned int i[2];
} u;
u.d = x;
if( sizeof(int) == 4 )
{
#ifdef IBMPC
if( ((u.i[1] & 0x7ff00000) == 0x7ff00000)
&& (((u.i[1] & 0x000fffff) != 0) || (u.i[0] != 0)))
return 1;
#endif
#ifdef DEC
if( (u.s[1] & 0x7fff) == 0)
{
if( (u.s[2] | u.s[1] | u.s[0]) != 0 )
return(1);
}
#endif
#ifdef MIEEE
if( ((u.i[0] & 0x7ff00000) == 0x7ff00000)
&& (((u.i[0] & 0x000fffff) != 0) || (u.i[1] != 0)))
return 1;
#endif
return(0);
}
else
{ /* size int not 4 */
#ifdef IBMPC
if( (u.s[3] & 0x7ff0) == 0x7ff0)
{
if( ((u.s[3] & 0x000f) | u.s[2] | u.s[1] | u.s[0]) != 0 )
return(1);
}
#endif
#ifdef DEC
if( (u.s[3] & 0x7fff) == 0)
{
if( (u.s[2] | u.s[1] | u.s[0]) != 0 )
return(1);
}
#endif
#ifdef MIEEE
if( (u.s[0] & 0x7ff0) == 0x7ff0)
{
if( ((u.s[0] & 0x000f) | u.s[1] | u.s[2] | u.s[3]) != 0 )
return(1);
}
#endif
return(0);
} /* size int not 4 */
#else
/* No NANS. */
return(0);
#endif
}
/* Return 1 if x is not infinite and is not a NaN. */
int isfinite(x)
double x;
{
#ifdef INFINITIES
union
{
double d;
unsigned short s[4];
unsigned int i[2];
} u;
u.d = x;
if( sizeof(int) == 4 )
{
#ifdef IBMPC
if( (u.i[1] & 0x7ff00000) != 0x7ff00000)
return 1;
#endif
#ifdef DEC
if( (u.s[3] & 0x7fff) != 0)
return 1;
#endif
#ifdef MIEEE
if( (u.i[0] & 0x7ff00000) != 0x7ff00000)
return 1;
#endif
return(0);
}
else
{
#ifdef IBMPC
if( (u.s[3] & 0x7ff0) != 0x7ff0)
return 1;
#endif
#ifdef DEC
if( (u.s[3] & 0x7fff) != 0)
return 1;
#endif
#ifdef MIEEE
if( (u.s[0] & 0x7ff0) != 0x7ff0)
return 1;
#endif
return(0);
}
#else
/* No INFINITY. */
return(1);
#endif
}

365
source/common/thirdparty/math/log.c vendored Normal file
View file

@ -0,0 +1,365 @@
/* log.c
*
* Natural logarithm
*
*
*
* SYNOPSIS:
*
* double x, y, log();
*
* y = log( x );
*
*
*
* DESCRIPTION:
*
* Returns the base e (2.718...) logarithm of x.
*
* The argument is separated into its exponent and fractional
* parts. If the exponent is between -1 and +1, the logarithm
* of the fraction is approximated by
*
* log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
*
* Otherwise, setting z = 2(x-1)/x+1),
*
* log(x) = z + z**3 P(z)/Q(z).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0.5, 2.0 150000 1.44e-16 5.06e-17
* IEEE +-MAXNUM 30000 1.20e-16 4.78e-17
* DEC 0, 10 170000 1.8e-17 6.3e-18
*
* In the tests over the interval [+-MAXNUM], the logarithms
* of the random arguments were uniformly distributed over
* [0, MAXLOG].
*
* ERROR MESSAGES:
*
* log singularity: x = 0; returns -INFINITY
* log domain: x < 0; returns NAN
*/
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. Neither the name of the <ORGANIZATION> nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
*/
#include "mconf.h"
static char fname[] = {"log"};
/* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
* 1/sqrt(2) <= x < sqrt(2)
*/
#ifdef UNK
static double P[] = {
1.01875663804580931796E-4,
4.97494994976747001425E-1,
4.70579119878881725854E0,
1.44989225341610930846E1,
1.79368678507819816313E1,
7.70838733755885391666E0,
};
static double Q[] = {
/* 1.00000000000000000000E0, */
1.12873587189167450590E1,
4.52279145837532221105E1,
8.29875266912776603211E1,
7.11544750618563894466E1,
2.31251620126765340583E1,
};
#endif
#ifdef DEC
static unsigned short P[] = {
0037777,0127270,0162547,0057274,
0041001,0054665,0164317,0005341,
0041451,0034104,0031640,0105773,
0041677,0011276,0123617,0160135,
0041701,0126603,0053215,0117250,
0041420,0115777,0135206,0030232,
};
static unsigned short Q[] = {
/*0040200,0000000,0000000,0000000,*/
0041220,0144332,0045272,0174241,
0041742,0164566,0035720,0130431,
0042246,0126327,0166065,0116357,
0042372,0033420,0157525,0124560,
0042271,0167002,0066537,0172303,
0041730,0164777,0113711,0044407,
};
#endif
#ifdef IBMPC
static unsigned short P[] = {
0x1bb0,0x93c3,0xb4c2,0x3f1a,
0x52f2,0x3f56,0xd6f5,0x3fdf,
0x6911,0xed92,0xd2ba,0x4012,
0xeb2e,0xc63e,0xff72,0x402c,
0xc84d,0x924b,0xefd6,0x4031,
0xdcf8,0x7d7e,0xd563,0x401e,
};
static unsigned short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0xef8e,0xae97,0x9320,0x4026,
0xc033,0x4e19,0x9d2c,0x4046,
0xbdbd,0xa326,0xbf33,0x4054,
0xae21,0xeb5e,0xc9e2,0x4051,
0x25b2,0x9e1f,0x200a,0x4037,
};
#endif
#ifdef MIEEE
static unsigned short P[] = {
0x3f1a,0xb4c2,0x93c3,0x1bb0,
0x3fdf,0xd6f5,0x3f56,0x52f2,
0x4012,0xd2ba,0xed92,0x6911,
0x402c,0xff72,0xc63e,0xeb2e,
0x4031,0xefd6,0x924b,0xc84d,
0x401e,0xd563,0x7d7e,0xdcf8,
};
static unsigned short Q[] = {
/*0x3ff0,0x0000,0x0000,0x0000,*/
0x4026,0x9320,0xae97,0xef8e,
0x4046,0x9d2c,0x4e19,0xc033,
0x4054,0xbf33,0xa326,0xbdbd,
0x4051,0xc9e2,0xeb5e,0xae21,
0x4037,0x200a,0x9e1f,0x25b2,
};
#endif
/* Coefficients for log(x) = z + z**3 P(z)/Q(z),
* where z = 2(x-1)/(x+1)
* 1/sqrt(2) <= x < sqrt(2)
*/
#ifdef UNK
static double R[3] = {
-7.89580278884799154124E-1,
1.63866645699558079767E1,
-6.41409952958715622951E1,
};
static double S[3] = {
/* 1.00000000000000000000E0,*/
-3.56722798256324312549E1,
3.12093766372244180303E2,
-7.69691943550460008604E2,
};
#endif
#ifdef DEC
static unsigned short R[12] = {
0140112,0020756,0161540,0072035,
0041203,0013743,0114023,0155527,
0141600,0044060,0104421,0050400,
};
static unsigned short S[12] = {
/*0040200,0000000,0000000,0000000,*/
0141416,0130152,0017543,0064122,
0042234,0006000,0104527,0020155,
0142500,0066110,0146631,0174731,
};
#endif
#ifdef IBMPC
static unsigned short R[12] = {
0x0e84,0xdc6c,0x443d,0xbfe9,
0x7b6b,0x7302,0x62fc,0x4030,
0x2a20,0x1122,0x0906,0xc050,
};
static unsigned short S[12] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0x6d0a,0x43ec,0xd60d,0xc041,
0xe40e,0x112a,0x8180,0x4073,
0x3f3b,0x19b3,0x0d89,0xc088,
};
#endif
#ifdef MIEEE
static unsigned short R[12] = {
0xbfe9,0x443d,0xdc6c,0x0e84,
0x4030,0x62fc,0x7302,0x7b6b,
0xc050,0x0906,0x1122,0x2a20,
};
static unsigned short S[12] = {
/*0x3ff0,0x0000,0x0000,0x0000,*/
0xc041,0xd60d,0x43ec,0x6d0a,
0x4073,0x8180,0x112a,0xe40e,
0xc088,0x0d89,0x19b3,0x3f3b,
};
#endif
#ifdef ANSIPROT
extern double frexp ( double, int * );
extern double ldexp ( double, int );
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern int isnan ( double );
extern int isfinite ( double );
#else
double frexp(), ldexp(), polevl(), p1evl();
int isnan(), isfinite();
#endif
#define SQRTH 0.70710678118654752440
extern double INFINITY, NAN;
double c_log(x)
double x;
{
int e;
#ifdef DEC
short *q;
#endif
double y, z;
#ifdef NANS
if( isnan(x) )
return(x);
#endif
#ifdef INFINITIES
if( x == INFINITY )
return(x);
#endif
/* Test for domain */
if( x <= 0.0 )
{
if( x == 0.0 )
{
mtherr( fname, SING );
return( -INFINITY );
}
else
{
mtherr( fname, DOMAIN );
return( NAN );
}
}
/* separate mantissa from exponent */
#ifdef DEC
q = (short *)&x;
e = *q; /* short containing exponent */
e = ((e >> 7) & 0377) - 0200; /* the exponent */
*q &= 0177; /* strip exponent from x */
*q |= 040000; /* x now between 0.5 and 1 */
#endif
/* Note, frexp is used so that denormal numbers
* will be handled properly.
*/
#ifdef IBMPC
x = frexp( x, &e );
/*
q = (short *)&x;
q += 3;
e = *q;
e = ((e >> 4) & 0x0fff) - 0x3fe;
*q &= 0x0f;
*q |= 0x3fe0;
*/
#endif
/* Equivalent C language standard library function: */
#ifdef UNK
x = frexp( x, &e );
#endif
#ifdef MIEEE
x = frexp( x, &e );
#endif
/* logarithm using log(x) = z + z**3 P(z)/Q(z),
* where z = 2(x-1)/x+1)
*/
if( (e > 2) || (e < -2) )
{
if( x < SQRTH )
{ /* 2( 2x-1 )/( 2x+1 ) */
e -= 1;
z = x - 0.5;
y = 0.5 * z + 0.5;
}
else
{ /* 2 (x-1)/(x+1) */
z = x - 0.5;
z -= 0.5;
y = 0.5 * x + 0.5;
}
x = z / y;
/* rational form */
z = x*x;
z = x * ( z * polevl( z, R, 2 ) / p1evl( z, S, 3 ) );
y = e;
z = z - y * 2.121944400546905827679e-4;
z = z + x;
z = z + e * 0.693359375;
goto ldone;
}
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
if( x < SQRTH )
{
e -= 1;
x = ldexp( x, 1 ) - 1.0; /* 2x - 1 */
}
else
{
x = x - 1.0;
}
/* rational form */
z = x*x;
#if DEC
y = x * ( z * polevl( x, P, 5 ) / p1evl( x, Q, 6 ) );
#else
y = x * ( z * polevl( x, P, 5 ) / p1evl( x, Q, 5 ) );
#endif
if( e )
y = y - e * 2.121944400546905827679e-4;
y = y - ldexp( z, -1 ); /* y - 0.5 * z */
z = x + y;
if( e )
z = z + e * 0.693359375;
ldone:
return( z );
}

274
source/common/thirdparty/math/log10.c vendored Normal file
View file

@ -0,0 +1,274 @@
/* log10.c
*
* Common logarithm
*
*
*
* SYNOPSIS:
*
* double x, y, log10();
*
* y = log10( x );
*
*
*
* DESCRIPTION:
*
* Returns logarithm to the base 10 of x.
*
* The argument is separated into its exponent and fractional
* parts. The logarithm of the fraction is approximated by
*
* log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0.5, 2.0 30000 1.5e-16 5.0e-17
* IEEE 0, MAXNUM 30000 1.4e-16 4.8e-17
* DEC 1, MAXNUM 50000 2.5e-17 6.0e-18
*
* In the tests over the interval [1, MAXNUM], the logarithms
* of the random arguments were uniformly distributed over
* [0, MAXLOG].
*
* ERROR MESSAGES:
*
* log10 singularity: x = 0; returns -INFINITY
* log10 domain: x < 0; returns NAN
*/
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. Neither the name of the <ORGANIZATION> nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
*/
#include "mconf.h"
static char fname[] = {"log10"};
/* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
* 1/sqrt(2) <= x < sqrt(2)
*/
#ifdef UNK
static double P[] = {
4.58482948458143443514E-5,
4.98531067254050724270E-1,
6.56312093769992875930E0,
2.97877425097986925891E1,
6.06127134467767258030E1,
5.67349287391754285487E1,
1.98892446572874072159E1
};
static double Q[] = {
/* 1.00000000000000000000E0, */
1.50314182634250003249E1,
8.27410449222435217021E1,
2.20664384982121929218E2,
3.07254189979530058263E2,
2.14955586696422947765E2,
5.96677339718622216300E1
};
#endif
#ifdef DEC
static unsigned short P[] = {
0034500,0046473,0051374,0135174,
0037777,0037566,0145712,0150321,
0040722,0002426,0031543,0123107,
0041356,0046513,0170752,0004346,
0041562,0071553,0023536,0163343,
0041542,0170221,0024316,0114216,
0041237,0016454,0046611,0104602
};
static unsigned short Q[] = {
/*0040200,0000000,0000000,0000000,*/
0041160,0100260,0067736,0102424,
0041645,0075552,0036563,0147072,
0042134,0125025,0021132,0025320,
0042231,0120211,0046030,0103271,
0042126,0172241,0052151,0120426,
0041556,0125702,0072116,0047103
};
#endif
#ifdef IBMPC
static unsigned short P[] = {
0x974f,0x6a5f,0x09a7,0x3f08,
0x5a1a,0xd979,0xe7ee,0x3fdf,
0x74c9,0xc66c,0x40a2,0x401a,
0x411d,0x7e3d,0xc9a9,0x403d,
0xdcdc,0x64eb,0x4e6d,0x404e,
0xd312,0x2519,0x5e12,0x404c,
0x3130,0x89b1,0xe3a5,0x4033
};
static unsigned short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0xd0a2,0x0dfb,0x1016,0x402e,
0x79c7,0x47ae,0xaf6d,0x4054,
0x455a,0xa44b,0x9542,0x406b,
0x10d7,0x2983,0x3411,0x4073,
0x3423,0x2a8d,0xde94,0x406a,
0xc9c8,0x4e89,0xd578,0x404d
};
#endif
#ifdef MIEEE
static unsigned short P[] = {
0x3f08,0x09a7,0x6a5f,0x974f,
0x3fdf,0xe7ee,0xd979,0x5a1a,
0x401a,0x40a2,0xc66c,0x74c9,
0x403d,0xc9a9,0x7e3d,0x411d,
0x404e,0x4e6d,0x64eb,0xdcdc,
0x404c,0x5e12,0x2519,0xd312,
0x4033,0xe3a5,0x89b1,0x3130
};
static unsigned short Q[] = {
0x402e,0x1016,0x0dfb,0xd0a2,
0x4054,0xaf6d,0x47ae,0x79c7,
0x406b,0x9542,0xa44b,0x455a,
0x4073,0x3411,0x2983,0x10d7,
0x406a,0xde94,0x2a8d,0x3423,
0x404d,0xd578,0x4e89,0xc9c8
};
#endif
#define SQRTH 0.70710678118654752440
#define L102A 3.0078125E-1
#define L102B 2.48745663981195213739E-4
#define L10EA 4.3359375E-1
#define L10EB 7.00731903251827651129E-4
#ifdef ANSIPROT
extern double frexp ( double, int * );
extern double ldexp ( double, int );
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern int isnan ( double );
extern int isfinite ( double );
#else
double frexp(), ldexp(), polevl(), p1evl();
int isnan(), isfinite();
#endif
extern double LOGE2, SQRT2, INFINITY, NAN;
double c_log10(x)
double x;
{
VOLATILE double z;
double y;
#ifdef DEC
short *q;
#endif
int e;
#ifdef NANS
if( isnan(x) )
return(x);
#endif
#ifdef INFINITIES
if( x == INFINITY )
return(x);
#endif
/* Test for domain */
if( x <= 0.0 )
{
if( x == 0.0 )
{
mtherr( fname, SING );
return( -INFINITY );
}
else
{
mtherr( fname, DOMAIN );
return( NAN );
}
}
/* separate mantissa from exponent */
#ifdef DEC
q = (short *)&x;
e = *q; /* short containing exponent */
e = ((e >> 7) & 0377) - 0200; /* the exponent */
*q &= 0177; /* strip exponent from x */
*q |= 040000; /* x now between 0.5 and 1 */
#endif
#ifdef IBMPC
x = frexp( x, &e );
/*
q = (short *)&x;
q += 3;
e = *q;
e = ((e >> 4) & 0x0fff) - 0x3fe;
*q &= 0x0f;
*q |= 0x3fe0;
*/
#endif
/* Equivalent C language standard library function: */
#ifdef UNK
x = frexp( x, &e );
#endif
#ifdef MIEEE
x = frexp( x, &e );
#endif
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
if( x < SQRTH )
{
e -= 1;
x = ldexp( x, 1 ) - 1.0; /* 2x - 1 */
}
else
{
x = x - 1.0;
}
/* rational form */
z = x*x;
y = x * ( z * polevl( x, P, 6 ) / p1evl( x, Q, 6 ) );
y = y - ldexp( z, -1 ); /* y - 0.5 * x**2 */
/* multiply log of fraction by log10(e)
* and base 2 exponent by log10(2)
*/
z = (x + y) * L10EB; /* accumulate terms in order of size */
z += y * L10EA;
z += x * L10EA;
z += e * L102B;
z += e * L102A;
return( z );
}

223
source/common/thirdparty/math/mconf.h vendored Normal file
View file

@ -0,0 +1,223 @@
/* mconf.h
*
* Common include file for math routines
*
*
*
* SYNOPSIS:
*
* #include "mconf.h"
*
*
*
* DESCRIPTION:
*
* This file contains definitions for error codes that are
* passed to the common error handling routine mtherr()
* (which see).
*
* The file also includes a conditional assembly definition
* for the type of computer arithmetic (IEEE, DEC, Motorola
* IEEE, or UNKnown).
*
* For Digital Equipment PDP-11 and VAX computers, certain
* IBM systems, and others that use numbers with a 56-bit
* significand, the symbol DEC should be defined. In this
* mode, most floating point constants are given as arrays
* of octal integers to eliminate decimal to binary conversion
* errors that might be introduced by the compiler.
*
* For little-endian computers, such as IBM PC, that follow the
* IEEE Standard for Binary Floating Point Arithmetic (ANSI/IEEE
* Std 754-1985), the symbol IBMPC should be defined. These
* numbers have 53-bit significands. In this mode, constants
* are provided as arrays of hexadecimal 16 bit integers.
*
* Big-endian IEEE format is denoted MIEEE. On some RISC
* systems such as Sun SPARC, double precision constants
* must be stored on 8-byte address boundaries. Since integer
* arrays may be aligned differently, the MIEEE configuration
* may fail on such machines.
*
* To accommodate other types of computer arithmetic, all
* constants are also provided in a normal decimal radix
* which one can hope are correctly converted to a suitable
* format by the available C language compiler. To invoke
* this mode, define the symbol UNK.
*
* An important difference among these modes is a predefined
* set of machine arithmetic constants for each. The numbers
* MACHEP (the machine roundoff error), MAXNUM (largest number
* represented), and several other parameters are preset by
* the configuration symbol. Check the file const.c to
* ensure that these values are correct for your computer.
*
* Configurations NANS, INFINITIES, MINUSZERO, and DENORMAL
* may fail on many systems. Verify that they are supposed
* to work on your computer.
*/
/*
Cephes Math Library Release 2.3: June, 1995
Copyright 1984, 1987, 1989, 1995 by Stephen L. Moshier
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. Neither the name of the <ORGANIZATION> nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
*/
/* Define if the `long double' type works. */
//#define HAVE_LONG_DOUBLE 0
/* Define as the return type of signal handlers (int or void). */
#define RETSIGTYPE void
/* Define if you have the ANSI C header files. */
#define STDC_HEADERS 1
/* Define if your processor stores words with the most significant
byte first (like Motorola and SPARC, unlike Intel and VAX). */
/* #undef WORDS_BIGENDIAN */
/* Define if floating point words are bigendian. */
/* #undef FLOAT_WORDS_BIGENDIAN */
/* The number of bytes in a int. */
#define SIZEOF_INT 4
/* Define if you have the <string.h> header file. */
#define HAVE_STRING_H 1
/* Name of package */
#define PACKAGE "cephes"
/* Version number of package */
#define VERSION "2.7"
/* Constant definitions for math error conditions
*/
#define DOMAIN 1 /* argument domain error */
#define SING 2 /* argument singularity */
#define OVERFLOW 3 /* overflow range error */
#define UNDERFLOW 4 /* underflow range error */
#define TLOSS 5 /* total loss of precision */
#define PLOSS 6 /* partial loss of precision */
#define EDOM 33
#define ERANGE 34
/* Complex numeral. */
typedef struct
{
double r;
double i;
} cmplx;
#ifdef HAVE_LONG_DOUBLE
/* Long double complex numeral. */
typedef struct
{
long double r;
long double i;
} cmplxl;
#endif
/* Type of computer arithmetic */
/* PDP-11, Pro350, VAX:
*/
/* #define DEC 1 */
/* Intel IEEE, low order words come first:
*/
//#define IBMPC 1
/* Motorola IEEE, high order words come first
* (Sun 680x0 workstation):
*/
/* #define MIEEE 1 */
/* UNKnown arithmetic, invokes coefficients given in
* normal decimal format. Beware of range boundary
* problems (MACHEP, MAXLOG, etc. in const.c) and
* roundoff problems in pow.c:
* (Sun SPARCstation)
*/
#define UNK 1
/* If you define UNK, then be sure to set BIGENDIAN properly. */
#ifdef FLOAT_WORDS_BIGENDIAN
#define BIGENDIAN 1
#else
#define BIGENDIAN 0
#endif
/* Define this `volatile' if your compiler thinks
* that floating point arithmetic obeys the associative
* and distributive laws. It will defeat some optimizations
* (but probably not enough of them).
*
* #define VOLATILE volatile
*/
#define VOLATILE
/* For 12-byte long doubles on an i386, pad a 16-bit short 0
* to the end of real constants initialized by integer arrays.
*
* #define XPD 0,
*
* Otherwise, the type is 10 bytes long and XPD should be
* defined blank (e.g., Microsoft C).
*
* #define XPD
*/
#define XPD 0,
/* Define to support tiny denormal numbers, else undefine. */
#define DENORMAL 1
/* Define to ask for infinity support, else undefine. */
#define INFINITIES 1
/* Define to ask for support of numbers that are Not-a-Number,
else undefine. This may automatically define INFINITIES in some files. */
#define NANS 1
/* Define to distinguish between -0.0 and +0.0. */
#define MINUSZERO 1
/* Define 1 for ANSI C atan2() function
See atan.c and clog.c. */
#define ANSIC 1
/* Get ANSI function prototypes, if you want them. */
#if 1
/* #ifdef __STDC__ */
#define ANSIPROT 1
int mtherr ( char *, int );
#else
int mtherr();
#endif
/* Variable for error reporting. See mtherr.c. */
extern int merror;

127
source/common/thirdparty/math/mtherr.c vendored Normal file
View file

@ -0,0 +1,127 @@
/* mtherr.c
*
* Library common error handling routine
*
*
*
* SYNOPSIS:
*
* char *fctnam;
* int code;
* int mtherr();
*
* mtherr( fctnam, code );
*
*
*
* DESCRIPTION:
*
* This routine may be called to report one of the following
* error conditions (in the include file mconf.h).
*
* Mnemonic Value Significance
*
* DOMAIN 1 argument domain error
* SING 2 function singularity
* OVERFLOW 3 overflow range error
* UNDERFLOW 4 underflow range error
* TLOSS 5 total loss of precision
* PLOSS 6 partial loss of precision
* EDOM 33 Unix domain error code
* ERANGE 34 Unix range error code
*
* The default version of the file prints the function name,
* passed to it by the pointer fctnam, followed by the
* error condition. The display is directed to the standard
* output device. The routine then returns to the calling
* program. Users may wish to modify the program to abort by
* calling exit() under severe error conditions such as domain
* errors.
*
* Since all error conditions pass control to this function,
* the display may be easily changed, eliminated, or directed
* to an error logging device.
*
* SEE ALSO:
*
* mconf.h
*
*/
/*
Cephes Math Library Release 2.0: April, 1987
Copyright 1984, 1987 by Stephen L. Moshier
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. Neither the name of the <ORGANIZATION> nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
#include <stdio.h>
#include "mconf.h"
int merror = 0;
/* Notice: the order of appearance of the following
* messages is bound to the error codes defined
* in mconf.h.
*/
static char *ermsg[7] = {
"unknown", /* error code 0 */
"domain", /* error code 1 */
"singularity", /* et seq. */
"overflow",
"underflow",
"total loss of precision",
"partial loss of precision"
};
int mtherr( name, code )
char *name;
int code;
{
/* Display string passed by calling program,
* which is supposed to be the name of the
* function in which the error occurred:
*/
printf( "\n%s ", name );
/* Set global error message word */
merror = code;
/* Display error message defined
* by the code argument.
*/
if( (code <= 0) || (code >= 7) )
code = 0;
printf( "%s error\n", ermsg[code] );
/* Return to calling
* program
*/
return( 0 );
}

122
source/common/thirdparty/math/polevl.c vendored Normal file
View file

@ -0,0 +1,122 @@
/* polevl.c
* p1evl.c
*
* Evaluate polynomial
*
*
*
* SYNOPSIS:
*
* int N;
* double x, y, coef[N+1], polevl[];
*
* y = polevl( x, coef, N );
*
*
*
* DESCRIPTION:
*
* Evaluates polynomial of degree N:
*
* 2 N
* y = C + C x + C x +...+ C x
* 0 1 2 N
*
* Coefficients are stored in reverse order:
*
* coef[0] = C , ..., coef[N] = C .
* N 0
*
* The function p1evl() assumes that coef[N] = 1.0 and is
* omitted from the array. Its calling arguments are
* otherwise the same as polevl().
*
*
* SPEED:
*
* In the interest of speed, there are no checks for out
* of bounds arithmetic. This routine is used by most of
* the functions in the library. Depending on available
* equipment features, the user may wish to rewrite the
* program in microcode or assembly language.
*
*/
/*
Cephes Math Library Release 2.1: December, 1988
Copyright 1984, 1987, 1988 by Stephen L. Moshier
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. Neither the name of the <ORGANIZATION> nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
double polevl( x, coef, N )
double x;
double coef[];
int N;
{
double ans;
int i;
double *p;
p = coef;
ans = *p++;
i = N;
do
ans = ans * x + *p++;
while( --i );
return( ans );
}
/* p1evl() */
/* N
* Evaluate polynomial when coefficient of x is 1.0.
* Otherwise same as polevl.
*/
double p1evl( x, coef, N )
double x;
double coef[];
int N;
{
double ans;
double *p;
int i;
p = coef;
ans = x + *p++;
i = N-1;
do
ans = ans * x + *p++;
while( --i );
return( ans );
}

780
source/common/thirdparty/math/pow.c vendored Normal file
View file

@ -0,0 +1,780 @@
/* pow.c
*
* Power function
*
*
*
* SYNOPSIS:
*
* double x, y, z, pow();
*
* z = pow( x, y );
*
*
*
* DESCRIPTION:
*
* Computes x raised to the yth power. Analytically,
*
* x**y = exp( y log(x) ).
*
* Following Cody and Waite, this program uses a lookup table
* of 2**-i/16 and pseudo extended precision arithmetic to
* obtain an extra three bits of accuracy in both the logarithm
* and the exponential.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -26,26 30000 4.2e-16 7.7e-17
* DEC -26,26 60000 4.8e-17 9.1e-18
* 1/26 < x < 26, with log(x) uniformly distributed.
* -26 < y < 26, y uniformly distributed.
* IEEE 0,8700 30000 1.5e-14 2.1e-15
* 0.99 < x < 1.01, 0 < y < 8700, uniformly distributed.
*
*
* ERROR MESSAGES:
*
* message condition value returned
* pow overflow x**y > MAXNUM INFINITY
* pow underflow x**y < 1/MAXNUM 0.0
* pow domain x<0 and y noninteger 0.0
*
*/
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. Neither the name of the <ORGANIZATION> nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
*/
#include "mconf.h"
static char fname[] = {"pow"};
#define SQRTH 0.70710678118654752440
#ifdef UNK
static double P[] = {
4.97778295871696322025E-1,
3.73336776063286838734E0,
7.69994162726912503298E0,
4.66651806774358464979E0
};
static double Q[] = {
/* 1.00000000000000000000E0, */
9.33340916416696166113E0,
2.79999886606328401649E1,
3.35994905342304405431E1,
1.39995542032307539578E1
};
/* 2^(-i/16), IEEE precision */
static double A[] = {
1.00000000000000000000E0,
9.57603280698573700036E-1,
9.17004043204671215328E-1,
8.78126080186649726755E-1,
8.40896415253714502036E-1,
8.05245165974627141736E-1,
7.71105412703970372057E-1,
7.38413072969749673113E-1,
7.07106781186547572737E-1,
6.77127773468446325644E-1,
6.48419777325504820276E-1,
6.20928906036742001007E-1,
5.94603557501360513449E-1,
5.69394317378345782288E-1,
5.45253866332628844837E-1,
5.22136891213706877402E-1,
5.00000000000000000000E-1
};
static double B[] = {
0.00000000000000000000E0,
1.64155361212281360176E-17,
4.09950501029074826006E-17,
3.97491740484881042808E-17,
-4.83364665672645672553E-17,
1.26912513974441574796E-17,
1.99100761573282305549E-17,
-1.52339103990623557348E-17,
0.00000000000000000000E0
};
static double R[] = {
1.49664108433729301083E-5,
1.54010762792771901396E-4,
1.33335476964097721140E-3,
9.61812908476554225149E-3,
5.55041086645832347466E-2,
2.40226506959099779976E-1,
6.93147180559945308821E-1
};
#define douba(k) A[k]
#define doubb(k) B[k]
#define MEXP 16383.0
#ifdef DENORMAL
#define MNEXP -17183.0
#else
#define MNEXP -16383.0
#endif
#endif
#ifdef DEC
static unsigned short P[] = {
0037776,0156313,0175332,0163602,
0040556,0167577,0052366,0174245,
0040766,0062753,0175707,0055564,
0040625,0052035,0131344,0155636,
};
static unsigned short Q[] = {
/*0040200,0000000,0000000,0000000,*/
0041025,0052644,0154404,0105155,
0041337,0177772,0007016,0047646,
0041406,0062740,0154273,0020020,
0041137,0177054,0106127,0044555,
};
static unsigned short A[] = {
0040200,0000000,0000000,0000000,
0040165,0022575,0012444,0103314,
0040152,0140306,0163735,0022071,
0040140,0146336,0166052,0112341,
0040127,0042374,0145326,0116553,
0040116,0022214,0012437,0102201,
0040105,0063452,0010525,0003333,
0040075,0004243,0117530,0006067,
0040065,0002363,0031771,0157145,
0040055,0054076,0165102,0120513,
0040045,0177326,0124661,0050471,
0040036,0172462,0060221,0120422,
0040030,0033760,0050615,0134251,
0040021,0141723,0071653,0010703,
0040013,0112701,0161752,0105727,
0040005,0125303,0063714,0044173,
0040000,0000000,0000000,0000000
};
static unsigned short B[] = {
0000000,0000000,0000000,0000000,
0021473,0040265,0153315,0140671,
0121074,0062627,0042146,0176454,
0121413,0003524,0136332,0066212,
0121767,0046404,0166231,0012553,
0121257,0015024,0002357,0043574,
0021736,0106532,0043060,0056206,
0121310,0020334,0165705,0035326,
0000000,0000000,0000000,0000000
};
static unsigned short R[] = {
0034173,0014076,0137624,0115771,
0035041,0076763,0003744,0111311,
0035656,0141766,0041127,0074351,
0036435,0112533,0073611,0116664,
0037143,0054106,0134040,0152223,
0037565,0176757,0176026,0025551,
0040061,0071027,0173721,0147572
};
/*
static double R[] = {
0.14928852680595608186e-4,
0.15400290440989764601e-3,
0.13333541313585784703e-2,
0.96181290595172416964e-2,
0.55504108664085595326e-1,
0.24022650695909537056e0,
0.69314718055994529629e0
};
*/
#define douba(k) (*(double *)&A[(k)<<2])
#define doubb(k) (*(double *)&B[(k)<<2])
#define MEXP 2031.0
#define MNEXP -2031.0
#endif
#ifdef IBMPC
static unsigned short P[] = {
0x5cf0,0x7f5b,0xdb99,0x3fdf,
0xdf15,0xea9e,0xddef,0x400d,
0xeb6f,0x7f78,0xccbd,0x401e,
0x9b74,0xb65c,0xaa83,0x4012,
};
static unsigned short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0x914e,0x9b20,0xaab4,0x4022,
0xc9f5,0x41c1,0xffff,0x403b,
0x6402,0x1b17,0xccbc,0x4040,
0xe92e,0x918a,0xffc5,0x402b,
};
static unsigned short A[] = {
0x0000,0x0000,0x0000,0x3ff0,
0x90da,0xa2a4,0xa4af,0x3fee,
0xa487,0xdcfb,0x5818,0x3fed,
0x529c,0xdd85,0x199b,0x3fec,
0xd3ad,0x995a,0xe89f,0x3fea,
0xf090,0x82a3,0xc491,0x3fe9,
0xa0db,0x422a,0xace5,0x3fe8,
0x0187,0x73eb,0xa114,0x3fe7,
0x3bcd,0x667f,0xa09e,0x3fe6,
0x5429,0xdd48,0xab07,0x3fe5,
0x2a27,0xd536,0xbfda,0x3fe4,
0x3422,0x4c12,0xdea6,0x3fe3,
0xb715,0x0a31,0x06fe,0x3fe3,
0x6238,0x6e75,0x387a,0x3fe2,
0x517b,0x3c7d,0x72b8,0x3fe1,
0x890f,0x6cf9,0xb558,0x3fe0,
0x0000,0x0000,0x0000,0x3fe0
};
static unsigned short B[] = {
0x0000,0x0000,0x0000,0x0000,
0x3707,0xd75b,0xed02,0x3c72,
0xcc81,0x345d,0xa1cd,0x3c87,
0x4b27,0x5686,0xe9f1,0x3c86,
0x6456,0x13b2,0xdd34,0xbc8b,
0x42e2,0xafec,0x4397,0x3c6d,
0x82e4,0xd231,0xf46a,0x3c76,
0x8a76,0xb9d7,0x9041,0xbc71,
0x0000,0x0000,0x0000,0x0000
};
static unsigned short R[] = {
0x937f,0xd7f2,0x6307,0x3eef,
0x9259,0x60fc,0x2fbe,0x3f24,
0xef1d,0xc84a,0xd87e,0x3f55,
0x33b7,0x6ef1,0xb2ab,0x3f83,
0x1a92,0xd704,0x6b08,0x3fac,
0xc56d,0xff82,0xbfbd,0x3fce,
0x39ef,0xfefa,0x2e42,0x3fe6
};
#define douba(k) (*(double *)&A[(k)<<2])
#define doubb(k) (*(double *)&B[(k)<<2])
#define MEXP 16383.0
#ifdef DENORMAL
#define MNEXP -17183.0
#else
#define MNEXP -16383.0
#endif
#endif
#ifdef MIEEE
static unsigned short P[] = {
0x3fdf,0xdb99,0x7f5b,0x5cf0,
0x400d,0xddef,0xea9e,0xdf15,
0x401e,0xccbd,0x7f78,0xeb6f,
0x4012,0xaa83,0xb65c,0x9b74
};
static unsigned short Q[] = {
0x4022,0xaab4,0x9b20,0x914e,
0x403b,0xffff,0x41c1,0xc9f5,
0x4040,0xccbc,0x1b17,0x6402,
0x402b,0xffc5,0x918a,0xe92e
};
static unsigned short A[] = {
0x3ff0,0x0000,0x0000,0x0000,
0x3fee,0xa4af,0xa2a4,0x90da,
0x3fed,0x5818,0xdcfb,0xa487,
0x3fec,0x199b,0xdd85,0x529c,
0x3fea,0xe89f,0x995a,0xd3ad,
0x3fe9,0xc491,0x82a3,0xf090,
0x3fe8,0xace5,0x422a,0xa0db,
0x3fe7,0xa114,0x73eb,0x0187,
0x3fe6,0xa09e,0x667f,0x3bcd,
0x3fe5,0xab07,0xdd48,0x5429,
0x3fe4,0xbfda,0xd536,0x2a27,
0x3fe3,0xdea6,0x4c12,0x3422,
0x3fe3,0x06fe,0x0a31,0xb715,
0x3fe2,0x387a,0x6e75,0x6238,
0x3fe1,0x72b8,0x3c7d,0x517b,
0x3fe0,0xb558,0x6cf9,0x890f,
0x3fe0,0x0000,0x0000,0x0000
};
static unsigned short B[] = {
0x0000,0x0000,0x0000,0x0000,
0x3c72,0xed02,0xd75b,0x3707,
0x3c87,0xa1cd,0x345d,0xcc81,
0x3c86,0xe9f1,0x5686,0x4b27,
0xbc8b,0xdd34,0x13b2,0x6456,
0x3c6d,0x4397,0xafec,0x42e2,
0x3c76,0xf46a,0xd231,0x82e4,
0xbc71,0x9041,0xb9d7,0x8a76,
0x0000,0x0000,0x0000,0x0000
};
static unsigned short R[] = {
0x3eef,0x6307,0xd7f2,0x937f,
0x3f24,0x2fbe,0x60fc,0x9259,
0x3f55,0xd87e,0xc84a,0xef1d,
0x3f83,0xb2ab,0x6ef1,0x33b7,
0x3fac,0x6b08,0xd704,0x1a92,
0x3fce,0xbfbd,0xff82,0xc56d,
0x3fe6,0x2e42,0xfefa,0x39ef
};
#define douba(k) (*(double *)&A[(k)<<2])
#define doubb(k) (*(double *)&B[(k)<<2])
#define MEXP 16383.0
#ifdef DENORMAL
#define MNEXP -17183.0
#else
#define MNEXP -16383.0
#endif
#endif
/* log2(e) - 1 */
#define LOG2EA 0.44269504088896340736
#define F W
#define Fa Wa
#define Fb Wb
#define G W
#define Ga Wa
#define Gb u
#define H W
#define Ha Wb
#define Hb Wb
#ifdef ANSIPROT
extern double floor ( double );
extern double fabs ( double );
extern double frexp ( double, int * );
extern double ldexp ( double, int );
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern double c_powi ( double, int );
extern int signbit ( double );
extern int isnan ( double );
extern int isfinite ( double );
static double reduc ( double );
#else
double floor(), fabs(), frexp(), ldexp();
double polevl(), p1evl(), c_powi();
int signbit(), isnan(), isfinite();
static double reduc();
#endif
extern double MAXNUM;
#ifdef INFINITIES
extern double INFINITY;
#endif
#ifdef NANS
extern double NAN;
#endif
#ifdef MINUSZERO
extern double NEGZERO;
#endif
double c_pow( x, y )
double x, y;
{
double w, z, W, Wa, Wb, ya, yb, u;
/* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */
double aw, ay, wy;
int e, i, nflg, iyflg, yoddint;
if( y == 0.0 )
return( 1.0 );
#ifdef NANS
if( isnan(x) )
return( x );
if( isnan(y) )
return( y );
#endif
if( y == 1.0 )
return( x );
#ifdef INFINITIES
if( !isfinite(y) && (x == 1.0 || x == -1.0) )
{
mtherr( "pow", DOMAIN );
#ifdef NANS
return( NAN );
#else
return( INFINITY );
#endif
}
#endif
if( x == 1.0 )
return( 1.0 );
if( y >= MAXNUM )
{
#ifdef INFINITIES
if( x > 1.0 )
return( INFINITY );
#else
if( x > 1.0 )
return( MAXNUM );
#endif
if( x > 0.0 && x < 1.0 )
return( 0.0);
if( x < -1.0 )
{
#ifdef INFINITIES
return( INFINITY );
#else
return( MAXNUM );
#endif
}
if( x > -1.0 && x < 0.0 )
return( 0.0 );
}
if( y <= -MAXNUM )
{
if( x > 1.0 )
return( 0.0 );
#ifdef INFINITIES
if( x > 0.0 && x < 1.0 )
return( INFINITY );
#else
if( x > 0.0 && x < 1.0 )
return( MAXNUM );
#endif
if( x < -1.0 )
return( 0.0 );
#ifdef INFINITIES
if( x > -1.0 && x < 0.0 )
return( INFINITY );
#else
if( x > -1.0 && x < 0.0 )
return( MAXNUM );
#endif
}
if( x >= MAXNUM )
{
#if INFINITIES
if( y > 0.0 )
return( INFINITY );
#else
if( y > 0.0 )
return( MAXNUM );
#endif
return(0.0);
}
/* Set iyflg to 1 if y is an integer. */
iyflg = 0;
w = floor(y);
if( w == y )
iyflg = 1;
/* Test for odd integer y. */
yoddint = 0;
if( iyflg )
{
ya = fabs(y);
ya = floor(0.5 * ya);
yb = 0.5 * fabs(w);
if( ya != yb )
yoddint = 1;
}
if( x <= -MAXNUM )
{
if( y > 0.0 )
{
#ifdef INFINITIES
if( yoddint )
return( -INFINITY );
return( INFINITY );
#else
if( yoddint )
return( -MAXNUM );
return( MAXNUM );
#endif
}
if( y < 0.0 )
{
#ifdef MINUSZERO
if( yoddint )
return( NEGZERO );
#endif
return( 0.0 );
}
}
nflg = 0; /* flag = 1 if x<0 raised to integer power */
if( x <= 0.0 )
{
if( x == 0.0 )
{
if( y < 0.0 )
{
#ifdef MINUSZERO
if( signbit(x) && yoddint )
return( -INFINITY );
#endif
#ifdef INFINITIES
return( INFINITY );
#else
return( MAXNUM );
#endif
}
if( y > 0.0 )
{
#ifdef MINUSZERO
if( signbit(x) && yoddint )
return( NEGZERO );
#endif
return( 0.0 );
}
return( 1.0 );
}
else
{
if( iyflg == 0 )
{ /* noninteger power of negative number */
mtherr( fname, DOMAIN );
#ifdef NANS
return(NAN);
#else
return(0.0L);
#endif
}
nflg = 1;
}
}
/* Integer power of an integer. */
if( iyflg )
{
i = (int)w;
w = floor(x);
if( (w == x) && (fabs(y) < 32768.0) )
{
w = c_powi( x, (int) y );
return( w );
}
}
if( nflg )
x = fabs(x);
/* For results close to 1, use a series expansion. */
w = x - 1.0;
aw = fabs(w);
ay = fabs(y);
wy = w * y;
ya = fabs(wy);
if((aw <= 1.0e-3 && ay <= 1.0)
|| (ya <= 1.0e-3 && ay >= 1.0))
{
z = (((((w*(y-5.)/720. + 1./120.)*w*(y-4.) + 1./24.)*w*(y-3.)
+ 1./6.)*w*(y-2.) + 0.5)*w*(y-1.) )*wy + wy + 1.;
goto done;
}
/* These are probably too much trouble. */
#if 0
w = y * log(x);
if (aw > 1.0e-3 && fabs(w) < 1.0e-3)
{
z = ((((((
w/7. + 1.)*w/6. + 1.)*w/5. + 1.)*w/4. + 1.)*w/3. + 1.)*w/2. + 1.)*w + 1.;
goto done;
}
if(ya <= 1.0e-3 && aw <= 1.0e-4)
{
z = (((((
wy*1./720.
+ (-w*1./48. + 1./120.) )*wy
+ ((w*17./144. - 1./12.)*w + 1./24.) )*wy
+ (((-w*5./16. + 7./24.)*w - 1./4.)*w + 1./6.) )*wy
+ ((((w*137./360. - 5./12.)*w + 11./24.)*w - 1./2.)*w + 1./2.) )*wy
+ (((((-w*1./6. + 1./5.)*w - 1./4)*w + 1./3.)*w -1./2.)*w ) )*wy
+ wy + 1.0;
goto done;
}
#endif
/* separate significand from exponent */
x = frexp( x, &e );
#if 0
/* For debugging, check for gross overflow. */
if( (e * y) > (MEXP + 1024) )
goto overflow;
#endif
/* Find significand of x in antilog table A[]. */
i = 1;
if( x <= douba(9) )
i = 9;
if( x <= douba(i+4) )
i += 4;
if( x <= douba(i+2) )
i += 2;
if( x >= douba(1) )
i = -1;
i += 1;
/* Find (x - A[i])/A[i]
* in order to compute log(x/A[i]):
*
* log(x) = log( a x/a ) = log(a) + log(x/a)
*
* log(x/a) = log(1+v), v = x/a - 1 = (x-a)/a
*/
x -= douba(i);
x -= doubb(i/2);
x /= douba(i);
/* rational approximation for log(1+v):
*
* log(1+v) = v - v**2/2 + v**3 P(v) / Q(v)
*/
z = x*x;
w = x * ( z * polevl( x, P, 3 ) / p1evl( x, Q, 4 ) );
w = w - ldexp( z, -1 ); /* w - 0.5 * z */
/* Convert to base 2 logarithm:
* multiply by log2(e)
*/
w = w + LOG2EA * w;
/* Note x was not yet added in
* to above rational approximation,
* so do it now, while multiplying
* by log2(e).
*/
z = w + LOG2EA * x;
z = z + x;
/* Compute exponent term of the base 2 logarithm. */
w = -i;
w = ldexp( w, -4 ); /* divide by 16 */
w += e;
/* Now base 2 log of x is w + z. */
/* Multiply base 2 log by y, in extended precision. */
/* separate y into large part ya
* and small part yb less than 1/16
*/
ya = reduc(y);
yb = y - ya;
F = z * y + w * yb;
Fa = reduc(F);
Fb = F - Fa;
G = Fa + w * ya;
Ga = reduc(G);
Gb = G - Ga;
H = Fb + Gb;
Ha = reduc(H);
w = ldexp( Ga+Ha, 4 );
/* Test the power of 2 for overflow */
if( w > MEXP )
{
#ifndef INFINITIES
mtherr( fname, OVERFLOW );
#endif
#ifdef INFINITIES
if( nflg && yoddint )
return( -INFINITY );
return( INFINITY );
#else
if( nflg && yoddint )
return( -MAXNUM );
return( MAXNUM );
#endif
}
if( w < (MNEXP - 1) )
{
#ifndef DENORMAL
mtherr( fname, UNDERFLOW );
#endif
#ifdef MINUSZERO
if( nflg && yoddint )
return( NEGZERO );
#endif
return( 0.0 );
}
e = (int)w;
Hb = H - Ha;
if( Hb > 0.0 )
{
e += 1;
Hb -= 0.0625;
}
/* Now the product y * log2(x) = Hb + e/16.0.
*
* Compute base 2 exponential of Hb,
* where -0.0625 <= Hb <= 0.
*/
z = Hb * polevl( Hb, R, 6 ); /* z = 2**Hb - 1 */
/* Express e/16 as an integer plus a negative number of 16ths.
* Find lookup table entry for the fractional power of 2.
*/
if( e < 0 )
i = 0;
else
i = 1;
i = e/16 + i;
e = 16*i - e;
w = douba( e );
z = w + w * z; /* 2**-e * ( 1 + (2**Hb-1) ) */
z = ldexp( z, i ); /* multiply by integer power of 2 */
done:
/* Negate if odd integer power of negative number */
if( nflg && yoddint )
{
#ifdef MINUSZERO
if( z == 0.0 )
z = NEGZERO;
else
#endif
z = -z;
}
return( z );
}
/* Find a multiple of 1/16 that is within 1/16 of x. */
static double reduc(x)
double x;
{
double t;
t = ldexp( x, 4 );
t = floor( t );
t = ldexp( t, -4 );
return(t);
}

210
source/common/thirdparty/math/powi.c vendored Normal file
View file

@ -0,0 +1,210 @@
/* powi.c
*
* Real raised to integer power
*
*
*
* SYNOPSIS:
*
* double x, y, powi();
* int n;
*
* y = powi( x, n );
*
*
*
* DESCRIPTION:
*
* Returns argument x raised to the nth power.
* The routine efficiently decomposes n as a sum of powers of
* two. The desired power is a product of two-to-the-kth
* powers of x. Thus to compute the 32767 power of x requires
* 28 multiplications instead of 32767 multiplications.
*
*
*
* ACCURACY:
*
*
* Relative error:
* arithmetic x domain n domain # trials peak rms
* DEC .04,26 -26,26 100000 2.7e-16 4.3e-17
* IEEE .04,26 -26,26 50000 2.0e-15 3.8e-16
* IEEE 1,2 -1022,1023 50000 8.6e-14 1.6e-14
*
* Returns MAXNUM on overflow, zero on underflow.
*
*/
/* powi.c */
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. Neither the name of the <ORGANIZATION> nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
*/
#include "mconf.h"
#ifdef ANSIPROT
extern double log ( double );
extern double frexp ( double, int * );
extern int signbit ( double );
#else
double log(), frexp();
int signbit();
#endif
extern double NEGZERO, INFINITY, MAXNUM, MAXLOG, MINLOG, LOGE2;
double c_powi( x, nn )
double x;
int nn;
{
int n, e, sign, asign, lx;
double w, y, s;
/* See pow.c for these tests. */
if( x == 0.0 )
{
if( nn == 0 )
return( 1.0 );
else if( nn < 0 )
return( INFINITY );
else
{
if( nn & 1 )
return( x );
else
return( 0.0 );
}
}
if( nn == 0 )
return( 1.0 );
if( nn == -1 )
return( 1.0/x );
if( x < 0.0 )
{
asign = -1;
x = -x;
}
else
asign = 0;
if( nn < 0 )
{
sign = -1;
n = -nn;
}
else
{
sign = 1;
n = nn;
}
/* Even power will be positive. */
if( (n & 1) == 0 )
asign = 0;
/* Overflow detection */
/* Calculate approximate logarithm of answer */
s = frexp( x, &lx );
e = (lx - 1)*n;
if( (e == 0) || (e > 64) || (e < -64) )
{
s = (s - 7.0710678118654752e-1) / (s + 7.0710678118654752e-1);
s = (2.9142135623730950 * s - 0.5 + lx) * nn * LOGE2;
}
else
{
s = LOGE2 * e;
}
if( s > MAXLOG )
{
mtherr( "powi", OVERFLOW );
y = INFINITY;
goto done;
}
#if DENORMAL
if( s < MINLOG )
{
y = 0.0;
goto done;
}
/* Handle tiny denormal answer, but with less accuracy
* since roundoff error in 1.0/x will be amplified.
* The precise demarcation should be the gradual underflow threshold.
*/
if( (s < (-MAXLOG+2.0)) && (sign < 0) )
{
x = 1.0/x;
sign = -sign;
}
#else
/* do not produce denormal answer */
if( s < -MAXLOG )
return(0.0);
#endif
/* First bit of the power */
if( n & 1 )
y = x;
else
y = 1.0;
w = x;
n >>= 1;
while( n )
{
w = w * w; /* arg to the 2-to-the-kth power */
if( n & 1 ) /* if that bit is set, then include in product */
y *= w;
n >>= 1;
}
if( sign < 0 )
y = 1.0/y;
done:
if( asign )
{
/* odd power of negative number */
if( y == 0.0 )
y = NEGZERO;
else
y = -y;
}
return(y);
}

411
source/common/thirdparty/math/sin.c vendored Normal file
View file

@ -0,0 +1,411 @@
/* sin.c
*
* Circular sine
*
*
*
* SYNOPSIS:
*
* double x, y, sin();
*
* y = sin( x );
*
*
*
* DESCRIPTION:
*
* Range reduction is into intervals of pi/4. The reduction
* error is nearly eliminated by contriving an extended precision
* modular arithmetic.
*
* Two polynomial approximating functions are employed.
* Between 0 and pi/4 the sine is approximated by
* x + x**3 P(x**2).
* Between pi/4 and pi/2 the cosine is represented as
* 1 - x**2 Q(x**2).
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC 0, 10 150000 3.0e-17 7.8e-18
* IEEE -1.07e9,+1.07e9 130000 2.1e-16 5.4e-17
*
* ERROR MESSAGES:
*
* message condition value returned
* sin total loss x > 1.073741824e9 0.0
*
* Partial loss of accuracy begins to occur at x = 2**30
* = 1.074e9. The loss is not gradual, but jumps suddenly to
* about 1 part in 10e7. Results may be meaningless for
* x > 2**49 = 5.6e14. The routine as implemented flags a
* TLOSS error for x > 2**30 and returns 0.0.
*/
/* cos.c
*
* Circular cosine
*
*
*
* SYNOPSIS:
*
* double x, y, cos();
*
* y = cos( x );
*
*
*
* DESCRIPTION:
*
* Range reduction is into intervals of pi/4. The reduction
* error is nearly eliminated by contriving an extended precision
* modular arithmetic.
*
* Two polynomial approximating functions are employed.
* Between 0 and pi/4 the cosine is approximated by
* 1 - x**2 Q(x**2).
* Between pi/4 and pi/2 the sine is represented as
* x + x**3 P(x**2).
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -1.07e9,+1.07e9 130000 2.1e-16 5.4e-17
* DEC 0,+1.07e9 17000 3.0e-17 7.2e-18
*/
/* sin.c */
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1985, 1995, 2000 by Stephen L. Moshier
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. Neither the name of the <ORGANIZATION> nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
*/
#include "mconf.h"
#ifdef UNK
static double sincof[] = {
1.58962301576546568060E-10,
-2.50507477628578072866E-8,
2.75573136213857245213E-6,
-1.98412698295895385996E-4,
8.33333333332211858878E-3,
-1.66666666666666307295E-1,
};
static double coscof[6] = {
-1.13585365213876817300E-11,
2.08757008419747316778E-9,
-2.75573141792967388112E-7,
2.48015872888517045348E-5,
-1.38888888888730564116E-3,
4.16666666666665929218E-2,
};
static double DP1 = 7.85398125648498535156E-1;
static double DP2 = 3.77489470793079817668E-8;
static double DP3 = 2.69515142907905952645E-15;
/* static double lossth = 1.073741824e9; */
#endif
#ifdef DEC
static unsigned short sincof[] = {
0030056,0143750,0177214,0163153,
0131727,0027455,0044510,0175352,
0033470,0167432,0131752,0042414,
0135120,0006400,0146776,0174027,
0036410,0104210,0104207,0137202,
0137452,0125252,0125252,0125103,
};
static unsigned short coscof[24] = {
0127107,0151115,0002060,0152325,
0031017,0072353,0155161,0174053,
0132623,0171173,0172542,0057056,
0034320,0006400,0147102,0023652,
0135666,0005540,0133012,0076213,
0037052,0125252,0125252,0125126,
};
/* 7.853981629014015197753906250000E-1 */
static unsigned short P1[] = {0040111,0007732,0120000,0000000,};
/* 4.960467869796758577649598009884E-10 */
static unsigned short P2[] = {0030410,0055060,0100000,0000000,};
/* 2.860594363054915898381331279295E-18 */
static unsigned short P3[] = {0021523,0011431,0105056,0001560,};
#define DP1 *(double *)P1
#define DP2 *(double *)P2
#define DP3 *(double *)P3
#endif
#ifdef IBMPC
static unsigned short sincof[] = {
0x9ccd,0x1fd1,0xd8fd,0x3de5,
0x1f5d,0xa929,0xe5e5,0xbe5a,
0x48a1,0x567d,0x1de3,0x3ec7,
0xdf03,0x19bf,0x01a0,0xbf2a,
0xf7d0,0x1110,0x1111,0x3f81,
0x5548,0x5555,0x5555,0xbfc5,
};
static unsigned short coscof[24] = {
0x1a9b,0xa086,0xfa49,0xbda8,
0x3f05,0x7b4e,0xee9d,0x3e21,
0x4bc6,0x7eac,0x7e4f,0xbe92,
0x44f5,0x19c8,0x01a0,0x3efa,
0x4f91,0x16c1,0xc16c,0xbf56,
0x554b,0x5555,0x5555,0x3fa5,
};
/*
7.85398125648498535156E-1,
3.77489470793079817668E-8,
2.69515142907905952645E-15,
*/
static unsigned short P1[] = {0x0000,0x4000,0x21fb,0x3fe9};
static unsigned short P2[] = {0x0000,0x0000,0x442d,0x3e64};
static unsigned short P3[] = {0x5170,0x98cc,0x4698,0x3ce8};
#define DP1 *(double *)P1
#define DP2 *(double *)P2
#define DP3 *(double *)P3
#endif
#ifdef MIEEE
static unsigned short sincof[] = {
0x3de5,0xd8fd,0x1fd1,0x9ccd,
0xbe5a,0xe5e5,0xa929,0x1f5d,
0x3ec7,0x1de3,0x567d,0x48a1,
0xbf2a,0x01a0,0x19bf,0xdf03,
0x3f81,0x1111,0x1110,0xf7d0,
0xbfc5,0x5555,0x5555,0x5548,
};
static unsigned short coscof[24] = {
0xbda8,0xfa49,0xa086,0x1a9b,
0x3e21,0xee9d,0x7b4e,0x3f05,
0xbe92,0x7e4f,0x7eac,0x4bc6,
0x3efa,0x01a0,0x19c8,0x44f5,
0xbf56,0xc16c,0x16c1,0x4f91,
0x3fa5,0x5555,0x5555,0x554b,
};
static unsigned short P1[] = {0x3fe9,0x21fb,0x4000,0x0000};
static unsigned short P2[] = {0x3e64,0x442d,0x0000,0x0000};
static unsigned short P3[] = {0x3ce8,0x4698,0x98cc,0x5170};
#define DP1 *(double *)P1
#define DP2 *(double *)P2
#define DP3 *(double *)P3
#endif
#ifdef ANSIPROT
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern double floor ( double );
extern double ldexp ( double, int );
extern int isnan ( double );
extern int isfinite ( double );
#else
double polevl(), floor(), ldexp();
int isnan(), isfinite();
#endif
extern double PIO4;
static double lossth = 1.073741824e9;
#ifdef NANS
extern double NAN;
#endif
#ifdef INFINITIES
extern double INFINITY;
#endif
double c_sin(x)
double x;
{
double y, z, zz;
int j, sign;
#ifdef MINUSZERO
if( x == 0.0 )
return(x);
#endif
#ifdef NANS
if( isnan(x) )
return(x);
if( !isfinite(x) )
{
mtherr( "sin", DOMAIN );
return(NAN);
}
#endif
/* make argument positive but save the sign */
sign = 1;
if( x < 0 )
{
x = -x;
sign = -1;
}
if( x > lossth )
{
mtherr( "sin", TLOSS );
return(0.0);
}
y = floor( x/PIO4 ); /* integer part of x/PIO4 */
/* strip high bits of integer part to prevent integer overflow */
z = ldexp( y, -4 );
z = floor(z); /* integer part of y/8 */
z = y - ldexp( z, 4 ); /* y - 16 * (y/16) */
j = (int)z; /* convert to integer for tests on the phase angle */
/* map zeros to origin */
if( j & 1 )
{
j += 1;
y += 1.0;
}
j = j & 07; /* octant modulo 360 degrees */
/* reflect in x axis */
if( j > 3)
{
sign = -sign;
j -= 4;
}
/* Extended precision modular arithmetic */
z = ((x - y * DP1) - y * DP2) - y * DP3;
zz = z * z;
if( (j==1) || (j==2) )
{
y = 1.0 - ldexp(zz,-1) + zz * zz * polevl( zz, coscof, 5 );
}
else
{
/* y = z + z * (zz * polevl( zz, sincof, 5 ));*/
y = z + z * z * z * polevl( zz, sincof, 5 );
}
if(sign < 0)
y = -y;
return(y);
}
double c_cos(x)
double x;
{
double y, z, zz;
int i;
int j, sign;
#ifdef NANS
if( isnan(x) )
return(x);
if( !isfinite(x) )
{
mtherr( "cos", DOMAIN );
return(NAN);
}
#endif
/* make argument positive */
sign = 1;
if( x < 0 )
x = -x;
if( x > lossth )
{
mtherr( "cos", TLOSS );
return(0.0);
}
y = floor( x/PIO4 );
z = ldexp( y, -4 );
z = floor(z); /* integer part of y/8 */
z = y - ldexp( z, 4 ); /* y - 16 * (y/16) */
/* integer and fractional part modulo one octant */
i = (int)z;
if( i & 1 ) /* map zeros to origin */
{
i += 1;
y += 1.0;
}
j = i & 07;
if( j > 3)
{
j -=4;
sign = -sign;
}
if( j > 1 )
sign = -sign;
/* Extended precision modular arithmetic */
z = ((x - y * DP1) - y * DP2) - y * DP3;
zz = z * z;
if( (j==1) || (j==2) )
{
/* y = z + z * (zz * polevl( zz, sincof, 5 ));*/
y = z + z * z * z * polevl( zz, sincof, 5 );
}
else
{
y = 1.0 - ldexp(zz,-1) + zz * zz * polevl( zz, coscof, 5 );
}
if(sign < 0)
y = -y;
return(y);
}
/* Degrees, minutes, seconds to radians: */
/* 1 arc second, in radians = 4.8481368110953599358991410e-5 */
#ifdef DEC
static unsigned short P648[] = {034513,054170,0176773,0116043,};
#define P64800 *(double *)P648
#else
static double P64800 = 4.8481368110953599358991410e-5;
#endif
double radian(d,m,s)
double d,m,s;
{
return( ((d*60.0 + m)*60.0 + s)*P64800 );
}

172
source/common/thirdparty/math/sinh.c vendored Normal file
View file

@ -0,0 +1,172 @@
/* sinh.c
*
* Hyperbolic sine
*
*
*
* SYNOPSIS:
*
* double x, y, sinh();
*
* y = sinh( x );
*
*
*
* DESCRIPTION:
*
* Returns hyperbolic sine of argument in the range MINLOG to
* MAXLOG.
*
* The range is partitioned into two segments. If |x| <= 1, a
* rational function of the form x + x**3 P(x)/Q(x) is employed.
* Otherwise the calculation is sinh(x) = ( exp(x) - exp(-x) )/2.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC +- 88 50000 4.0e-17 7.7e-18
* IEEE +-MAXLOG 30000 2.6e-16 5.7e-17
*
*/
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. Neither the name of the <ORGANIZATION> nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
*/
#include "mconf.h"
#ifdef UNK
static double P[] = {
-7.89474443963537015605E-1,
-1.63725857525983828727E2,
-1.15614435765005216044E4,
-3.51754964808151394800E5
};
static double Q[] = {
/* 1.00000000000000000000E0,*/
-2.77711081420602794433E2,
3.61578279834431989373E4,
-2.11052978884890840399E6
};
#endif
#ifdef DEC
static unsigned short P[] = {
0140112,0015377,0042731,0163255,
0142043,0134721,0146177,0123761,
0143464,0122706,0034353,0006017,
0144653,0140536,0157665,0054045
};
static unsigned short Q[] = {
/*0040200,0000000,0000000,0000000,*/
0142212,0155404,0133513,0022040,
0044015,0036723,0173271,0011053,
0145400,0150407,0023710,0001034
};
#endif
#ifdef IBMPC
static unsigned short P[] = {
0x3cd6,0xe8bb,0x435f,0xbfe9,
0xf4fe,0x398f,0x773a,0xc064,
0x6182,0xc71d,0x94b8,0xc0c6,
0xab05,0xdbf6,0x782b,0xc115
};
static unsigned short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0x6484,0x96e9,0x5b60,0xc071,
0x2245,0x7ed7,0xa7ba,0x40e1,
0x0044,0xe4f9,0x1a20,0xc140
};
#endif
#ifdef MIEEE
static unsigned short P[] = {
0xbfe9,0x435f,0xe8bb,0x3cd6,
0xc064,0x773a,0x398f,0xf4fe,
0xc0c6,0x94b8,0xc71d,0x6182,
0xc115,0x782b,0xdbf6,0xab05
};
static unsigned short Q[] = {
0xc071,0x5b60,0x96e9,0x6484,
0x40e1,0xa7ba,0x7ed7,0x2245,
0xc140,0x1a20,0xe4f9,0x0044
};
#endif
#ifdef ANSIPROT
extern double fabs ( double );
extern double c_exp ( double );
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
#else
double fabs(), c_exp(), polevl(), p1evl();
#endif
extern double INFINITY, MINLOG, MAXLOG, LOGE2;
double c_sinh(x)
double x;
{
double a;
#ifdef MINUSZERO
if( x == 0.0 )
return(x);
#endif
a = fabs(x);
if( (x > (MAXLOG + LOGE2)) || (x > -(MINLOG-LOGE2) ) )
{
mtherr( "sinh", DOMAIN );
if( x > 0 )
return( INFINITY );
else
return( -INFINITY );
}
if( a > 1.0 )
{
if( a >= (MAXLOG - LOGE2) )
{
a = c_exp(0.5*a);
a = (0.5 * a) * a;
if( x < 0 )
a = -a;
return(a);
}
a = c_exp(a);
a = 0.5*a - (0.5/a);
if( x < 0 )
a = -a;
return(a);
}
a *= a;
return( x + x * a * (polevl(a,P,3)/p1evl(a,Q,3)) );
}

202
source/common/thirdparty/math/sqrt.c vendored Normal file
View file

@ -0,0 +1,202 @@
/* _sqrt.c
*
* Square root
*
*
*
* SYNOPSIS:
*
* double x, y, _sqrt();
*
* y = _sqrt( x );
*
*
*
* DESCRIPTION:
*
* Returns the square root of x.
*
* Range reduction involves isolating the power of two of the
* argument and using a polynomial approximation to obtain
* a rough value for the square root. Then Heron's iteration
* is used three times to converge to an accurate value.
*
*
*
* ACCURACY:
*
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC 0, 10 60000 2.1e-17 7.9e-18
* IEEE 0,1.7e308 30000 1.7e-16 6.3e-17
*
*
* ERROR MESSAGES:
*
* message condition value returned
* _sqrt domain x < 0 0.0
*
*/
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. Neither the name of the <ORGANIZATION> nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
*/
#include "mconf.h"
#ifdef ANSIPROT
extern double frexp ( double, int * );
extern double ldexp ( double, int );
#else
double frexp(), ldexp();
#endif
extern double SQRT2; /* _sqrt2 = 1.41421356237309504880 */
double c_sqrt(x)
double x;
{
int e;
#ifndef UNK
short *q;
#endif
double z, w;
if( x <= 0.0 )
{
if( x < 0.0 )
mtherr( "_sqrt", DOMAIN );
return( 0.0 );
}
w = x;
/* separate exponent and significand */
#ifdef UNK
z = frexp( x, &e );
#endif
#ifdef DEC
q = (short *)&x;
e = ((*q >> 7) & 0377) - 0200;
*q &= 0177;
*q |= 040000;
z = x;
#endif
/* Note, frexp and ldexp are used in order to
* handle denormal numbers properly.
*/
#ifdef IBMPC
z = frexp( x, &e );
q = (short *)&x;
q += 3;
/*
e = ((*q >> 4) & 0x0fff) - 0x3fe;
*q &= 0x000f;
*q |= 0x3fe0;
z = x;
*/
#endif
#ifdef MIEEE
z = frexp( x, &e );
q = (short *)&x;
/*
e = ((*q >> 4) & 0x0fff) - 0x3fe;
*q &= 0x000f;
*q |= 0x3fe0;
z = x;
*/
#endif
/* approximate square root of number between 0.5 and 1
* relative error of approximation = 7.47e-3
*/
x = 4.173075996388649989089E-1 + 5.9016206709064458299663E-1 * z;
/* adjust for odd powers of 2 */
if( (e & 1) != 0 )
x *= SQRT2;
/* re-insert exponent */
#ifdef UNK
x = ldexp( x, (e >> 1) );
#endif
#ifdef DEC
*q += ((e >> 1) & 0377) << 7;
*q &= 077777;
#endif
#ifdef IBMPC
x = ldexp( x, (e >> 1) );
/*
*q += ((e >>1) & 0x7ff) << 4;
*q &= 077777;
*/
#endif
#ifdef MIEEE
x = ldexp( x, (e >> 1) );
/*
*q += ((e >>1) & 0x7ff) << 4;
*q &= 077777;
*/
#endif
/* Newton iterations: */
#ifdef UNK
x = 0.5*(x + w/x);
x = 0.5*(x + w/x);
x = 0.5*(x + w/x);
#endif
/* Note, assume the square root cannot be denormal,
* so it is safe to use integer exponent operations here.
*/
#ifdef DEC
x += w/x;
*q -= 0200;
x += w/x;
*q -= 0200;
x += w/x;
*q -= 0200;
#endif
#ifdef IBMPC
x += w/x;
*q -= 0x10;
x += w/x;
*q -= 0x10;
x += w/x;
*q -= 0x10;
#endif
#ifdef MIEEE
x += w/x;
*q -= 0x10;
x += w/x;
*q -= 0x10;
x += w/x;
*q -= 0x10;
#endif
return(x);
}

328
source/common/thirdparty/math/tan.c vendored Normal file
View file

@ -0,0 +1,328 @@
/* tan.c
*
* Circular tangent
*
*
*
* SYNOPSIS:
*
* double x, y, tan();
*
* y = tan( x );
*
*
*
* DESCRIPTION:
*
* Returns the circular tangent of the radian argument x.
*
* Range reduction is modulo pi/4. A rational function
* x + x**3 P(x**2)/Q(x**2)
* is employed in the basic interval [0, pi/4].
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC +-1.07e9 44000 4.1e-17 1.0e-17
* IEEE +-1.07e9 30000 2.9e-16 8.1e-17
*
* ERROR MESSAGES:
*
* message condition value returned
* tan total loss x > 1.073741824e9 0.0
*
*/
/* cot.c
*
* Circular cotangent
*
*
*
* SYNOPSIS:
*
* double x, y, cot();
*
* y = cot( x );
*
*
*
* DESCRIPTION:
*
* Returns the circular cotangent of the radian argument x.
*
* Range reduction is modulo pi/4. A rational function
* x + x**3 P(x**2)/Q(x**2)
* is employed in the basic interval [0, pi/4].
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE +-1.07e9 30000 2.9e-16 8.2e-17
*
*
* ERROR MESSAGES:
*
* message condition value returned
* cot total loss x > 1.073741824e9 0.0
* cot singularity x = 0 INFINITY
*
*/
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. Neither the name of the <ORGANIZATION> nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
*/
#include "mconf.h"
#ifdef UNK
static double P[] = {
-1.30936939181383777646E4,
1.15351664838587416140E6,
-1.79565251976484877988E7
};
static double Q[] = {
/* 1.00000000000000000000E0,*/
1.36812963470692954678E4,
-1.32089234440210967447E6,
2.50083801823357915839E7,
-5.38695755929454629881E7
};
static double DP1 = 7.853981554508209228515625E-1;
static double DP2 = 7.94662735614792836714E-9;
static double DP3 = 3.06161699786838294307E-17;
static double lossth = 1.073741824e9;
#endif
#ifdef DEC
static unsigned short P[] = {
0143514,0113306,0111171,0174674,
0045214,0147545,0027744,0167346,
0146210,0177526,0114514,0105660
};
static unsigned short Q[] = {
/*0040200,0000000,0000000,0000000,*/
0043525,0142457,0072633,0025617,
0145241,0036742,0140525,0162256,
0046276,0146176,0013526,0143573,
0146515,0077401,0162762,0150607
};
/* 7.853981629014015197753906250000E-1 */
static unsigned short P1[] = {0040111,0007732,0120000,0000000,};
/* 4.960467869796758577649598009884E-10 */
static unsigned short P2[] = {0030410,0055060,0100000,0000000,};
/* 2.860594363054915898381331279295E-18 */
static unsigned short P3[] = {0021523,0011431,0105056,0001560,};
#define DP1 *(double *)P1
#define DP2 *(double *)P2
#define DP3 *(double *)P3
static double lossth = 1.073741824e9;
#endif
#ifdef IBMPC
static unsigned short P[] = {
0x3f38,0xd24f,0x92d8,0xc0c9,
0x9ddd,0xa5fc,0x99ec,0x4131,
0x9176,0xd329,0x1fea,0xc171
};
static unsigned short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0x6572,0xeeb3,0xb8a5,0x40ca,
0xbc96,0x582a,0x27bc,0xc134,
0xd8ef,0xc2ea,0xd98f,0x4177,
0x5a31,0x3cbe,0xafe0,0xc189
};
/*
7.85398125648498535156E-1,
3.77489470793079817668E-8,
2.69515142907905952645E-15,
*/
static unsigned short P1[] = {0x0000,0x4000,0x21fb,0x3fe9};
static unsigned short P2[] = {0x0000,0x0000,0x442d,0x3e64};
static unsigned short P3[] = {0x5170,0x98cc,0x4698,0x3ce8};
#define DP1 *(double *)P1
#define DP2 *(double *)P2
#define DP3 *(double *)P3
static double lossth = 1.073741824e9;
#endif
#ifdef MIEEE
static unsigned short P[] = {
0xc0c9,0x92d8,0xd24f,0x3f38,
0x4131,0x99ec,0xa5fc,0x9ddd,
0xc171,0x1fea,0xd329,0x9176
};
static unsigned short Q[] = {
0x40ca,0xb8a5,0xeeb3,0x6572,
0xc134,0x27bc,0x582a,0xbc96,
0x4177,0xd98f,0xc2ea,0xd8ef,
0xc189,0xafe0,0x3cbe,0x5a31
};
static unsigned short P1[] = {
0x3fe9,0x21fb,0x4000,0x0000
};
static unsigned short P2[] = {
0x3e64,0x442d,0x0000,0x0000
};
static unsigned short P3[] = {
0x3ce8,0x4698,0x98cc,0x5170,
};
#define DP1 *(double *)P1
#define DP2 *(double *)P2
#define DP3 *(double *)P3
static double lossth = 1.073741824e9;
#endif
#ifdef ANSIPROT
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern double floor ( double );
extern double ldexp ( double, int );
extern int isnan ( double );
extern int isfinite ( double );
static double tancot(double, int);
#else
double polevl(), p1evl(), floor(), ldexp();
static double tancot();
int isnan(), isfinite();
#endif
extern double PIO4;
extern double INFINITY;
extern double NAN;
double c_tan(x)
double x;
{
#ifdef MINUSZERO
if( x == 0.0 )
return(x);
#endif
#ifdef NANS
if( isnan(x) )
return(x);
if( !isfinite(x) )
{
mtherr( "tan", DOMAIN );
return(NAN);
}
#endif
return( tancot(x,0) );
}
double c_cot(x)
double x;
{
if( x == 0.0 )
{
mtherr( "cot", SING );
return( INFINITY );
}
return( tancot(x,1) );
}
static double tancot( xx, cotflg )
double xx;
int cotflg;
{
double x, y, z, zz;
int j, sign;
/* make argument positive but save the sign */
if( xx < 0 )
{
x = -xx;
sign = -1;
}
else
{
x = xx;
sign = 1;
}
if( x > lossth )
{
if( cotflg )
mtherr( "cot", TLOSS );
else
mtherr( "tan", TLOSS );
return(0.0);
}
/* compute x mod PIO4 */
y = floor( x/PIO4 );
/* strip high bits of integer part */
z = ldexp( y, -3 );
z = floor(z); /* integer part of y/8 */
z = y - ldexp( z, 3 ); /* y - 16 * (y/16) */
/* integer and fractional part modulo one octant */
j = (int)z;
/* map zeros and singularities to origin */
if( j & 1 )
{
j += 1;
y += 1.0;
}
z = ((x - y * DP1) - y * DP2) - y * DP3;
zz = z * z;
if( zz > 1.0e-14 )
y = z + z * (zz * polevl( zz, P, 2 )/p1evl(zz, Q, 4));
else
y = z;
if( j & 2 )
{
if( cotflg )
y = -y;
else
y = -1.0/y;
}
else
{
if( cotflg )
y = 1.0/y;
}
if( sign < 0 )
y = -y;
return( y );
}

165
source/common/thirdparty/math/tanh.c vendored Normal file
View file

@ -0,0 +1,165 @@
/* tanh.c
*
* Hyperbolic tangent
*
*
*
* SYNOPSIS:
*
* double x, y, tanh();
*
* y = tanh( x );
*
*
*
* DESCRIPTION:
*
* Returns hyperbolic tangent of argument in the range MINLOG to
* MAXLOG.
*
* A rational function is used for |x| < 0.625. The form
* x + x**3 P(x)/Q(x) of Cody _& Waite is employed.
* Otherwise,
* tanh(x) = sinh(x)/cosh(x) = 1 - 2/(exp(2x) + 1).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC -2,2 50000 3.3e-17 6.4e-18
* IEEE -2,2 30000 2.5e-16 5.8e-17
*
*/
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. Neither the name of the <ORGANIZATION> nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
*/
#include "mconf.h"
#ifdef UNK
static double P[] = {
-9.64399179425052238628E-1,
-9.92877231001918586564E1,
-1.61468768441708447952E3
};
static double Q[] = {
/* 1.00000000000000000000E0,*/
1.12811678491632931402E2,
2.23548839060100448583E3,
4.84406305325125486048E3
};
#endif
#ifdef DEC
static unsigned short P[] = {
0140166,0161335,0053753,0075126,
0141706,0111520,0070463,0040552,
0142711,0153001,0101300,0025430
};
static unsigned short Q[] = {
/*0040200,0000000,0000000,0000000,*/
0041741,0117624,0051300,0156060,
0043013,0133720,0071251,0127717,
0043227,0060201,0021020,0020136
};
#endif
#ifdef IBMPC
static unsigned short P[] = {
0x6f4b,0xaafd,0xdc5b,0xbfee,
0x682d,0x0e26,0xd26a,0xc058,
0x0563,0x3058,0x3ac0,0xc099
};
static unsigned short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0x1b86,0x8a58,0x33f2,0x405c,
0x35fa,0x0e55,0x76fa,0x40a1,
0x040c,0x2442,0xec10,0x40b2
};
#endif
#ifdef MIEEE
static unsigned short P[] = {
0xbfee,0xdc5b,0xaafd,0x6f4b,
0xc058,0xd26a,0x0e26,0x682d,
0xc099,0x3ac0,0x3058,0x0563
};
static unsigned short Q[] = {
0x405c,0x33f2,0x8a58,0x1b86,
0x40a1,0x76fa,0x0e55,0x35fa,
0x40b2,0xec10,0x2442,0x040c
};
#endif
#ifdef ANSIPROT
extern double fabs ( double );
extern double c_exp ( double );
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
#else
double fabs(), c_exp(), polevl(), p1evl();
#endif
extern double MAXLOG;
double c_tanh(x)
double x;
{
double s, z;
#ifdef MINUSZERO
if( x == 0.0 )
return(x);
#endif
z = fabs(x);
if( z > 0.5 * MAXLOG )
{
if( x > 0 )
return( 1.0 );
else
return( -1.0 );
}
if( z >= 0.625 )
{
s = c_exp(2.0*z);
z = 1.0 - 2.0/(s + 1.0);
if( x < 0 )
z = -z;
}
else
{
if( x == 0.0 )
return(x);
s = x * x;
z = polevl( s, P, 2 )/p1evl(s, Q, 3);
z = x * s * z;
z = x + z;
}
return( z );
}

236
source/common/thirdparty/md5.cpp vendored Normal file
View file

@ -0,0 +1,236 @@
/*
* This code implements the MD5 message-digest algorithm.
* The algorithm is due to Ron Rivest. This code was
* written by Colin Plumb in 1993, no copyright is claimed.
* This code is in the public domain; do with it what you wish.
*
* Equivalent code is available from RSA Data Security, Inc.
* This code has been tested against that, and is equivalent,
* except that you don't need to include two pages of legalese
* with every copy.
*
* To compute the message digest of a chunk of bytes, declare an
* MD5Context structure, pass it to MD5Init, call MD5Update as
* needed on buffers full of bytes, and then call MD5Final, which
* will fill a supplied 16-byte array with the digest.
*/
#include <string.h> /* for memcpy() */
#include <algorithm>
#include "md5.h"
#ifdef __BIG_ENDIAN__
void byteSwap(uint32_t *buf, unsigned words)
{
uint8_t *p = (uint8_t *)buf;
do {
*buf++ = (uint32_t)((unsigned)p[3] << 8 | p[2]) << 16 |
((unsigned)p[1] << 8 | p[0]);
p += 4;
} while (--words);
}
#else
#define byteSwap(buf,words)
#endif
/*
* Start MD5 accumulation. Set bit count to 0 and buffer to mysterious
* initialization constants.
*/
void MD5Context::Init()
{
buf[0] = 0x67452301;
buf[1] = 0xefcdab89;
buf[2] = 0x98badcfe;
buf[3] = 0x10325476;
bytes[0] = 0;
bytes[1] = 0;
}
/*
* Update context to reflect the concatenation of another buffer full
* of bytes.
*/
void MD5Context::Update(const uint8_t *buf, unsigned len)
{
uint32_t t;
/* Update byte count */
t = bytes[0];
if ((bytes[0] = t + len) < t)
bytes[1]++; /* Carry from low to high */
t = 64 - (t & 0x3f); /* Space available in ctx->in (at least 1) */
if (t > len) {
memcpy((uint8_t *)in + 64 - t, buf, len);
return;
}
/* First chunk is an odd size */
if (t < 64)
{
memcpy((uint8_t *)in + 64 - t, buf, t);
byteSwap(in, 16);
MD5Transform(this->buf, in);
buf += t;
len -= t;
}
/* Process data in 64-byte chunks */
while (len >= 64)
{
memcpy(in, buf, 64);
byteSwap(in, 16);
MD5Transform(this->buf, in);
buf += 64;
len -= 64;
}
/* Handle any remaining bytes of data. */
memcpy(in, buf, len);
}
/*
* Final wrapup - pad to 64-byte boundary with the bit pattern
* 1 0* (64-bit count of bits processed, MSB-first)
*/
void MD5Context::Final(uint8_t digest[16])
{
int count = bytes[0] & 0x3f; /* Number of bytes in ctx->in */
uint8_t *p = (uint8_t *)in + count;
/* Set the first char of padding to 0x80. There is always room. */
*p++ = 0x80;
/* Bytes of padding needed to make 56 bytes (-8..55) */
count = 56 - 1 - count;
if (count < 0) /* Padding forces an extra block */
{
memset(p, 0, count + 8);
byteSwap(in, 16);
MD5Transform(buf, in);
p = (uint8_t *)in;
count = 56;
}
memset(p, 0, count);
byteSwap(in, 14);
/* Append length in bits and transform */
in[14] = bytes[0] << 3;
in[15] = (bytes[1] << 3) | (bytes[0] >> 29);
MD5Transform(buf, in);
byteSwap(buf, 4);
memcpy(digest, buf, 16);
memset(this, 0, sizeof(*this)); /* In case it's sensitive */
}
#ifndef ASM_MD5
/* The four core functions - F1 is optimized somewhat */
/* #define F1(x, y, z) (x & y | ~x & z) */
#define F1(x, y, z) (z ^ (x & (y ^ z)))
#define F2(x, y, z) F1(z, x, y)
#define F3(x, y, z) (x ^ y ^ z)
#define F4(x, y, z) (y ^ (x | ~z))
/* This is the central step in the MD5 algorithm. */
#define MD5STEP(f,w,x,y,z,in,s) \
(w += f(x,y,z) + in, w = (w<<s | w>>(32-s)) + x)
/*
* The core of the MD5 algorithm, this alters an existing MD5 hash to
* reflect the addition of 16 longwords of new data. MD5Update blocks
* the data and converts bytes into longwords for this routine.
*/
void
MD5Transform(uint32_t buf[4], const uint32_t in[16])
{
uint32_t a, b, c, d;
a = buf[0];
b = buf[1];
c = buf[2];
d = buf[3];
MD5STEP(F1, a, b, c, d, in[0] + 0xd76aa478, 7);
MD5STEP(F1, d, a, b, c, in[1] + 0xe8c7b756, 12);
MD5STEP(F1, c, d, a, b, in[2] + 0x242070db, 17);
MD5STEP(F1, b, c, d, a, in[3] + 0xc1bdceee, 22);
MD5STEP(F1, a, b, c, d, in[4] + 0xf57c0faf, 7);
MD5STEP(F1, d, a, b, c, in[5] + 0x4787c62a, 12);
MD5STEP(F1, c, d, a, b, in[6] + 0xa8304613, 17);
MD5STEP(F1, b, c, d, a, in[7] + 0xfd469501, 22);
MD5STEP(F1, a, b, c, d, in[8] + 0x698098d8, 7);
MD5STEP(F1, d, a, b, c, in[9] + 0x8b44f7af, 12);
MD5STEP(F1, c, d, a, b, in[10] + 0xffff5bb1, 17);
MD5STEP(F1, b, c, d, a, in[11] + 0x895cd7be, 22);
MD5STEP(F1, a, b, c, d, in[12] + 0x6b901122, 7);
MD5STEP(F1, d, a, b, c, in[13] + 0xfd987193, 12);
MD5STEP(F1, c, d, a, b, in[14] + 0xa679438e, 17);
MD5STEP(F1, b, c, d, a, in[15] + 0x49b40821, 22);
MD5STEP(F2, a, b, c, d, in[1] + 0xf61e2562, 5);
MD5STEP(F2, d, a, b, c, in[6] + 0xc040b340, 9);
MD5STEP(F2, c, d, a, b, in[11] + 0x265e5a51, 14);
MD5STEP(F2, b, c, d, a, in[0] + 0xe9b6c7aa, 20);
MD5STEP(F2, a, b, c, d, in[5] + 0xd62f105d, 5);
MD5STEP(F2, d, a, b, c, in[10] + 0x02441453, 9);
MD5STEP(F2, c, d, a, b, in[15] + 0xd8a1e681, 14);
MD5STEP(F2, b, c, d, a, in[4] + 0xe7d3fbc8, 20);
MD5STEP(F2, a, b, c, d, in[9] + 0x21e1cde6, 5);
MD5STEP(F2, d, a, b, c, in[14] + 0xc33707d6, 9);
MD5STEP(F2, c, d, a, b, in[3] + 0xf4d50d87, 14);
MD5STEP(F2, b, c, d, a, in[8] + 0x455a14ed, 20);
MD5STEP(F2, a, b, c, d, in[13] + 0xa9e3e905, 5);
MD5STEP(F2, d, a, b, c, in[2] + 0xfcefa3f8, 9);
MD5STEP(F2, c, d, a, b, in[7] + 0x676f02d9, 14);
MD5STEP(F2, b, c, d, a, in[12] + 0x8d2a4c8a, 20);
MD5STEP(F3, a, b, c, d, in[5] + 0xfffa3942, 4);
MD5STEP(F3, d, a, b, c, in[8] + 0x8771f681, 11);
MD5STEP(F3, c, d, a, b, in[11] + 0x6d9d6122, 16);
MD5STEP(F3, b, c, d, a, in[14] + 0xfde5380c, 23);
MD5STEP(F3, a, b, c, d, in[1] + 0xa4beea44, 4);
MD5STEP(F3, d, a, b, c, in[4] + 0x4bdecfa9, 11);
MD5STEP(F3, c, d, a, b, in[7] + 0xf6bb4b60, 16);
MD5STEP(F3, b, c, d, a, in[10] + 0xbebfbc70, 23);
MD5STEP(F3, a, b, c, d, in[13] + 0x289b7ec6, 4);
MD5STEP(F3, d, a, b, c, in[0] + 0xeaa127fa, 11);
MD5STEP(F3, c, d, a, b, in[3] + 0xd4ef3085, 16);
MD5STEP(F3, b, c, d, a, in[6] + 0x04881d05, 23);
MD5STEP(F3, a, b, c, d, in[9] + 0xd9d4d039, 4);
MD5STEP(F3, d, a, b, c, in[12] + 0xe6db99e5, 11);
MD5STEP(F3, c, d, a, b, in[15] + 0x1fa27cf8, 16);
MD5STEP(F3, b, c, d, a, in[2] + 0xc4ac5665, 23);
MD5STEP(F4, a, b, c, d, in[0] + 0xf4292244, 6);
MD5STEP(F4, d, a, b, c, in[7] + 0x432aff97, 10);
MD5STEP(F4, c, d, a, b, in[14] + 0xab9423a7, 15);
MD5STEP(F4, b, c, d, a, in[5] + 0xfc93a039, 21);
MD5STEP(F4, a, b, c, d, in[12] + 0x655b59c3, 6);
MD5STEP(F4, d, a, b, c, in[3] + 0x8f0ccc92, 10);
MD5STEP(F4, c, d, a, b, in[10] + 0xffeff47d, 15);
MD5STEP(F4, b, c, d, a, in[1] + 0x85845dd1, 21);
MD5STEP(F4, a, b, c, d, in[8] + 0x6fa87e4f, 6);
MD5STEP(F4, d, a, b, c, in[15] + 0xfe2ce6e0, 10);
MD5STEP(F4, c, d, a, b, in[6] + 0xa3014314, 15);
MD5STEP(F4, b, c, d, a, in[13] + 0x4e0811a1, 21);
MD5STEP(F4, a, b, c, d, in[4] + 0xf7537e82, 6);
MD5STEP(F4, d, a, b, c, in[11] + 0xbd3af235, 10);
MD5STEP(F4, c, d, a, b, in[2] + 0x2ad7d2bb, 15);
MD5STEP(F4, b, c, d, a, in[9] + 0xeb86d391, 21);
buf[0] += a;
buf[1] += b;
buf[2] += c;
buf[3] += d;
}
#endif

38
source/common/thirdparty/md5.h vendored Normal file
View file

@ -0,0 +1,38 @@
/*
* This is the header file for the MD5 message-digest algorithm.
* The algorithm is due to Ron Rivest. This code was
* written by Colin Plumb in 1993, no copyright is claimed.
* This code is in the public domain; do with it what you wish.
*
* Equivalent code is available from RSA Data Security, Inc.
* This code has been tested against that, and is equivalent,
* except that you don't need to include two pages of legalese
* with every copy.
*
* To compute the message digest of a chunk of bytes, declare an
* MD5Context structure, pass it to MD5Init, call MD5Update as
* needed on buffers full of bytes, and then call MD5Final, which
* will fill a supplied 16-byte array with the digest.
*/
#ifndef MD5_H
#define MD5_H
struct MD5Context
{
MD5Context() { Init(); }
void Init();
void Update(const uint8_t *buf, unsigned len);
void Final(uint8_t digest[16]);
private:
uint32_t buf[4];
uint32_t bytes[2];
uint32_t in[16];
};
void MD5Transform(uint32_t buf[4], uint32_t const in[16]);
#endif /* !MD5_H */

View file

@ -0,0 +1,29 @@
Copyright (c) 2006,2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima
University. All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above
copyright notice, this list of conditions and the following
disclaimer in the documentation and/or other materials provided
with the distribution.
* Neither the name of the Hiroshima University nor the names of
its contributors may be used to endorse or promote products
derived from this software without specific prior written
permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

View file

@ -0,0 +1,156 @@
/**
* @file SFMT-alti.h
*
* @brief SIMD oriented Fast Mersenne Twister(SFMT)
* pseudorandom number generator
*
* @author Mutsuo Saito (Hiroshima University)
* @author Makoto Matsumoto (Hiroshima University)
*
* Copyright (C) 2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima
* University. All rights reserved.
*
* The new BSD License is applied to this software.
* see LICENSE.txt
*/
#ifndef SFMT_ALTI_H
#define SFMT_ALTI_H
inline static vector unsigned int vec_recursion(vector unsigned int a,
vector unsigned int b,
vector unsigned int c,
vector unsigned int d)
ALWAYSINLINE;
/**
* This function represents the recursion formula in AltiVec and BIG ENDIAN.
* @param a a 128-bit part of the interal state array
* @param b a 128-bit part of the interal state array
* @param c a 128-bit part of the interal state array
* @param d a 128-bit part of the interal state array
* @return output
*/
inline static vector unsigned int vec_recursion(vector unsigned int a,
vector unsigned int b,
vector unsigned int c,
vector unsigned int d) {
const vector unsigned int sl1 = ALTI_SL1;
const vector unsigned int sr1 = ALTI_SR1;
#ifdef ONLY64
const vector unsigned int mask = ALTI_MSK64;
const vector unsigned char perm_sl = ALTI_SL2_PERM64;
const vector unsigned char perm_sr = ALTI_SR2_PERM64;
#else
const vector unsigned int mask = ALTI_MSK;
const vector unsigned char perm_sl = ALTI_SL2_PERM;
const vector unsigned char perm_sr = ALTI_SR2_PERM;
#endif
vector unsigned int v, w, x, y, z;
x = vec_perm(a, (vector unsigned int)perm_sl, perm_sl);
v = a;
y = vec_sr(b, sr1);
z = vec_perm(c, (vector unsigned int)perm_sr, perm_sr);
w = vec_sl(d, sl1);
z = vec_xor(z, w);
y = vec_and(y, mask);
v = vec_xor(v, x);
z = vec_xor(z, y);
z = vec_xor(z, v);
return z;
}
/**
* This function fills the internal state array with pseudorandom
* integers.
*/
inline static void gen_rand_all(void) {
int i;
vector unsigned int r, r1, r2;
r1 = sfmt[N - 2].s;
r2 = sfmt[N - 1].s;
for (i = 0; i < N - POS1; i++) {
r = vec_recursion(sfmt[i].s, sfmt[i + POS1].s, r1, r2);
sfmt[i].s = r;
r1 = r2;
r2 = r;
}
for (; i < N; i++) {
r = vec_recursion(sfmt[i].s, sfmt[i + POS1 - N].s, r1, r2);
sfmt[i].s = r;
r1 = r2;
r2 = r;
}
}
/**
* This function fills the user-specified array with pseudorandom
* integers.
*
* @param array an 128-bit array to be filled by pseudorandom numbers.
* @param size number of 128-bit pesudorandom numbers to be generated.
*/
inline static void gen_rand_array(w128_t *array, int size) {
int i, j;
vector unsigned int r, r1, r2;
r1 = sfmt[N - 2].s;
r2 = sfmt[N - 1].s;
for (i = 0; i < N - POS1; i++) {
r = vec_recursion(sfmt[i].s, sfmt[i + POS1].s, r1, r2);
array[i].s = r;
r1 = r2;
r2 = r;
}
for (; i < N; i++) {
r = vec_recursion(sfmt[i].s, array[i + POS1 - N].s, r1, r2);
array[i].s = r;
r1 = r2;
r2 = r;
}
/* main loop */
for (; i < size - N; i++) {
r = vec_recursion(array[i - N].s, array[i + POS1 - N].s, r1, r2);
array[i].s = r;
r1 = r2;
r2 = r;
}
for (j = 0; j < 2 * N - size; j++) {
sfmt[j].s = array[j + size - N].s;
}
for (; i < size; i++) {
r = vec_recursion(array[i - N].s, array[i + POS1 - N].s, r1, r2);
array[i].s = r;
sfmt[j++].s = r;
r1 = r2;
r2 = r;
}
}
#ifndef ONLY64
#if defined(__APPLE__)
#define ALTI_SWAP (vector unsigned char) \
(4, 5, 6, 7, 0, 1, 2, 3, 12, 13, 14, 15, 8, 9, 10, 11)
#else
#define ALTI_SWAP {4, 5, 6, 7, 0, 1, 2, 3, 12, 13, 14, 15, 8, 9, 10, 11}
#endif
/**
* This function swaps high and low 32-bit of 64-bit integers in user
* specified array.
*
* @param array an 128-bit array to be swaped.
* @param size size of 128-bit array.
*/
inline static void swap(w128_t *array, int size) {
int i;
const vector unsigned char perm = ALTI_SWAP;
for (i = 0; i < size; i++) {
array[i].s = vec_perm(array[i].s, (vector unsigned int)perm, perm);
}
}
#endif
#endif

View file

@ -0,0 +1,70 @@
#ifndef SFMT_PARAMS_H
#define SFMT_PARAMS_H
/*----------------------
the parameters of SFMT
following definitions are in paramsXXXX.h file.
----------------------*/
/** the pick up position of the array.
#define POS1 122
*/
/** the parameter of shift left as four 32-bit registers.
#define SL1 18
*/
/** the parameter of shift left as one 128-bit register.
* The 128-bit integer is shifted by (SL2 * 8) bits.
#define SL2 1
*/
/** the parameter of shift right as four 32-bit registers.
#define SR1 11
*/
/** the parameter of shift right as one 128-bit register.
* The 128-bit integer is shifted by (SL2 * 8) bits.
#define SR2 1
*/
/** A bitmask, used in the recursion. These parameters are introduced
* to break symmetry of SIMD.
#define MSK1 0xdfffffefU
#define MSK2 0xddfecb7fU
#define MSK3 0xbffaffffU
#define MSK4 0xbffffff6U
*/
/** These definitions are part of a 128-bit period certification vector.
#define PARITY1 0x00000001U
#define PARITY2 0x00000000U
#define PARITY3 0x00000000U
#define PARITY4 0xc98e126aU
*/
#if MEXP == 607
#include "SFMT-params607.h"
#elif MEXP == 1279
#include "SFMT-params1279.h"
#elif MEXP == 2281
#include "SFMT-params2281.h"
#elif MEXP == 4253
#include "SFMT-params4253.h"
#elif MEXP == 11213
#include "SFMT-params11213.h"
#elif MEXP == 19937
#include "SFMT-params19937.h"
#elif MEXP == 44497
#include "SFMT-params44497.h"
#elif MEXP == 86243
#include "SFMT-params86243.h"
#elif MEXP == 132049
#include "SFMT-params132049.h"
#elif MEXP == 216091
#include "SFMT-params216091.h"
#else
#error "MEXP is not valid."
#undef MEXP
#endif
#endif /* SFMT_PARAMS_H */

View file

@ -0,0 +1,46 @@
#ifndef SFMT_PARAMS11213_H
#define SFMT_PARAMS11213_H
#define POS1 68
#define SL1 14
#define SL2 3
#define SR1 7
#define SR2 3
#define MSK1 0xeffff7fbU
#define MSK2 0xffffffefU
#define MSK3 0xdfdfbfffU
#define MSK4 0x7fffdbfdU
#define PARITY1 0x00000001U
#define PARITY2 0x00000000U
#define PARITY3 0xe8148000U
#define PARITY4 0xd0c7afa3U
/* PARAMETERS FOR ALTIVEC */
#if defined(__APPLE__) /* For OSX */
#define ALTI_SL1 (vector unsigned int)(SL1, SL1, SL1, SL1)
#define ALTI_SR1 (vector unsigned int)(SR1, SR1, SR1, SR1)
#define ALTI_MSK (vector unsigned int)(MSK1, MSK2, MSK3, MSK4)
#define ALTI_MSK64 \
(vector unsigned int)(MSK2, MSK1, MSK4, MSK3)
#define ALTI_SL2_PERM \
(vector unsigned char)(3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10)
#define ALTI_SL2_PERM64 \
(vector unsigned char)(3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2)
#define ALTI_SR2_PERM \
(vector unsigned char)(5,6,7,0,9,10,11,4,13,14,15,8,19,19,19,12)
#define ALTI_SR2_PERM64 \
(vector unsigned char)(13,14,15,0,1,2,3,4,19,19,19,8,9,10,11,12)
#else /* For OTHER OSs(Linux?) */
#define ALTI_SL1 {SL1, SL1, SL1, SL1}
#define ALTI_SR1 {SR1, SR1, SR1, SR1}
#define ALTI_MSK {MSK1, MSK2, MSK3, MSK4}
#define ALTI_MSK64 {MSK2, MSK1, MSK4, MSK3}
#define ALTI_SL2_PERM {3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10}
#define ALTI_SL2_PERM64 {3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2}
#define ALTI_SR2_PERM {5,6,7,0,9,10,11,4,13,14,15,8,19,19,19,12}
#define ALTI_SR2_PERM64 {13,14,15,0,1,2,3,4,19,19,19,8,9,10,11,12}
#endif /* For OSX */
#define IDSTR "SFMT-11213:68-14-3-7-3:effff7fb-ffffffef-dfdfbfff-7fffdbfd"
#endif /* SFMT_PARAMS11213_H */

View file

@ -0,0 +1,46 @@
#ifndef SFMT_PARAMS1279_H
#define SFMT_PARAMS1279_H
#define POS1 7
#define SL1 14
#define SL2 3
#define SR1 5
#define SR2 1
#define MSK1 0xf7fefffdU
#define MSK2 0x7fefcfffU
#define MSK3 0xaff3ef3fU
#define MSK4 0xb5ffff7fU
#define PARITY1 0x00000001U
#define PARITY2 0x00000000U
#define PARITY3 0x00000000U
#define PARITY4 0x20000000U
/* PARAMETERS FOR ALTIVEC */
#if defined(__APPLE__) /* For OSX */
#define ALTI_SL1 (vector unsigned int)(SL1, SL1, SL1, SL1)
#define ALTI_SR1 (vector unsigned int)(SR1, SR1, SR1, SR1)
#define ALTI_MSK (vector unsigned int)(MSK1, MSK2, MSK3, MSK4)
#define ALTI_MSK64 \
(vector unsigned int)(MSK2, MSK1, MSK4, MSK3)
#define ALTI_SL2_PERM \
(vector unsigned char)(3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10)
#define ALTI_SL2_PERM64 \
(vector unsigned char)(3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2)
#define ALTI_SR2_PERM \
(vector unsigned char)(7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14)
#define ALTI_SR2_PERM64 \
(vector unsigned char)(15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14)
#else /* For OTHER OSs(Linux?) */
#define ALTI_SL1 {SL1, SL1, SL1, SL1}
#define ALTI_SR1 {SR1, SR1, SR1, SR1}
#define ALTI_MSK {MSK1, MSK2, MSK3, MSK4}
#define ALTI_MSK64 {MSK2, MSK1, MSK4, MSK3}
#define ALTI_SL2_PERM {3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10}
#define ALTI_SL2_PERM64 {3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2}
#define ALTI_SR2_PERM {7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14}
#define ALTI_SR2_PERM64 {15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14}
#endif /* For OSX */
#define IDSTR "SFMT-1279:7-14-3-5-1:f7fefffd-7fefcfff-aff3ef3f-b5ffff7f"
#endif /* SFMT_PARAMS1279_H */

View file

@ -0,0 +1,46 @@
#ifndef SFMT_PARAMS132049_H
#define SFMT_PARAMS132049_H
#define POS1 110
#define SL1 19
#define SL2 1
#define SR1 21
#define SR2 1
#define MSK1 0xffffbb5fU
#define MSK2 0xfb6ebf95U
#define MSK3 0xfffefffaU
#define MSK4 0xcff77fffU
#define PARITY1 0x00000001U
#define PARITY2 0x00000000U
#define PARITY3 0xcb520000U
#define PARITY4 0xc7e91c7dU
/* PARAMETERS FOR ALTIVEC */
#if defined(__APPLE__) /* For OSX */
#define ALTI_SL1 (vector unsigned int)(SL1, SL1, SL1, SL1)
#define ALTI_SR1 (vector unsigned int)(SR1, SR1, SR1, SR1)
#define ALTI_MSK (vector unsigned int)(MSK1, MSK2, MSK3, MSK4)
#define ALTI_MSK64 \
(vector unsigned int)(MSK2, MSK1, MSK4, MSK3)
#define ALTI_SL2_PERM \
(vector unsigned char)(1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8)
#define ALTI_SL2_PERM64 \
(vector unsigned char)(1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0)
#define ALTI_SR2_PERM \
(vector unsigned char)(7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14)
#define ALTI_SR2_PERM64 \
(vector unsigned char)(15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14)
#else /* For OTHER OSs(Linux?) */
#define ALTI_SL1 {SL1, SL1, SL1, SL1}
#define ALTI_SR1 {SR1, SR1, SR1, SR1}
#define ALTI_MSK {MSK1, MSK2, MSK3, MSK4}
#define ALTI_MSK64 {MSK2, MSK1, MSK4, MSK3}
#define ALTI_SL2_PERM {1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8}
#define ALTI_SL2_PERM64 {1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0}
#define ALTI_SR2_PERM {7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14}
#define ALTI_SR2_PERM64 {15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14}
#endif /* For OSX */
#define IDSTR "SFMT-132049:110-19-1-21-1:ffffbb5f-fb6ebf95-fffefffa-cff77fff"
#endif /* SFMT_PARAMS132049_H */

View file

@ -0,0 +1,46 @@
#ifndef SFMT_PARAMS19937_H
#define SFMT_PARAMS19937_H
#define POS1 122
#define SL1 18
#define SL2 1
#define SR1 11
#define SR2 1
#define MSK1 0xdfffffefU
#define MSK2 0xddfecb7fU
#define MSK3 0xbffaffffU
#define MSK4 0xbffffff6U
#define PARITY1 0x00000001U
#define PARITY2 0x00000000U
#define PARITY3 0x00000000U
#define PARITY4 0x13c9e684U
/* PARAMETERS FOR ALTIVEC */
#if defined(__APPLE__) /* For OSX */
#define ALTI_SL1 (vector unsigned int)(SL1, SL1, SL1, SL1)
#define ALTI_SR1 (vector unsigned int)(SR1, SR1, SR1, SR1)
#define ALTI_MSK (vector unsigned int)(MSK1, MSK2, MSK3, MSK4)
#define ALTI_MSK64 \
(vector unsigned int)(MSK2, MSK1, MSK4, MSK3)
#define ALTI_SL2_PERM \
(vector unsigned char)(1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8)
#define ALTI_SL2_PERM64 \
(vector unsigned char)(1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0)
#define ALTI_SR2_PERM \
(vector unsigned char)(7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14)
#define ALTI_SR2_PERM64 \
(vector unsigned char)(15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14)
#else /* For OTHER OSs(Linux?) */
#define ALTI_SL1 {SL1, SL1, SL1, SL1}
#define ALTI_SR1 {SR1, SR1, SR1, SR1}
#define ALTI_MSK {MSK1, MSK2, MSK3, MSK4}
#define ALTI_MSK64 {MSK2, MSK1, MSK4, MSK3}
#define ALTI_SL2_PERM {1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8}
#define ALTI_SL2_PERM64 {1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0}
#define ALTI_SR2_PERM {7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14}
#define ALTI_SR2_PERM64 {15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14}
#endif /* For OSX */
#define IDSTR "SFMT-19937:122-18-1-11-1:dfffffef-ddfecb7f-bffaffff-bffffff6"
#endif /* SFMT_PARAMS19937_H */

View file

@ -0,0 +1,46 @@
#ifndef SFMT_PARAMS216091_H
#define SFMT_PARAMS216091_H
#define POS1 627
#define SL1 11
#define SL2 3
#define SR1 10
#define SR2 1
#define MSK1 0xbff7bff7U
#define MSK2 0xbfffffffU
#define MSK3 0xbffffa7fU
#define MSK4 0xffddfbfbU
#define PARITY1 0xf8000001U
#define PARITY2 0x89e80709U
#define PARITY3 0x3bd2b64bU
#define PARITY4 0x0c64b1e4U
/* PARAMETERS FOR ALTIVEC */
#if defined(__APPLE__) /* For OSX */
#define ALTI_SL1 (vector unsigned int)(SL1, SL1, SL1, SL1)
#define ALTI_SR1 (vector unsigned int)(SR1, SR1, SR1, SR1)
#define ALTI_MSK (vector unsigned int)(MSK1, MSK2, MSK3, MSK4)
#define ALTI_MSK64 \
(vector unsigned int)(MSK2, MSK1, MSK4, MSK3)
#define ALTI_SL2_PERM \
(vector unsigned char)(3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10)
#define ALTI_SL2_PERM64 \
(vector unsigned char)(3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2)
#define ALTI_SR2_PERM \
(vector unsigned char)(7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14)
#define ALTI_SR2_PERM64 \
(vector unsigned char)(15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14)
#else /* For OTHER OSs(Linux?) */
#define ALTI_SL1 {SL1, SL1, SL1, SL1}
#define ALTI_SR1 {SR1, SR1, SR1, SR1}
#define ALTI_MSK {MSK1, MSK2, MSK3, MSK4}
#define ALTI_MSK64 {MSK2, MSK1, MSK4, MSK3}
#define ALTI_SL2_PERM {3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10}
#define ALTI_SL2_PERM64 {3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2}
#define ALTI_SR2_PERM {7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14}
#define ALTI_SR2_PERM64 {15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14}
#endif /* For OSX */
#define IDSTR "SFMT-216091:627-11-3-10-1:bff7bff7-bfffffff-bffffa7f-ffddfbfb"
#endif /* SFMT_PARAMS216091_H */

View file

@ -0,0 +1,46 @@
#ifndef SFMT_PARAMS2281_H
#define SFMT_PARAMS2281_H
#define POS1 12
#define SL1 19
#define SL2 1
#define SR1 5
#define SR2 1
#define MSK1 0xbff7ffbfU
#define MSK2 0xfdfffffeU
#define MSK3 0xf7ffef7fU
#define MSK4 0xf2f7cbbfU
#define PARITY1 0x00000001U
#define PARITY2 0x00000000U
#define PARITY3 0x00000000U
#define PARITY4 0x41dfa600U
/* PARAMETERS FOR ALTIVEC */
#if defined(__APPLE__) /* For OSX */
#define ALTI_SL1 (vector unsigned int)(SL1, SL1, SL1, SL1)
#define ALTI_SR1 (vector unsigned int)(SR1, SR1, SR1, SR1)
#define ALTI_MSK (vector unsigned int)(MSK1, MSK2, MSK3, MSK4)
#define ALTI_MSK64 \
(vector unsigned int)(MSK2, MSK1, MSK4, MSK3)
#define ALTI_SL2_PERM \
(vector unsigned char)(1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8)
#define ALTI_SL2_PERM64 \
(vector unsigned char)(1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0)
#define ALTI_SR2_PERM \
(vector unsigned char)(7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14)
#define ALTI_SR2_PERM64 \
(vector unsigned char)(15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14)
#else /* For OTHER OSs(Linux?) */
#define ALTI_SL1 {SL1, SL1, SL1, SL1}
#define ALTI_SR1 {SR1, SR1, SR1, SR1}
#define ALTI_MSK {MSK1, MSK2, MSK3, MSK4}
#define ALTI_MSK64 {MSK2, MSK1, MSK4, MSK3}
#define ALTI_SL2_PERM {1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8}
#define ALTI_SL2_PERM64 {1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0}
#define ALTI_SR2_PERM {7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14}
#define ALTI_SR2_PERM64 {15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14}
#endif /* For OSX */
#define IDSTR "SFMT-2281:12-19-1-5-1:bff7ffbf-fdfffffe-f7ffef7f-f2f7cbbf"
#endif /* SFMT_PARAMS2281_H */

View file

@ -0,0 +1,46 @@
#ifndef SFMT_PARAMS4253_H
#define SFMT_PARAMS4253_H
#define POS1 17
#define SL1 20
#define SL2 1
#define SR1 7
#define SR2 1
#define MSK1 0x9f7bffffU
#define MSK2 0x9fffff5fU
#define MSK3 0x3efffffbU
#define MSK4 0xfffff7bbU
#define PARITY1 0xa8000001U
#define PARITY2 0xaf5390a3U
#define PARITY3 0xb740b3f8U
#define PARITY4 0x6c11486dU
/* PARAMETERS FOR ALTIVEC */
#if defined(__APPLE__) /* For OSX */
#define ALTI_SL1 (vector unsigned int)(SL1, SL1, SL1, SL1)
#define ALTI_SR1 (vector unsigned int)(SR1, SR1, SR1, SR1)
#define ALTI_MSK (vector unsigned int)(MSK1, MSK2, MSK3, MSK4)
#define ALTI_MSK64 \
(vector unsigned int)(MSK2, MSK1, MSK4, MSK3)
#define ALTI_SL2_PERM \
(vector unsigned char)(1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8)
#define ALTI_SL2_PERM64 \
(vector unsigned char)(1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0)
#define ALTI_SR2_PERM \
(vector unsigned char)(7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14)
#define ALTI_SR2_PERM64 \
(vector unsigned char)(15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14)
#else /* For OTHER OSs(Linux?) */
#define ALTI_SL1 {SL1, SL1, SL1, SL1}
#define ALTI_SR1 {SR1, SR1, SR1, SR1}
#define ALTI_MSK {MSK1, MSK2, MSK3, MSK4}
#define ALTI_MSK64 {MSK2, MSK1, MSK4, MSK3}
#define ALTI_SL2_PERM {1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8}
#define ALTI_SL2_PERM64 {1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0}
#define ALTI_SR2_PERM {7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14}
#define ALTI_SR2_PERM64 {15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14}
#endif /* For OSX */
#define IDSTR "SFMT-4253:17-20-1-7-1:9f7bffff-9fffff5f-3efffffb-fffff7bb"
#endif /* SFMT_PARAMS4253_H */

View file

@ -0,0 +1,46 @@
#ifndef SFMT_PARAMS44497_H
#define SFMT_PARAMS44497_H
#define POS1 330
#define SL1 5
#define SL2 3
#define SR1 9
#define SR2 3
#define MSK1 0xeffffffbU
#define MSK2 0xdfbebfffU
#define MSK3 0xbfbf7befU
#define MSK4 0x9ffd7bffU
#define PARITY1 0x00000001U
#define PARITY2 0x00000000U
#define PARITY3 0xa3ac4000U
#define PARITY4 0xecc1327aU
/* PARAMETERS FOR ALTIVEC */
#if defined(__APPLE__) /* For OSX */
#define ALTI_SL1 (vector unsigned int)(SL1, SL1, SL1, SL1)
#define ALTI_SR1 (vector unsigned int)(SR1, SR1, SR1, SR1)
#define ALTI_MSK (vector unsigned int)(MSK1, MSK2, MSK3, MSK4)
#define ALTI_MSK64 \
(vector unsigned int)(MSK2, MSK1, MSK4, MSK3)
#define ALTI_SL2_PERM \
(vector unsigned char)(3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10)
#define ALTI_SL2_PERM64 \
(vector unsigned char)(3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2)
#define ALTI_SR2_PERM \
(vector unsigned char)(5,6,7,0,9,10,11,4,13,14,15,8,19,19,19,12)
#define ALTI_SR2_PERM64 \
(vector unsigned char)(13,14,15,0,1,2,3,4,19,19,19,8,9,10,11,12)
#else /* For OTHER OSs(Linux?) */
#define ALTI_SL1 {SL1, SL1, SL1, SL1}
#define ALTI_SR1 {SR1, SR1, SR1, SR1}
#define ALTI_MSK {MSK1, MSK2, MSK3, MSK4}
#define ALTI_MSK64 {MSK2, MSK1, MSK4, MSK3}
#define ALTI_SL2_PERM {3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10}
#define ALTI_SL2_PERM64 {3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2}
#define ALTI_SR2_PERM {5,6,7,0,9,10,11,4,13,14,15,8,19,19,19,12}
#define ALTI_SR2_PERM64 {13,14,15,0,1,2,3,4,19,19,19,8,9,10,11,12}
#endif /* For OSX */
#define IDSTR "SFMT-44497:330-5-3-9-3:effffffb-dfbebfff-bfbf7bef-9ffd7bff"
#endif /* SFMT_PARAMS44497_H */

View file

@ -0,0 +1,46 @@
#ifndef SFMT_PARAMS607_H
#define SFMT_PARAMS607_H
#define POS1 2
#define SL1 15
#define SL2 3
#define SR1 13
#define SR2 3
#define MSK1 0xfdff37ffU
#define MSK2 0xef7f3f7dU
#define MSK3 0xff777b7dU
#define MSK4 0x7ff7fb2fU
#define PARITY1 0x00000001U
#define PARITY2 0x00000000U
#define PARITY3 0x00000000U
#define PARITY4 0x5986f054U
/* PARAMETERS FOR ALTIVEC */
#if defined(__APPLE__) /* For OSX */
#define ALTI_SL1 (vector unsigned int)(SL1, SL1, SL1, SL1)
#define ALTI_SR1 (vector unsigned int)(SR1, SR1, SR1, SR1)
#define ALTI_MSK (vector unsigned int)(MSK1, MSK2, MSK3, MSK4)
#define ALTI_MSK64 \
(vector unsigned int)(MSK2, MSK1, MSK4, MSK3)
#define ALTI_SL2_PERM \
(vector unsigned char)(3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10)
#define ALTI_SL2_PERM64 \
(vector unsigned char)(3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2)
#define ALTI_SR2_PERM \
(vector unsigned char)(5,6,7,0,9,10,11,4,13,14,15,8,19,19,19,12)
#define ALTI_SR2_PERM64 \
(vector unsigned char)(13,14,15,0,1,2,3,4,19,19,19,8,9,10,11,12)
#else /* For OTHER OSs(Linux?) */
#define ALTI_SL1 {SL1, SL1, SL1, SL1}
#define ALTI_SR1 {SR1, SR1, SR1, SR1}
#define ALTI_MSK {MSK1, MSK2, MSK3, MSK4}
#define ALTI_MSK64 {MSK2, MSK1, MSK4, MSK3}
#define ALTI_SL2_PERM {3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10}
#define ALTI_SL2_PERM64 {3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2}
#define ALTI_SR2_PERM {5,6,7,0,9,10,11,4,13,14,15,8,19,19,19,12}
#define ALTI_SR2_PERM64 {13,14,15,0,1,2,3,4,19,19,19,8,9,10,11,12}
#endif /* For OSX */
#define IDSTR "SFMT-607:2-15-3-13-3:fdff37ff-ef7f3f7d-ff777b7d-7ff7fb2f"
#endif /* SFMT_PARAMS607_H */

View file

@ -0,0 +1,46 @@
#ifndef SFMT_PARAMS86243_H
#define SFMT_PARAMS86243_H
#define POS1 366
#define SL1 6
#define SL2 7
#define SR1 19
#define SR2 1
#define MSK1 0xfdbffbffU
#define MSK2 0xbff7ff3fU
#define MSK3 0xfd77efffU
#define MSK4 0xbf9ff3ffU
#define PARITY1 0x00000001U
#define PARITY2 0x00000000U
#define PARITY3 0x00000000U
#define PARITY4 0xe9528d85U
/* PARAMETERS FOR ALTIVEC */
#if defined(__APPLE__) /* For OSX */
#define ALTI_SL1 (vector unsigned int)(SL1, SL1, SL1, SL1)
#define ALTI_SR1 (vector unsigned int)(SR1, SR1, SR1, SR1)
#define ALTI_MSK (vector unsigned int)(MSK1, MSK2, MSK3, MSK4)
#define ALTI_MSK64 \
(vector unsigned int)(MSK2, MSK1, MSK4, MSK3)
#define ALTI_SL2_PERM \
(vector unsigned char)(25,25,25,25,3,25,25,25,7,0,1,2,11,4,5,6)
#define ALTI_SL2_PERM64 \
(vector unsigned char)(7,25,25,25,25,25,25,25,15,0,1,2,3,4,5,6)
#define ALTI_SR2_PERM \
(vector unsigned char)(7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14)
#define ALTI_SR2_PERM64 \
(vector unsigned char)(15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14)
#else /* For OTHER OSs(Linux?) */
#define ALTI_SL1 {SL1, SL1, SL1, SL1}
#define ALTI_SR1 {SR1, SR1, SR1, SR1}
#define ALTI_MSK {MSK1, MSK2, MSK3, MSK4}
#define ALTI_MSK64 {MSK2, MSK1, MSK4, MSK3}
#define ALTI_SL2_PERM {25,25,25,25,3,25,25,25,7,0,1,2,11,4,5,6}
#define ALTI_SL2_PERM64 {7,25,25,25,25,25,25,25,15,0,1,2,3,4,5,6}
#define ALTI_SR2_PERM {7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14}
#define ALTI_SR2_PERM64 {15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14}
#endif /* For OSX */
#define IDSTR "SFMT-86243:366-6-7-19-1:fdbffbff-bff7ff3f-fd77efff-bf9ff3ff"
#endif /* SFMT_PARAMS86243_H */

View file

@ -0,0 +1,121 @@
/**
* @file SFMT-sse2.h
* @brief SIMD oriented Fast Mersenne Twister(SFMT) for Intel SSE2
*
* @author Mutsuo Saito (Hiroshima University)
* @author Makoto Matsumoto (Hiroshima University)
*
* @note We assume LITTLE ENDIAN in this file
*
* Copyright (C) 2006, 2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima
* University. All rights reserved.
*
* The new BSD License is applied to this software, see LICENSE.txt
*/
#ifndef SFMT_SSE2_H
#define SFMT_SSE2_H
PRE_ALWAYS static __m128i mm_recursion(__m128i *a, __m128i *b, __m128i c,
__m128i d, __m128i mask) ALWAYSINLINE;
/**
* This function represents the recursion formula.
* @param a a 128-bit part of the interal state array
* @param b a 128-bit part of the interal state array
* @param c a 128-bit part of the interal state array
* @param d a 128-bit part of the interal state array
* @param mask 128-bit mask
* @return output
*/
PRE_ALWAYS static __m128i mm_recursion(__m128i *a, __m128i *b,
__m128i c, __m128i d, __m128i mask) {
__m128i v, x, y, z;
x = _mm_load_si128(a);
y = _mm_srli_epi32(*b, SR1);
z = _mm_srli_si128(c, SR2);
v = _mm_slli_epi32(d, SL1);
z = _mm_xor_si128(z, x);
z = _mm_xor_si128(z, v);
x = _mm_slli_si128(x, SL2);
y = _mm_and_si128(y, mask);
z = _mm_xor_si128(z, x);
z = _mm_xor_si128(z, y);
return z;
}
/**
* This function fills the internal state array with pseudorandom
* integers.
*/
inline static void gen_rand_all(void) {
int i;
__m128i r, r1, r2, mask;
mask = _mm_set_epi32(MSK4, MSK3, MSK2, MSK1);
r1 = _mm_load_si128(&sfmt[N - 2].si);
r2 = _mm_load_si128(&sfmt[N - 1].si);
for (i = 0; i < N - POS1; i++) {
r = mm_recursion(&sfmt[i].si, &sfmt[i + POS1].si, r1, r2, mask);
_mm_store_si128(&sfmt[i].si, r);
r1 = r2;
r2 = r;
}
for (; i < N; i++) {
r = mm_recursion(&sfmt[i].si, &sfmt[i + POS1 - N].si, r1, r2, mask);
_mm_store_si128(&sfmt[i].si, r);
r1 = r2;
r2 = r;
}
}
/**
* This function fills the user-specified array with pseudorandom
* integers.
*
* @param array an 128-bit array to be filled by pseudorandom numbers.
* @param size number of 128-bit pesudorandom numbers to be generated.
*/
inline static void gen_rand_array(w128_t *array, int size) {
int i, j;
__m128i r, r1, r2, mask;
mask = _mm_set_epi32(MSK4, MSK3, MSK2, MSK1);
r1 = _mm_load_si128(&sfmt[N - 2].si);
r2 = _mm_load_si128(&sfmt[N - 1].si);
for (i = 0; i < N - POS1; i++) {
r = mm_recursion(&sfmt[i].si, &sfmt[i + POS1].si, r1, r2, mask);
_mm_store_si128(&array[i].si, r);
r1 = r2;
r2 = r;
}
for (; i < N; i++) {
r = mm_recursion(&sfmt[i].si, &array[i + POS1 - N].si, r1, r2, mask);
_mm_store_si128(&array[i].si, r);
r1 = r2;
r2 = r;
}
/* main loop */
for (; i < size - N; i++) {
r = mm_recursion(&array[i - N].si, &array[i + POS1 - N].si, r1, r2,
mask);
_mm_store_si128(&array[i].si, r);
r1 = r2;
r2 = r;
}
for (j = 0; j < 2 * N - size; j++) {
r = _mm_load_si128(&array[j + size - N].si);
_mm_store_si128(&sfmt[j].si, r);
}
for (; i < size; i++) {
r = mm_recursion(&array[i - N].si, &array[i + POS1 - N].si, r1, r2,
mask);
_mm_store_si128(&array[i].si, r);
_mm_store_si128(&sfmt[j++].si, r);
r1 = r2;
r2 = r;
}
}
#endif

583
source/common/thirdparty/sfmt/SFMT.cpp vendored Normal file
View file

@ -0,0 +1,583 @@
/**
* @file SFMT.c
* @brief SIMD oriented Fast Mersenne Twister(SFMT)
*
* @author Mutsuo Saito (Hiroshima University)
* @author Makoto Matsumoto (Hiroshima University)
*
* Copyright (C) 2006,2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima
* University. All rights reserved.
*
* The new BSD License is applied to this software, see LICENSE.txt
*/
#include <string.h>
#include <assert.h>
#include "SFMTObj.h"
#include "SFMT-params.h"
#if defined(__BIG_ENDIAN__) && !defined(__amd64) && !defined(BIG_ENDIAN64)
#define BIG_ENDIAN64 1
#endif
#if defined(HAVE_ALTIVEC) && !defined(BIG_ENDIAN64)
#define BIG_ENDIAN64 1
#endif
#if defined(ONLY64) && !defined(BIG_ENDIAN64)
#if defined(__GNUC__)
#error "-DONLY64 must be specified with -DBIG_ENDIAN64"
#endif
#undef ONLY64
#endif
/** a parity check vector which certificate the period of 2^{MEXP} */
static const uint32_t parity[4] = { PARITY1, PARITY2, PARITY3, PARITY4 };
/*----------------
STATIC FUNCTIONS
----------------*/
inline static int idxof(int i);
inline static void rshift128(w128_t *out, w128_t const *in, int shift);
inline static void lshift128(w128_t *out, w128_t const *in, int shift);
inline static uint32_t func1(uint32_t x);
inline static uint32_t func2(uint32_t x);
#if defined(BIG_ENDIAN64) && !defined(ONLY64)
inline static void swap(w128_t *array, int size);
#endif
// These SIMD versions WILL NOT work as-is. I'm not even sure SSE2 is
// safe to provide as a runtime option without significant changes to
// how the state is stored, since the VC++ docs warn that:
// Using variables of type __m128i will cause the compiler to generate
// the SSE2 movdqa instruction. This instruction does not cause a fault
// on Pentium III processors but will result in silent failure, with
// possible side effects caused by whatever instructions movdqa
// translates into on Pentium III processors.
#if defined(HAVE_ALTIVEC)
#include "SFMT-alti.h"
#elif defined(HAVE_SSE2)
#include "SFMT-sse2.h"
#endif
/**
* This function simulate a 64-bit index of LITTLE ENDIAN
* in BIG ENDIAN machine.
*/
#ifdef ONLY64
inline static int idxof(int i) {
return i ^ 1;
}
#else
inline static int idxof(int i) {
return i;
}
#endif
/**
* This function simulates SIMD 128-bit right shift by the standard C.
* The 128-bit integer given in in is shifted by (shift * 8) bits.
* This function simulates the LITTLE ENDIAN SIMD.
* @param out the output of this function
* @param in the 128-bit data to be shifted
* @param shift the shift value
*/
#ifdef ONLY64
inline static void rshift128(w128_t *out, w128_t const *in, int shift) {
uint64_t th, tl, oh, ol;
th = ((uint64_t)in->u[2] << 32) | ((uint64_t)in->u[3]);
tl = ((uint64_t)in->u[0] << 32) | ((uint64_t)in->u[1]);
oh = th >> (shift * 8);
ol = tl >> (shift * 8);
ol |= th << (64 - shift * 8);
out->u[0] = (uint32_t)(ol >> 32);
out->u[1] = (uint32_t)ol;
out->u[2] = (uint32_t)(oh >> 32);
out->u[3] = (uint32_t)oh;
}
#else
inline static void rshift128(w128_t *out, w128_t const *in, int shift) {
uint64_t th, tl, oh, ol;
th = ((uint64_t)in->u[3] << 32) | ((uint64_t)in->u[2]);
tl = ((uint64_t)in->u[1] << 32) | ((uint64_t)in->u[0]);
oh = th >> (shift * 8);
ol = tl >> (shift * 8);
ol |= th << (64 - shift * 8);
out->u[1] = (uint32_t)(ol >> 32);
out->u[0] = (uint32_t)ol;
out->u[3] = (uint32_t)(oh >> 32);
out->u[2] = (uint32_t)oh;
}
#endif
/**
* This function simulates SIMD 128-bit left shift by the standard C.
* The 128-bit integer given in in is shifted by (shift * 8) bits.
* This function simulates the LITTLE ENDIAN SIMD.
* @param out the output of this function
* @param in the 128-bit data to be shifted
* @param shift the shift value
*/
#ifdef ONLY64
inline static void lshift128(w128_t *out, w128_t const *in, int shift) {
uint64_t th, tl, oh, ol;
th = ((uint64_t)in->u[2] << 32) | ((uint64_t)in->u[3]);
tl = ((uint64_t)in->u[0] << 32) | ((uint64_t)in->u[1]);
oh = th << (shift * 8);
ol = tl << (shift * 8);
oh |= tl >> (64 - shift * 8);
out->u[0] = (uint32_t)(ol >> 32);
out->u[1] = (uint32_t)ol;
out->u[2] = (uint32_t)(oh >> 32);
out->u[3] = (uint32_t)oh;
}
#else
inline static void lshift128(w128_t *out, w128_t const *in, int shift) {
uint64_t th, tl, oh, ol;
th = ((uint64_t)in->u[3] << 32) | ((uint64_t)in->u[2]);
tl = ((uint64_t)in->u[1] << 32) | ((uint64_t)in->u[0]);
oh = th << (shift * 8);
ol = tl << (shift * 8);
oh |= tl >> (64 - shift * 8);
out->u[1] = (uint32_t)(ol >> 32);
out->u[0] = (uint32_t)ol;
out->u[3] = (uint32_t)(oh >> 32);
out->u[2] = (uint32_t)oh;
}
#endif
/**
* This function represents the recursion formula.
* @param r output
* @param a a 128-bit part of the internal state array
* @param b a 128-bit part of the internal state array
* @param c a 128-bit part of the internal state array
* @param d a 128-bit part of the internal state array
*/
#if (!defined(HAVE_ALTIVEC)) && (!defined(HAVE_SSE2))
#ifdef ONLY64
inline static void do_recursion(w128_t *r, w128_t *a, w128_t *b, w128_t *c,
w128_t *d) {
w128_t x;
w128_t y;
lshift128(&x, a, SL2);
rshift128(&y, c, SR2);
r->u[0] = a->u[0] ^ x.u[0] ^ ((b->u[0] >> SR1) & MSK2) ^ y.u[0]
^ (d->u[0] << SL1);
r->u[1] = a->u[1] ^ x.u[1] ^ ((b->u[1] >> SR1) & MSK1) ^ y.u[1]
^ (d->u[1] << SL1);
r->u[2] = a->u[2] ^ x.u[2] ^ ((b->u[2] >> SR1) & MSK4) ^ y.u[2]
^ (d->u[2] << SL1);
r->u[3] = a->u[3] ^ x.u[3] ^ ((b->u[3] >> SR1) & MSK3) ^ y.u[3]
^ (d->u[3] << SL1);
}
#else
inline static void do_recursion(w128_t *r, w128_t *a, w128_t *b, w128_t *c,
w128_t *d) {
w128_t x;
w128_t y;
lshift128(&x, a, SL2);
rshift128(&y, c, SR2);
r->u[0] = a->u[0] ^ x.u[0] ^ ((b->u[0] >> SR1) & MSK1) ^ y.u[0]
^ (d->u[0] << SL1);
r->u[1] = a->u[1] ^ x.u[1] ^ ((b->u[1] >> SR1) & MSK2) ^ y.u[1]
^ (d->u[1] << SL1);
r->u[2] = a->u[2] ^ x.u[2] ^ ((b->u[2] >> SR1) & MSK3) ^ y.u[2]
^ (d->u[2] << SL1);
r->u[3] = a->u[3] ^ x.u[3] ^ ((b->u[3] >> SR1) & MSK4) ^ y.u[3]
^ (d->u[3] << SL1);
}
#endif
#endif
#if (!defined(HAVE_ALTIVEC)) && (!defined(HAVE_SSE2))
/**
* This function fills the internal state array with pseudorandom
* integers.
*/
void SFMTObj::GenRandAll()
{
int i;
w128_t *r1, *r2;
r1 = &sfmt.w128[SFMT::N - 2];
r2 = &sfmt.w128[SFMT::N - 1];
for (i = 0; i < SFMT::N - POS1; i++) {
do_recursion(&sfmt.w128[i], &sfmt.w128[i], &sfmt.w128[i + POS1], r1, r2);
r1 = r2;
r2 = &sfmt.w128[i];
}
for (; i < SFMT::N; i++) {
do_recursion(&sfmt.w128[i], &sfmt.w128[i], &sfmt.w128[i + POS1 - SFMT::N], r1, r2);
r1 = r2;
r2 = &sfmt.w128[i];
}
}
/**
* This function fills the user-specified array with pseudorandom
* integers.
*
* @param array an 128-bit array to be filled by pseudorandom numbers.
* @param size number of 128-bit pseudorandom numbers to be generated.
*/
void SFMTObj::GenRandArray(w128_t *array, int size)
{
int i, j;
w128_t *r1, *r2;
r1 = &sfmt.w128[SFMT::N - 2];
r2 = &sfmt.w128[SFMT::N - 1];
for (i = 0; i < SFMT::N - POS1; i++) {
do_recursion(&array[i], &sfmt.w128[i], &sfmt.w128[i + POS1], r1, r2);
r1 = r2;
r2 = &array[i];
}
for (; i < SFMT::N; i++) {
do_recursion(&array[i], &sfmt.w128[i], &array[i + POS1 - SFMT::N], r1, r2);
r1 = r2;
r2 = &array[i];
}
for (; i < size - SFMT::N; i++) {
do_recursion(&array[i], &array[i - SFMT::N], &array[i + POS1 - SFMT::N], r1, r2);
r1 = r2;
r2 = &array[i];
}
for (j = 0; j < 2 * SFMT::N - size; j++) {
sfmt.w128[j] = array[j + size - SFMT::N];
}
for (; i < size; i++, j++) {
do_recursion(&array[i], &array[i - SFMT::N], &array[i + POS1 - SFMT::N], r1, r2);
r1 = r2;
r2 = &array[i];
sfmt.w128[j] = array[i];
}
}
#endif
#if defined(BIG_ENDIAN64) && !defined(ONLY64) && !defined(HAVE_ALTIVEC)
inline static void swap(w128_t *array, int size) {
int i;
uint32_t x, y;
for (i = 0; i < size; i++) {
x = array[i].u[0];
y = array[i].u[2];
array[i].u[0] = array[i].u[1];
array[i].u[2] = array[i].u[3];
array[i].u[1] = x;
array[i].u[3] = y;
}
}
#endif
/**
* This function represents a function used in the initialization
* by init_by_array
* @param x 32-bit integer
* @return 32-bit integer
*/
static uint32_t func1(uint32_t x)
{
return (x ^ (x >> 27)) * (uint32_t)1664525UL;
}
/**
* This function represents a function used in the initialization
* by init_by_array
* @param x 32-bit integer
* @return 32-bit integer
*/
static uint32_t func2(uint32_t x)
{
return (x ^ (x >> 27)) * (uint32_t)1566083941UL;
}
/**
* This function certificate the period of 2^{MEXP}
*/
void SFMTObj::PeriodCertification()
{
int inner = 0;
int i, j;
uint32_t work;
for (i = 0; i < 4; i++)
inner ^= sfmt.u[idxof(i)] & parity[i];
for (i = 16; i > 0; i >>= 1)
inner ^= inner >> i;
inner &= 1;
/* check OK */
if (inner == 1) {
return;
}
/* check NG, and modification */
for (i = 0; i < 4; i++) {
work = 1;
for (j = 0; j < 32; j++) {
if ((work & parity[i]) != 0) {
sfmt.u[idxof(i)] ^= work;
return;
}
work = work << 1;
}
}
}
/*----------------
PUBLIC FUNCTIONS
----------------*/
/**
* This function returns the minimum size of array used for \b
* fill_array32() function.
* @return minimum size of array used for FillArray32() function.
*/
int SFMTObj::GetMinArraySize32()
{
return SFMT::N32;
}
/**
* This function returns the minimum size of array used for \b
* fill_array64() function.
* @return minimum size of array used for FillArray64() function.
*/
int SFMTObj::GetMinArraySize64()
{
return SFMT::N64;
}
#ifndef ONLY64
/**
* This function generates and returns 32-bit pseudorandom number.
* init_gen_rand or init_by_array must be called before this function.
* @return 32-bit pseudorandom number
*/
unsigned int SFMTObj::GenRand32()
{
uint32_t r;
assert(initialized);
if (idx >= SFMT::N32)
{
GenRandAll();
idx = 0;
}
r = sfmt.u[idx++];
return r;
}
#endif
/**
* This function generates and returns 64-bit pseudorandom number.
* init_gen_rand or init_by_array must be called before this function.
* The function gen_rand64 should not be called after gen_rand32,
* unless an initialization is again executed.
* @return 64-bit pseudorandom number
*/
uint64_t SFMTObj::GenRand64()
{
#if defined(BIG_ENDIAN64) && !defined(ONLY64)
uint32_t r1, r2;
#else
uint64_t r;
#endif
assert(initialized);
assert(idx % 2 == 0);
if (idx >= SFMT::N32)
{
GenRandAll();
idx = 0;
}
#if defined(BIG_ENDIAN64) && !defined(ONLY64)
r1 = sfmt.u[idx];
r2 = sfmt.u[idx + 1];
idx += 2;
return ((uint64_t)r2 << 32) | r1;
#else
r = sfmt.u64[idx / 2];
idx += 2;
return r;
#endif
}
#ifndef ONLY64
/**
* This function generates pseudorandom 32-bit integers in the
* specified array[] by one call. The number of pseudorandom integers
* is specified by the argument size, which must be at least 624 and a
* multiple of four. The generation by this function is much faster
* than the following gen_rand function.
*
* For initialization, init_gen_rand or init_by_array must be called
* before the first call of this function. This function can not be
* used after calling gen_rand function, without initialization.
*
* @param array an array where pseudorandom 32-bit integers are filled
* by this function. The pointer to the array must be \b "aligned"
* (namely, must be a multiple of 16) in the SIMD version, since it
* refers to the address of a 128-bit integer. In the standard C
* version, the pointer is arbitrary.
*
* @param size the number of 32-bit pseudorandom integers to be
* generated. size must be a multiple of 4, and greater than or equal
* to (MEXP / 128 + 1) * 4.
*
* @note \b memalign or \b posix_memalign is available to get aligned
* memory. Mac OSX doesn't have these functions, but \b malloc of OSX
* returns the pointer to the aligned memory block.
*/
void SFMTObj::FillArray32(uint32_t *array, int size)
{
assert(initialized);
assert(idx == SFMT::N32);
assert(size % 4 == 0);
assert(size >= SFMT::N32);
GenRandArray((w128_t *)array, size / 4);
idx = SFMT::N32;
}
#endif
/**
* This function generates pseudorandom 64-bit integers in the
* specified array[] by one call. The number of pseudorandom integers
* is specified by the argument size, which must be at least 312 and a
* multiple of two. The generation by this function is much faster
* than the following gen_rand function.
*
* For initialization, init_gen_rand or init_by_array must be called
* before the first call of this function. This function can not be
* used after calling gen_rand function, without initialization.
*
* @param array an array where pseudorandom 64-bit integers are filled
* by this function. The pointer to the array must be "aligned"
* (namely, must be a multiple of 16) in the SIMD version, since it
* refers to the address of a 128-bit integer. In the standard C
* version, the pointer is arbitrary.
*
* @param size the number of 64-bit pseudorandom integers to be
* generated. size must be a multiple of 2, and greater than or equal
* to (MEXP / 128 + 1) * 2
*
* @note \b memalign or \b posix_memalign is available to get aligned
* memory. Mac OSX doesn't have these functions, but \b malloc of OSX
* returns the pointer to the aligned memory block.
*/
void SFMTObj::FillArray64(uint64_t *array, int size)
{
assert(initialized);
assert(idx == SFMT::N32);
assert(size % 2 == 0);
assert(size >= SFMT::N64);
GenRandArray((w128_t *)array, size / 2);
idx = SFMT::N32;
#if defined(BIG_ENDIAN64) && !defined(ONLY64)
swap((w128_t *)array, size / 2);
#endif
}
/**
* This function initializes the internal state array with a 32-bit
* integer seed.
*
* @param seed a 32-bit integer used as the seed.
*/
void SFMTObj::InitGenRand(uint32_t seed)
{
int i;
sfmt.u[idxof(0)] = seed;
for (i = 1; i < SFMT::N32; i++)
{
sfmt.u[idxof(i)] = 1812433253UL * (sfmt.u[idxof(i - 1)]
^ (sfmt.u[idxof(i - 1)] >> 30))
+ i;
}
idx = SFMT::N32;
PeriodCertification();
#ifndef NDEBUG
initialized = 1;
#endif
}
/**
* This function initializes the internal state array,
* with an array of 32-bit integers used as the seeds
* @param init_key the array of 32-bit integers, used as a seed.
* @param key_length the length of init_key.
*/
void SFMTObj::InitByArray(uint32_t *init_key, int key_length)
{
int i, j, count;
uint32_t r;
int lag;
int mid;
int size = SFMT::N * 4;
if (size >= 623) {
lag = 11;
} else if (size >= 68) {
lag = 7;
} else if (size >= 39) {
lag = 5;
} else {
lag = 3;
}
mid = (size - lag) / 2;
memset(&sfmt, 0x8b, sizeof(sfmt));
if (key_length + 1 > SFMT::N32) {
count = key_length + 1;
} else {
count = SFMT::N32;
}
r = func1(sfmt.u[idxof(0)] ^ sfmt.u[idxof(mid)] ^ sfmt.u[idxof(SFMT::N32 - 1)]);
sfmt.u[idxof(mid)] += r;
r += key_length;
sfmt.u[idxof(mid + lag)] += r;
sfmt.u[idxof(0)] = r;
count--;
for (i = 1, j = 0; (j < count) && (j < key_length); j++)
{
r = func1(sfmt.u[idxof(i)] ^ sfmt.u[idxof((i + mid) % SFMT::N32)] ^ sfmt.u[idxof((i + SFMT::N32 - 1) % SFMT::N32)]);
sfmt.u[idxof((i + mid) % SFMT::N32)] += r;
r += init_key[j] + i;
sfmt.u[idxof((i + mid + lag) % SFMT::N32)] += r;
sfmt.u[idxof(i)] = r;
i = (i + 1) % SFMT::N32;
}
for (; j < count; j++)
{
r = func1(sfmt.u[idxof(i)] ^ sfmt.u[idxof((i + mid) % SFMT::N32)] ^ sfmt.u[idxof((i + SFMT::N32 - 1) % SFMT::N32)]);
sfmt.u[idxof((i + mid) % SFMT::N32)] += r;
r += i;
sfmt.u[idxof((i + mid + lag) % SFMT::N32)] += r;
sfmt.u[idxof(i)] = r;
i = (i + 1) % SFMT::N32;
}
for (j = 0; j < SFMT::N32; j++)
{
r = func2(sfmt.u[idxof(i)] + sfmt.u[idxof((i + mid) % SFMT::N32)] + sfmt.u[idxof((i + SFMT::N32 - 1) % SFMT::N32)]);
sfmt.u[idxof((i + mid) % SFMT::N32)] ^= r;
r -= i;
sfmt.u[idxof((i + mid + lag) % SFMT::N32)] ^= r;
sfmt.u[idxof(i)] = r;
i = (i + 1) % SFMT::N32;
}
idx = SFMT::N32;
PeriodCertification();
#ifndef NDEBUG
initialized = 1;
#endif
}

120
source/common/thirdparty/sfmt/SFMT.h vendored Normal file
View file

@ -0,0 +1,120 @@
/**
* @file SFMT.h
*
* @brief SIMD oriented Fast Mersenne Twister(SFMT) pseudorandom
* number generator
*
* @author Mutsuo Saito (Hiroshima University)
* @author Makoto Matsumoto (Hiroshima University)
*
* Copyright (C) 2006, 2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima
* University. All rights reserved.
*
* The new BSD License is applied to this software.
* see LICENSE.txt
*
* @note We assume that your system has inttypes.h. If your system
* doesn't have inttypes.h, you have to typedef uint32_t and uint64_t,
* and you have to define PRIu64 and PRIx64 in this file as follows:
* @verbatim
typedef unsigned int uint32_t
typedef unsigned long long uint64_t
#define PRIu64 "llu"
#define PRIx64 "llx"
@endverbatim
* uint32_t must be exactly 32-bit unsigned integer type (no more, no
* less), and uint64_t must be exactly 64-bit unsigned integer type.
* PRIu64 and PRIx64 are used for printf function to print 64-bit
* unsigned int and 64-bit unsigned int in hexadecimal format.
*/
#ifndef SFMT_H
#define SFMT_H
#ifndef PRIu64
#if defined(_MSC_VER) || defined(__BORLANDC__)
#define PRIu64 "I64u"
#define PRIx64 "I64x"
#else
#define PRIu64 "llu"
#define PRIx64 "llx"
#endif
#endif
#if defined(__GNUC__)
#define ALWAYSINLINE __attribute__((always_inline))
#else
#define ALWAYSINLINE
#endif
#if defined(_MSC_VER)
#if _MSC_VER >= 1200
#define PRE_ALWAYS __forceinline
#else
#define PRE_ALWAYS inline
#endif
#else
#define PRE_ALWAYS inline
#endif
/*------------------------------------------------------
128-bit SIMD data type for Altivec, SSE2 or standard C
------------------------------------------------------*/
#if defined(HAVE_ALTIVEC)
#if !defined(__APPLE__)
#include <altivec.h>
#endif
/** 128-bit data structure */
union w128_t {
vector unsigned int s;
uint32_t u[4];
uint64_t u64[2];
};
#elif defined(HAVE_SSE2)
#include <emmintrin.h>
/** 128-bit data structure */
union w128_t {
__m128i si;
uint32_t u[4];
uint64_t u64[2];
};
#else
/** 128-bit data structure */
union w128_t {
uint32_t u[4];
uint64_t u64[2];
};
#endif
/*-----------------
BASIC DEFINITIONS
-----------------*/
/** Mersenne Exponent. The period of the sequence
* is a multiple of 2^MEXP-1. */
#if !defined(MEXP)
// [RH] The default MEXP for SFMT is 19937, but since that consumes
// quite a lot of space for state, and we're using lots of different
// RNGs, default to something smaller.
#define MEXP 607
#endif
namespace SFMT
{
/** SFMT generator has an internal state array of 128-bit integers,
* and N is its size. */
enum { N = MEXP / 128 + 1 };
/** N32 is the size of internal state array when regarded as an array
* of 32-bit integers.*/
enum { N32 = N * 4 };
/** N64 is the size of internal state array when regarded as an array
* of 64-bit integers.*/
enum { N64 = N * 2 };
};
#endif

43
source/common/thirdparty/sfmt/SFMTObj.h vendored Normal file
View file

@ -0,0 +1,43 @@
#pragma once
#include <stdint.h>
#include "SFMT.h"
struct SFMTObj
{
void Init(uint32_t seed1, uint32_t seed2)
{
uint32_t seeds[2] = { seed1, seed2 };
InitByArray(seeds, 2);
}
void GenRandAll();
void GenRandArray(w128_t *array, int size);
void PeriodCertification();
int GetMinArraySize32();
int GetMinArraySize64();
unsigned int GenRand32();
uint64_t GenRand64();
void FillArray32(uint32_t *array, int size);
void FillArray64(uint64_t *array, int size);
void InitGenRand(uint32_t seed);
void InitByArray(uint32_t *init_key, int key_length);
protected:
/** index counter to the 32-bit internal state array */
int idx;
/** the 128-bit internal state array */
union
{
w128_t w128[SFMT::N];
unsigned int u[SFMT::N32];
uint64_t u64[SFMT::N64];
} sfmt;
/** a flag: it is 0 if and only if the internal state is not yet
* initialized. */
#ifndef NDEBUG
bool initialized = false;
#endif
};

176
source/common/thirdparty/strnatcmp.c vendored Normal file
View file

@ -0,0 +1,176 @@
/* -*- mode: c; c-file-style: "k&r" -*-
strnatcmp.c -- Perform 'natural order' comparisons of strings in C.
Copyright (C) 2000, 2004 by Martin Pool <mbp sourcefrog net>
This software is provided 'as-is', without any express or implied
warranty. In no event will the authors be held liable for any damages
arising from the use of this software.
Permission is granted to anyone to use this software for any purpose,
including commercial applications, and to alter it and redistribute it
freely, subject to the following restrictions:
1. The origin of this software must not be misrepresented; you must not
claim that you wrote the original software. If you use this software
in a product, an acknowledgment in the product documentation would be
appreciated but is not required.
2. Altered source versions must be plainly marked as such, and must not be
misrepresented as being the original software.
3. This notice may not be removed or altered from any source distribution.
*/
/* partial change history:
*
* 2004-10-10 mbp: Lift out character type dependencies into macros.
*
* Eric Sosman pointed out that ctype functions take a parameter whose
* value must be that of an unsigned int, even on platforms that have
* negative chars in their default char type.
*/
#include <ctype.h>
#include <assert.h>
#include "strnatcmp.h"
/* These are defined as macros to make it easier to adapt this code to
* different characters types or comparison functions. */
static __inline int
nat_isdigit(nat_char a)
{
return isdigit((unsigned char) a);
}
static __inline int
nat_isspace(nat_char a)
{
return isspace((unsigned char) a);
}
static __inline nat_char
nat_toupper(nat_char a)
{
return toupper((unsigned char) a);
}
static int
compare_right(nat_char const *a, nat_char const *b)
{
int bias = 0;
/* The longest run of digits wins. That aside, the greatest
value wins, but we can't know that it will until we've scanned
both numbers to know that they have the same magnitude, so we
remember it in BIAS. */
for (;; a++, b++) {
if (!nat_isdigit(*a) && !nat_isdigit(*b))
return bias;
else if (!nat_isdigit(*a))
return -1;
else if (!nat_isdigit(*b))
return +1;
else if (*a < *b) {
if (!bias)
bias = -1;
} else if (*a > *b) {
if (!bias)
bias = +1;
} else if (!*a && !*b)
return bias;
}
return 0;
}
static int
compare_left(nat_char const *a, nat_char const *b)
{
/* Compare two left-aligned numbers: the first to have a
different value wins. */
for (;; a++, b++) {
if (!nat_isdigit(*a) && !nat_isdigit(*b))
return 0;
else if (!nat_isdigit(*a))
return -1;
else if (!nat_isdigit(*b))
return +1;
else if (*a < *b)
return -1;
else if (*a > *b)
return +1;
}
return 0;
}
static int strnatcmp0(nat_char const *a, nat_char const *b, int fold_case)
{
int ai, bi;
nat_char ca, cb;
int fractional, result;
assert(a && b);
ai = bi = 0;
while (1) {
ca = a[ai]; cb = b[bi];
/* skip over leading spaces or zeros */
while (nat_isspace(ca))
ca = a[++ai];
while (nat_isspace(cb))
cb = b[++bi];
/* process run of digits */
if (nat_isdigit(ca) && nat_isdigit(cb)) {
fractional = (ca == '0' || cb == '0');
if (fractional) {
if ((result = compare_left(a+ai, b+bi)) != 0)
return result;
} else {
if ((result = compare_right(a+ai, b+bi)) != 0)
return result;
}
}
if (!ca && !cb) {
/* The strings compare the same. Perhaps the caller
will want to call strcmp to break the tie. */
return 0;
}
if (fold_case) {
ca = nat_toupper(ca);
cb = nat_toupper(cb);
}
if (ca < cb)
return -1;
else if (ca > cb)
return +1;
++ai; ++bi;
}
}
int strnatcmp(nat_char const *a, nat_char const *b) {
return strnatcmp0(a, b, 0);
}
/* Compare, recognizing numeric string and ignoring case. */
int strnatcasecmp(nat_char const *a, nat_char const *b) {
return strnatcmp0(a, b, 1);
}

39
source/common/thirdparty/strnatcmp.h vendored Normal file
View file

@ -0,0 +1,39 @@
/* -*- mode: c; c-file-style: "k&r" -*-
strnatcmp.c -- Perform 'natural order' comparisons of strings in C.
Copyright (C) 2000, 2004 by Martin Pool <mbp sourcefrog net>
This software is provided 'as-is', without any express or implied
warranty. In no event will the authors be held liable for any damages
arising from the use of this software.
Permission is granted to anyone to use this software for any purpose,
including commercial applications, and to alter it and redistribute it
freely, subject to the following restrictions:
1. The origin of this software must not be misrepresented; you must not
claim that you wrote the original software. If you use this software
in a product, an acknowledgment in the product documentation would be
appreciated but is not required.
2. Altered source versions must be plainly marked as such, and must not be
misrepresented as being the original software.
3. This notice may not be removed or altered from any source distribution.
*/
#ifdef __cplusplus
extern "C"
{
#endif
/* CUSTOMIZATION SECTION
*
* You can change this typedef, but must then also change the inline
* functions in strnatcmp.c */
typedef char nat_char;
int strnatcmp(nat_char const *a, nat_char const *b);
int strnatcasecmp(nat_char const *a, nat_char const *b);
#ifdef __cplusplus
}
#endif

View file

@ -0,0 +1,127 @@
#include <stdint.h>
#include <ctype.h>
#include <string.h>
/* ======================================================================== */
/* By Paul Hsieh (C) 2004, 2005. Covered under the Paul Hsieh derivative
license. See:
http://www.azillionmonkeys.com/qed/weblicense.html for license details.
http://www.azillionmonkeys.com/qed/hash.html */
#undef get16bits
#if (defined(__GNUC__) && defined(__i386__)) || defined(__WATCOMC__) \
|| defined(_MSC_VER) || defined (__BORLANDC__) || defined (__TURBOC__)
#define get16bits(d) (*((const uint16_t *) (d)))
#endif
#if !defined (get16bits)
#define get16bits(d) ((((uint32_t)(((const uint8_t *)(d))[1])) << 8)\
+(uint32_t)(((const uint8_t *)(d))[0]) )
#endif
uint32_t SuperFastHash (const char *data, size_t len)
{
uint32_t hash = 0, tmp;
size_t rem;
if (len == 0 || data == NULL) return 0;
rem = len & 3;
len >>= 2;
/* Main loop */
for (;len > 0; len--)
{
hash += get16bits (data);
tmp = (get16bits (data+2) << 11) ^ hash;
hash = (hash << 16) ^ tmp;
data += 2*sizeof (uint16_t);
hash += hash >> 11;
}
/* Handle end cases */
switch (rem)
{
case 3: hash += get16bits (data);
hash ^= hash << 16;
hash ^= data[sizeof (uint16_t)] << 18;
hash += hash >> 11;
break;
case 2: hash += get16bits (data);
hash ^= hash << 11;
hash += hash >> 17;
break;
case 1: hash += *data;
hash ^= hash << 10;
hash += hash >> 1;
}
/* Force "avalanching" of final 127 bits */
hash ^= hash << 3;
hash += hash >> 5;
hash ^= hash << 4;
hash += hash >> 17;
hash ^= hash << 25;
hash += hash >> 6;
return hash;
}
/* A modified version to do a case-insensitive hash */
#undef get16bits
#define get16bits(d) ((((uint32_t)tolower(((const uint8_t *)(d))[1])) << 8)\
+(uint32_t)tolower(((const uint8_t *)(d))[0]) )
uint32_t SuperFastHashI (const char *data, size_t len)
{
uint32_t hash = 0, tmp;
size_t rem;
if (len <= 0 || data == NULL) return 0;
rem = len & 3;
len >>= 2;
/* Main loop */
for (;len > 0; len--)
{
hash += get16bits (data);
tmp = (get16bits (data+2) << 11) ^ hash;
hash = (hash << 16) ^ tmp;
data += 2*sizeof (uint16_t);
hash += hash >> 11;
}
/* Handle end cases */
switch (rem)
{
case 3: hash += get16bits (data);
hash ^= hash << 16;
hash ^= tolower(data[sizeof (uint16_t)]) << 18;
hash += hash >> 11;
break;
case 2: hash += get16bits (data);
hash ^= hash << 11;
hash += hash >> 17;
break;
case 1: hash += tolower(*data);
hash ^= hash << 10;
hash += hash >> 1;
}
/* Force "avalanching" of final 127 bits */
hash ^= hash << 3;
hash += hash >> 5;
hash ^= hash << 4;
hash += hash >> 17;
hash ^= hash << 25;
hash += hash >> 6;
return hash;
}
/* ======================================================================== */

View file

@ -0,0 +1,20 @@
#pragma once
#include <stdint.h>
uint32_t SuperFastHash (const char *data, size_t len);
uint32_t SuperFastHashI (const char *data, size_t len);
inline unsigned int MakeKey(const char* s)
{
if (s == NULL)
{
return 0;
}
return SuperFastHashI(s, strlen(s));
}
inline unsigned int MakeKey(const char* s, size_t len)
{
return SuperFastHashI(s, len);
}

View file

@ -35,7 +35,7 @@
#include "cmdlib.h"
#include "i_system.h"
#include "findfile.h"
#include <sys/types.h>
#include <sys/stat.h>

View file

@ -0,0 +1,77 @@
/*
** colormatcher.h
**
**---------------------------------------------------------------------------
** Copyright 1998-2006 Randy Heit
** All rights reserved.
**
** Redistribution and use in source and binary forms, with or without
** modification, are permitted provided that the following conditions
** are met:
**
** 1. Redistributions of source code must retain the above copyright
** notice, this list of conditions and the following disclaimer.
** 2. Redistributions in binary form must reproduce the above copyright
** notice, this list of conditions and the following disclaimer in the
** documentation and/or other materials provided with the distribution.
** 3. The name of the author may not be used to endorse or promote products
** derived from this software without specific prior written permission.
**
** THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
** IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
** OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
** IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
** INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
** NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
** DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
** THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
** (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
** THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
**---------------------------------------------------------------------------
**
** Once upon a time, this tried to be a fast closest color finding system.
** It was, but the results were not as good as I would like, so I didn't
** actually use it. But I did keep the code around in case I ever felt like
** revisiting the problem. I never did, so now it's relegated to the mists
** of SVN history, and this is just a thin wrapper around BestColor().
**
*/
#ifndef __COLORMATCHER_H__
#define __COLORMATCHER_H__
#include "palutil.h"
int BestColor (const uint32_t *pal_in, int r, int g, int b, int first, int num);
class FColorMatcher
{
public:
FColorMatcher () = default;
FColorMatcher (const uint32_t *palette) { Pal = reinterpret_cast<const PalEntry*>(palette); }
FColorMatcher (const FColorMatcher &other) = default;
void SetPalette (const uint32_t *palette) { Pal = reinterpret_cast<const PalEntry*>(palette); }
uint8_t Pick (int r, int g, int b)
{
if (Pal == nullptr)
return 1;
return (uint8_t)BestColor ((uint32_t *)Pal, r, g, b, 1, 255);
}
uint8_t Pick (PalEntry pe)
{
return Pick(pe.r, pe.g, pe.b);
}
FColorMatcher &operator= (const FColorMatcher &other) = default;
private:
const PalEntry *Pal;
};
extern FColorMatcher ColorMatcher;
#endif //__COLORMATCHER_H__

View file

@ -0,0 +1,115 @@
/*
** engineerrors.cpp
** Contains error classes that can be thrown around
**
**---------------------------------------------------------------------------
** Copyright 1998-2016 Randy Heit
** Copyright 2005-2020 Christoph Oelckers
** All rights reserved.
**
** Redistribution and use in source and binary forms, with or without
** modification, are permitted provided that the following conditions
** are met:
**
** 1. Redistributions of source code must retain the above copyright
** notice, this list of conditions and the following disclaimer.
** 2. Redistributions in binary form must reproduce the above copyright
** notice, this list of conditions and the following disclaimer in the
** documentation and/or other materials provided with the distribution.
** 3. The name of the author may not be used to endorse or promote products
** derived from this software without specific prior written permission.
**
** THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
** IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
** OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
** IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
** INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
** NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
** DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
** THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
** (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
** THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
**---------------------------------------------------------------------------
**
*/
bool gameisdead;
#ifdef _WIN32
#include <windows.h>
#include "zstring.h"
void I_DebugPrint(const char *cp)
{
if (IsDebuggerPresent())
{
auto wstr = WideString(cp);
OutputDebugStringW(wstr.c_str());
}
}
#else
void I_DebugPrint(const char *cp)
{
}
#endif
#include "engineerrors.h"
#include "printf.h"
//==========================================================================
//
// I_Error
//
// Throw an error that will send us to the console if we are far enough
// along in the startup process.
//
//==========================================================================
void I_Error(const char *error, ...)
{
va_list argptr;
char errortext[MAX_ERRORTEXT];
va_start(argptr, error);
myvsnprintf(errortext, MAX_ERRORTEXT, error, argptr);
va_end(argptr);
I_DebugPrint(errortext);
throw CRecoverableError(errortext);
}
//==========================================================================
//
// I_FatalError
//
// Throw an error that will end the game.
//
//==========================================================================
extern FILE *Logfile;
void I_FatalError(const char *error, ...)
{
static bool alreadyThrown = false;
gameisdead = true;
if (!alreadyThrown) // ignore all but the first message -- killough
{
alreadyThrown = true;
char errortext[MAX_ERRORTEXT];
va_list argptr;
va_start(argptr, error);
myvsnprintf(errortext, MAX_ERRORTEXT, error, argptr);
va_end(argptr);
I_DebugPrint(errortext);
// Record error to log (if logging)
if (Logfile)
{
fprintf(Logfile, "\n**** DIED WITH FATAL ERROR:\n%s\n", errortext);
fflush(Logfile);
}
throw CFatalError(errortext);
}
std::terminate(); // recursive I_FatalErrors must immediately terminate.
}

View file

@ -0,0 +1,116 @@
/*
** doomerrors.h
** Contains error classes that can be thrown around
**
**---------------------------------------------------------------------------
** Copyright 1998-2006 Randy Heit
** All rights reserved.
**
** Redistribution and use in source and binary forms, with or without
** modification, are permitted provided that the following conditions
** are met:
**
** 1. Redistributions of source code must retain the above copyright
** notice, this list of conditions and the following disclaimer.
** 2. Redistributions in binary form must reproduce the above copyright
** notice, this list of conditions and the following disclaimer in the
** documentation and/or other materials provided with the distribution.
** 3. The name of the author may not be used to endorse or promote products
** derived from this software without specific prior written permission.
**
** THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
** IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
** OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
** IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
** INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
** NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
** DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
** THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
** (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
** THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
**---------------------------------------------------------------------------
**
*/
#ifndef __ERRORS_H__
#define __ERRORS_H__
#include <string.h>
#include <stdio.h>
#include <exception>
#include <stdexcept>
#include "basics.h"
#define MAX_ERRORTEXT 1024
class CEngineError : public std::exception
{
public:
CEngineError ()
{
m_Message[0] = '\0';
}
CEngineError (const char *message)
{
SetMessage (message);
}
void SetMessage (const char *message)
{
strncpy (m_Message, message, MAX_ERRORTEXT-1);
m_Message[MAX_ERRORTEXT-1] = '\0';
}
void AppendMessage(const char *message)
{
size_t len = strlen(m_Message);
strncpy(m_Message + len, message, MAX_ERRORTEXT - 1 - len);
m_Message[MAX_ERRORTEXT - 1] = '\0';
}
const char *GetMessage (void) const
{
if (m_Message[0] != '\0')
return (const char *)m_Message;
else
return NULL;
}
char const *what() const noexcept override
{
return m_Message;
}
protected:
char m_Message[MAX_ERRORTEXT];
};
class CRecoverableError : public CEngineError
{
public:
CRecoverableError() : CEngineError() {}
CRecoverableError(const char *message) : CEngineError(message) {}
};
class CFatalError : public CEngineError
{
public:
CFatalError() : CEngineError() {}
CFatalError(const char *message) : CEngineError(message) {}
};
class CExitEvent : public std::exception
{
int m_reason;
public:
CExitEvent(int reason) { m_reason = reason; }
char const *what() const noexcept override
{
return "The game wants to exit";
}
int Reason() const { return m_reason; }
};
void I_ShowFatalError(const char *message);
void I_Error (const char *error, ...) GCCPRINTF(1,2);
void I_FatalError (const char *error, ...) GCCPRINTF(1,2);
#endif //__ERRORS_H__

View file

@ -0,0 +1,190 @@
/*
** findfile.cpp
** Warpper around the native directory scanning API
**
**---------------------------------------------------------------------------
** Copyright 1998-2016 Randy Heit
** Copyright 2005-2020 Christoph Oelckers
**
** Redistribution and use in source and binary forms, with or without
** modification, are permitted provided that the following conditions
** are met:
**
** 1. Redistributions of source code must retain the above copyright
** notice, this list of conditions and the following disclaimer.
** 2. Redistributions in binary form must reproduce the above copyright
** notice, this list of conditions and the following disclaimer in the
** documentation and/or other materials provided with the distribution.
** 3. The name of the author may not be used to endorse or promote products
** derived from this software without specific prior written permission.
**
** THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
** IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
** OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
** IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
** INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
** NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
** DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
** THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
** (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
** THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
**---------------------------------------------------------------------------
**
*/
#include "findfile.h"
#include "zstring.h"
#ifndef _WIN32
#include <fnmatch.h>
static const char *pattern;
#if defined(__APPLE__) && MAC_OS_X_VERSION_MAX_ALLOWED < 1080
static int matchfile(struct dirent *ent)
#else
static int matchfile(const struct dirent *ent)
#endif
{
return fnmatch(pattern, ent->d_name, FNM_NOESCAPE) == 0;
}
void *I_FindFirst(const char *const filespec, findstate_t *const fileinfo)
{
FString dir;
const char *const slash = strrchr(filespec, '/');
if (slash)
{
pattern = slash + 1;
dir = FString(filespec, slash - filespec + 1);
fileinfo->path = dir;
}
else
{
pattern = filespec;
dir = ".";
}
fileinfo->current = 0;
fileinfo->count = scandir(dir.GetChars(), &fileinfo->namelist, matchfile, alphasort);
if (fileinfo->count > 0)
{
return fileinfo;
}
return (void *)-1;
}
int I_FindNext(void *const handle, findstate_t *const fileinfo)
{
findstate_t *const state = static_cast<findstate_t *>(handle);
if (state->current < fileinfo->count)
{
return ++state->current < fileinfo->count ? 0 : -1;
}
return -1;
}
int I_FindClose(void *const handle)
{
findstate_t *const state = static_cast<findstate_t *>(handle);
if (handle != (void *)-1 && state->count > 0)
{
for (int i = 0; i < state->count; ++i)
{
free(state->namelist[i]);
}
free(state->namelist);
state->namelist = nullptr;
state->count = 0;
}
return 0;
}
int I_FindAttr(findstate_t *const fileinfo)
{
dirent *const ent = fileinfo->namelist[fileinfo->current];
const FString path = fileinfo->path + ent->d_name;
bool isdir;
if (DirEntryExists(path, &isdir))
{
return isdir ? FA_DIREC : 0;
}
return 0;
}
#else
#include <windows.h>
//==========================================================================
//
// I_FindFirst
//
// Start a pattern matching sequence.
//
//==========================================================================
void *I_FindFirst(const char *filespec, findstate_t *fileinfo)
{
static_assert(sizeof(WIN32_FIND_DATAW) == sizeof(fileinfo->FindData), "FindData size mismatch");
auto widespec = WideString(filespec);
fileinfo->UTF8Name = "";
return FindFirstFileW(widespec.c_str(), (LPWIN32_FIND_DATAW)&fileinfo->FindData);
}
//==========================================================================
//
// I_FindNext
//
// Return the next file in a pattern matching sequence.
//
//==========================================================================
int I_FindNext(void *handle, findstate_t *fileinfo)
{
fileinfo->UTF8Name = "";
return !FindNextFileW((HANDLE)handle, (LPWIN32_FIND_DATAW)&fileinfo->FindData);
}
//==========================================================================
//
// I_FindClose
//
// Finish a pattern matching sequence.
//
//==========================================================================
int I_FindClose(void *handle)
{
return FindClose((HANDLE)handle);
}
//==========================================================================
//
// I_FindName
//
// Returns the name for an entry
//
//==========================================================================
const char *I_FindName(findstate_t *fileinfo)
{
if (fileinfo->UTF8Name.IsEmpty()) fileinfo->UTF8Name = fileinfo->FindData.Name;
return fileinfo->UTF8Name.GetChars();
}
#endif

View file

@ -1,16 +1,56 @@
#pragma once
// Directory searching routines
#include <stdint.h>
#include "zstring.h"
#ifndef _WIN32
struct findstate_t
{
private:
FString path;
struct dirent **namelist;
int current;
int count;
friend void *I_FindFirst(const char *filespec, findstate_t *fileinfo);
friend int I_FindNext(void *handle, findstate_t *fileinfo);
friend const char *I_FindName(findstate_t *fileinfo);
friend int I_FindAttr(findstate_t *fileinfo);
friend int I_FindClose(void *handle);
};
int I_FindAttr (findstate_t *fileinfo);
inline const char *I_FindName(findstate_t *fileinfo)
{
return (fileinfo->namelist[fileinfo->current]->d_name);
}
#define FA_RDONLY 1
#define FA_HIDDEN 2
#define FA_SYSTEM 4
#define FA_DIREC 8
#define FA_ARCH 16
#else
// Mirror WIN32_FIND_DATAW in <winbase.h>
struct findstate_t
{
private:
struct FileTime
{
uint32_t lo, hi;
};
struct WinData
{
uint32_t Attribs;
uint32_t Times[3 * 2];
FileTime Times[3];
uint32_t Size[2];
uint32_t Reserved[2];
wchar_t Name[260];
@ -25,9 +65,6 @@ private:
friend int I_FindAttr(findstate_t *fileinfo);
};
void *I_FindFirst (const char *filespec, findstate_t *fileinfo);
int I_FindNext (void *handle, findstate_t *fileinfo);
int I_FindClose (void *handle);
const char *I_FindName(findstate_t *fileinfo);
inline int I_FindAttr(findstate_t *fileinfo)
@ -40,3 +77,9 @@ inline int I_FindAttr(findstate_t *fileinfo)
#define FA_SYSTEM 0x00000004
#define FA_DIREC 0x00000010
#define FA_ARCH 0x00000020
#endif
void *I_FindFirst (const char *filespec, findstate_t *fileinfo);
int I_FindNext (void *handle, findstate_t *fileinfo);
int I_FindClose (void *handle);

Some files were not shown because too many files have changed in this diff Show more