Logarithms

[GabrielRhodes]

Created a function for calculating logarithms, including a natural logarithm function. Euler's number is calculated using the Taylor series for e, which could possibly be improved with binary splitting if needed. The logarithm of a number is calculated using the algorithm published by Daniel Shanks here:

https://www.ams.org/journals/mcom/1954-08-046/S0025-5718-1954-0061464-9/S0025-5718-1954-0061464-9.pdf

Example driver code would be:
let a = new BigNumber('9')
console.log(a.logBase(3))
console.log(BigNumber.e().naturalLog())
// output:
// 2
// 1

Returns NaN in the event that the user passes Infinity or a number <= 0 as a parameter. Similarly, if the base of the log is set to 1, the result is also NaN since log(1, n) = undefined

[magmel48]

@MikeMcl could you please consider this PR?

[MikeMcl]

This PR has not been tested properly by its author and is not likely to be merged.

It's a decent effort but the implementations of the log and ln functions are not reliable or performant enough to be included here.

Here is a zip file of a folder containing some test files that compare the PR against decimal.js.

test.zip

The PR is in the bignumber.mjs file and I have made some corrections to it in bignumber_amended.mjs. With the corrections the implementations seem to be reliable although I have done limited testing.

Unzip the zip file, and from a command line in the resulting directory use Node.js to run the test files.

Run multiple logarithm tests:
$ node logs.mjs {count} {base}

Run multiple natural logarithm tests:
$ node lns.mjs {count}

Run a single logarithm test:
$ node log.mjs {n} {base} {decimal places} {rounding mode}

Run a single natural logarithm test:
$ node ln.mjs {n} {decimal places} {rounding mode}

To test the amended version, amend the import in each test file so that bignumber_amended.mjs is used instead of bignumber.mjs.

[MikeMcl]

@magmel48

Review complete. This PR (amended) could be published separately as extension methods and if the author wants to do that I could add a link to it in the README here but I am not going to include these implementations internally.

[ronknott]

Hi Mike    Sorry I have not emailed you earlier about a log function for BigNumber.  I emailed a few years ago about my online library of maths functions to go with your  BigNumbers. I am a huge admirer of that library of yours and, as a former maths lecturer at University of Surrey, have found it so useful as I use JavaScript to make my maths pages and Calculators.   A LG and exp function are now part of  an online Calculator sitting as a user interface which evaluates an input expressions to any number of dps and can converts them to continued fractions too. The number of dps required is also an input number on the Calculator and BNdps() looks it up to get the result accurate to that precision.You can try it, with full details of functions available atExact, Fractions and Quadratics Calculatorr-knott.surrey.ac.ukIf you want to use the code, copied below, or want to adapt it, please feel free. I have tested it but there may be pathological cases where it fails but I've not found any yet. With thanks for your all your work on BigNumber!Ron KnottThe code for log and exp is based on well-known maths formulas and methods (see reference in the code) for arbitrary precision computations and is more modern than Shanks method.  I calculate the log to base e then use one division to convert to base b.  gt and lt are the greater than and less than operations as functions, x.gt(y) means x>y. X.sqr() Is X^2const BN1=new BigNumber(1);function BNlog(x){    //using Maclaurin series if x>0:  y=(x-1)/(x+1); Log(x)=2(y+y^3/3+y^5/5+y^7/7...)  //see Omondi Computer Arithmetic Systems page 420  if(typeof x=="number"||typeof x=="string")x=new BigNumber(x);  var xIN=x,mant=0,e=BN1.exp(),epow,pow=0,z;  if(x.isNegative())HALT("!! Cannot take the log of a negative number: "+x)  else if(x.isZero())L=BN1    else { if(x.gt(e)){pow=-1;epow=BN1;while((z=epow.times(e)).lt(x)){epow=z;pow++;};                    x=x.div(epow);pow++;}; //putmsg("BNloga: "+x+" pow="+pow);         if(x.lt(BN1)){pow=0;epow=e;while(x.times(epow).lt(1)){pow--;epow=epow.times(e)};                 x=x.times(epow);pow--};         var y=x.minus(1).div(x.plus(1));         var ysqr=y.sqr(),L=y,term=y,i=1,zero=new BigNumber("1e-"+(BNdps()+9));         while(term.gt(zero)){i+=2;term=term.times(ysqr);L=L.plus(term.div(i))};         L=L.times(2).plus(pow); // putmsg("BNlog max term was "+i+" lib val="+Math.log(xIN));  };  return L};BigNumber.prototype.exp=function(){ //e^this for any value of this    if(this.isInteger()) return BNe.pow(this)     var sum=BN1,x=this,term=BN1,zero=new BigNumber("1e-"+(BNdps()+5)),den=0,sgn1;    if(this.isZero())return BN1;    if(this.isNegative())x=this.negated();    // ASSERT x is +ve, use e^x = 1+x+x^2/2!+x^3/3!+...    while(term.abs().gt(zero)&&den<500)    {den++;term=term.times(x).div(den);    sum=sum.plus(term)};    if(this.isNegative())sum=BN1.div(sum);    return sum };Sent from my iPadOn 4 May 2023, at 19:09, Michael M ***@***.***> wrote: @MikeMcl commented on this pull request. In bignumber.js:

