mirror of
https://repo.dactyloidae.xyz/Dactyloidae/UXP.git
synced 2026-09-25 18:07:31 +09:00
Issue #1258 - Part 1: Import mailnews, ldap, and mork from comm-esr52.9.1
This commit is contained in:
parent
23e0d82436
commit
e400f4130a
1564 changed files with 510348 additions and 0 deletions
6
mailnews/extensions/bayesian-spam-filter/moz.build
Normal file
6
mailnews/extensions/bayesian-spam-filter/moz.build
Normal file
|
|
@ -0,0 +1,6 @@
|
|||
# vim: set filetype=python:
|
||||
# This Source Code Form is subject to the terms of the Mozilla Public
|
||||
# License, v. 2.0. If a copy of the MPL was not distributed with this
|
||||
# file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
||||
|
||||
DIRS += ['src']
|
||||
11
mailnews/extensions/bayesian-spam-filter/src/moz.build
Normal file
11
mailnews/extensions/bayesian-spam-filter/src/moz.build
Normal file
|
|
@ -0,0 +1,11 @@
|
|||
# vim: set filetype=python:
|
||||
# This Source Code Form is subject to the terms of the Mozilla Public
|
||||
# License, v. 2.0. If a copy of the MPL was not distributed with this
|
||||
# file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
||||
|
||||
SOURCES += [
|
||||
'nsBayesianFilter.cpp',
|
||||
]
|
||||
|
||||
FINAL_LIBRARY = 'mail'
|
||||
|
||||
2758
mailnews/extensions/bayesian-spam-filter/src/nsBayesianFilter.cpp
Normal file
2758
mailnews/extensions/bayesian-spam-filter/src/nsBayesianFilter.cpp
Normal file
File diff suppressed because it is too large
Load diff
404
mailnews/extensions/bayesian-spam-filter/src/nsBayesianFilter.h
Normal file
404
mailnews/extensions/bayesian-spam-filter/src/nsBayesianFilter.h
Normal file
|
|
@ -0,0 +1,404 @@
|
|||
/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
|
||||
/* This Source Code Form is subject to the terms of the Mozilla Public
|
||||
* License, v. 2.0. If a copy of the MPL was not distributed with this
|
||||
* file, You can obtain one at http://mozilla.org/MPL/2.0/. */
|
||||
|
||||
#ifndef nsBayesianFilter_h__
|
||||
#define nsBayesianFilter_h__
|
||||
|
||||
#include <stdio.h>
|
||||
#include "nsCOMPtr.h"
|
||||
#include "nsIMsgFilterPlugin.h"
|
||||
#include "nsISemanticUnitScanner.h"
|
||||
#include "PLDHashTable.h"
|
||||
#include "nsITimer.h"
|
||||
#include "nsTArray.h"
|
||||
#include "nsStringGlue.h"
|
||||
#include "nsWeakReference.h"
|
||||
#include "nsIObserver.h"
|
||||
|
||||
// XXX can't simply byte align arenas, must at least 2-byte align.
|
||||
#define PL_ARENA_CONST_ALIGN_MASK 1
|
||||
#include "plarena.h"
|
||||
|
||||
#define DEFAULT_MIN_INTERVAL_BETWEEN_WRITES 15*60*1000
|
||||
|
||||
struct Token;
|
||||
class TokenEnumeration;
|
||||
class TokenAnalyzer;
|
||||
class nsIMsgWindow;
|
||||
class nsIMimeHeaders;
|
||||
class nsIUTF8StringEnumerator;
|
||||
struct BaseToken;
|
||||
struct CorpusToken;
|
||||
|
||||
/**
|
||||
* Helper class to enumerate Token objects in a PLDHashTable
|
||||
* safely and without copying (see bugzilla #174859). The
|
||||
* enumeration is safe to use until an Add()
|
||||
* or Remove() is performed on the table.
|
||||
*/
|
||||
class TokenEnumeration {
|
||||
public:
|
||||
TokenEnumeration(PLDHashTable* table);
|
||||
bool hasMoreTokens();
|
||||
BaseToken* nextToken();
|
||||
|
||||
private:
|
||||
PLDHashTable::Iterator mIterator;
|
||||
};
|
||||
|
||||
// A trait is some aspect of a message, like being junk or tagged as
|
||||
// Personal, that the statistical classifier should track. The Trait
|
||||
// structure is a per-token representation of information pertaining to
|
||||
// a message trait.
|
||||
//
|
||||
// Traits per token are maintained as a linked list.
|
||||
//
|
||||
struct TraitPerToken
|
||||
{
|
||||
uint32_t mId; // identifying number for a trait
|
||||
uint32_t mCount; // count of messages with this token and trait
|
||||
uint32_t mNextLink; // index in mTraitStore for the next trait, or 0
|
||||
// for none
|
||||
TraitPerToken(uint32_t aId, uint32_t aCount); // inititializer
|
||||
};
|
||||
|
||||
// An Analysis is the statistical results for a particular message, a
|
||||
// particular token, and for a particular pair of trait/antitrait, that
|
||||
// is then used in subsequent analysis to score the message.
|
||||
//
|
||||
// Analyses per token are maintained as a linked list.
|
||||
//
|
||||
struct AnalysisPerToken
|
||||
{
|
||||
uint32_t mTraitIndex; // index representing a protrait/antitrait pair.
|
||||
// So if we are analyzing 3 different traits, then
|
||||
// the first trait is 0, the second 1, etc.
|
||||
double mDistance; // absolute value of mProbability - 0.5
|
||||
double mProbability; // relative indicator of match of trait to token
|
||||
uint32_t mNextLink; // index in mAnalysisStore for the Analysis object
|
||||
// for the next trait index, or 0 for none.
|
||||
// initializer
|
||||
AnalysisPerToken(uint32_t aTraitIndex, double aDistance, double aProbability);
|
||||
};
|
||||
|
||||
class TokenHash {
|
||||
public:
|
||||
|
||||
virtual ~TokenHash();
|
||||
/**
|
||||
* Clears out the previous message tokens.
|
||||
*/
|
||||
nsresult clearTokens();
|
||||
uint32_t countTokens();
|
||||
TokenEnumeration getTokens();
|
||||
BaseToken* add(const char* word);
|
||||
|
||||
protected:
|
||||
TokenHash(uint32_t entrySize);
|
||||
PLArenaPool mWordPool;
|
||||
uint32_t mEntrySize;
|
||||
PLDHashTable mTokenTable;
|
||||
char* copyWord(const char* word, uint32_t len);
|
||||
BaseToken* get(const char* word);
|
||||
};
|
||||
|
||||
class Tokenizer: public TokenHash {
|
||||
public:
|
||||
Tokenizer();
|
||||
~Tokenizer();
|
||||
|
||||
Token* get(const char* word);
|
||||
|
||||
// The training set keeps an occurrence count on each word. This count
|
||||
// is supposed to count the # of messsages it occurs in.
|
||||
// When add/remove is called while tokenizing a message and NOT the training set,
|
||||
//
|
||||
Token* add(const char* word, uint32_t count = 1);
|
||||
|
||||
Token* copyTokens();
|
||||
|
||||
void tokenize(const char* text);
|
||||
|
||||
/**
|
||||
* Creates specific tokens based on the mime headers for the message being tokenized
|
||||
*/
|
||||
void tokenizeHeaders(nsIUTF8StringEnumerator * aHeaderNames, nsIUTF8StringEnumerator * aHeaderValues);
|
||||
|
||||
void tokenizeAttachment(const char * aContentType, const char * aFileName);
|
||||
|
||||
nsCString mBodyDelimiters; // delimiters for body tokenization
|
||||
nsCString mHeaderDelimiters; // delimiters for header tokenization
|
||||
|
||||
// arrays of extra headers to tokenize / to not tokenize
|
||||
nsTArray<nsCString> mEnabledHeaders;
|
||||
nsTArray<nsCString> mDisabledHeaders;
|
||||
// Delimiters used in tokenizing a particular header.
|
||||
// Parallel array to mEnabledHeaders
|
||||
nsTArray<nsCString> mEnabledHeadersDelimiters;
|
||||
bool mCustomHeaderTokenization; // Are there any preference-set tokenization customizations?
|
||||
uint32_t mMaxLengthForToken; // maximum length of a token
|
||||
// should we convert iframe to div during tokenization?
|
||||
bool mIframeToDiv;
|
||||
|
||||
private:
|
||||
|
||||
void tokenize_ascii_word(char * word);
|
||||
void tokenize_japanese_word(char* chunk);
|
||||
inline void addTokenForHeader(const char * aTokenPrefix, nsACString& aValue,
|
||||
bool aTokenizeValue = false, const char* aDelimiters = nullptr);
|
||||
nsresult stripHTML(const nsAString& inString, nsAString& outString);
|
||||
// helper function to escape \n, \t, etc from a CString
|
||||
void UnescapeCString(nsCString& aCString);
|
||||
|
||||
private:
|
||||
nsCOMPtr<nsISemanticUnitScanner> mScanner;
|
||||
};
|
||||
|
||||
/**
|
||||
* Implements storage of a collection of message tokens and counts for
|
||||
* a corpus of classified messages
|
||||
*/
|
||||
|
||||
class CorpusStore: public TokenHash {
|
||||
public:
|
||||
CorpusStore();
|
||||
~CorpusStore();
|
||||
|
||||
/**
|
||||
* retrieve the token structure for a particular string
|
||||
*
|
||||
* @param word the character representation of the token
|
||||
*
|
||||
* @return token structure containing counts, null if not found
|
||||
*/
|
||||
CorpusToken* get(const char* word);
|
||||
|
||||
/**
|
||||
* add tokens to the storage, or increment counts if already exists.
|
||||
*
|
||||
* @param aTokenizer tokenizer for the list of tokens to remember
|
||||
* @param aTraitId id for the trait whose counts will be remembered
|
||||
* @param aCount number of new messages represented by the token list
|
||||
*/
|
||||
void rememberTokens(Tokenizer& aTokenizer, uint32_t aTraitId, uint32_t aCount);
|
||||
|
||||
/**
|
||||
* decrement counts for tokens in the storage, removing if all counts
|
||||
* are zero
|
||||
*
|
||||
* @param aTokenizer tokenizer for the list of tokens to forget
|
||||
* @param aTraitId id for the trait whose counts will be removed
|
||||
* @param aCount number of messages represented by the token list
|
||||
*/
|
||||
void forgetTokens(Tokenizer& aTokenizer, uint32_t aTraitId, uint32_t aCount);
|
||||
|
||||
/**
|
||||
* write the corpus information to file storage
|
||||
*
|
||||
* @param aMaximumTokenCount prune tokens if number of tokens exceeds
|
||||
* this value. == 0 for no pruning
|
||||
*/
|
||||
void writeTrainingData(uint32_t aMaximumTokenCount);
|
||||
|
||||
/**
|
||||
* read the corpus information from file storage
|
||||
*/
|
||||
void readTrainingData();
|
||||
|
||||
/**
|
||||
* delete the local corpus storage file and data
|
||||
*/
|
||||
nsresult resetTrainingData();
|
||||
|
||||
/**
|
||||
* get the count of messages whose tokens are stored that are associated
|
||||
* with a trait
|
||||
*
|
||||
* @param aTraitId identifier for the trait
|
||||
* @return number of messages for that trait
|
||||
*/
|
||||
uint32_t getMessageCount(uint32_t aTraitId);
|
||||
|
||||
/**
|
||||
* set the count of messages whose tokens are stored that are associated
|
||||
* with a trait
|
||||
*
|
||||
* @param aTraitId identifier for the trait
|
||||
* @param aCount number of messages for that trait
|
||||
*/
|
||||
void setMessageCount(uint32_t aTraitId, uint32_t aCount);
|
||||
|
||||
/**
|
||||
* get the count of messages associated with a particular token and trait
|
||||
*
|
||||
* @param token the token string and associated counts
|
||||
* @param aTraitId identifier for the trait
|
||||
*/
|
||||
uint32_t getTraitCount(CorpusToken *token, uint32_t aTraitId);
|
||||
|
||||
/**
|
||||
* Add (or remove) data from a particular file to the corpus data.
|
||||
*
|
||||
* @param aFile the file with the data, in the format:
|
||||
*
|
||||
* Format of the trait file for version 1:
|
||||
* [0xFCA93601] (the 01 is the version)
|
||||
* for each trait to write:
|
||||
* [id of trait to write] (0 means end of list)
|
||||
* [number of messages per trait]
|
||||
* for each token with non-zero count
|
||||
* [count]
|
||||
* [length of word]word
|
||||
*
|
||||
* @param aIsAdd should the data be added, or removed? true if adding,
|
||||
* else removing.
|
||||
*
|
||||
* @param aRemapCount number of items in the parallel arrays aFromTraits,
|
||||
* aToTraits. These arrays allow conversion of the
|
||||
* trait id stored in the file (which may be originated
|
||||
* externally) to the trait id used in the local corpus
|
||||
* (which is defined locally using nsIMsgTraitService).
|
||||
*
|
||||
* @param aFromTraits array of trait ids used in aFile. If aFile contains
|
||||
* trait ids that are not in this array, they are not
|
||||
* remapped, but assummed to be local trait ids.
|
||||
*
|
||||
* @param aToTraits array of trait ids, corresponding to elements of
|
||||
* aFromTraits, that represent the local trait ids to be
|
||||
* used in storing data from aFile into the local corpus.
|
||||
*
|
||||
*/
|
||||
nsresult UpdateData(nsIFile *aFile, bool aIsAdd,
|
||||
uint32_t aRemapCount, uint32_t *aFromTraits,
|
||||
uint32_t *aToTraits);
|
||||
|
||||
/**
|
||||
* remove all counts (message and tokens) for a trait id
|
||||
*
|
||||
* @param aTrait trait id for the trait to remove
|
||||
*/
|
||||
nsresult ClearTrait(uint32_t aTrait);
|
||||
|
||||
protected:
|
||||
|
||||
/**
|
||||
* return the local corpus storage file for junk traits
|
||||
*/
|
||||
nsresult getTrainingFile(nsIFile ** aFile);
|
||||
|
||||
/**
|
||||
* return the local corpus storage file for non-junk traits
|
||||
*/
|
||||
nsresult getTraitFile(nsIFile ** aFile);
|
||||
|
||||
/**
|
||||
* read token strings from the data file
|
||||
*
|
||||
* @param stream file stream with token data
|
||||
* @param fileSize file size
|
||||
* @param aTraitId id for the trait whose counts will be read
|
||||
* @param aIsAdd true to add the counts, false to remove them
|
||||
*
|
||||
* @return true if successful, false if error
|
||||
*/
|
||||
bool readTokens(FILE* stream, int64_t fileSize, uint32_t aTraitId,
|
||||
bool aIsAdd);
|
||||
|
||||
/**
|
||||
* write token strings to the data file
|
||||
*/
|
||||
bool writeTokens(FILE* stream, bool shrink, uint32_t aTraitId);
|
||||
|
||||
/**
|
||||
* remove counts for a token string
|
||||
*/
|
||||
void remove(const char* word, uint32_t aTraitId, uint32_t aCount);
|
||||
|
||||
/**
|
||||
* add counts for a token string, adding the token string if new
|
||||
*/
|
||||
CorpusToken* add(const char* word, uint32_t aTraitId, uint32_t aCount);
|
||||
|
||||
/**
|
||||
* change counts in a trait in the traits array, adding the trait if needed
|
||||
*/
|
||||
nsresult updateTrait(CorpusToken* token, uint32_t aTraitId,
|
||||
int32_t aCountChange);
|
||||
nsCOMPtr<nsIFile> mTrainingFile; // file used to store junk training data
|
||||
nsCOMPtr<nsIFile> mTraitFile; // file used to store non-junk
|
||||
// training data
|
||||
nsTArray<TraitPerToken> mTraitStore; // memory for linked-list of counts
|
||||
uint32_t mNextTraitIndex; // index in mTraitStore to first empty
|
||||
// TraitPerToken
|
||||
nsTArray<uint32_t> mMessageCounts; // count of messages per trait
|
||||
// represented in the store
|
||||
nsTArray<uint32_t> mMessageCountsId; // Parallel array to mMessageCounts, with
|
||||
// the corresponding trait ID
|
||||
};
|
||||
|
||||
class nsBayesianFilter : public nsIJunkMailPlugin, nsIMsgCorpus,
|
||||
nsIObserver, nsSupportsWeakReference {
|
||||
public:
|
||||
NS_DECL_ISUPPORTS
|
||||
NS_DECL_NSIMSGFILTERPLUGIN
|
||||
NS_DECL_NSIJUNKMAILPLUGIN
|
||||
NS_DECL_NSIMSGCORPUS
|
||||
NS_DECL_NSIOBSERVER
|
||||
|
||||
nsBayesianFilter();
|
||||
|
||||
nsresult Init();
|
||||
|
||||
nsresult tokenizeMessage(const char* messageURI, nsIMsgWindow *aMsgWindow, TokenAnalyzer* analyzer);
|
||||
void classifyMessage(Tokenizer& tokens, const char* messageURI,
|
||||
nsIJunkMailClassificationListener* listener);
|
||||
|
||||
void classifyMessage(
|
||||
Tokenizer& tokenizer,
|
||||
const char* messageURI,
|
||||
nsTArray<uint32_t>& aProTraits,
|
||||
nsTArray<uint32_t>& aAntiTraits,
|
||||
nsIJunkMailClassificationListener* listener,
|
||||
nsIMsgTraitClassificationListener* aTraitListener,
|
||||
nsIMsgTraitDetailListener* aDetailListener);
|
||||
|
||||
void observeMessage(Tokenizer& tokens, const char* messageURI,
|
||||
nsTArray<uint32_t>& oldClassifications,
|
||||
nsTArray<uint32_t>& newClassifications,
|
||||
nsIJunkMailClassificationListener* listener,
|
||||
nsIMsgTraitClassificationListener* aTraitListener);
|
||||
|
||||
|
||||
protected:
|
||||
virtual ~nsBayesianFilter();
|
||||
|
||||
static void TimerCallback(nsITimer* aTimer, void* aClosure);
|
||||
|
||||
CorpusStore mCorpus;
|
||||
double mJunkProbabilityThreshold;
|
||||
int32_t mMaximumTokenCount;
|
||||
bool mTrainingDataDirty;
|
||||
int32_t mMinFlushInterval; // in milliseconds, must be positive
|
||||
//and not too close to 0
|
||||
nsCOMPtr<nsITimer> mTimer;
|
||||
|
||||
// index in mAnalysisStore for first empty AnalysisPerToken
|
||||
uint32_t mNextAnalysisIndex;
|
||||
// memory for linked list of AnalysisPerToken objects
|
||||
nsTArray<AnalysisPerToken> mAnalysisStore;
|
||||
/**
|
||||
* Determine the location in mAnalysisStore where the AnalysisPerToken
|
||||
* object for a particular token and trait is stored
|
||||
*/
|
||||
uint32_t getAnalysisIndex(Token& token, uint32_t aTraitIndex);
|
||||
/**
|
||||
* Set the value of the AnalysisPerToken object for a particular
|
||||
* token and trait
|
||||
*/
|
||||
nsresult setAnalysis(Token& token, uint32_t aTraitIndex,
|
||||
double aDistance, double aProbability);
|
||||
};
|
||||
|
||||
#endif // _nsBayesianFilter_h__
|
||||
|
|
@ -0,0 +1,22 @@
|
|||
/* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- */
|
||||
/* This Source Code Form is subject to the terms of the Mozilla Public
|
||||
* License, v. 2.0. If a copy of the MPL was not distributed with this
|
||||
* file, You can obtain one at http://mozilla.org/MPL/2.0/. */
|
||||
|
||||
#ifndef nsBayesianFilterCID_h__
|
||||
#define nsBayesianFilterCID_h__
|
||||
|
||||
#include "nsISupports.h"
|
||||
#include "nsIFactory.h"
|
||||
#include "nsIComponentManager.h"
|
||||
|
||||
#include "nsIMsgMdnGenerator.h"
|
||||
|
||||
#define NS_BAYESIANFILTER_CONTRACTID \
|
||||
"@mozilla.org/messenger/filter-plugin;1?name=bayesianfilter"
|
||||
#define NS_BAYESIANFILTER_CID \
|
||||
{ /* F1070BFA-D539-11D6-90CA-00039310A47A */ \
|
||||
0xF1070BFA, 0xD539, 0x11D6, \
|
||||
{ 0x90, 0xCA, 0x00, 0x03, 0x93, 0x10, 0xA4, 0x7A }}
|
||||
|
||||
#endif /* nsBayesianFilterCID_h__ */
|
||||
259
mailnews/extensions/bayesian-spam-filter/src/nsIncompleteGamma.h
Normal file
259
mailnews/extensions/bayesian-spam-filter/src/nsIncompleteGamma.h
Normal file
|
|
@ -0,0 +1,259 @@
|
|||
/* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
|
||||
/* This Source Code Form is subject to the terms of the Mozilla Public
|
||||
* License, v. 2.0. If a copy of the MPL was not distributed with this
|
||||
* file, You can obtain one at http://mozilla.org/MPL/2.0/. */
|
||||
|
||||
#ifndef nsIncompleteGamma_h__
|
||||
#define nsIncompleteGamma_h__
|
||||
|
||||
/* An implementation of the incomplete gamma functions for real
|
||||
arguments. P is defined as
|
||||
|
||||
x
|
||||
/
|
||||
1 [ a - 1 - t
|
||||
P(a, x) = -------- I t e dt
|
||||
Gamma(a) ]
|
||||
/
|
||||
0
|
||||
|
||||
and
|
||||
|
||||
infinity
|
||||
/
|
||||
1 [ a - 1 - t
|
||||
Q(a, x) = -------- I t e dt
|
||||
Gamma(a) ]
|
||||
/
|
||||
x
|
||||
|
||||
so that P(a,x) + Q(a,x) = 1.
|
||||
|
||||
Both a series expansion and a continued fraction exist. This
|
||||
implementation uses the more efficient method based on the arguments.
|
||||
|
||||
Either case involves calculating a multiplicative term:
|
||||
e^(-x)*x^a/Gamma(a).
|
||||
Here we calculate the log of this term. Most math libraries have a
|
||||
"lgamma" function but it is not re-entrant. Some libraries have a
|
||||
"lgamma_r" which is re-entrant. Use it if possible. I have included a
|
||||
simple replacement but it is certainly not as accurate.
|
||||
|
||||
Relative errors are almost always < 1e-10 and usually < 1e-14. Very
|
||||
small and very large arguments cause trouble.
|
||||
|
||||
The region where a < 0.5 and x < 0.5 has poor error properties and is
|
||||
not too stable. Get a better routine if you need results in this
|
||||
region.
|
||||
|
||||
The error argument will be set negative if there is a domain error or
|
||||
positive for an internal calculation error, currently lack of
|
||||
convergence. A value is always returned, though.
|
||||
|
||||
*/
|
||||
|
||||
#include <math.h>
|
||||
#include <float.h>
|
||||
|
||||
// the main routine
|
||||
static double nsIncompleteGammaP (double a, double x, int *error);
|
||||
|
||||
// nsLnGamma(z): either a wrapper around lgamma_r or the internal function.
|
||||
// C_m = B[2*m]/(2*m*(2*m-1)) where B is a Bernoulli number
|
||||
static const double C_1 = 1.0 / 12.0;
|
||||
static const double C_2 = -1.0 / 360.0;
|
||||
static const double C_3 = 1.0 / 1260.0;
|
||||
static const double C_4 = -1.0 / 1680.0;
|
||||
static const double C_5 = 1.0 / 1188.0;
|
||||
static const double C_6 = -691.0 / 360360.0;
|
||||
static const double C_7 = 1.0 / 156.0;
|
||||
static const double C_8 = -3617.0 / 122400.0;
|
||||
static const double C_9 = 43867.0 / 244188.0;
|
||||
static const double C_10 = -174611.0 / 125400.0;
|
||||
static const double C_11 = 77683.0 / 5796.0;
|
||||
|
||||
// truncated asymptotic series in 1/z
|
||||
static inline double lngamma_asymp (double z)
|
||||
{
|
||||
double w, w2, sum;
|
||||
w = 1.0 / z;
|
||||
w2 = w * w;
|
||||
sum = w * (w2 * (w2 * (w2 * (w2 * (w2 * (w2 * (w2 * (w2 * (w2
|
||||
* (C_11 * w2 + C_10) + C_9) + C_8) + C_7) + C_6)
|
||||
+ C_5) + C_4) + C_3) + C_2) + C_1);
|
||||
|
||||
return sum;
|
||||
}
|
||||
|
||||
struct fact_table_s
|
||||
{
|
||||
double fact;
|
||||
double lnfact;
|
||||
};
|
||||
|
||||
// for speed and accuracy
|
||||
static const struct fact_table_s FactTable[] = {
|
||||
{1.000000000000000, 0.0000000000000000000000e+00},
|
||||
{1.000000000000000, 0.0000000000000000000000e+00},
|
||||
{2.000000000000000, 6.9314718055994530942869e-01},
|
||||
{6.000000000000000, 1.7917594692280550007892e+00},
|
||||
{24.00000000000000, 3.1780538303479456197550e+00},
|
||||
{120.0000000000000, 4.7874917427820459941458e+00},
|
||||
{720.0000000000000, 6.5792512120101009952602e+00},
|
||||
{5040.000000000000, 8.5251613610654142999881e+00},
|
||||
{40320.00000000000, 1.0604602902745250228925e+01},
|
||||
{362880.0000000000, 1.2801827480081469610995e+01},
|
||||
{3628800.000000000, 1.5104412573075515295248e+01},
|
||||
{39916800.00000000, 1.7502307845873885839769e+01},
|
||||
{479001600.0000000, 1.9987214495661886149228e+01},
|
||||
{6227020800.000000, 2.2552163853123422886104e+01},
|
||||
{87178291200.00000, 2.5191221182738681499610e+01},
|
||||
{1307674368000.000, 2.7899271383840891566988e+01},
|
||||
{20922789888000.00, 3.0671860106080672803835e+01},
|
||||
{355687428096000.0, 3.3505073450136888885825e+01},
|
||||
{6402373705728000., 3.6395445208033053576674e+01}
|
||||
};
|
||||
#define FactTableLength (int)(sizeof(FactTable)/sizeof(FactTable[0]))
|
||||
|
||||
// for speed
|
||||
static const double ln_2pi_2 = 0.918938533204672741803; // log(2*PI)/2
|
||||
|
||||
/* A simple lgamma function, not very robust.
|
||||
|
||||
Valid for z_in > 0 ONLY.
|
||||
|
||||
For z_in > 8 precision is quite good, relative errors < 1e-14 and
|
||||
usually better. For z_in < 8 relative errors increase but are usually
|
||||
< 1e-10. In two small regions, 1 +/- .001 and 2 +/- .001 errors
|
||||
increase quickly.
|
||||
*/
|
||||
static double nsLnGamma (double z_in, int *gsign)
|
||||
{
|
||||
double scale, z, sum, result;
|
||||
*gsign = 1;
|
||||
|
||||
int zi = (int) z_in;
|
||||
if (z_in == (double) zi)
|
||||
{
|
||||
if (0 < zi && zi <= FactTableLength)
|
||||
return FactTable[zi - 1].lnfact; // gamma(z) = (z-1)!
|
||||
}
|
||||
|
||||
for (scale = 1.0, z = z_in; z < 8.0; ++z)
|
||||
scale *= z;
|
||||
|
||||
sum = lngamma_asymp (z);
|
||||
result = (z - 0.5) * log (z) - z + ln_2pi_2 - log (scale);
|
||||
result += sum;
|
||||
return result;
|
||||
}
|
||||
|
||||
// log( e^(-x)*x^a/Gamma(a) )
|
||||
static inline double lnPQfactor (double a, double x)
|
||||
{
|
||||
int gsign; // ignored because a > 0
|
||||
return a * log (x) - x - nsLnGamma (a, &gsign);
|
||||
}
|
||||
|
||||
static double Pseries (double a, double x, int *error)
|
||||
{
|
||||
double sum, term;
|
||||
const double eps = 2.0 * DBL_EPSILON;
|
||||
const int imax = 5000;
|
||||
int i;
|
||||
|
||||
sum = term = 1.0 / a;
|
||||
for (i = 1; i < imax; ++i)
|
||||
{
|
||||
term *= x / (a + i);
|
||||
sum += term;
|
||||
if (fabs (term) < eps * fabs (sum))
|
||||
break;
|
||||
}
|
||||
|
||||
if (i >= imax)
|
||||
*error = 1;
|
||||
|
||||
return sum;
|
||||
}
|
||||
|
||||
static double Qcontfrac (double a, double x, int *error)
|
||||
{
|
||||
double result, D, C, e, f, term;
|
||||
const double eps = 2.0 * DBL_EPSILON;
|
||||
const double small =
|
||||
DBL_EPSILON * DBL_EPSILON * DBL_EPSILON * DBL_EPSILON;
|
||||
const int imax = 5000;
|
||||
int i;
|
||||
|
||||
// modified Lentz method
|
||||
f = x - a + 1.0;
|
||||
if (fabs (f) < small)
|
||||
f = small;
|
||||
C = f + 1.0 / small;
|
||||
D = 1.0 / f;
|
||||
result = D;
|
||||
for (i = 1; i < imax; ++i)
|
||||
{
|
||||
e = i * (a - i);
|
||||
f += 2.0;
|
||||
D = f + e * D;
|
||||
if (fabs (D) < small)
|
||||
D = small;
|
||||
D = 1.0 / D;
|
||||
C = f + e / C;
|
||||
if (fabs (C) < small)
|
||||
C = small;
|
||||
term = C * D;
|
||||
result *= term;
|
||||
if (fabs (term - 1.0) < eps)
|
||||
break;
|
||||
}
|
||||
|
||||
if (i >= imax)
|
||||
*error = 1;
|
||||
return result;
|
||||
}
|
||||
|
||||
static double nsIncompleteGammaP (double a, double x, int *error)
|
||||
{
|
||||
double result, dom, ldom;
|
||||
// domain errors. the return values are meaningless but have
|
||||
// to return something.
|
||||
*error = -1;
|
||||
if (a <= 0.0)
|
||||
return 1.0;
|
||||
if (x < 0.0)
|
||||
return 0.0;
|
||||
*error = 0;
|
||||
if (x == 0.0)
|
||||
return 0.0;
|
||||
|
||||
ldom = lnPQfactor (a, x);
|
||||
dom = exp (ldom);
|
||||
// might need to adjust the crossover point
|
||||
if (a <= 0.5)
|
||||
{
|
||||
if (x < a + 1.0)
|
||||
result = dom * Pseries (a, x, error);
|
||||
else
|
||||
result = 1.0 - dom * Qcontfrac (a, x, error);
|
||||
}
|
||||
else
|
||||
{
|
||||
if (x < a)
|
||||
result = dom * Pseries (a, x, error);
|
||||
else
|
||||
result = 1.0 - dom * Qcontfrac (a, x, error);
|
||||
}
|
||||
|
||||
// not clear if this can ever happen
|
||||
if (result > 1.0)
|
||||
result = 1.0;
|
||||
if (result < 0.0)
|
||||
result = 0.0;
|
||||
return result;
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
Loading…
Add table
Add a link
Reference in a new issue