doctor/rappaurio-sae501_502
1
0
Fork 0
mirror of https://scm.univ-tours.fr/22107988t/rappaurio-sae501_502.git synced 2025-08-29 01:55:58 +02:00
Code Issues Packages Projects Releases Wiki Activity
Files
223677a8b44638baec500ab318dd6c4bb99f1efa
rappaurio-sae501_502/app/node_modules/cejs/data/math.js
22107988t a94f68f22a v1.0 du site web
2023-09-25 13:27:24 +02:00

4315 lines
120 KiB
JavaScript
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

/**
* @name CeL function for mathematics.
* @fileoverview 本檔案包含了數學演算相關的 functions。
*
* TODO: 方程式圖形顯示 by SVG
*
* @see http://www.wolframalpha.com/<br />
* http://www.numberempire.com/
*/
// More examples: see /_test suite/test.js
'use strict';
// 'use asm';
// --------------------------------------------------------------------------------------------
// 不採用 if 陳述式，可以避免 Eclipse JSDoc 與 format 多縮排一層。
typeof CeL === 'function' && CeL.run({
// module name
name : 'data.math',
require : 'data.code.compatibility.|data.native.set_bind',
// 設定不匯出的子函式。
// no_extend : '*',
// 為了方便格式化程式碼，因此將 module 函式主體另外抽出。
code : module_code
});
function module_code(library_namespace) {
// requiring
var set_bind = this.r('set_bind');
var has_bigint = library_namespace.env.has_bigint;
/**
* null module constructor
*
* @class 數學相關的 functions
*/
var _// JSDT:_module_
= function() {
// null module constructor
};
/**
* for JSDT: 有 prototype 才會將之當作 Class
*/
_// JSDT:_module_
.prototype = {};
/**
* <code>
數位
十分位 tenths digit
整數 whole number
</code>
*/
// ---------------------------------------------------------------------//
// basic constants. 定義基本常數。
var
/**
* empty product, or nullary product, 乘法單位元素.<br />
* number * MULTIPLICATIVE_IDENTITY === number.<br />
* 2/2, 3/3, ..
*
* MULTIPLICATIVE_IDENTITY = 1
*
* @type {Number}
* @constant
*
* @see https://en.wikipedia.org/wiki/Identity_element
* https://en.wikipedia.org/wiki/Empty_product
*/
MULTIPLICATIVE_IDENTITY = 1 / 1,
/**
* Any nonzero number raised by the exponent 0 is 1.<br />
* (any number) ^ 0 === Math.pow(number, 0) === ZERO_EXPONENT<br />
* Math.pow(2, 0), Math.pow(3, 0), ..
*
* ZERO_EXPONENT = 1
*
* @type {Number}
* @constant
*
* @see https://en.wikipedia.org/wiki/Exponentiation
*/
ZERO_EXPONENT = Math.pow(1, 0),
/**
* absorbing element, zero element.<br />
* number * ABSORBING_ELEMENT === ABSORBING_ELEMENT<br />
* Math.pow(2, 0), Math.pow(3, 0), ..
*
* @type {Number}
* @constant
*
* @see https://en.wikipedia.org/wiki/Absorbing_element
*/
ABSORBING_ELEMENT = 0,
/**
* multiplication sign. e.g., '⋅', '*', '×'.
*
* @type {String}
* @constant
*
* @see https://en.wikipedia.org/wiki/Multiplication_sign
* https://en.wikipedia.org/wiki/Interpunct
*/
MULTIPLICATION_SIGN = '⋅',
/**
* default base = 10.<br />
* 內定：10位數。應與 parseInt() 一致。
*
* @type {Natural}
* @constant
*/
DEFAULT_BASE = parseInt('10'),
/**
* The biggest integer we can square. 超過此數則無法安全操作平方。
*
* @type {Natural}
* @constant
*/
sqrt_max_integer = Math.sqrt(Number.MAX_SAFE_INTEGER) | 0;
// ---------------------------------------------------------------------//
/**
* use Horner's method to calculate the value of polynomial.
*
* @param {Array}coefficients
* coefficients of polynomial.<br />
* coefficients: [ degree 0, degree 1, degree 2, ... ]
* @param {Number}variable
* value of (x)
*
* @returns {Number} the value of polynomial
*
* @see https://en.wikipedia.org/wiki/Horner%27s_method
*/
function polynomial_value(coefficients, variable) {
return coefficients.reduceRight(function(value, coefficient) {
return value * variable + coefficient;
});
}
_.polynomial_value = polynomial_value;
_// JSDT:_module_
.
/**
* 輾轉相除 n1/n2 或 小數 n1/1 轉成 整數/整數。
*
* @param {Natural}n1
* number 1
* @param {Natural}[n2]
* number 2
* @param {Natural}times
* maximum times 次數, 1,2,..
*
* @return {Array} 連分數序列 (continued fraction) ** 負數視 _.mutual_division.done
* 而定!
*/
mutual_division = function mutual_division(n1, n2, times) {
var q = [], c;
if (isNaN(times) || times <= 0)
times = 80;
if (!n2 || isNaN(n2))
n2 = 1;
if (!Number.isInteger(n1)) {
c = n1;
var i = 9, f = n2;
while (i--) {
// 以整數運算比較快！這樣會造成整數多4%，浮點數多1/3倍的時間，但仍值得。
f *= DEFAULT_BASE;
c *= DEFAULT_BASE;
if (Number.isInteger(c)) {
n1 = c;
n2 = f;
break;
}
}
}
// 連分數負數之處理。更沒問題的: (n1 < 0?1:0) ^ (n2 < 0?1:0)
if (_.mutual_division.mode && ((n1 < 0) ^ (n2 < 0))) {
// 使兩數皆為正
if (n2 < 0)
n2 = -n2;
else
n1 = -n1;
q.push(-(1 + (n1 - (c = n1 % n2)) / n2));
n1 = n2;
n2 -= c;
}
// old:
if (false) {
while (b && n--) {
// 2.08s@10000
// 可能因為少設定（=）一次c所以較快。但（若輸入不為整數）不確保d為整數？用Math.floor((a-(c=a%b))/b)可確保，速度與下式一樣快。
c = a % b;
d.push((a - c) / b);
a = b;
b = c;
// 2.14s@10000:mutual_division(.142857)
// d.push(c=Math.floor(a/b)),c=a-b*c,a=b,b=c;
// 2.2s@10000
// d.push(Math.floor(a/b)),b=a%(c=b),a=c;
}
if (n)
d.push(0);
}
// 2.4s@10000
// 可能因為少設定（=）一次c所以較快。但（若輸入不為整數）不確保d為整數？用Math.floor((a-(c=a%b))/b)可確保，速度與下式一樣快。
while (times--)
if (n2) {
c = n1 % n2;
q.push((n1 - c) / n2);
n1 = n2;
n2 = c;
} else {
// [ ... , done mark, (最後非零的餘數。若原 n1, n2 皆為整數，則此值為
// GCD。但請注意:這邊是已經經過前面為了以整數運算，增加倍率過的數值!!) ]
q.push(_.mutual_division.done, n1);
// library_namespace.debug('done: ' + q);
break;
}
/**
* <code>
// 2.26s@10000
while(b&&n--)if(d.push((a-(c=a%b))/b),a=b,!(b=c)){d.push(0);break;}
var m=1;c=1;while(m&&n--)d.push(m=++c%2?b?(a-(a%=b))/b:0:a?(b-(b%=a))/a:0);//buggy
</code>
*/
return q;
};
_// JSDT:_module_
.mutual_division.done = -7;// ''
_// JSDT:_module_
.
/**
* !!mode:連分數處理，對負數僅有最初一數為負。
*/
mutual_division.mode = 0;
_// JSDT:_module_
.
/**
* 取得連分數序列的數值。
*
* @param {Array}sequence
* 序列
* @param {Natural}[max_no]
* maximum no. 取至第 max_no 個
*
* @return {Array}連分數序列的數值
*
* @requires mutual_division.done
*/
continued_fraction = function(sequence, max_no) {
if (!Array.isArray(sequence) || !sequence.length)
return sequence;
if (sequence.at(-2) === _.mutual_division.done)
sequence.length -= 2;
if (sequence.length < 1)
return sequence;
if (!max_no/* || max_no < 2 */|| max_no > sequence.length)
max_no = sequence.length;
var a, b;
if (max_no % 2) {
b = 1;
a = 0;
} else {
a = 1;
b = 0;
}
if (false) {
sequence[max_no++] = 1;
if (--max_no % 2) {
b = sequence[max_no];
a = s[--max_no];
} else {
a = sequence[max_no];
b = sequence[--max_no];
}
}
if (false)
library_namespace.debug('a=' + a + ', b=' + b + ', max_no='
+ max_no);
while (max_no--)
if (max_no % 2)
b += a * sequence[max_no];
else
a += b * sequence[max_no];
if (false)
library_namespace.debug('a=' + a + ', b=' + b);
return [ a, b ];
};
// quadratic (m√r + i) / D → continued fraction [... , [period ...]]
// Rosen, Kenneth H. (2011). Elementary Number Theory and its Applications
// (6th edition). Boston: Pearson Addison-Wesley. pp. 508–511.
// https://en.wikipedia.org/wiki/Periodic_continued_fraction
// https://en.wikipedia.org/wiki/Square_root_of_2
// https://en.wikipedia.org/wiki/Square_root#As_periodic_continued_fractions
// https://en.wikipedia.org/wiki/Generalized_continued_fraction#Roots_of_positive_numbers
function quadratic_to_continued_fraction(r, m, i, D) {
if (r < 0) {
throw 'The root is negative!';
}
if (!i)
i = 0;
if (!D)
D = 1;
if (!m)
m = 1;
else if (m < 0) {
m = -m;
i = -i;
D = -D;
}
// (m√r + i) / D
// = (√(r m^2) + i) / D
// = (√(d in book) + P0) / Q0
var d = m * m * r,
//
P = i, Q = D,
// A: α in book.
A, a, sequence = [], ptr = sequence, start_PQ;
// Be sure Q0 | (d - P0^2)
if ((d - P * P) % Q !== 0)
P *= Q, d *= Q * Q, Q *= Q;
// assert: now: Q0 | (d - P0^2)
for (var sqrt = Math.sqrt(d), t;;) {
if (start_PQ) {
if (P === start_PQ[0] && Q === start_PQ[1])
return sequence;
} else if (0 < (t = sqrt - P) && t < Q) {
// test if α is purely periodic.
start_PQ = [ P, Q ];
sequence.push(ptr = []);
}
ptr.push(a = Math.floor(A = (sqrt + P) / Q));
library_namespace.debug(((sequence === ptr ? 0
: sequence.length - 1)
+ ptr.length - 1)
+ ': P='
+ P
+ ', Q='
+ Q
+ ', α≈'
+ (DEFAULT_BASE * A | 0)
/ DEFAULT_BASE + ', a=' + a, 3);
// set next Pn = a(n-1)Q(n-1) - P(n-1), Qn = (d - Pn^2) / Q(n-1).
P = a * Q - P;
Q = (d - P * P) / Q;
if (Q === 0)
// is not a quadratic irrationality?
return sequence;
// assert: Pn, Qn are both integers.
}
}
_.quadratic_to_continued_fraction = quadratic_to_continued_fraction;
// get the first solution of Pell's equation: x^2 - d y^2 = 1 or -1.
// https://en.wikipedia.org/wiki/Pell%27s_equation
// Rosen, Kenneth H. (2005). Elementary Number Theory and its Applications
// (5th edition). Boston: Pearson Addison-Wesley. pp. 542-545.
function solve_Pell(d, n, NO) {
// TODO
// use CeL.data.math.quadratic.solve_Pell instead
;
}
// _.solve_Pell = solve_Pell;
_// JSDT:_module_
.
/**
* The best rational approximation. 取得值最接近之有理數 (use 連分數 continued fraction),
* 取近似值. c.f., 調日法 在分子或分母小於下一個漸進分數的分數中，其值是最接近精準值的近似值。
*
* @param {Number}number
* number
* @param {Number}[rate]
* 比例在 rate 以上
* @param {Natural}[max_no]
* maximum no. 最多取至序列第 max_no 個 TODO : 並小於 l: limit
*
* @return {Array}[分子, 分母, 誤差]
*
* @requires mutual_division,continued_fraction
* @see https://en.wikipedia.org/wiki/Continued_fraction#Best_to_rational_numbers
*/
to_rational_number = function(number, rate, max_no) {
if (!rate)
// This is a magic number: 我們無法準確得知其界限為何。
rate = 65536;
var d = _
.mutual_division(number, 1, max_no && max_no > 0 ? max_no : 20), i = 0, a, b = d[0], done = _.mutual_division.done;
if (!b)
b = d[++i];
while (++i < d.length && (a = d[i]) !== done)
if (a / b < rate)
b = a;
else
break;
if (false)
library_namespace.debug(number
+ ' '
+
// 連分數表示 (continued fraction)
(d.length > 1 && d.at(-2) === _.mutual_division.done ? '='
+ ' [<em>'
+ d[0]
+ ';'
+ d.slice(1, i).join(', ')
+ '</em>'
+ (i < d.length - 2 ? ', '
+ d.slice(i, -2).join(', ') : '')
+ '] ... ' + d.slice(-1)
:
// 約等於的符號是≈或≒，不等於的符號是≠。
// https://zh.wikipedia.org/wiki/%E7%AD%89%E4%BA%8E
'≈'
+ ' [<em>'
+ d[0]
+ ';'
+ d.slice(1, i).join(', ')
+ '</em>'
+ (i < d.length ? ', '
+ d.slice(i).join(', ') : '')
+ ']: ' + d.length + ',' + i + ',' + d[i]));
d = _.continued_fraction(d, i);
if (d[1] < 0) {
d[0] = -d[0];
d[1] = -d[1];
}
if (false)
library_namespace.debug('→ ' + d[0] + '/' + d[1]);
// [ {Integer}±numerator, {Natural}denominator ]
return [ d[0], d[1], d[0] / d[1] - number ];
};
// 正規化帶分數。 to_mixed_fraction
// 2019/7/10 14:22:6
// mixed_fraction =
// [ {Integer}±whole, {Integer}±numerator, {Natural|Undefined}denominator ]
function normalize_mixed_fraction(mixed_fraction) {
if (typeof mixed_fraction === 'number') {
// treat as float
mixed_fraction = [ mixed_fraction ];
}
var whole = mixed_fraction[0] || ABSORBING_ELEMENT;
var numerator = mixed_fraction[1] || ABSORBING_ELEMENT;
var denominator = mixed_fraction[2] || MULTIPLICATIVE_IDENTITY;
// {Natural}denominator >= 1
if (denominator < 0) {
denominator = -denominator;
numerator = -numerator;
}
// {Natural}denominator ∈ ℤ
if (!Number.isInteger(denominator)) {
// assert: is float
denominator = _.to_rational_number(denominator);
numerator *= typeof numerator === 'bigint' ? BigInt(denominator[1])
: denominator[1];
denominator = denominator[0];
}
// {Integer}±numerator ∈ ℤ
if (!Number.isInteger(numerator)) {
// assert: is float
numerator = _.to_rational_number(numerator);
denominator *= typeof denominator === 'bigint' ? BigInt(numerator[1])
: numerator[1];
numerator = numerator[0];
}
// assert: {Natural}denominator, {Integer}±numerator
// TODO: 約分here。
var using_bigint;
// {Integer}±whole ∈ ℤ
if (typeof whole === 'number' && !Number.isInteger(whole)) {
// assert: whole is float
whole = _.to_rational_number(whole);
var LCM = _.LCM(denominator, whole[1]);
if (typeof LCM === 'bigint') {
using_bigint = true;
// convert all numbers to the same type.
numerator = BigInt(numerator);
denominator = BigInt(denominator);
whole = whole.map(BigInt);
}
numerator = numerator
// (LCM / denominator) === GCD * whole[1]
* (LCM / denominator)
+ (whole[0] < 0 ? -(-whole[0] % whole[1]) : whole[0]
% whole[1]) * (LCM / whole[1]);
denominator = LCM;
if (using_bigint) {
whole = whole[0] / whole[1];
} else {
whole = whole[0] < 0 ? -Math.floor(-whole[0] / whole[1]) : Math
.floor(whole[0] / whole[1]);
}
}
// TODO: convert all numbers to the same type.
// whole, numerator 必須同符號。
if (whole * numerator < 0) {
if (whole < 0) {
whole++;
// numerator > 0
numerator -= denominator;
} else {
whole--;
// numerator < 0
numerator += denominator;
}
}
// 處理假分數。同時會處理絕對值為整數之問題。
if (Math.absolute(numerator) >= denominator) {
whole += using_bigint ? numerator / denominator : Math
.floor(numerator / denominator);
numerator %= denominator;
}
// 約分。
if (numerator == ABSORBING_ELEMENT) {
// normalize
denominator = MULTIPLICATIVE_IDENTITY;
} else {
var GCD = _.GCD(numerator, denominator);
if (GCD >= 2) {
if (using_bigint)
GCD = BigInt(GCD);
numerator /= GCD;
denominator /= GCD;
}
}
// export
mixed_fraction = Object.assign([ whole, numerator, denominator ], {
valueOf : mixed_fraction_valueOf,
toString : mixed_fraction_toString
});
return mixed_fraction;
}
function mixed_fraction_valueOf() {
if (!this[1])
return this[0];
return this[0] + this[1] / this[2];
}
function mixed_fraction_toString() {
if (!this[1])
return String(this[0]);
if (!this[0])
return this[1] + '/' + this[2];
return this[0] + (this[1] < 0 ? '' : '+') + this[1] + '/' + this[2];
}
_.normalize_mixed_fraction = normalize_mixed_fraction;
// ------------------------------------------------------------------------
// 正規化數字成 integer 或 bigint
// 在大量計算前，盡可能先轉換成普通 {Number} 以加快速度。
// cohandler(may convert to number)
function to_int_or_bigint(value, cohandler) {
var number;
if (typeof value === 'bigint') {
number = Number(value);
if (Number.isSafeInteger(number)) {
cohandler && cohandler(true);
return number;
} else {
cohandler && cohandler(false);
return value;
}
}
// 這方法無法準確處理像 `1e38/7`, `10/7` 這樣的情況。
if (typeof value === 'number') {
number = Math.round(value);
if (!Number.isSafeInteger(number)) {
throw new RangeError('Cannot convert number ' + value
+ ' to safe integer!');
}
cohandler && cohandler(true);
return Math.round(number);
}
number = parseInt(value);
if (Number.isSafeInteger(number)) {
cohandler && cohandler(true);
return number;
}
if (!has_bigint)
throw new RangeError('Cannot convert ' + number
+ ' to safe integer!');
cohandler && cohandler(false);
return BigInt(value);
}
// Let all elements of {Array}this the same type: int, else bigint.
// 可能的話應該將絕對值最大的數字放在前面，早點判別出是否需要用 {BigInt}。
function array_to_int_or_bigint() {
// assert: {Array}this
// cache int values
if (this.some(function(value, index) {
value = to_int_or_bigint(value);
this[index] = value;
return typeof value === 'bigint';
}, this)) {
// must using bigint
this.forEach(function(value, index) {
this[index] = BigInt(value);
}, this);
}
// assert: all elements of `this` is in the same type.
// typeof this[0] === typeof this[1]
return this;
}
// 可用於 {BigInt} 之 Math.abs
// https://en.wikipedia.org/wiki/Absolute_value
function absolute(value) {
return value < 0 ? -value : value;
}
Math.absolute = absolute;
/**
* 求多個數之 GCD(Greatest Common Divisor, 最大公因數/公約數).<br />
* Using Euclidean algorithm(輾轉相除法).<br />
*
* TODO: 判斷互質.
*
* @param {Integers}number_array
* number array
*
* @returns {Natural} GCD of the numbers specified
*/
function GCD(number_array) {
if (arguments.length > 1) {
// Array.from()
number_array = Array.prototype.slice.call(arguments);
}
// 正規化數字。
number_array = number_array.map(function(value) {
return Math.absolute(to_int_or_bigint(value));
})
// 由小至大排序可以減少計算次數?? 最起碼能夠延後使用 {BigInt} 的時機。
.sort(library_namespace.general_ascending)
// .unique_sorted()
;
// console.log(number_array);
// 不在此先設定 gcd = number_array[0]，是為了讓每個數字通過資格檢驗。
var index = 0, length = number_array.length, gcd = 0, remainder, number;
// assert: 所有數字皆已先轉換成數字，並已轉為絕對值。
while (index < length) {
number = number_array[index++];
if (number >= 1) {
gcd = number;
break;
}
}
// console.log(gcd);
while (index < length && 2 <= gcd) {
number = number_array[index++];
if (!(number >= 1))
continue;
if (typeof number === 'bigint') {
number %= BigInt(gcd);
// [ gcd, number ] = [ gcd, number ].to_int_or_bigint();
remainder = [ gcd, number ].to_int_or_bigint();
gcd = remainder[0];
number = remainder[1];
}
// assert: typeof gcd === typeof number
// console.log([ gcd, number ]);
// Euclidean algorithm 輾轉相除法。
while ((remainder = number % gcd) >= 1) {
number = gcd;
// 使用絕對值最小的餘數。為了要處理 {BigInt}，因此不採用 Math.min()。
// gcd = Math.min(remainder, gcd - remainder);
gcd = gcd - remainder < remainder ? gcd - remainder : remainder;
}
}
if (typeof gcd === 'bigint'
&& Number.isSafeInteger(number = Number(gcd))) {
gcd = number;
}
return gcd;
}
_// JSDT:_module_
.GCD = GCD;
_// JSDT:_module_
.
/**
* 求多個數之 LCM(Least Common Multiple, 最小公倍數): method 1.<br />
* Using 類輾轉相除法.<br />
*
* @param {Integers}number_array
* number array
*
* @returns {Natural} LCM of the numbers specified
*/
LCM = function LCM(number_array) {
if (arguments.length > 1) {
// Array.from()
number_array = Array.prototype.slice.call(arguments);
}
// 正規化數字。
number_array = number_array.map(function(value) {
return Math.absolute(to_int_or_bigint(value));
})
// .sort().reverse()
;
if (number_array.some(function(number) {
return number == 0;
})) {
// 允許 0:
return 0;
}
var lcm = number_array[0];
for (var index = 1, length = number_array.length; index < length; index++) {
var number = number_array[index];
// assert: {Integer}number
var gcd = _.GCD(number, lcm);
if (typeof number === typeof gcd) {
number /= gcd;
if (typeof number === typeof lcm) {
gcd = lcm * number;
if (Number.isSafeInteger(gcd)) {
lcm = gcd;
} else if (has_bigint) {
lcm = BigInt(lcm) * BigInt(number);
} else {
throw new RangeError('LCM is not safe integer!');
}
} else {
// assert: {BigInt}number or {BigInt}lcm
lcm = BigInt(lcm) * BigInt(number);
}
} else {
// assert: {BigInt}number, {Number}gcd
lcm = BigInt(lcm) * (number / BigInt(gcd));
}
}
return lcm;
};
_// JSDT:_module_
.
/**
* 求多個數之 LCM(Least Common Multiple, 最小公倍數): method 1.<br />
* Using 類輾轉相除法.<br />
*
* TODO: 更快的方法： 短除法? 一次算出 GCD, LCM?
*
* @param {Integers}number_array
* number array
*
* @returns {Natural} LCM of the numbers specified
*/
LCM3 = function LCM3(number_array) {
if (arguments.length > 1) {
// Array.from()
number_array = Array.prototype.slice.call(arguments);
}
// 正規化數字。
number_array = number_array.map(function(value) {
return Math.absolute(to_int_or_bigint(value));
})
// .sort().reverse()
;
if (number_array.some(function(number) {
return number == 0;
})) {
// 允許 0:
return 0;
}
var lcm = number_array[0], using_bigint;
for (var index = 1, length = number_array.length; index < length; index++) {
var number = number_array[index];
// assert: {Integer}number
if (typeof number !== typeof lcm) {
// assert: number, lcm 有一個是 bigint。
using_bigint = true;
lcm = BigInt(lcm);
number = BigInt(number);
}
// console.log([ lcm, number ]);
var number0 = number;
var lcm0 = lcm;
// 倒反版的 Euclidean algorithm 輾轉相除法.
// 反覆讓兩方各自加到比對方大的倍數，當兩者相同時，即為 lcm。
while (lcm !== number) {
// console.log([ lcm0, number0, lcm, number ]);
if (lcm > number) {
var remainder = -lcm % number0;
if (remainder) {
number = lcm + remainder + number0;
if (!using_bigint && has_bigint
&& !Number.isSafeInteger(number)) {
using_bigint = true;
number0 = BigInt(number0);
lcm0 = BigInt(lcm0);
number = BigInt(lcm + remainder) + number0;
lcm = BigInt(lcm);
}
} else {
// number0 整除 lcm: 取 lcm 即可.
break;
}
} else {
var remainder = -number % lcm0;
if (remainder) {
lcm = number + remainder + lcm0;
if (!using_bigint && has_bigint
&& !Number.isSafeInteger(lcm)) {
using_bigint = true;
number0 = BigInt(number0);
lcm0 = BigInt(lcm0);
lcm = BigInt(number + remainder) + BigInt(lcm0);
number = BigInt(number);
}
} else {
// lcm0 整除 number: 取 number 即可.
lcm = number;
break;
}
}
}
}
return lcm;
};
_// JSDT:_module_
.
/**
* 求多個數之 LCM(Least Common Multiple, 最小公倍數): method 2.<br />
* Using 類輾轉相除法.<br />
*
* @param {Integers}number_array
* number array
*
* @returns {Integer} LCM of the numbers specified
*/
LCM2 = function LCM2(number_array) {
if (arguments.length > 1) {
// Array.from()
number_array = Array.prototype.slice.call(arguments);
}
var i = 0, l = number_array.length, lcm = 1, r, n, num, gcd;
for (; i < l && lcm; i++) {
// 每個數字都要做運算，雖可確保正確，但沒有效率!
if (!isNaN(num = n = Math.abs(parseInt(number_array[i])))) {
gcd = lcm;
// Euclidean algorithm.
while (r = n % gcd)
n = gcd, gcd = r;
lcm = num / gcd * lcm;
}
}
return lcm;
};
/**
* Get <a href="https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm"
* accessdate="2013/8/3 19:45">Extended Euclidean algorithm</a><br />
*
* @param {Integer}n1
* number 1
* @param {Integer}n2
* number 2
* @returns [ GCD, m1, m2 ]: GCD = m1 * n1 + m2 * n2
*
* @see division_with_remainder() @ data.math
* @since 2013/8/3 20:24:30
*/
function extended_GCD(n1, n2) {
var remainder, quotient, using_g1 = false, using_bigint,
// 前一group [dividend 應乘的倍數, divisor 應乘的倍數]
m1g1 = 1, m2g1 = 0;
if (typeof n1 === 'bigint' || typeof n2 === 'bigint') {
// convert all numbers to the same type.
n1 = BigInt(n1);
n2 = BigInt(n2);
m1g1 = BigInt(m1g1);
m2g1 = BigInt(m2g1);
using_bigint = true;
}
// 前前group [dividend 應乘的倍數, divisor 應乘的倍數]
var m1g2 = /* 0 */m2g1, m2g2 = /* 1 */m1g1;
while (remainder = n1 % n2) {
quotient = (n1 - remainder) / n2;
if (!using_bigint) {
// assert: typeof quotient === 'number'
quotient = Math.floor(quotient);
}
// 現 group = remainder = 前前group - quotient * 前一group
if (using_g1 = !using_g1)
m1g1 -= quotient * m1g2, m2g1 -= quotient * m2g2;
else
m1g2 -= quotient * m1g1, m2g2 -= quotient * m2g1;
// swap numbers
n1 = n2;
n2 = remainder;
}
return using_g1 ? [ n2, m1g1, m2g1 ] : [ n2, m1g2, m2g2 ];
}
// extended GCD algorithm
_.EGCD = extended_GCD;
/**
* 帶餘除法 Euclidean division。<br />
* 除非設定 closest，否則預設 remainder ≥ 0.
*
* @param {Number}dividend
* 被除數。
* @param {Number}divisor
* 除數。
* @param {Boolean}[closest]
* get the closest quotient
*
* @returns {Array} [ {Integer}quotient 商, {Number}remainder 餘數 ]
*
* @see http://stackoverflow.com/questions/14997165/fastest-way-to-get-a-positive-modulo-in-c-c
* @see extended_GCD() @ data.math
*
* @since 2015/10/31 10:4:45
*/
function division_with_remainder(dividend, divisor, closest) {
if (false)
return [ Math.floor(dividend / divisor),
// 轉正。保證餘數值非負數。
(dividend % divisor + divisor) % divisor ];
var remainder = dividend % divisor;
if (closest) {
if (remainder != 0
// 0 !== 0n
&& Math.absolute(remainder + remainder) > Math.absolute(divisor))
if (remainder < 0)
remainder += Math.absolute(divisor);
else
remainder -= Math.absolute(divisor);
} else if (remainder < 0) {
// assert: (-0 < 0) === false
remainder += Math.absolute(divisor);
}
dividend = (dividend - remainder) / divisor;
if (typeof dividend === 'number') {
dividend = Math.round(dividend);
} else {
// assert: typeof dividend === 'bigint'
}
return [ dividend, remainder ];
}
// 帶餘數除法 division with remainder
_.division = division_with_remainder;
/**
* 取得所有分母為 denominator，分子分母互質的循環小數的循環節位數。<br />
* Repeating decimal: get period (repetend length)
*
* @param {Natural}denominator
* 分母
* @param {Boolean}with_transient
* 亦取得非循環節部分位數
* @param {Natural}min
* 必須最小長度，在測試大量數字時使用。若發現長度必小於 min 則即時跳出。效果不俗 (test
* Euler_26(1e7))。
*
* @returns {Array}[{Number}period length 循環節位數 < denominator,
* {Number}transient 非循環節部分位數 ]
*
* @see https://en.wikipedia.org/wiki/Repeating_decimal#Reciprocals_of_composite_integers_coprime_to_10
*/
function period_length(denominator, with_transient, min) {
// 去除所有 2 或 5 的因子。
var non_repeating = 0, non_repeating_5 = 0;
while (denominator % 5 === 0)
denominator /= 5, non_repeating_5++;
while (denominator % 2 === 0)
denominator /= 2, non_repeating++;
if (non_repeating < non_repeating_5)
non_repeating = non_repeating_5;
if (denominator === 1 || denominator <= min)
return with_transient ? [ 0, non_repeating ] : 0;
for (var length = 1, remainder = 1;; length++) {
remainder = remainder * DEFAULT_BASE % denominator;
if (remainder === 1)
return with_transient ? [ length, non_repeating ] : length;
}
}
_.period_length = period_length;
// ---------------------------------------------------------------------//
/**
* 從數集 set 中挑出某些數，使其積最接近指定的數 target。<br />
* To picks some numbers from set, so the product is approximately the
* target number.
*
* TODO: improve/optimize
*
* @param {Array}set
* number set of {Natural}
* @param {Natural}target
* target number
* @param {Object}[options]
* 附加參數/設定特殊功能與選項
*
* @returns {Array}某些數，其積最接近 target。
*
* @see http://stackoverflow.com/questions/19572043/given-a-target-sum-and-a-set-of-integers-find-the-closest-subset-of-numbers-tha
*/
function closest_product(set, target, options) {
var status, minor_data;
if (Array.isArray(options)) {
status = options;
minor_data = status[0];
options = minor_data.options;
} else {
// 初始化
if (!options)
options = new Boolean;
else if (typeof options === 'number')
options = {
direction : options
};
else if (typeof options === 'boolean')
options = {
sorted : options
};
minor_data = [ Infinity ];
minor_data.options = options;
// status = [ [minor, set of minor], product, set of product ]
status = [ minor_data, ZERO_EXPONENT, [] ];
if (!options.sorted)
set = set.clone()
// 由小至大排序。
.sort(library_namespace.ascending);
}
// direction = -1: 僅接受小於 target 的積。
// direction = +1: 僅接受大於 target 的積。
var direction = options.direction,
//
product = status[1], selected = status[2];
set.some(function(natural, index) {
if (selected[index])
// 已經處理過，跳過。
return;
var _product = product * natural, _selected,
/** {Number}差 ≥ 0 */
difference = Math.abs(target - _product),
// 是否發現新極小值。採用 minor_data 而不 cache 是因為此間 minor_data 可能已經改變。
_status = difference <= minor_data[0];
library_namespace.debug(target + '=' + (target - _product) + '+'
+ natural + '×' + product + ', ' + product + '='
+ (set.filter(function(n, index) {
return selected[index];
}).join('⋅') || 1), 6, 'closest_product');
if (target < _product) {
library_namespace.debug('target < _product, direction: '
+ direction, 6, 'closest_product');
if (!_status || direction < 0) {
library_namespace.debug('積已經過大，之後不會有合適的。', 5,
'closest_product');
return true;
}
}
_selected = selected.clone();
_selected[index] = true;
if (_status && (!(direction > 0) || target <= _product)) {
_status = set.filter(function(n, index) {
return _selected[index];
}).join(closest_product.separator);
if (difference === minor_data[0]) {
if (minor_data.includes(_status)) {
library_namespace.debug('已經處理過相同的，跳過。', 5,
'closest_product');
return;
}
minor_data.push(_status);
} else {
minor_data.clear();
minor_data.push(difference, _status);
}
library_namespace.debug('發現極小值:' + target + '=' + difference
+ '+' + natural + '×' + product + ', ' + product + '='
+ (set.filter(function(n, index) {
return selected[index];
}).join('⋅') || 1), 3, 'closest_product');
}
_status = [ minor_data, _product, _selected ];
library_namespace.debug('繼續探究是否有更小的差:' + _status.join(';'), 4,
'closest_product');
closest_product(set, target, _status);
});
return minor_data.length > 1 && minor_data;
}
closest_product.separator = MULTIPLICATION_SIGN;
_.closest_product = closest_product;
// TODO:將數列分為積最接近的兩組。
/**
* Get <a
* href="https://en.wikipedia.org/wiki/Modular_multiplicative_inverse"
* accessdate="2013/8/3 20:10">modular multiplicative inverse</a> (模反元素)
*
* TODO:<br />
* untested!
*
* @param {Integer}number
* number
* @param {Integer}modulo
* modulo
*
* @returns {Integer} modular multiplicative inverse
*
* @since 2013/8/3 20:24:30
*/
function modular_inverse(number, modulo) {
number = extended_GCD(number, modulo);
if (number[0] == 1)
return (number = number[1]) < 0 ? number + modulo : number;
}
_.modular_inverse = modular_inverse;
// factorial_cache[ n ] = n!
// factorial_cache = [ 0! = 1, 1!, 2!, ... ]
var factorial_cache = [ 1 ], factorial_cache_to;
/**
* Get the factorial (階乘) of (natural).<br />
*
* @param {ℕ⁰:Natural+0}natural
* safe integer. 0–18
*
* @returns {Natural}natural的階乘.
*
* @see https://en.wikipedia.org/wiki/Factorial
*/
function factorial(natural) {
var length = factorial_cache.length;
if (length <= natural && !factorial_cache_to) {
var f = factorial_cache[--length];
while (length++ < natural)
if (isFinite(f *= length))
factorial_cache.push(f);
else {
factorial_cache_to = length - 1;
break;
}
}
return natural < length ? factorial_cache[natural] : Infinity;
}
// var factorial_map = CeL.math.factorial.map(9);
// generate factorial map
factorial.map = function(natural) {
if (!natural)
natural = 9;
if (factorial_cache.length <= natural && !factorial_cache_to)
factorial(natural);
return factorial_cache.slice(0, natural + 1);
};
_.factorial = factorial;
// ---------------------------------------------------------------------//
/**
* http://www.math.umbc.edu/~campbell/NumbThy/Class/Programming/JavaScript.html
* http://aoki2.si.gunma-u.ac.jp/JavaScript/
*/
/**
* 得到開方數，相當於 Math.floor(Math.sqrt(number)) === Math.sqrt(number) | 0. get
* integer square root. TODO: use 牛頓法
*
* @param {Number}
* positive number
*
* @return r, r^2 ≤ number < (r+1)^2
*
* @see <a href="http://www.azillionmonkeys.com/qed/sqroot.html"
* accessdate="2010/3/11 18:37">Paul Hsieh's Square Root page</a><br />
* <a
* href="http://www.embeddedrelated.com/usenet/embedded/show/114789-1.php"
* accessdate="2010/3/11 18:34">Suitable Integer Square Root Algorithm
* for 32-64-Bit Integers on Inexpensive Microcontroller? |
* Comp.Arch.Embedded | EmbeddedRelated.com</a>
*/
function floor_sqrt(number) {
// return Math.sqrt(number) | 0;
if (!Number.isFinite(number = Math.floor(number)))
return;
var g = 0, v, h, t;
while ((t = g << 1) < (v = number - g * g)) {
// library_namespace.debug(t + ', ' + v);
h = 1;
while (h * (h + t) <= v)
// 因為型別轉關係，還是保留 << 而不用 *2
h <<= 1;// h *= 2;
g += h >> 1;// h / 2;//
}
if (false)
library_namespace.debug('end: ' + t + ', ' + v);
return g;
}
_.floor_sqrt = floor_sqrt;
// count digits of integer: using .digit_length()
function ceil_log(number, base) {
if (!number)
return 0;
if (!base)
base = DEFAULT_BASE;
// assert: base >= 2, base === (base | 0)
number = Math.abs(number);
// ideal
return Math.ceil(base === 10 ? Math.log10(number) : base === 2 ? Math
.log2(number)
// TODO: base = 2^n
: Math.log(number) / Math.log(base));
// slow... should use multiply by exponents
// Logarithm
var log = 0;
if (number < ZERO_EXPONENT) {
while (number < ZERO_EXPONENT) {
number *= base;
if (false)
library_namespace.debug(number);
log--;
}
if (number !== ZERO_EXPONENT)
// 修正。
log++;
} else {
while (number > ZERO_EXPONENT) {
// 因為可能損失 base^exp + (...) 之剩餘部分，因此不能僅採用 Math.floor(number /
// base)
// 但如此較費時。
number /= base;
if (false)
library_namespace.log(number);
log++;
}
}
return log;
}
// add binding
_.ceil_log = ceil_log;
/** {Object}all possible last 2 digits of square number */
var square_ending = Object.create(null);
[ 0, 1, 4, 9, 16, 21, 24, 25, 29, 36, 41, 44, 49, 56, 61, 64, 69, 76, 81,
84, 89, 96 ].forEach(function(n) {
square_ending[n] = null;
});
// Squarity testing
// 檢測 ({Natural}number) 是否為完全平方數
// a square number or perfect square. TODO: use 牛頓法
// is square number, n²
function is_square(number) {
// 快速判定 possible_square(number)
// https://www.johndcook.com/blog/2008/11/17/fast-way-to-test-whether-a-number-is-a-square/
// 0x0213 = parseInt('1111110111101100', 2).toString(0x10)
if (0xFDEC & (1 << (number & 0xF))) {
return false;
}
// TRUE only if number % 16 === 0, 1, 4, 9
// %16 有4個: http://oeis.org/A023105
if (!((number % 100) in square_ending)) {
return false;
}
number = Math.sqrt(number);
return number === (number | 0) && number;
// another method
// https://en.wikipedia.org/wiki/Methods_of_computing_square_roots
// https://gmplib.org/manual/Perfect-Square-Algorithm.html
var sqrt = floor_sqrt(number);
return sqrt * sqrt === number && sqrt;
}
_.is_square = is_square;
// 檢測 ({Natural}f1 * {Natural}f2) 是否為完全平方數
function product_is_square(f1, f2) {
if (f1 === f2) {
return true;
// e.g., r * p^2, r * p^2
}
// e.g., p^2 * r, q^2 * r
var product = f1 * f2;
if (Number.isSafeInteger(product)) {
return is_square(product);
}
// 除法不比較快。
if (f1 > f2) {
// swap
var tmp = f1;
f1 = f2;
f2 = tmp;
}
if (!Number.isSafeInteger(f2)) {
library_namespace.error('The number ' + f2
+ ' is NOT a safe number!');
}
// assert: f1 < f2
if (f2 % f1 === 0) {
// e.g., r, r * p^2
return is_square(f2 / f1);
}
if (is_square(f1)) {
// e.g., p^2, q^2
return is_square(f2);
}
if (is_square(f2)) {
return false;
}
for (var index = 0, length = Math.min(10, primes.length); index < length; index++) {
var p = primes[index], p2 = p * p;
tmp = false;
while (f1 % p2 === 0) {
f1 /= p2;
tmp = true;
}
while (f2 % p2 === 0) {
f2 /= p2;
tmp = true;
}
if (f1 % p === 0) {
while (f1 % p === 0 && f2 % p === 0) {
f1 /= p;
f2 /= p;
tmp = true;
}
}
if (tmp) {
product = f1 * f2;
if (Number.isSafeInteger(product)) {
return is_square(product);
}
}
}
// 找GCD。較慢，但沒辦法。
var gcd = GCD(f1, f2);
return is_square(f1 / gcd)
//
&& is_square(f2 / gcd);
}
_.product_is_square = product_is_square;
/**
* <code>
n(n+1)/2=T, n∈ℕ, n=?
n(n+1)/2=T, n=?
n = 1/2 (sqrt(8 T+1)-1)
Reduce[(n (1 + n))/2 == T, n]
n = 1/2 (sqrt(8 T+1)-1)
Reduce[n(3n−1)/2==P, n]
n = 1/6 (sqrt(24 P+1)+1)
// hexagonal
Reduce[n(2n−1)==H, n]
n = 1/4 (sqrt(8 H+1)+1)
</code>
*/
function is_triangular(natural) {
// https://en.wikipedia.org/wiki/Triangular_number
var sqrt = is_square(8 * natural + 1);
return sqrt && sqrt % 2 === 1;
}
_.is_triangular = is_triangular;
function is_generalized_pentagonal(generalized) {
// https://en.wikipedia.org/wiki/Pentagonal_number
return is_square(24 * generalized + 1);
}
_.is_generalized_pentagonal = is_generalized_pentagonal;
function is_pentagonal(natural) {
// https://en.wikipedia.org/wiki/Pentagonal_number
var sqrt = is_square(24 * natural + 1);
return sqrt && sqrt % 6 === 5;
}
_.is_pentagonal = is_pentagonal;
// 素勾股數 primitive Pythagorean triple
var 素勾股數 = [], last_Pythagorean_m = 2;
/**
* primitive Pythagorean triples 素勾股數組/素商高數組/素畢氏三元數
*
* @param {Natural}limit
* limit of m. 若欲改成 limit of 斜邊，請輸入斜邊長後自行 filter。
*
* @returns {Array}primitive Pythagorean triple list
*
* @see https://en.wikipedia.org/wiki/Pythagorean_triple#Generating_a_triple
* The triple generated by Euclid's formula is primitive if and only if
* m and n are coprime and m − n is odd.
*/
function Pythagorean_list(limit) {
if (last_Pythagorean_m < limit) {
for (var m = last_Pythagorean_m; m < limit; m++) {
for (var m2 = m * m, m_2 = 2 * m, n, n2 = n = m % 2 === 0 ? 1
: 0;
// 設 m > n 互質且均是正整數，m 和 n 有一個是偶數，
// 計算出來的 (a, b, c) 就是素勾股數。所有素勾股數可用列式找出
n < m; n += 2, n2 = n * n) {
if (GCD(m, n) === 1) {
var a = m2 - n2, b = m_2 * n, c = m2 + n2;
// let a < b < c
素勾股數.push(a < b ? [ a, b, c ] : [ b, a, c ]);
}
}
}
last_Pythagorean_m = limit;
}
return 素勾股數;
}
_.Pythagorean_list = Pythagorean_list;
// ---------------------------------------------------------------------//
// Catalan_number[0] = 1
var Catalan_number_list = [ 1 ];
// Catalan numbers
// @see https://en.wikipedia.org/wiki/Catalan_number
function Catalan_number(NO) {
if (NO < Catalan_number_list.length) {
// use cache
return Catalan_number_list[NO];
}
var n = Catalan_number_list.length - 1,
//
this_Catalan_number = Catalan_number_list[n];
for (; n < NO; n++) {
this_Catalan_number = this_Catalan_number * (4 * n + 2) / (n + 2);
Catalan_number_list.push(this_Catalan_number);
}
return this_Catalan_number;
}
_.Catalan_number_list = Catalan_number_list;
_.Catalan_number = Catalan_number;
// ---------------------------------------------------------------------//
/** {Array}Collatz_conjecture_steps[number] = steps. cache 以加快速度。 */
var Collatz_conjecture_steps_cache = [ , 1 ];
if (false) {
// 此法費時 1.5 倍， 12s → 19s
Collatz_conjecture_steps_cache = new Array(1000001);
Collatz_conjecture_steps_cache[1] = 1;
}
// assert: Collatz_conjecture_steps_cache[1] === 1 (因程式判別方法需要此項)
// Collatz conjecture
// https://en.wikipedia.org/wiki/Collatz_conjecture
function Collatz_conjecture(natural) {
if (!(natural > 0))
return;
var chain = [ natural ];
while (natural > 1) {
chain.push(natural % 2 === 0 ? natural /= 2
: (natural = natural * 3 + 1));
}
// Collatz_conjecture_steps_cache[natural] = chain.length;
// return all terms
return chain;
}
// 為計算 steps 特殊化。
// assert: CeL.Collatz_conjecture.steps(natural) ===
// CeL.Collatz_conjecture(natural).length
function Collatz_conjecture_steps(natural) {
if (!(natural > 0))
return;
var chain = [];
while (!(natural in Collatz_conjecture_steps_cache)) {
chain.push(natural);
if (natural % 2 === 0)
natural /= 2;
else
natural = natural * 3 + 1;
}
var steps = Collatz_conjecture_steps_cache[natural] + chain.length, s = steps;
// 紀錄 steps。
chain.forEach(function(natural) {
Collatz_conjecture_steps_cache[natural] = s--;
});
return steps;
}
/**
* <code>
backwards 反向:
1000000: 153 steps
999999: 259 steps
999667: 290 steps
999295: 396 steps
997823: 440 steps
970599: 458 steps
939497: 507 steps
837799: 525 steps
</code>
*/
// search the longest chain / sequence below ((natural))
function Collatz_conjecture_longest(natural) {
if (!(natural > 0))
return;
// maximum steps
var max_steps = 0, max_steps_natural;
// brute force
for (var n = 1, steps, _n; n <= natural; n++) {
if (n in Collatz_conjecture_steps_cache)
steps = Collatz_conjecture_steps_cache[_n = n];
else {
steps = Collatz_conjecture_steps(_n = n);
// 預先快速處理所有 2倍數字。採用此方法，約可增加 5% 速度。不採用此方法，n 正反向速度差不多。
while (_n * 2 <= natural) {
Collatz_conjecture_steps_cache[_n *= 2] = ++steps;
}
}
if (max_steps < steps) {
library_namespace.debug(natural + ': ' + steps + ' steps', 3,
'Collatz_conjecture.longest');
max_steps = steps;
max_steps_natural = _n;
}
}
return [ max_steps_natural, max_steps ];
}
_.Collatz_conjecture = Collatz_conjecture;
Collatz_conjecture.steps = Collatz_conjecture_steps;
Collatz_conjecture.longest = Collatz_conjecture_longest;
// ---------------------------------------------------------------------//
// https://en.wikipedia.org/wiki/Memoization
/** {Array}質數列表。 cache / memoization 以加快速度。 */
var primes = [ 2, 3, 5 ],
/**
* last prime tested.<br />
* assert: last_prime_tested is ((6n ± 1)). 因此最起碼應該從 5 開始。
*
* @type {Natural}
*/
last_prime_tested = primes.at(-1);
// https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
// the sieve of Eratosthenes 篩法
function prime_sieve(limit, limit_index) {
// var list = _.number_array(limit + 1, 0, Int8Array);
var list = new Array(limit + 1);
// 重建 re-build list (table)
primes.forEach(function(prime) {
for (var number = prime; number <= limit;) {
// list[number += prime] = 1;
list[number += prime] = true;
}
});
for (var n = last_prime_tested; n <= limit;) {
if (list[++n])
continue;
// n is prime
// library_namespace.debug(n + ' is prime');
primes.push(n);
if (limit_index && primes.length > limit_index)
break;
// 登記所有倍數。
for (var number = n; number <= limit;) {
// list[number += prime] = 1;
list[number += prime] = true;
}
}
last_prime_tested = primes.at(-1);
return primes;
}
_.prime_sieve = prime_sieve;
// integer: number to test
function test_is_prime(integer, index, sqrt) {
// assert: Number.isInteger(integer), integer ≥ 0
index |= 0;
if (!sqrt)
sqrt = floor_sqrt(integer);
// 採用試除法, use trial division。
// 從第一個質數一直除到 ≤ sqrt(integer) 之質數
for (var prime, length = primes.length; index < length;) {
if (integer % (prime = primes[index++]) === 0)
// return: prime factor found
return integer === prime ? false : prime;
if (sqrt < prime)
return false;
}
// 質數列表中的質數尚無法檢測 integer。
}
/**
* Get the prime[index] or prime list.
*
* @param {Natural}[index]
* prime index starts from 1
* @param {Natural}[limit]
* the upper boundary of prime value
*
* @returns {Natural}prime value
*/
function prime(index, limit) {
if (!(index > 0)) {
if (limit > 0) {
index = prime_pi(limit);
return primes.slice(0, index);
}
return primes;
}
if (primes.length < index) {
if (false && index - primes.length > 1e6) {
// using the sieve of Eratosthenes 篩法
// 沒比較快。
// 詳細數量應採 prime π(x)。
// https://zh.wikipedia.org/wiki/%E8%B3%AA%E6%95%B8%E5%AE%9A%E7%90%86
prime_sieve(limit || index * 10, index);
} else {
// assert: last_prime_tested is ((6n ± 1))
/**
* {Boolean}p1 === true: last_prime_tested is 6n+1.<br />
* else: last_prime_tested is 6n-1
*/
var p1 = last_prime_tested % 6 === 1;
for (; primes.length < index
&& Number.isSafeInteger(last_prime_tested);) {
last_prime_tested += (p1 = !p1) ? 2 : 4;
// 實質為 https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
if (!test_is_prime(last_prime_tested, 2))
primes.push(last_prime_tested);
if (limit && limit <= last_prime_tested)
break;
}
library_namespace.debug('last prime tested = '
+ last_prime_tested, 2, 'prime');
}
}
return primes[index - 1];
}
_.prime = prime;
// prime #5484598 = 94906249, the biggest prime <
// Math.sqrt(Number.MAX_SAFE_INTEGER) - 1.
// the 2nd biggest prime is 94906247.
// CeL.prime(CeL.prime_pi(Number.MAX_SAFE_INTEGER = 2^53 - 1)) =
// 9007199254740881
function prime_pi(value) {
value = Math.floor(Math.abs(value));
if (last_prime_tested < value)
prime(value, value);
// +1: index of function prime() starts from 1!
return primes.search_sorted(value, true) + 1;
}
_.prime_pi = prime_pi;
/**
* Get the primorial (質數階乘, p_n#) of (NO).<br />
*
* @param {Natural}NO
* safe integer. 1–13
*
* @returns {Natural}p_NO的質數階乘.
*
* @see https://en.wikipedia.org/wiki/Primorial
*/
function primorial(NO) {
if (!(NO >= 1))
return MULTIPLICATIVE_IDENTITY;
prime(NO);
var index = 0, product = MULTIPLICATIVE_IDENTITY;
while (index < NO)
product *= primes[index++];
return product;
}
/**
* Get the primorial (質數階乘, n#) of (natural).<br />
*
* @param {Natural}natural
* safe integer. 2–42
*
* @returns {Natural}natural的質數階乘.
*
* @see https://en.wikipedia.org/wiki/Primorial
*/
function primorial_natural(natural) {
// 2: primes[0]
if (!(natural >= 2))
return MULTIPLICATIVE_IDENTITY;
var index = 0, length = prime_pi(natural), product = MULTIPLICATIVE_IDENTITY;
while (index < length)
product *= primes[index++];
return product;
}
_.primorial = primorial;
primorial.natural = primorial_natural;
// return multiplicand × multiplier % modulus
// assert: 三者皆為 natural number, and Number.isSafeInteger() is OK.
// max(multiplicand, multiplier) < modulus. 否則會出現錯誤!
function multiply_modulo(multiplicand, multiplier, modulus) {
var quotient = multiplicand * multiplier;
if (Number.isSafeInteger(quotient))
return quotient % modulus;
// 避免 overflow
if (multiplicand > multiplier)
quotient = multiplicand, multiplicand = multiplier,
multiplier = quotient;
if (quotient === 1)
throw new Error('Please use data.math.integer instead!');
quotient = Math.floor(modulus / multiplicand);
quotient = (multiplicand * (multiplier % quotient) - Math
.floor(multiplier / quotient)
* (modulus % multiplicand))
% modulus;
return quotient;
}
_.multiply_modulo = multiply_modulo;
// return integer ^ exponent % modulus
// assert: 三者皆為 natural number, and Number.isSafeInteger() is OK. 否則會出現錯誤!
function power_modulo(integer, exponent, modulus) {
for (var remainder = 1, power = integer % modulus;;) {
if (exponent % 2 === 1)
remainder = multiply_modulo(remainder, power, modulus);
if ((exponent >>= 1) === 0)
return remainder;
if ((power = multiply_modulo(power, power, modulus)) === 1)
return remainder;
}
}
_.power_modulo = power_modulo;
function power_modulo(natural, exponent, modulus) {
var remainder = 1;
for (natural %= modulus; exponent > 0; natural = natural * natural
% modulus, exponent >>= 1)
if (exponent % 2 === 1)
remainder = remainder * natural % modulus;
return remainder;
}
_.power_modulo = power_modulo;
// Miller–Rabin primality test
// return true: is composite, undefined: probable prime (PRP) / invalid
// number
// https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
function Miller_Rabin(natural, times) {
if (natural % 2 === 0)
return natural !== 2;
if (!(natural < sqrt_max_integer) || natural < 2)
return;
var n_1 = natural - 1, d = n_1, s = 0, a, x;
do {
s++, d /= 2;
} while (d % 2 === 0);
// assert: s > 0
if (!(times |= 0))
times = 3;
for (var prime_index = 0; prime_index < times;) {
// 3rd 起用大數( > 307) 偵測。
x = power_modulo(prime(++prime_index + (prime_index < 3 ? 0 : 60)),
d, natural);
if (x === 1 || x === n_1)
continue;
var i = 1, j = 1;
for (; i < s; i++) {
x = x * x % natural;
if (x === 1)
return true;
if (x === n_1) {
j = 0;
break;
}
}
if (j)
// composite
return true;
}
// probable prime to base 2, 3.
// probable prime (PRP), probably prime
// https://en.wikipedia.org/wiki/Probable_prime
}
_.Miller_Rabin = Miller_Rabin;
/**
* Test if ((natural)) is not prime. 是否不為質數。
*
* 對大數，僅能確定為合數，不是質數；不能保證是質數。
*
* @param {Natural}natural
* natural number to test
*
* @returns true: is composite.<br />
* false: is prime.<br />
* prime: min prime factor. The least prime divisor of ((natural)).<br />
* undefined: probable prime (PRP) / invalid number.
*/
function not_prime(natural) {
if (!Number.isSafeInteger(natural) || natural < 2)
return true;
var result;
if (false) {
var sqrt = floor_sqrt(natural = result);
result = 0;
while ((result = test_is_prime(natural, result, sqrt)) === undefined) {
// 多取一些質數。
prime((result = primes.length) + 1);
}
}
// warming up. 為 Miller_Rabin() 暖身。
prime(70);
// 先從耗費少的檢測開始。
// 先檢測此數是否在質數列表中。
if (natural <= last_prime_tested)
// -1: NOT_FOUND
return primes.search_sorted(natural) === -1;
result = primes.length < 1e3
&& last_prime_tested * last_prime_tested < natural
// ↑ 1e3: 當有太多質數要測，test_is_prime()就不划算了。
? undefined : test_is_prime(natural);
if (result === undefined)
result = Miller_Rabin(natural);
if (result === undefined) {
if (primes.length < 1e5)
// 多取一些質數。一般說來，產生這表的速度頗快。
prime(1e5);
result = test_is_prime(natural);
}
return result;
}
_.not_prime = not_prime;
function Pollards_rho_1980(natural) {
if (natural % 2 === 0)
return 2;
if (natural % 3 === 0)
return 3;
if (!(natural < sqrt_max_integer) || natural < 2)
return;
// reset initial value
// assert: natural > max(y, m, c)
var x = natural - 4, y = 2 + (Math.random() * x | 0), m = 2 + (Math
.random()
* x | 0), c = 2 + (Math.random() * x | 0), r = 1, q = 1, i, k, ys, G;
do {
for (x = y, i = r; i--;)
// f(x) = (x^2 + c) % N
y = (y * y % natural + c) % natural;
k = 0;
do {
for (ys = y, i = Math.min(m, r - k); i--;) {
// f(x) = (x^2 + c) % N
y = (y * y % natural + c) % natural;
q = q * Math.abs(x - y) % natural;
}
G = GCD(q, natural);
k += m;
} while (k < r && G === 1);
r *= 2;
// TODO: 當 r 過大，例如為十位數以上之質數時，過於消耗時間。
} while (G === 1);
if (natural === G)
do {
// f(x) = (x^2 + c) % N
ys = (ys * ys % natural + c) % natural;
G = GCD(Math.abs(x - ys), natural);
} while (G === 1);
return natural === G && G;
}
_.Pollards_rho = Pollards_rho_1980;
// ---------------------------------------------------------------------//
// CeL.factorize(natural > 1).toString(exponentiation_sign,
// multiplication_sign)
function factors_toString(exponentiation_sign, multiplication_sign) {
if (!exponentiation_sign && !Number.prototype.to_super)
exponentiation_sign = true;
if (!multiplication_sign)
// https://en.wikipedia.org/wiki/Multiplication_sign
multiplication_sign = '⋅';
else if (multiplication_sign === true)
// https://en.wikipedia.org/wiki/Interpunct
multiplication_sign = '×';
// others: '*'
var list = [];
for ( var factor in this) {
if (this[factor] > ZERO_EXPONENT)
// expand exponentiation
factor += exponentiation_sign === true ? (multiplication_sign + factor)
.repeat(this[factor] - ZERO_EXPONENT)
// https://en.wikipedia.org/wiki/Exponentiation#In_programming_languages
// exponentiation_sign: ^, **, ↑, ^^, ⋆
: exponentiation_sign ? exponentiation_sign
+ this[factor]
// https://en.wikipedia.org/wiki/Prime_factor
// To shorten prime factorizations, factors are often
// expressed in powers (multiplicities).
: this[factor].to_super();
list.push(factor);
}
return list.join(multiplication_sign);
}
// 計算所有因數個數。
// 計算所有質因數個數: Object.keys(factors).length
// 質因數列表: Object.keys(factors)
function count_all_factors() {
var count = MULTIPLICATIVE_IDENTITY;
for ( var prime in this) {
// exponent
count *= this[prime] + 1;
}
// re-define count
Object.defineProperty(this, 'count', {
enumerable : false,
value : count
});
return count;
}
// 歐拉函數 φ(n), Euler's totient function, Euler's phi function
// 是小於或等於n的正整數中與n互質的數的數目。
// https://en.wikipedia.org/wiki/Euler's_totient_function
function coprime() {
var count = this.natural;
for ( var prime in this) {
count = count / prime * (prime - 1);
}
// re-define coprime
Object.defineProperty(this, 'coprime', {
enumerable : false,
value : count
});
return count;
}
/**
* 取得某數的質因數分解，整數分解/因式分解/素因子分解, prime factorization, get floor factor.<br />
* 唯一分解定理(The Unique Factorization Theorem)告訴我們素因子分解是唯一的，這即是稱為算術基本定理 (The
* Fundamental Theorem of Arithmeric) 的數學金科玉律。<br />
*
* use Object.keys(factors) to get primes
*
* @param {Natural}natural
* integer number ≥ 2
* @param {Natural}radix
* output radix
* @param {Natural}index
* start prime index
*
* @return {Object}prime factors { prime1:power1, prime2:power2, ... }
*
* @requires floor_sqrt
*
* @see <a href="http://homepage2.nifty.com/m_kamada/math/10001.htm"
* accessdate="2010/3/11 18:7">Factorizations of 100...001</a>
*/
function factorize(natural, radix, index, factors) {
if (!Number.isSafeInteger(natural) || natural < 2
/**
* javascript 可以表示的最大整數值 = 10^21-2^16-1 = 999999999999999934463
*
* @see http://www.highdots.com/forums/javascript/how-js-numbers-represented-internally-166538-4.html
*/
// && !(1 < (natural = Math.floor(Math.abs(natural))) && natural <
// 999999999999999934469)
)
return;
if (!radix)
radix = undefined;
try {
radix.toString(radix);
} catch (e) {
// IE8?
radix = DEFAULT_BASE;
}
if (!factors)
factors = Object.create(null);
Object.defineProperties(factors, {
natural : {
enumerable : false,
value : natural
},
toString : {
enumerable : false,
value : factors_toString
},
count : {
enumerable : false,
configurable : true,
get : count_all_factors
},
coprime : {
enumerable : false,
configurable : true,
get : coprime
}
});
index |= 0;
var p = 1, sqrt = floor_sqrt(natural);
for (var power, length = primes.length; p <= sqrt;)
// 採用試除法, use trial division。
if (natural % (p = index < length ? primes[index++]
// find enough primes
: prime(++index)) === 0) {
for (power = 1; (natural /= p) % p === 0;)
power++;
factors[p.toString(radix)] = power;
sqrt = floor_sqrt(natural);
}
if (1 < natural)
factors[natural.toString(radix)] = 1;
return factors;
// 為了獲得確實結果，在 (Number.MAX_SAFE_INTEGER) 範圍內不採用 Pollard's rho。
// 事實上，若加上 (natural < sqrt_max_integer) 的限制，
// 一般說來在此範圍內使用 Pollard's rho 亦不切實際。
if (sqrt < p) {
// assert: natural is now prime.
if (1 < natural)
factors[natural.toString(radix)] = 1;
return factors;
}
var fA = [], fac = function(i) {
if (sqrt_max_integer <= natural || not_prime(i)) {
var p, count = 3;
while (count-- && !(p = Pollards_rho_1980(i)))
;
if (p) {
fac(p);
fac(i / p);
return;
} else
library_namespace
.warn('factorize: 無法分解'
+ (library_namespace.is_debug()
&& Miller_Rabin(i) ? '合數' : '')
+ '因子 [' + i.toString()
+ ']；您或許有必要自行質因數分解此數！');
}
fA.push(i);
};
fac(natural);
if (Array.isArray(factors)) {
// TODO
}
fA.sort(library_namespace.ascending);
fA.forEach(function(p) {
p = p.toString(radix);
if (p in factors)
factors[p]++;
else
factors[p] = 1;
});
return factors;
}
factorize._toString = factors_toString;
_.factorize = factorize;
// test
(function() {
function count(n) {
var a = factorize(n), s = '', v = 1;
if (a) {
for ( var i in a) {
s += '*' + i + (a[i] > 1 ? '^' + a[i] : '');
v *= Math.pow(i, a[i]);
}
s = s.substr(1) + '=' + v + '=' + n;
} else
s = 'error! ' + n;
document.getElementById('result').value += s
+ '\n-------------------------------------------\n';
}
});
/**
* 得到第一個質數因子。
*
* @param {Natural}natural
* natural number ≥ 2
*
* @returns{Natural}the first prime factor of ((natural))
*/
function first_factor(natural) {
for (var p = 1, sqrt = floor_sqrt(natural), index = 0, length = primes.length; p <= sqrt;)
// 採用試除法, use trial division。
if (natural % (p = index < length ? primes[index++]
// find enough primes
: prime(++index)) === 0)
return p;
return natural;
}
_.first_factor = first_factor;
/**
* Get the summation map 1–limit of proper factors.
*
* A proper factor of a positive integer n is a factor of n other than 1 or
* n (Derbyshire 2004, p. 32).<br />
* A positive divisor of n which is different from n is called a proper
* divisor or an aliquot part of n.
*
* @param {Natural}limit
* 處理到哪個數字。include limit itself.
* @param {Object}[options]
* 附加參數/設定特殊功能與選項
*
* @returns {Natural}summation of proper factors
*
* @see http://mathworld.wolfram.com/ProperFactor.html
*/
function factor_sum_map(limit, options) {
var add_1, add_self,
// default: {Natural}∑ summation of proper factors
get_sum = true, processor, list;
if (options) {
if (typeof options === 'function')
processor = options.processor;
add_1 = options.add_1;
add_self = options.add_self;
list = options.list;
if (options.all_factors)
add_1 = add_self = true;
processor = options.processor;
if (typeof processor === 'string' && factor_sum_map[processor])
processor = factor_sum_map[processor];
if (typeof processor !== 'function') {
processor = undefined;
// options.list: get factor list instead of summation.
get_sum = !options.get_list;
}
}
// assert: limit ≥ 1
// ++limit: number up to ((limit)), but need ((limit+1)) elements.
++limit;
var
// ((index)) starts from 2.
// skip 0: needless, natural numbers starts from 1.
// skip 1: already precessed by .fill(1).
index = list ? 0 : 2,
// options.add_1: every number has factor 1,
// set this if you want include 1 into sum.
factor_map = get_sum ? _.number_array(limit, add_1 ? 1 : 0,
options.type)
//
: add_1 ? [ , [ 1 ] ] : [];
if (false)
// factor_map[0] is nonsense 無意義，預設成 0。
factor_map[0] = 0;
if (options && typeof options.preprocessor === 'function')
factor_map = options.preprocessor(factor_map) || factor_map;
// generate factor map: a kind of sieve method 篩法.
// https://en.wikipedia.org/wiki/Sieve_theory
for (;; index++) {
var number = list ? list[index] : index;
if (!(number < limit))
break;
for (var n = add_self ? number : 2 * number; n < limit; n += number) {
// 處理所有 ((number)) 之倍數。
if (processor)
// processor : function(factor_map, factor, natural) {;}
processor(factor_map, number, n);
else if (get_sum)
// 將所有 ((number)) 之倍數都加上 ((number))。
// Append ((number)) to every multiple of ((number)).
//
// factor_map[0] is nonsense 無意義
// factor_map[number>0]
// = summation of the proper factors of number.
factor_map[n] += number;
else if (n in factor_map)
factor_map[n].push(number);
else
factor_map[n] = add_1 ? [ 1, number ] : [ number ];
}
}
library_namespace.debug('factor map: [' + factor_map.length + '] '
+ factor_map.slice(0, 30).join(';') + '...', 1,
'factor_sum_map');
return factor_map;
}
_.factor_sum_map = factor_sum_map;
// count factors
factor_sum_map.count = function(factor_map, factor, natural) {
factor_map[natural]++;
};
/**
* count coprime numbers below ((limit)).<br />
* 計算所有比各數字小，並與各數字互質的數。<br />
* 歐拉函數 φ(n), Euler's totient function, Euler's phi function
* 是小於或等於n的正整數中與n互質的數的數目。
*
* @param {Natural}limit
* natural number ≥ 2
*
* @returns {Array}coprime map
*
* @see function coprime()<br />
* https://en.wikipedia.org/wiki/Euler's_totient_function
*
*/
function coprime_map(limit, options) {
return factor_sum_map(limit, {
add_self : true,
// 初始化。
preprocessor : function(factor_map) {
factor_map.forEach(function(v, i) {
factor_map[i] = i;
});
},
list : prime(0, limit),
processor : function(factor_map, prime, number) {
factor_map[number] = factor_map[number] / prime * (prime - 1);
}
});
}
_.coprime_map = coprime_map;
/**
* Get perfect number list.
*
* @param {Natural}limit
* 處理到哪個數字。include limit itself.
* @param {Number}type
* type>0: abundant number.<br />
* type<0: deficient number.<br />
* default: perfect number
*
* @returns {Array}number list
*/
function perfect_numbers(limit, type) {
var numbers = [], factor_map = factor_sum_map(limit, {
add_1 : true
});
// skip 0: needless, natural numbers starts from 1.
// 僅對過剩數才需要做此處置。
if (type > 0)
factor_map[0] = 0;
// 雖然設定 add_1，但對 1 本身，應該為 0。
factor_map[1] = 0;
factor_map.forEach(type > 0 ? function(factor_sum, index) {
// abundant number or excessive number. 過剩數又稱作豐數或盈數
// https://en.wikipedia.org/wiki/Abundant_number
if (factor_sum > index)
numbers.push(index);
} : type < 0 ? function(factor_sum, index) {
// deficient or deficient number. 虧數又稱作缺數
// https://en.wikipedia.org/wiki/Deficient_number
if (factor_sum < index)
numbers.push(index);
} : function(factor_sum, index) {
// perfect numbers. 完全數，又稱完美數或完備數
// https://en.wikipedia.org/wiki/Perfect_number
if (factor_sum === index)
numbers.push(index);
});
library_namespace.debug('numbers: [' + numbers.length + '] '
+ numbers.slice(0, 30) + '...', 1, 'perfect_numbers');
return numbers;
}
_.perfect_numbers = perfect_numbers;
// ---------------------------------------------------------------------//
// 回文數 palindromic number or numeral palindrome
// http://www.csie.ntnu.edu.tw/~u91029/Palindrome.html
function palindrome_list(limit, base) {
if (!base)
base = DEFAULT_BASE;
// 個位數皆為回文數。
var list = new Array(Math.min(base, limit)).fill(ABSORBING_ELEMENT)
.map(function(v, i) {
return i;
});
if (limit <= base)
return list;
for (var power = 1, next_power = base, n;;) {
// 2e 位數(e), e.g., 1001
// left: 10^e–10^(e+1)-1
for (var l = power; l < next_power; l++) {
var left = l.toString(base),
// right side
right = left.split('').reverse().join('');
n = parseInt(left + right, base);
if (n >= limit)
break;
list.push(n);
}
if (n >= limit)
break;
// 2e+1 位數, e.g., 10201
// left: 10^e–10^(e+1)-1
for (var l = power; l < next_power; l++) {
var left = l.toString(base),
// right side
right = left.split('').reverse().join('');
for (var middle = 0; middle < base; middle++) {
n = parseInt(left + middle + right, base);
if (n >= limit)
break;
list.push(n);
}
}
if (n >= limit)
break;
power = next_power;
next_power *= base;
}
return list;
}
_.palindrome_list = palindrome_list;
function is_palindrome(natural, base) {
if (typeof natural !== 'srting')
// assert: typeof natural !== 'number'
natural = natural.toString(base);
// 將這個數的數字按相反的順序重新排列後，所得到的數和原來的數一樣。
return natural === natural.split('').reverse().join('');
}
_.is_palindrome = is_palindrome;
// ---------------------------------------------------------------------//
_// JSDT:_module_
.
/**
* 猜測一個數可能的次方數。
*
* @param {Number}
* number 數字
* @param {Boolean}
* type false: base 為整數, true: base 為有理數
*
* @returns [{Integer} base 分子, {Integer} base 分母, {Integer} exponent 分子,
* {Integer} exponent 分母]
*
* @since 2005/2/18 19:20 未完成
*/
guess_exponent = function(number, type) {
var bn, bd, en = 1, ed, sq = [ 1, number ], t, q,
// default error, accuracy, stopping tolerance, 容許誤差
error = Number.EPSILON,
//
g = function(n) {
q = _.to_rational_number(n, 99999);
if ((!type || q[1] === 1) && !(q[0] > 99999 && q[1] > 99999)
&& q[2] / n < error)
bn = q[0], bd = q[1], ed = t;
};
if (!ed)
g(sq[t = 1]);
if (!ed)
g(sq[t = 2] = sq[1] * sq[1]);
if (!ed)
g(sq[t = 3] = sq[1] * sq[2]);
if (!ed)
g(sq[t = 4] = sq[2] * sq[2]);
if (!ed)
g(sq[t = 5] = sq[2] * sq[3]);
if (!ed)
bn = number, bd = ed = 1;
return [ bn, bd, en, ed ];
};
/**
* get random prime(s)
*
* @param {Integer}count
* 個數
* @param {Array}exclude
* 排除
* @param {Boolean}all_different
*
* @returns random prime / random prime array
*
* @since 2009/10/21 11:57:47
*/
function get_random_prime(count, exclude, all_different) {
var i, j, p = [], l;
if (!count || count < 1)
count = 1;
if (!get_random_prime.excluded)
get_random_prime.excluded = [];
if (exclude)
exclude = [];
// 先行準備好足夠的 primes。
prime(2 * count, 2 * count);
for (j = 0; j < count; j++) {
// timeout
l = 80;
do {
i = Math.round(10 * Math.tan(Math.random() * 1.5));
if (!--l)
// timeout
return;
} while (get_random_prime.excluded[i]);
p.push(primes[i]);
if (exclude)
exclude.push(i);
}
// 選完才排除本次選的
if (exclude)
for (j = 0, l = exclude.length; j < l; j++) {
i = exclude[j];
if (get_random_prime.excluded[i])
get_random_prime.excluded[i]++;
else
get_random_prime.excluded[i] = 1;
}
return count === 1 ? p[0] : p;
}
// return [GCD, n1, n2, ..]
get_random_prime.get_different_number_set = function(count, till, GCD_till) {
delete this.excluded;
if (!GCD_till)
GCD_till = 1e5;
if (!till)
till = 1e5;
/**
* 求乘積, 乘到比till小就回傳.
*
* @param nums
* num array
* @param till
* @returns {Number}
*/
function get_product(nums, till) {
var p = 1, i = 0, l = nums.length;
for (; i < l; i++) {
if (till && p * nums[i] > till)
break;
p *= nums[i];
}
return p;
}
var GCD = get_product(this(20, 1), GCD_till), na = [], n_e = [], n, i = 0, out;
n_e[GCD] = 1;
for (; i < count; i++) {
out = 80; // timeout
do {
n = this(20);
n.unshift(GCD);
n = get_product(n, till);
} while (n_e[n] && --out);
n_e[n] = 1;
na.push(n);
}
if (typeof lcm == 'function')
na.LCM = lcm(na);
na.GCD = GCD;
return na;
};
// ------------------------------------------------------------------------------------------------------//
/**
* 求取反函數 caculator[-1](result)
*
* @deprecated
*/
function get_boundary(caculator, result, down, up, limit) {
if (up - down === 0)
return up;
var boundary, value, increase;
// assert: caculator(down) – caculator(up) 為嚴格遞增/嚴格遞減函數。
if (caculator(up) - caculator(down) < 0) {
// swap.
boundary = up;
up = down;
down = boundary;
}
// assert: caculator(down) < caculator(up)
increase = down < up;
if (!(limit > 0))
limit = 800;
do {
boundary = (up + down) / 2;
// console.log(down + ' – ' + boundary + ' – ' + up);
if (boundary === down || boundary === up)
return boundary;
value = result - caculator(boundary);
if (value === 0) {
if (result - caculator(down) === 0) {
down = boundary;
value = true;
}
if (result - caculator(up) === 0) {
up = boundary;
value = true;
}
if (value && (increase ? up - down > 0 : up - down < 0))
continue;
return boundary;
}
if (value > 0)
down = boundary;
else
up = boundary;
} while (--limit > 0 && (increase ? up - down > 0 : up - down < 0));
throw 'get_boundary: caculator is not either strictly increasing or decreasing?';
}
/**
* 求根/求取反函數 equation^-1(y)。 using Secant method.
*
* @param {Function}equation
* 演算式, mapping function
* @param {Number}x0
* 內插法(線性插值 Interpolation)求值之自變數 variable 下限，設定初始近似值。
* @param {Number}x1
* 內插法(線性插值 Interpolation)求值之自變數 variable 上限，設定初始近似值。
* @param {Number}[y]
* 目標值。default: 0. get (equation^-1)(y)
* @param {Object}[options]
* 附加參數/設定特殊功能與選項
*
* @returns {Number}root: equation(root)≈y
*
* @see Interpolation 以內插值替換
* https://en.wikipedia.org/wiki/Root-finding_algorithm
* https://en.wikipedia.org/wiki/Secant_method
*/
function secant_method(equation, x0, x1, y, options) {
// default error, accuracy, stopping tolerance, 容許誤差
var error = Number.EPSILON;
if (!options)
options = Object.create(null);
else if (options > 0)
error = options;
else if (options.error > 0)
error = Math.abs(options.error);
y = +y || 0;
var count = (options.count || 40) | 0,
// assert: y0 = equation(x0)
y0 = 'y0' in options ? options.y0 : equation(x0),
// assert: y1 = equation(x1)
y1 = 'y1' in options ? options.y1 : equation(x1),
//
x2 = x1, y2 = y1;
if (typeof options.start_OK === 'function'
// 初始測試: Invalid initial value, 不合理的初始值，因此毋須繼續。
&& !options.start_OK(y0, y1))
return;
// main loop
while (error < Math.abs(y2 - y) && count-- > 0
// 分母不應為 0 或 NaN。
&& (y0 -= y1)
// 測試已達極限，已經得到相當好的效果。無法取得更精準值。
// assert: else: x0===x, 可能是因為誤差已過小。
&& ((x2 = x1 - (x1 - x0) * (y - y1) / y0) !== x1 || x1 !== x0)) {
// evaluate result
y2 = equation(x2);
if (false)
library_namespace.debug(count + ': ' + x2 + ',' + y2 + ' → '
+ (y2 - y));
// shift items
x0 = x1, y0 = y1;
x1 = x2, y1 = y2;
}
return x2;
}
_.secant_method = secant_method;
/**
* 求根/求取反函數 equation^-1(y)。 using Sidi's generalized secant method.
*
* @param {Function}equation
* 演算式, mapping function
* @param {Number}x0
* 求值之自變數 variable 下限，設定初始近似值。
* @param {Number}x1
* 求值之自變數 variable 上限，設定初始近似值。
* @param {Number}[y]
* 目標值。default: 0. get (equation^-1)(y)
* @param {Object}[options]
* 附加參數/設定特殊功能與選項
*
* @returns {Number}root: equation(root)≈y
*
* @see https://en.wikipedia.org/wiki/Root-finding_algorithm
* @see https://en.wikipedia.org/wiki/Sidi's_generalized_secant_method
*/
function Sidi_method(equation, x0, x1, y, options) {
// default error, accuracy, stopping tolerance, 容許誤差
var error = Number.EPSILON;
if (!options)
options = Object.create(null);
else if (typeof options === 'number' && options > 0)
// ↑ @Firefox/44.0:
// isNaN(Object.create(null)): TypeError: can't convert v to number
// Number.isNaN(Object.create(null)) === false
error = options;
else if (options.error > 0)
error = Math.abs(options.error);
y = +y || 0;
var count = (options.count || 40) | 0,
// assert: y0 = equation(x0)
y0 = 'y0' in options ? options.y0 : equation(x0),
// assert: y1 = equation(x1)
y1 = 'y1' in options ? options.y1 : equation(x1);
if (typeof options.start_OK === 'function'
// 初始測試: Invalid initial value, 不合理的初始值，因此毋須繼續。
&& !options.start_OK(y0, y1))
return;
// initialization
var x2 = x1 - (x1 - x0) * (y1 - y) / (y1 - y0),
//
y2 = equation(x2), x3 = x2, y3 = y2,
// divided differences, 1階差商
y10 = (y1 - y0) / (x1 - x0), y21 = (y2 - y1) / (x2 - x1),
// 2階差商
y210 = (y21 - y10) / (x2 - x0),
// 暫時使用。
denominator;
// main loop of Sidi's generalized secant method (take k = 2)
while (error < Math.abs(y3 - y)
&& count-- > 0
// 檢查是否兩個差距極小的不同輸入，獲得相同輸出。
&& y21 !== 0
// 分母不應為 0 或 NaN。
&& (denominator = y21 + y210 * (x2 - x1))
// Avram Sidi (2008), "Generalization Of The Secant Method For
// Nonlinear Equations"
// 可能需要考量會不會有循環的問題。
&& ((x3 = x2 - (y2 - y) / denominator) !== x2 || x2 !== x1 || x1 !== x0)) {
// evaluate result
y3 = equation(x3);
if (false)
console.log(count + ': ' + x3 + ',' + y3 + ' → error '
+ (y3 - y));
// shift items
x0 = x1, y0 = y1;
x1 = x2, y1 = y2;
x2 = x3, y2 = y3;
// reckon divided differences
y10 = y21;
y21 = (y2 - y1) / (x2 - x1);
// y210 = (y21 - y10) / (x2 - x0);
// incase y21 === y10
if (y210 = y21 - y10)
y210 /= x2 - x0;
if (false)
console.log('divided differences: ' + [ y10, y21, y210 ]);
}
return x3;
}
_.Sidi_method = Sidi_method;
/**
* 以 Brent's method 求根/求取反函數 equation^-1(y)。
*
* TODO: 牛頓法, options.derivative
*
* @param {Function}equation
* 演算式, mapping function
* @param {Number}x0
* 求值之自變數 variable 下限，設定初始近似值。
* @param {Number}x1
* 求值之自變數 variable 上限，設定初始近似值。
* @param {Number}[y]
* 目標值。default: 0. get (equation^-1)(y)
* @param {Object}[options]
* 附加參數/設定特殊功能與選項
*
* @returns {Number}root: equation(root)≈y
*
* @see https://en.wikipedia.org/wiki/Root-finding_algorithm
* @see https://en.wikipedia.org/wiki/Brent%27s_method
* @see http://www.boost.org/doc/libs/1_58_0/boost/math/tools/minima.hpp
* @see http://people.sc.fsu.edu/~jburkardt/cpp_src/brent/brent.cpp
* @see http://www.cscjournals.org/manuscript/Journals/IJEA/volume4/Issue1/IJEA-33.pdf
* @see http://www.cscjournals.org/manuscript/Journals/IJEA/volume2/Issue1/IJEA-7.pdf
*/
function Brent_method(equation, x0, x1, y, options) {
// default error, accuracy, stopping tolerance, 容許誤差
// @see Number.EPSILON
var error = 0;
if (!options)
options = Object.create(null);
else if (typeof options === 'number' && options > 0)
error = options;
else if (options.error > 0)
error = Math.abs(options.error);
y = +y || 0;
var count = (options.count || 40) | 0,
// assert: y0 = equation(x0)
y0 = 'y0' in options ? options.y0 : equation(x0),
// assert: y1 = equation(x1)
y1 = 'y1' in options ? options.y1 : equation(x1);
if (typeof options.start_OK === 'function'
// 初始測試: Invalid initial value, 不合理的初始值，因此毋須繼續。
&& !options.start_OK(y0, y1))
return;
if (y0 === y)
return x0;
if (y1 === y)
return x1;
if (!((y0 - y) * (y1 - y) < 0))
// the root is not bracketed.
// 但亦有可能只是所取範圍過大，若中間多挑幾點，或許能找到合適的值。
// 回退至使用 Sidi's generalized secant method。
return Sidi_method(equation, x0, x1, y, options);
// copy from double zero () @
// http://people.sc.fsu.edu/~jburkardt/cpp_src/brent/brent.cpp
// rewrite when I have time...
var a = x0, b = x1, c, d, e, sa, sb, sc, fa, fb, fc, m, tol, p, q, r, s;
//
// Make local copies of A and B.
//
c = sa = a;
fc = fa = y0 - y;
sb = b;
fb = y1 - y;
d = e = b - a;
while (true) {
if (Math.abs(fc) < Math.abs(fb)) {
sa = sb, sb = c, c = sa;
fa = fb, fb = fc, fc = fa;
}
// tol = 2 * Number.EPSILON * Math.abs(sb) + error;
tol = 0;
if (sa === sb || Math.abs(m = (c - sb) / 2) <= tol || fb === 0) {
break;
}
if (Math.abs(e) < tol || Math.abs(fa) <= Math.abs(fb)) {
d = e = m;
} else {
s = fb / fa;
if (sa === c) {
p = 2 * m * s;
q = 1 - s;
} else {
q = fa / fc;
r = fb / fc;
p = s * (2 * m * q * (q - r) - (sb - sa) * (r - 1));
q = (q - 1) * (r - 1) * (s - 1);
}
if (0 < p) {
q = -q;
} else {
p = -p;
}
s = e;
e = d;
if (2 * p < 3 * m * q - Math.abs(tol * q)
&& p < Math.abs(0.5 * s * q)) {
d = p / q;
} else {
d = e = m;
}
}
sa = sb;
fa = fb;
if (tol < Math.abs(d)) {
sb = sb + d;
} else if (0 < m) {
sb = sb + tol;
} else {
sb = sb - tol;
}
fb = equation(sb) - y;
if ((0 < fb && 0 < fc) || (fb <= 0 && fc <= 0)) {
c = sa;
fc = fa;
d = e = sb - sa;
}
}
return sb;
}
_.Brent_method = Brent_method;
_.find_root = Brent_method;
/**
* 不用微分求取局部極小值 get local minimum using Brent's golden section search
*
* 注意: 不同的輸入，即使對同一極值，亦可能會得出不同的輸出值。
*
* @param {Function}equation
* 演算式, mapping function
* @param {Number}min
* 求值之自變數 variable 下限，設定初始近似值。
* @param {Number}max
* maximum. 求值之自變數 variable 上限，設定初始近似值。
* @param {Object}[options]
* 附加參數/設定特殊功能與選項
*
* @returns {Number}[x,fx]: equation(x)≈fx≈minimum
*
* @see https://en.wikipedia.org/wiki/Maxima_and_minima
* https://en.wikipedia.org/wiki/Golden_section_search
* https://zh.wikipedia.org/wiki/%E6%9E%81%E5%80%BC#.E6.B1.82.E6.9E.81.E5.80.BC.E7.9A.84.E6.96.B9.E6.B3.95
* https://zh.wikipedia.org/wiki/%E5%8F%98%E5%88%86%E6%B3%95
* http://maths-people.anu.edu.au/~brent/pub/pub011.html
* http://math.stackexchange.com/questions/58787/looking-for-numerical-methods-for-finding-local-maxima-and-minima-of-a-function
*/
function Brent_minima(equation, min, max, options) {
// default error, accuracy, stopping tolerance, 容許誤差
// @see Number.EPSILON
var error = 0;
if (!options)
options = Object.create(null);
else if (options > 0)
error = options;
else if (options.error > 0)
error = Math.abs(options.error);
var count = (options.count || 40) | 0, data, data_u;
// copy from std::pair<T, T> brent_find_minima @
// http://www.boost.org/doc/libs/1_58_0/boost/math/tools/minima.hpp
// rewrite when I have time...
var tolerance = error,
// minima so far
x,
// second best point
w,
// previous value of w
v,
// most recent evaluation point
u,
// The distance moved in the last step
delta,
// The distance moved in the step before last
delta2,
// function evaluations at u, v, w, x
fu, fv, fw, fx,
// midpoint of ((min)) and ((max))
mid,
// minimal relative movement in x
fract1, fract2,
// golden ratio, don't need too much precision here!
golden = 0.3819660;
x = w = v = max;
fx = equation(x);
if (Array.isArray(fx))
data = fx[1], fx = fx[0];
fw = fv = fx;
delta2 = delta = 0;
do {
// get midpoint
mid = (min + max) / 2;
// work out if we're done already:
fract1 = tolerance * Math.abs(x) + tolerance / 4;
fract2 = 2 * fract1;
if (Math.abs(x - mid) <= (fract2 - (max - min) / 2))
break;
if (Math.abs(delta2) > fract1) {
// try and construct a parabolic fit:
var r = (x - w) * (fx - fv), q = (x - v) * (fx - fw), p = (x - v)
* q - (x - w) * r;
q = 2 * (q - r);
if (q > 0)
p = -p;
q = Math.abs(q);
var td = delta2;
delta2 = delta;
// determine whether a parabolic step is acceptible or not:
if ((Math.abs(p) >= Math.abs(q * td / 2))
|| (p <= q * (min - x)) || (p >= q * (max - x))) {
// nope, try golden section instead
delta2 = (x >= mid) ? min - x : max - x;
delta = golden * delta2;
} else {
// whew, parabolic fit:
delta = p / q;
u = x + delta;
if (((u - min) < fract2) || ((max - u) < fract2))
delta = (mid - x) < 0 ? -Math.abs(fract1) : Math
.abs(fract1);
}
} else {
// golden section:
delta2 = (x >= mid) ? min - x : max - x;
delta = golden * delta2;
}
// update current position:
u = (Math.abs(delta) >= fract1) ? (x + delta)
: (delta > 0 ? (x + Math.abs(fract1)) : (x - Math
.abs(fract1)));
fu = equation(u);
if (Array.isArray(fu))
data_u = fu[1], fu = fu[0];
if (fu <= fx) {
// good new point is an improvement!
// update brackets:
if (u >= x)
min = x;
else
max = x;
if (x === u)
// 無改變。
break;
// update control points:
v = w, w = x, x = u;
fv = fw, fw = fx, fx = fu;
data = data_u;
} else {
// Oh dear, point u is worse than what we have already,
// even so it *must* be better than one of our endpoints:
if (u < x)
min = u;
else
max = u;
if ((fu <= fw) || (w == x)) {
// however it is at least second best:
v = w, w = u;
fv = fw, fw = fu;
} else if ((fu <= fv) || (v == x) || (v == w)) {
// third best:
v = u;
fv = fu;
}
}
} while (--count);
(x = new Number(x)).y = fx;
if (data !== undefined)
if (library_namespace.is_Object(data))
Object.assign(x, data);
else
x.data = data;
return x;
}
_.Brent_minima = Brent_minima;
// 不用微分求取求取局部極小值 get local minimum
_.find_minima = Brent_minima;
_.find_maxima = function(equation, min, max, options) {
return Brent_minima(function(x) {
return -equation(x);
}, min, max, options);
};
// 不用微分求取局部極值 get local maximum and minimum
// _.find_extremum = Brent_minima;
// ------------------------------------------------------------------------------------------------------//
// 組合數學反向思考: 有重複的=全-沒有重複的
// 解法之所以錯誤往往是因為重複計數。
// 裝載問題
// http://codex.wiki/post/117994-555
// 正整數拆分/數字拆解演算法：將數字拆分成最大元素不大於 max 的組合
// http://www.nowamagic.net/algorithm/algorithm_IntegerDivisionDynamicProgramming.php
// https://openhome.cc/Gossip/AlgorithmGossip/SeparateNumber.htm
// https://en.wikipedia.org/wiki/Knapsack_problem
/**
* Get the count of integer partitions. 整數分拆: 將正整數 sum 拆分，表達成一些正整數的和。
*
* TODO: part
*
* @param {Natural}sum
* integer to be apart.
* @param {Natural|Array}[part_count]
* TODO: 給出恰好劃分成 part_count 個整數的劃分。
* @param {Array}[summands]
* 給出只包括 summands 的劃分。
* @param {Array}cache
* cache[sum][summands] = count
*
* @returns {Natural}組合方法數
*
* @see 貪心算法,貪心法 https://en.wikipedia.org/wiki/Greedy_algorithm
*
* @inner
*/
function count_partitions(sum, part_count, summands, cache) {
var key = summands.join(''), _c = cache[sum];
// https://en.wikipedia.org/wiki/Memoization
if (!_c)
_c = cache[sum] = [];
else if (key in _c)
return _c[key];
library_namespace.debug([ sum, summands ], 3);
// 不更動 summands
summands = summands.slice();
var summand = summands.pop(), count = sum / summand | 0;
sum %= summand;
if (summands.length === 0) {
// 檢查當前解是否是可行解。
// 若有餘數，表示此法不通。
// e.g., 以2元分3元
return sum === 0 ? 1 : 0;
}
var counter = 0;
for (; count >= 0; count--, sum += summand) {
library_namespace.debug(summand + '⋅' + count + '+' + sum + '; '
+ counter + '; ' + summands, 3);
if (sum === 0) {
// e.g., 以 5元分100元，當count===20時，此時也算一次。
counter++;
} else {
counter += count_partitions(sum, part_count, summands, cache);
}
}
return _c[key] = counter;
}
/**
* Get the count of integer partitions. 整數分拆: 將正整數 sum 拆分，表達成一些正整數的和。
* 兌換/分桶/分配問題
*
* TODO: part, count of summands, options
*
* TODO: 部分應該有 O(1) 的方法。 see
* http://www.mobile01.com/topicdetail.php?f=37&t=2195318&p=3
*
* @param {Natural}sum
* integer to be apart.
* @param {Natural|Array}[part_count]
* TODO: 給出恰好劃分成 part_count 個整數的劃分。
* @param {Array}[summands]
* 給出只包括 summands 的劃分。
*
* @returns {Natural}組合方法數
*
* @see https://en.wikipedia.org/wiki/Partition_%28number_theory%29
* 組合數學/總價值固定之錢幣排列組合方法數<br />
* http://www.cnblogs.com/python27/p/3303721.html
* http://mathworld.wolfram.com/Partition.html
* http://mathworld.wolfram.com/PartitionFunctionP.html
* http://www.zhihu.com/question/21075235
* http://blog.csdn.net/iheng_scau/article/details/8170669
*/
function integer_partitions(sum, part_count, summands) {
// assert: sum≥0
if (!sum)
return 0;
if (!summands || !summands.length) {
if (part_count)
throw 'integer_partitions: NYI';
if (sum < 0)
// p(n) = 0 for n negative.
return 0;
// the partition function p(n)
// https://en.wikipedia.org/wiki/Partition_%28number_theory%29#Partitions_in_a_rectangle_and_Gaussian_binomial_coefficients
var count = _.number_array(sum + 1);
// p(0) = 1
count[0] = 1;
for (var i = 1; i <= sum; i++)
for (var j = i; j <= sum; j++)
count[j] += count[j - i];
library_namespace.debug(count, 3);
return count[sum];
}
// 檢查是否有解。
if (sum % GCD.apply(null, summands) !== 0)
return;
summands = summands.slice().sort(library_namespace.ascending);
return count_partitions(sum, part_count, summands, []);
}
_.integer_partitions = integer_partitions;
// ------------------------------------------------------------------------------------------------------//
_// JSDT:_module_
.
/**
* VBScript has a Hex() function but JScript does not.
*
* @param {Number}number
* number
*
* @return {String} number in hex
*
* @example alert('0x'+CeL.hex(16725))
*/
hex = function(number) {
return ((number = Number(number)) < 0 ? number + 0x100000000 : number)
.toString(16);
};
_// JSDT:_module_
.
/**
* 補數計算。 正數的補數即為自身。若要求得互補之後的數字，請設成負數。
*
* @param {Number}number
* number
*
* @return {Number} base 1: 1's Complement, 2: 2's Complement, (TODO: 3, 4,
* ..)
* @example alert(complement())
* @see http://www.tomzap.com/notes/DigitalSystemsEngEE316/1sAnd2sComplement.pdf
* https://en.wikipedia.org/wiki/Method_of_complements
* https://en.wikipedia.org/wiki/Signed_number_representations
* @since 2010/3/12 23:47:52
*/
complement = function() {
return this.from.apply(this, arguments);
};
_// JSDT:_module_
.complement.prototype = {
base : 2,
// 1,2,..
bits : 8,
// radix complement or diminished radix complement.
// https://en.wikipedia.org/wiki/Method_of_complements
diminished : 0,
/**
* 正負符號. 正: 0/false, 負 negative value:!=0 / true
*/
sign : 0,
// get the value
valueOf : function() {
return this.sign ? -this.value : this.value;
},
/**
* set value
*/
set : function(value) {
var m = Number(value), a = Math.abs(m);
if (isNaN(m) || m && a < 1e-8 || a > 1e12) {
library_namespace.debug('complement.set: error number: '
+ value);
return;
}
this.sign = m < 0;
// this.value 僅有正值
this.value = a;
return this;
},
/**
* input
*/
from : function(number, base, diminished) {
// 正規化
number = ('' + (number || 0)).replace(/\s+$|^[\s0]+/g, '') || '0';
if (false)
library_namespace.debug(number + ':' + number.length + ','
+ this.bits);
// 整數部分位數
var value = number.indexOf('.'), tmp;
// -1: NOT_FOUND
if (value == -1)
value = number.length;
// TODO: not good, need to optimize
if (value > this.bits)
// throw 'overflow';
library_namespace.error('complement.from: overflow: ' + value);
if (typeof diminished === 'undefined')
// illegal setup
diminished = this.diminished;
else
this.diminished = diminished;
if ((base = Math.floor(base)) && base > 0) {
if (base === 1)
base = 2, this.diminished = 1;
this.base = base;
} else
// illegal base
base = this.base;
if (false)
library_namespace.debug(base + "'s Complement");
// TODO: 僅對 integer 有效
value = parseInt(number, base);
tmp = Math.pow(base, this.bits - 1);
if (value >= tmp * base)
// throw 'overflow';
library_namespace.error('complement.from: overflow: ' + value);
// library_namespace.debug('compare ' + value + ',' + tmp);
if (value < tmp)
this.sign = 0;
else {
library_namespace.debug('負數 '
+ (tmp * base - (diminished ? 1 : 0)) + '-' + value
+ '=' + (tmp * base - (diminished ? 1 : 0) - value), 3);
this.sign = 1;
value = tmp * base - (diminished ? 1 : 0) - value;
}
this.value = value;
// library_namespace.debug(number + ' → '+this.valueOf());
return this;
},
/**
* output
*/
to : function(base, diminished) {
if (!(base = Math.floor(base)) || base < 1)
base = this.base;
else if (base === 1)
base = 2, diminished = 1;
if (typeof diminished === 'undefined')
diminished = this.diminished;
var value = this.value, tmp = Math.pow(base, this.bits - 1);
if (value > tmp || value === tmp && (diminished || !this.sign))
// throw 'overflow';
library_namespace.error('complement.to: overflow: '
+ (this.sign ? '-' : '+') + value);
if (this.sign) {
tmp *= base;
if (diminished)
// TODO: 僅對 integer 有效
tmp--;
// library_namespace.debug('負數 ' + value + '，summation=' + tmp);
// 負數，添上兩補數之和
value = tmp - value;
}
if (false)
library_namespace.debug('value: ' + (this.sign ? '-' : '+')
+ value);
value = value.toString(Math.max(2, this.base));
return value;
}
};
_// JSDT:_module_
.complement.prototype.toString = _.complement.prototype.to;
/**
* haversines
*
* @param {Number}θ
* angle (in radians)
*
* @returns {Number}haversines(θ)
*
* @see https://en.wikipedia.org/wiki/Haversine_formula
*/
function hav(θ) {
return (1 - Math.cos(θ)) / 2;
// hav(θ) = sin^2(θ/2)
}
// ---------------------------------------------------------------------//
/**
* 取差集。<br />
* 警告: 若有成員為 {Object}，或許先自行處理會比較有效率。
*
* set_1 → set_1 - set_2 <br />
* set_2 → set_2 - set_1
*
* @param {Array|Object}set_1
* set / hash 1
* @param {Array|Object}set_2
* set / hash 2
* @param {Boolean}[clone]
* clone === true 表示先 clone，不直接在原物件上寫入結果。但若為 {Array}，則*必定*為
* clone!
*
* @returns {Array}[ complement of set_1, complement of set_2 ]
*/
function get_set_complement(set_1, set_2, clone) {
if (!set_1 || !set_2 || typeof set_1 !== 'object'
|| typeof set_2 !== 'object') {
throw new Error('get_set_complement: Invalid type');
}
var hash_1, hash_2;
// 確保 keys 準備好，並把為 hash 的全部轉到 hash_*。
if (!Array.isArray(set_1)) {
hash_1 = clone ? Object.clone(set_1) : set_1;
set_1 = Object.keys(hash_1);
}
if (!Array.isArray(set_2)) {
hash_2 = clone ? Object.clone(set_2) : set_2;
set_2 = Object.keys(hash_2);
}
// assert: set_1, set_2 are {Array}
// 從比較小的來處理比較快。
// 若 set_1 比較短，則看 hash_2 是否存在；若沒 hash_2 則將key指到比較短的 _2。
// 依之後的演算法，hash_2 必須存在，又需建造，因此務必為比較短的，因為消耗高。
// 因此若是另一方已經有 hash，則直接用之。
var key_is_2 = !(set_1.length > set_2.length ? hash_1 : hash_2);
library_namespace.debug('key_is_2: ' + key_is_2, 3,
'get_set_complement');
if (key_is_2) {
// swap 1, 2
var tmp = set_1;
set_1 = set_2;
set_2 = tmp;
tmp = hash_1;
hash_1 = hash_2;
hash_2 = tmp;
}
var no_hash_2 = !hash_2;
if (no_hash_2) {
library_namespace.debug('建造 hash: 依之後的演算法，hash_2 必須存在。', 3,
'get_set_complement');
hash_2 = set_2.to_hash();
}
var resort_1 = [];
set_1.forEach(function(item) {
// assert: item is in _1
// 處理 set_1 之成員為 {Object} 的情況。
// 警告:這邊需要與 Array.prototype.product() 採用相同的 to string 方法。
var key = typeof item === 'object' ? JSON.stringify(item) : item;
if (key in hash_2) {
library_namespace.debug('key 在 _1 + 在 _2: 在兩方都刪掉: ' + key, 5,
'get_set_complement');
delete hash_2[key];
if (hash_1) {
delete hash_1[key];
}
} else {
library_namespace.debug('key 在 _1 不在 _2: 留下 _1 的 key: ' + key,
5, 'get_set_complement');
if (!hash_1) {
resort_1.push(item);
}
// _2 本來就沒有，不動。
}
});
set_1 = resort_1;
if (no_hash_2) {
// 得造出 set_2。
// 維持 set_2 的順序，並避免去掉重複。
var keep_order = true;
if (keep_order) {
var resort_2 = [];
set_2.forEach(function(item) {
// 處理 set_2 之成員為 {Object} 的情況。
// 警告:這邊需要與前面採用相同的 to string 方法。
var key = typeof item === 'object' ? JSON.stringify(item)
: item;
if (key in hash_2) {
resort_2.push(item);
}
});
set_2 = resort_2;
} else {
// 這會比較快，但實際應用上會造成不確定性：不能確定截掉的是哪一個。
set_2 = Object.keys(hash_2);
}
hash_2 = undefined;
}
if (key_is_2) {
// swap 1, 2
var tmp = set_1;
set_1 = set_2;
set_2 = tmp;
tmp = hash_1;
hash_1 = hash_2;
hash_2 = tmp;
}
return [ hash_1 || set_1, hash_2 || set_2 ];
}
_.get_set_complement = get_set_complement;
// ---------------------------------------------------------------------//
// CeL.math.number_array()
if (library_namespace.typed_arrays) {
// https://developer.mozilla.org/en-US/docs/Web/JavaScript/Typed_arrays
(_.number_array = function number_array_ArrayBuffer(size, fill, type) {
if (false) {
var buffer = new ArrayBuffer(type.BYTES_PER_ELEMENT * size);
}
// DataView
var array = new (type || _.number_array.default_type)(
// buffer
size);
if (fill)
array.fill(fill);
return array;
})
// should TypedArray. e.g., Int32Array, Uint32Array
.default_type = Int32Array;
} else {
_.number_array = function number_array_Array(size, fill) {
// 經過 .fill() 以定義每個元素，這樣在 .forEach() 時才會遍歷到。
return new Array(size).fill(fill || 0);
};
}
/**
* Create digit value table. 建構出位數值表。
*
* 對一般問題，所要求的，即是以 greedy algorithm （貪心算法,貪心法）遞歸搜索，從
* digit_table[0–末位數]各選出一位數值，使其總合為0。<br />
* 因為各位數值有其特性，因此可能存有些技巧以降低所需處理之數據量。
*
* @param {Array}initial_value
* 當位數每一位數值應當減去的初始值。
* @param {Object}[options]
* 附加參數/設定特殊功能與選項
*
* @returns {Array}digit value table
*/
function digit_table(initial_value, options) {
var base = options && options.base || DEFAULT_BASE;
if (typeof initial_value === 'string') {
if (initial_value === 'factorial') {
initial_value = factorial.map(base - 1);
} else {
var matched = initial_value.match(/power[\s:=]*(\d+)/i);
if (matched) {
matched = +matched[1];
// initial_value[digit] = digit^power
initial_value = _.number_array(base).map(
function(id, index) {
return Math.pow(index, matched);
});
}
}
}
// assert: initial_value 之元素皆 {Natural} > 0
library_namespace.debug(initial_value, 3);
/**
* value of each digit.
*
* table[exponent=0–(maximum exponent)][digit=0–9] =<br />
* {Number} digit*base^exponent - digit^exponent
*
* @type {Array} [][]
*/
var table = [],
/**
* accumulated min. 自個位數起累積的最小值。
*
* sum_min[exponent=0–(maximum exponent)] =<br />
* ∑自0至(exponent-1)位累積的(digit value之最小值)
*
* @type {Array}
*/
sum_min = [];
table.min = sum_min;
// 準備好 digit_value table, min/maximum value。
for (var exponent = 0,
/** {Natural}power = base^exponent */
power = 1;
// 不需要此項限制，照理來說應該在其之前即已跳出。
// exponent < base
; exponent++) {
/** {Array}value_array[digit=0–9] = {Number}位數值(digit value) */
var value_array = _.number_array(base),
/** {Number}本位數各數字位數值的最小值，不包含0 */
min = Infinity;
// 計算 (base^exponent) 之位數值(digit value)，並記錄最小值。
for (var digit = 0; digit < base; digit++) {
var digit_value = value_array[digit] = digit * power
- initial_value[digit];
// 因為0不能當數字頭，[0] 不設定極值。
if (digit > 0 && digit_value < min)
min = digit_value;
}
/** {Number}本位數各數字位數值的最小值，包含0 */
var min_0 = Math.min(min, value_array[0]);
// console.log(value_array);
if (exponent > 0) {
// 記錄自個位數起之累積之最小值。
var last_min = sum_min[exponent - 1];
min += last_min;
min_0 += last_min;
}
if (min > 0)
// 對 n 位數，數值範圍為 base^(n-1)–base^n-1。但位數值和若已經過大，
// 代表此位數以上，如第 (n+1) 位數，就算每個位數值和都取最小值，總和也不可能為0。
break;
sum_min.push(min_0);
table.push(value_array);
power *= base;
}
if (!options || !options.find)
return table;
// ------------------------------------------
// process: 以 greedy algorithm （貪心算法,貪心法）遞歸搜索
/**
* calculate sum of digit values
*/
function calculate_sum(sum, exponent, digits) {
// 測試是否應終結。
if (exponent < 0) {
// 0,1 為當然結果。
if (sum === 0 && (digits = +digits) > min)
result.push(digits);
return;
}
if (sum > 0 && sum + sum_min[exponent] > 0)
// 接下來的總和也不可能為0。
return;
// 遞歸搜索
table[exponent--].forEach(function(v, d) {
calculate_sum(digits || d ? sum + v : sum, exponent,
// 對於 0 開頭者，視做少一位數，同時不計算本位數之 [0]。
digits || d ? digits + d : digits);
});
}
var result = [],
/** {Natural}所得值需要大於此值。 */
min;
if (options && ('min' in options))
min = options.min;
else if (!initial_value.some(function(v, d) {
if (d && v !== d) {
min = d;
return true;
}
}))
min = 0;
calculate_sum(0, table.length - 1, '');
library_namespace.debug(result, 2);
return result;
}
_.digit_table = digit_table;
// ---------------------------------------------------------------------//
// assert: {ℕ⁰:Natural+0}((this))
// fast than String(natural).chars(), natural.toString().chars() or
// natural.toString(base).chars()
// TODO: 處理小數/負數/大數
function Number_digits(base) {
if (!((base |= 0) >= 2))
base = DEFAULT_BASE;
var natural = Math.floor(Math.abs(this)), digits = [];
do {
digits.unshift(natural % base);
} while ((natural = Math.floor(natural / base)) > 0);
return digits;
}
// 數字和, 位數和
// natural.digits(base).sum()
// https://en.wikipedia.org/wiki/Digit_sum
// to get digit root: using natural.digit_sum(base) % base
function Number_digit_sum(base) {
if (!((base |= 0) >= 2))
base = DEFAULT_BASE;
var natural = Math.floor(Math.abs(this)), sum = ABSORBING_ELEMENT;
do {
sum += natural % base;
} while ((natural = Math.floor(natural / base)) > 0);
return sum;
}
// count digits of integer
// = floor(log_10(base)) + 1
// TODO: 測試二分搜尋 [base,base^2,base^3,...] 方法之效率。 e.g., [10,100,100,...]
function Number_digit_length(base) {
if (!((base |= 0) >= 2))
base = DEFAULT_BASE;
return (base === 10 ? Math.log10(this) | 0
//
: base === 2 ? Math.log2(this) | 0
// TODO: base = 2^n
: Math.log(this) / Math.log(base) | 0) + 1;
// slow... should use multiply by exponents, or Math.clz32()
var natural = Math.floor(Math.abs(this)), digits = 0;
do {
digits++;
} while ((natural = Math.floor(natural / base)) > 0);
return digits;
}
// reverse the digits
// assert: {ℕ⁰:Natural+0}((this))
// TODO: 處理小數/負數/大數
function Number_reverse(base) {
if (!base)
base = DEFAULT_BASE;
// ABSORBING_ELEMENT
var natural = +this, reversed = 0;
while (natural > 0) {
reversed = reversed * base + (natural % base);
natural = Math.floor(natural / base);
}
return reversed;
}
function Number_is_palindromic(base) {
return this === this.reverse(base);
}
// for palindromic number or numeral palindrome 迴文數, 回文數
// http://articles.leetcode.com/2012/01/palindrome-number.html
function String_is_palindromic(chars) {
if (!chars)
return false;
for (var index = 0, l_index = chars.length - 1; index < l_index; index++, l_index--)
if (chars.charAt(index) !== chars.charAt(l_index))
return false;
return true;
}
function Number_tail(natural, base) {
// TODO
;
}
// ---------------------------------------------------------------------//
/**
* Test if the array is an arithmetic progression.<br />
* 判斷 ((this)) 是否為等差數列/連續整數。<br />
* O(n)
*
* @param {String}type
* 'integer': arithmetic integers,<br />
* 'consecutive': consecutive integers 連續整數型別, if the array
* contains only consecutive integers / consecutive values;<br />
* 'odd': odd consecutive integers,<br />
* 'even': even consecutive integers
*
* @returns {Boolean}is AP
*
* @see https://simple.wikipedia.org/wiki/Consecutive_integer
*/
function Array_is_AP(type) {
var length = this.length;
if (length <= 1)
return length === 1;
var number = this[1],
/** {Boolean}為奇數或偶數型別。 */
parity = type === 'odd' || type === 'even';
if (parity && number % 2 !== (type === 'odd' ? 1 : 0))
return false;
var difference = number - this[0];
if (type && difference !== (parity ? 2 : type === 'consecutive' ? 1
// 當前只要設定 type，皆為整數型別。
: Math.floor(difference)))
return false;
for (var index = 2; index < length; index++)
if (this[index] !== (number += difference))
return false;
return true;
}
/**
* <code>
Sum[m + n ((M - m)/(l - 1)), {n, 0, -1 + l}]
=
l(m+M)/2
Sum[(m + n (M - m)/(l - 1) - b)^2, {n, 0, l - 1}]
=
(l (6 b^2 (l-1)-6 b (l-1) (m+M)+(2 l-1) m^2+2 (l-2) m M+(2 l-1) M^2))/(6 (l-1))
</code>
*/
/**
* Test if the array combines an arithmetic progression.<br />
* 判斷 ((this)) 是否可組成等差數列/連續整數，不計較次序。<br />
* O(n)
*
* @param {String}type
* 'integer': arithmetic integers,<br />
* 'consecutive': consecutive integers 連續整數型別, if the array
* contains only consecutive integers / consecutive values;<br />
* 'odd': odd consecutive integers,<br />
* 'even': even consecutive integers
* @param {Integer}MIN
* acceptable minimum
*
* @returns {Number}difference or {Boolean}false
*
* @see https://simple.wikipedia.org/wiki/Consecutive_integer
*/
function Array_combines_AP(type, MIN) {
var length = this.length;
if (length <= 1)
return length === 1;
var min = Infinity,
// maximum
max = -Infinity,
// ABSORBING_ELEMENT
sum = 0,
// 不能僅由 min/max/sum 即定奪是否等差。
// 但尚未確認如此條件即已充分!!
square_sum = 0, square_base = this[0],
/** {Boolean}為奇數或偶數型別。 */
parity = type === 'odd' || type === 'even';
if (this.some(function(number) {
if (number < min) {
min = number;
if (min < MIN)
return true;
}
if (max < number)
max = number;
if (parity ? max - min > 2 * (length - 1) : type === 'consecutive'
&& max - min > length - 1)
return true;
sum += number;
number -= square_base;
square_sum += number * number;
}))
return false;
if (2 * sum !== (max + min) * length)
return false;
// check sum of square
var min_p_max = min + max, min_m_max = min * max;
if (square_sum !== ((2 * length * (min_p_max * min_p_max - min_m_max)
- min_p_max * min_p_max - 2 * min_m_max)
/ (6 * (length - 1)) - square_base * (min_p_max - square_base))
* length)
return false;
if (parity && min % 2 !== (type === 'odd' ? 1 : 0))
return false;
var difference = (max - min) / (length - 1);
if (type && difference !== (parity ? 2 : type === 'consecutive' ? 1
// 當前只要設定 type，皆為整數型別。
: Math.floor(difference)))
return false;
return difference;
}
// ---------------------------------------------------------------------//
// 檢查是否有重複的數字。
// use CeL.PATTERN_duplicated(string) to test if string contains duplicated
// chars.
_.PATTERN_duplicated = /(.).*?\1/;
// all the same characters
_.PATTERN_all_same = /^(.)\1*$/;
// http://en.cppreference.com/w/cpp/algorithm/is_permutation
function Array_is_permutation(sequence_2) {
var sequence_1 = this, last = sequence_1.length;
if (last !== sequence_2.length)
return false;
if (sequence_1 === sequence_2)
return true;
while (last-- > 0 && sequence_1[last] === sequence_2[last])
;
if (last < 0)
return true;
last++;
// assert: sequence_1, sequence_2 有不同。
var start = 0;
while (sequence_1[start] === sequence_2[start])
start++;
// count elements
var processed = new Set();
for (; start < last; start++) {
var element = sequence_1[start];
if (processed.has(element))
// skip element counted
continue;
processed.add(element);
// O(n^2)
var index = start, difference = sequence_2[index++] === element ? 0
: 1;
for (; index < last; index++) {
if (sequence_1[index] === element)
difference++;
if (sequence_2[index] === element)
difference--;
}
if (difference !== 0)
return false;
}
return true;
}
// 純量變數 (scalar variable)
function String_is_permutation(sequence_2) {
var sequence_1 = this, last = sequence_1.length;
if (last !== sequence_2.length)
return false;
// 直接比較較快。
if (sequence_1 === sequence_2)
return true;
// processed elements
var processed = Object.create(null);
for (var start = 0; start < last; start++) {
var element = sequence_1.charAt(start);
if (element in processed)
// skip element counted
continue;
processed[element] = null;
if (sequence_1.count_of(element, start + 1) + 1
// O(n^2)
!== sequence_2.count_of(element)) {
return false;
}
}
return true;
}
/**
* 關於尋找相同排列的數字，亦可採用紀錄各數字和的方法。<br />
* 另外，若數字之數量遠小於計算 .is_permutation() 之工作量，則 cache .sort() 反而會快很多。
*/
function Number_is_permutation(sequence_2) {
return this.digit_sum() === (+sequence_2).digit_sum()
// ↑ 先測試數字和是否相同。
&& String(this).is_permutation(String(sequence_2));
}
// ---------------------------------------------------------------------//
/**
* 按升序/降序排列處理每一序列至最後。<br />
* contain exactly the same digits, but in a different order.
*
* 注意: 不會先做排序!
*
* @param {Function}handler
* 處理 function
* @param {Boolean}descending
* default: ascending (small→big 升序序列為最小排列，降序序列為最大的排列), or will
* be descending (big→small 降序)
* @param {Boolean}inplace
* no clone, do not clone array.
*
* @returns {Array}the last array processed
*/
function Array_for_permutation(handler, descending, inplace) {
var array = inplace ? this : this.clone(), last_index = array.length - 1;
if (last_index >= 0)
handler(array);
if (last_index < 1)
return;
var index;
do {
// 求出下一個按升序排列序列。
// http://en.cppreference.com/w/cpp/algorithm/next_permutation
// http://leonard1853.iteye.com/blog/1450085
// http://www.cplusplus.com/reference/algorithm/next_permutation/
index = last_index;
for (var now = array[index], _next; index > 0;) {
// search [index]=now < [index+1]=_next
_next = now;
now = array[--index];
if (descending ? now > _next : now < _next) {
var later_index = last_index;
// search [index]=now < [later_index]=_next
while (true) {
_next = array[later_index];
if (descending ? _next < now : _next > now)
break;
later_index--;
}
// swap [index]=now, [later_index]=_next
array[later_index] = now;
array[index] = _next;
// reverse elements 元素: [index+1] to [last_index]
for (index++, later_index = last_index; index < later_index; index++, later_index--) {
_next = array[index];
array[index] = array[later_index];
array[later_index] = _next;
}
if (handler(array))
index = 0;
break;
}
}
} while (index > 0);
return array;
}
function String_for_permutation(handler, descending, sort) {
if (this.length < 2) {
handler(this);
return;
}
var array = this.split('');
if (sort)
typeof sort === 'function' ? array.sort(sort) : array.sort();
return array.for_permutation(function(array) {
return handler(array.join(''));
}, descending, true);
}
function Number_for_permutation(handler, descending, sort, base) {
if (!base)
base = 10;
if (this < base) {
handler(this);
return;
}
var array = this.digits();
if (sort) {
if (typeof sort === 'function')
array.sort(sort);
else
array.sort();
}
return +array.for_permutation(function(array) {
return handler(+array.join(''));
}, descending, true).join('');
}
// ---------------------------------------------------------------------//
// combinatorics 組合數學
/**
* <code>
[1,2,4,8,16,32].for_combination(3,function(s){console.log(s);})
CeL.for_combination(6,3,function(s){console.log(s);})
TODO:
CeL.for_combination(6,3,function(s){console.log(s);},true)
</code>
*/
// 自 elements 中提取出 select 個元素的組合方法。
// 共 C(elements, select) = elements! / select! / (elements-select)! 種組合數量。
// next_combination
// select ((select)) elements, ((select))-selection
function for_combination(elements, select, handler, descending) {
if (!((select |= 0) > 0))
// nothing select
return;
var map;
if (Array.isArray(elements)) {
// elements as map array
map = elements;
elements = map.length;
}
var index = 0,
/** {Array}index array */
selected = [];
// initialization
for (; index < select; index++)
selected.push(index);
while (true) {
// TODO: descending
handler(map ? selected.map(function(index) {
return descending ? map[elements - index] : map[index];
}) : descending ? selected.map(function(index) {
return elements - index - 1;
}).reverse() : selected);
index = 1;
for (; index < select
&& selected[index] === selected[index - 1] + 1; index++)
;
if (++selected[--index] === elements)
break;
--index;
if (selected[index] !== index)
for (; index >= 0; index--)
selected[index] = index;
}
}
_.for_combination = for_combination;
function Array_for_combination(select, handler, descending, inplace) {
return for_combination(inplace ? this : this.clone(), select, handler,
descending);
}
// ---------------------------------------------------------------------//
// export 導出.
_.MULTIPLICATIVE_IDENTITY = MULTIPLICATIVE_IDENTITY;
_.ZERO_EXPONENT = ZERO_EXPONENT;
_.ABSORBING_ELEMENT = ABSORBING_ELEMENT;
_.MULTIPLICATION_SIGN = MULTIPLICATION_SIGN;
library_namespace.set_method(String.prototype, {
is_palindromic : set_bind(String_is_palindromic),
is_permutation : String_is_permutation,
for_permutation : String_for_permutation
});
library_namespace.set_method(String, {
is_palindromic : String_is_palindromic
});
library_namespace.set_method(Number.prototype, {
// division, divided_by
divided : set_bind(division_with_remainder, true),
floor_sqrt : set_bind(floor_sqrt),
ceil_log : set_bind(ceil_log),
digits : Number_digits,
digit_sum : Number_digit_sum,
digit_length : Number_digit_length,
reverse : Number_reverse,
is_palindromic : Number_is_palindromic,
is_permutation : Number_is_permutation,
for_permutation : Number_for_permutation
});
library_namespace.set_method(Array.prototype, {
to_int_or_bigint : array_to_int_or_bigint,
is_AP : Array_is_AP,
combines_AP : Array_combines_AP,
is_permutation : Array_is_permutation,
for_permutation : Array_for_permutation,
for_combination : Array_for_combination
});
library_namespace.set_method(Math, {
hav : hav
});
return (_// JSDT:_module_
);
}
Reference in New Issue View Git Blame Copy Permalink