@@ -1223,6 +1246,44 @@

})(); + // Perform logorithm using the method by Daniel Shanks https://www.ams.org/journals/mcom/1954-08-046/S0025-5718-1954-0061464-9/S0025-5718-1954-0061464-9.pdf + log = (function () { There is no reason to use an IIFE (immediately invoked function expression) here when you are only returning a function. It should just be function log(x, y, dp, rm) { In bignumber.js:

@@ -1223,6 +1246,44 @@

})(); + // Perform logorithm using the method by Daniel Shanks https://www.ams.org/journals/mcom/1954-08-046/S0025-5718-1954-0061464-9/S0025-5718-1954-0061464-9.pdf + log = (function () { + // x: base, y: number + return function (x, y, dp, rm) { + if (x.comparedTo(0) + y.comparedTo(0) !== 2 || !x.isFinite() || !y.isFinite() || y.valueOf() === '1') { There is no reason to call valueOf here and get the string value of a bignumber to see if it is 1. Use the built-in bignumber ONE and one of the comparison methods, e.g. y.eq(ONE), or ideally the compare function as it avoids unnecessary object creation. compare(y, ONE) === 0 In bignumber.js:

@@ -1223,6 +1246,44 @@

})(); + // Perform logorithm using the method by Daniel Shanks https://www.ams.org/journals/mcom/1954-08-046/S0025-5718-1954-0061464-9/S0025-5718-1954-0061464-9.pdf + log = (function () { + // x: base, y: number + return function (x, y, dp, rm) { + if (x.comparedTo(0) + y.comparedTo(0) !== 2 || !x.isFinite() || !y.isFinite() || y.valueOf() === '1') { + return BigNumber(NaN) + } + var A = y + var B = x + var C = BigNumber(1) + var D = BigNumber(0) + var E = BigNumber(0) + var F = BigNumber(1) + var s = 1 + var one = BigNumber('1') As stated above, use bignumber ONE internally. In bignumber.js:

+ // Perform logorithm using the method by Daniel Shanks https://www.ams.org/journals/mcom/1954-08-046/S0025-5718-1954-0061464-9/S0025-5718-1954-0061464-9.pdf

+ log = (function () { + // x: base, y: number + return function (x, y, dp, rm) { + if (x.comparedTo(0) + y.comparedTo(0) !== 2 || !x.isFinite() || !y.isFinite() || y.valueOf() === '1') { + return BigNumber(NaN) + } + var A = y + var B = x + var C = BigNumber(1) + var D = BigNumber(0) + var E = BigNumber(0) + var F = BigNumber(1) + var s = 1 + var one = BigNumber('1') + if (A.comparedTo(one) === -1) { if (A.lt(ONE)) { // or if (compare(A, ONE) === -1) { In bignumber.js:

+ log = (function () {

+ // x: base, y: number + return function (x, y, dp, rm) { + if (x.comparedTo(0) + y.comparedTo(0) !== 2 || !x.isFinite() || !y.isFinite() || y.valueOf() === '1') { + return BigNumber(NaN) + } + var A = y + var B = x + var C = BigNumber(1) + var D = BigNumber(0) + var E = BigNumber(0) + var F = BigNumber(1) + var s = 1 + var one = BigNumber('1') + if (A.comparedTo(one) === -1) { + A = one.div(A) Better to use the div function when performing division internally. As well as avoiding unnecessary object creation, it allows decimal places and rounding mode arguments to be passed in. I think some extra guard digits will be needed here also. // Add 10 guard digits and use rounding mode 1, which is truncation. A = div(ONE, A, dp + 10, 1); In bignumber.js:

+ var B = x

+ var C = BigNumber(1) + var D = BigNumber(0) + var E = BigNumber(0) + var F = BigNumber(1) + var s = 1 + var one = BigNumber('1') + if (A.comparedTo(one) === -1) { + A = one.div(A) + s *= -1 + } + if (B.comparedTo(one) === -1) { + B = one.div(B) + s *= -1 + } + for (let i = 0; i < dp * 1.1 + 5; i) { This loop could be better refactored to: const p = dp * 1.1 + 5; for (let i = 0; i < p && !B.eq(ONE);) { if (A.gte(B)) { [A, C, D] = [div(A, B, p, 1), C.plus(E), D.plus(F)]; } else { [A, B, C, D, E, F] = [B, A, E, F, C, D]; i++; } } In bignumber.js:

+ if (B.comparedTo(one) === -1) {

+ B = one.div(B) + s *= -1 + } + for (let i = 0; i < dp * 1.1 + 5; i) { + if (A.comparedTo(B) > -1 && B.valueOf() !== '1') { + [A, C, D] = [A.div(B), C.plus(E), D.plus(F)] + } else if (B.valueOf() === '1') { + break + } else { + [A, B, C, D, E, F] = [B, A, E, F, C, D] + i++ + } + + } + return round(E.div(F).multipliedBy(s), dp, rm) Use the internal div function, and I think that the sign can just be set directly rather than by multiplying. E.s = s; return div(E, F, dp, rm); In bignumber.js:

@@ -808,6 +810,27 @@

return sum; }; + /* + * Returns e to dp decimal places + * Runtime could be improved with binary splitting at expense of storage + * Runtime of 5,000 digits was ~3s + */ + + BigNumber.e = BigNumber.euler = function (dp, rm) { I'd prefer this function to be named exp, following the lead of Math.exp and decimal.js. In bignumber.js:

@@ -808,6 +810,27 @@

return sum; }; + /* + * Returns e to dp decimal places + * Runtime could be improved with binary splitting at expense of storage + * Runtime of 5,000 digits was ~3s + */ + + BigNumber.e = BigNumber.euler = function (dp, rm) { + let oldDp = DECIMAL_PLACES + dp = dp | oldDp I assume you meant dp = dp || oldDp here. Anyway, neither would handle a dp of 0 correctly. Proper argument handling would be: if (dp == null) { dp = DECIMAL_PLACES; rm = ROUNDING_MODE; } else { intCheck(dp, 0, MAX); if (rm == null) rm = ROUNDING_MODE; else intCheck(rm, 0, 8); } In bignumber.js:

@@ -808,6 +810,27 @@

return sum; }; + /* + * Returns e to dp decimal places + * Runtime could be improved with binary splitting at expense of storage + * Runtime of 5,000 digits was ~3s + */ + + BigNumber.e = BigNumber.euler = function (dp, rm) { + let oldDp = DECIMAL_PLACES + dp = dp | oldDp + BigNumber.set({DECIMAL_PLACES: dp + 1}) There is no reason to ever use the BigNumber.set method internally. You can just set the value of DECIMAL_PLACES directly (as long as you are within the clone function, which is most of the file). DECIMAL_PLACES = dp + 1; Although, there is need to change the value of DECIMAL_PLACES here anyway if you use the internal div function and just pass in the dp you want to use. e = e.plus(div(ONE, divisor, dp + 1, 1)); In bignumber.js:

@@ -808,6 +810,27 @@

return sum; }; + /* + * Returns e to dp decimal places + * Runtime could be improved with binary splitting at expense of storage + * Runtime of 5,000 digits was ~3s + */ + + BigNumber.e = BigNumber.euler = function (dp, rm) { + let oldDp = DECIMAL_PLACES + dp = dp | oldDp + BigNumber.set({DECIMAL_PLACES: dp + 1}) + let e = BigNumber('1') let e = BigNumber(ONE); In bignumber.js:

+ /*

+ * Returns e to dp decimal places + * Runtime could be improved with binary splitting at expense of storage + * Runtime of 5,000 digits was ~3s + */ + + BigNumber.e = BigNumber.euler = function (dp, rm) { + let oldDp = DECIMAL_PLACES + dp = dp | oldDp + BigNumber.set({DECIMAL_PLACES: dp + 1}) + let e = BigNumber('1') + let divisor = BigNumber('1') + let one = BigNumber('1') + for (let i = 2; i < dp * 10; i++) { + e = e.plus(one.div(divisor)) + divisor = divisor.multipliedBy(BigNumber(i)) Methods such as multipliedBy and the BigNumber constructor, accept a number, string or bignumber so it is just divisor = divisor.multipliedBy(i). In bignumber.js:

+ * Runtime could be improved with binary splitting at expense of storage

+ * Runtime of 5,000 digits was ~3s + */ + + BigNumber.e = BigNumber.euler = function (dp, rm) { + let oldDp = DECIMAL_PLACES + dp = dp | oldDp + BigNumber.set({DECIMAL_PLACES: dp + 1}) + let e = BigNumber('1') + let divisor = BigNumber('1') + let one = BigNumber('1') + for (let i = 2; i < dp * 10; i++) { + e = e.plus(one.div(divisor)) + divisor = divisor.multipliedBy(BigNumber(i)) + } + e = BigNumber(round(e, dp, rm, true).toString()) There is no need to call toString on the bignumber returned by round, and it is inefficient to do so. In bignumber.js:

+ */

+ + BigNumber.e = BigNumber.euler = function (dp, rm) { + let oldDp = DECIMAL_PLACES + dp = dp | oldDp + BigNumber.set({DECIMAL_PLACES: dp + 1}) + let e = BigNumber('1') + let divisor = BigNumber('1') + let one = BigNumber('1') + for (let i = 2; i < dp * 10; i++) { + e = e.plus(one.div(divisor)) + divisor = divisor.multipliedBy(BigNumber(i)) + } + e = BigNumber(round(e, dp, rm, true).toString()) + BigNumber.set({DECIMAL_PLACES: oldDp}) + return BigNumber(round(e, dp, rm, true).toString()) It should be just return round(e, dp, rm, true) The BigNumber constructor accepts bignumbers so you do not need to call toString to create a new bignumber from an existing one. It is also better to use new BigNumber rather than just use BigNumber which just recalls itself using new anyway. Actually, it should be return round(e, dp + e.e + 1, rm, true) because the round function's second argument is the required number of significant digits not the required number of decimal places. In bignumber.js:

@@ -808,6 +810,27 @@

return sum; }; + /* + * Returns e to dp decimal places + * Runtime could be improved with binary splitting at expense of storage + * Runtime of 5,000 digits was ~3s + */ + + BigNumber.e = BigNumber.euler = function (dp, rm) { + let oldDp = DECIMAL_PLACES + dp = dp | oldDp + BigNumber.set({DECIMAL_PLACES: dp + 1}) + let e = BigNumber('1') + let divisor = BigNumber('1') + let one = BigNumber('1') + for (let i = 2; i < dp * 10; i++) { dp might be 0, so it's not right to use dp * 10 here. In bignumber.js:

@@ -1669,6 +1730,18 @@

}; + /* + * Return a new BigNumber whose value is the value of log with base this BigNumber of the value + * BigNumber(y, b), rounded according to DECIMAL_PLACES and ROUNDING_MODE + */ + P.logBase = P.log = function (y, b) { + return log(this, new BigNumber(y, b), DECIMAL_PLACES, ROUNDING_MODE) + } + + P.naturalLog = P.ln = function () { + return log(this, BigNumber.euler(DECIMAL_PLACES), DECIMAL_PLACES, ROUNDING_MODE) The call to BigNumber.euler is going to need some guard digits. For example return log(this, BigNumber.euler(DECIMAL_PLACES + 5), DECIMAL_PLACES, ROUNDING_MODE) —Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>