/** Functions for performing calculations to arbitrary precision. A good reference is Ronald W. Potter, "Arbitrary Precision Calculation of Selected Higher Functions." References to "Potter" and equation numbers are from this book. */ use pi2.frink // See http://en.literateprograms.org/Arbitrary-precision_elementary_mathematical_functions_(Python) /** This is a somewhat naive implementation of e^x which has a Maclaurin series of exp[x] = 1 + x + x^2/2! + x^3/3! + ... It could probably be done faster with a binary splitting algorithm. See eBinarySplitting.frink for a sample. */ arbitraryExp[x, digits=getPrecision[], debug=false] := { if debug println["in arbitraryExp[$x, $digits]"] origPrec = getPrecision[] setPrecision[digits+3] var s try { if x < 0 { s = 1.0 / arbitraryExp[abs[x],digits,debug] } else { xpower = 1. ns = 0. s = 1 n = 0. factorial = 1. while s != ns { s = ns term = xpower / factorial if debug println["term is $term"] ns = s + term xpower = xpower * x n = n + 1. factorial = factorial * n } if debug println["s is $s"] } setPrecision[digits] retval = 1.0 s if debug println["arbitraryExp[$x, $digits] returning $retval"] return retval } finally setPrecision[origPrec] } /** Natural log to arbitrary precision. This uses a cubic convergence algorithm (that is, the number of correct digits in the result approximately triple with each iteration) with adaptive precision using equation 3.47 in Potter. It is significantly faster than the previous algorithm that did not have adaptive precision and used a quadratic Newton's method algorithm. */ arbitraryLn[x, digits=getPrecision[], debug=false] := { if debug println["in arbitraryLn2[$x]"] if x <= 0 return "Error: Arbitrary logs of negative numbers not yet implemented." origPrec = getPrecision[] try { setPrecision[10] eps = 10.^-(digits+1) // A good initial estimate is needed if (x < 10^308) { y = ln[x] prec = 15 } else // TODO: Store ln[2] somewhere. { y = approxLog2[x] * ln[2] prec = 10 // Is this reasonable? approxLog2 has a lot of latitude. } if debug println["Epsilon is $eps"] // Use Newton's method do { setPrecision[prec] y2 = y le = arbitraryExp[-y, prec] zn = 1 - x le y = y - zn(1 + zn/2) if debug println["y is $y"] prec = prec * 3 if (prec > digits + 3) prec = digits + 3 } while (prec < digits) or (abs[y2-y] > eps) setPrecision[digits] retval = 1. y if debug println["arbitraryLn about to return $retval"] return retval } finally setPrecision[origPrec] } // Arbitrary-precision power x^p // This uses the relationship that x^p = exp[p * ln[x]] arbitraryPow[x, p, digits = getPrecision[], debug=false ] := { // TODO: Make this work much faster for integer and rational powers. prec = getPrecision[] try { workingdigits = digits + 2 if digits <= 12 workingdigits = digits + 4 setPrecision[workingdigits] ret = arbitraryExp[p * arbitraryLn[x, workingdigits, debug], workingdigits, debug] setPrecision[digits] return 1. * ret } finally setPrecision[prec] } // Arbitrary log to the base 10. arbitraryLog[x, digits=getPrecision[], debug=false] := { origPrec = getPrecision[] try { setPrecision[digits+2] // TODO: Store ln[10] somewhere. retval = arbitraryLn[x, digits+2] / arbitraryLn[10, digits+2, debug] setPrecision[digits] return 1. retval } finally setPrecision[origPrec] } // This method computes sine of a number to an arbitrary precision. // This method is actually a dispatcher function which conditions the values // and tries to dispatch to the appropriate method which will be most likely // to converge rapidly. arbitrarySin[x, digits=getPrecision[], debug=false] := { origPrec = getPrecision[] try { // If x is something like 10^50, we actually need to work with // 50+digits at this point to get a meaningful result. extradigits = max[0, ceil[approxLog2[abs[x]]/ 3.219]] // Approx. log 10 if debug println["Extradigits is " + extradigits] setPrecision[digits+extradigits+4] if debug println["Dividing by pi to " + (digits + extradigits + 4) + " digits"] pi = Pi.getPi[digits+extradigits+4] // Break up one period of a sinewave into octants, each with width pi/4 octant = floor[(x / (pi/4)) mod 8] // Adjust x into [0, 2 pi] x = x mod (2 pi) if debug println["Octant is $octant"] if debug println["Adjusted value is $x"] if octant == 0 val = arbitrarySinTaylor[x, digits] else if octant == 1 val = arbitraryCosTaylor[-(x - pi/2), digits] else if octant == 2 val = arbitraryCosTaylor[x - pi/2, digits] else if octant == 3 or octant == 4 val = -arbitrarySinTaylor[x-pi, digits] else if octant == 5 val = -arbitraryCosTaylor[-(x - 3/2 pi), digits] else if octant == 6 val = -arbitraryCosTaylor[x - 3/2 pi, digits] else val = arbitrarySinTaylor[x - 2 pi, digits] setPrecision[digits] return 1. * val } finally setPrecision[origPrec] } /* This method computes cosine of a number to an arbitrary precision. This method actually just calls arbitrarySin[x + pi/2] */ arbitraryCos[x, digits=getPrecision[]] := { origPrec = getPrecision[] // If x is something like 10^50, we actually need to work with // 50+digits at this point to get a meaningful result. extradigits = max[0, ceil[approxLog2[abs[x]]/ 3.219]] // Approx. log 10 if debug println["Extradigits is " + extradigits] if debug println["Dividing by pi to " + (digits + extradigits + 4) + " digits"] pi = Pi.getPi[digits+extradigits+4] try { setPrecision[digits+extradigits+4] pi = Pi.getPi[digits+extradigits+4] arg = x+pi/2 return arbitrarySin[arg, digits] } finally setPrecision[origPrec] } // Arbitrary-precision sine arbitrarySinTaylor[x, digits=getPrecision[], returnInterval = false] := { origPrec = getPrecision[] eps = 10.^-(digits+3) terms = new array try { setPrecision[digits+5] x = x / radian // Factor out radians if we use them pi = Pi.getPi[digits+5] x = x mod (2 pi) if x > pi x = x - 2 pi num = x sum = x term = 3 denom = 1 factor = -x*x terms.push[sum] do { prevSum = sum num = num * factor denom = denom * (term-1) * term part = num/denom sum = sum + part term = term + 2 terms.push[part] } while prevSum != sum // println["terms for sin is $term"] sum = sum[reverse[terms]] setPrecision[digits] return 1. * sum } finally setPrecision[origPrec] } // Cosine for arbitrary digits. We could write this in terms of the sine // function (cos[x] = sin[x + pi/2]) but it's faster and more accurate // (especially around pi/2) to write it as the Taylor series expansion. arbitraryCosTaylor[x, digits=getPrecision[]] := { origPrec = getPrecision[] eps = 10.^-(digits+4) terms = new array try { setPrecision[digits+4] x = x / radian // Factor out radians if we use them pi = Pi.getPi[digits+4] x = x mod (2 pi) // println["Effective x is $x"] num = 1 sum = 1 term = 2 denom = 1 factor = -x*x terms.push[sum] do { prevSum = sum num = num * factor denom = denom * (term-1) * term part = num/denom sum = sum + part // println["$term $part $sum"] term = term + 2 terms.push[part] } while prevSum != sum // println["terms for cos is $term"] sum = sum[reverse[terms]] setPrecision[digits] return 1. * sum } finally setPrecision[origPrec] } // Tangent for arbitrary digits. This is written in terms of // sin[x]/cos[x] but it seems to behave badly around pi/2 where cos[x] goes to // zero. // // TODO: Make this a series expansion with the tangent numbers. This might // be more efficient. // See: http://mathworld.wolfram.com/TangentNumber.html // also TangentNumbers.frink // which calculate these numbers directly and efficiently. // // We could also try using Newton's method to invert arctan[x] which // has a simple series expansion, // arctan[x] = sum[(-1)^k x^(2k+1) / (2k + 1), {k, 0, infinity}] // but this only converges for abs[x] <= 1, x != +/- i // // See: // Fast Algorithms for High-Precision Computation of Elementary Functions, // Richard P. Brent, 2006 // https://pdfs.semanticscholar.org/bf5a/ce09214f071251bfae3a09a91100e77d7ff6.pdf // tangentnumbers.frink arbitraryTan[x, digits=getPrecision[], debug=false] := { // If x is something like 10^50, we actually need to work with // 50+digits at this point to get a meaningful result. extradigits = max[0, ceil[approxLog2[abs[x]]/ 3.219]] // Approx. log 10 if debug println["Extradigits is " + extradigits] origPrec = getPrecision[] try { setPrecision[digits+extradigits+4] retval = arbitrarySin[x, digits+4] / arbitraryCos[x, digits+4] setPrecision[digits] return 1. * retval } finally setPrecision[origPrec] } // Polylogarithm. See: // https://en.wikipedia.org/wiki/Polylogarithm polylog[s, z, digits = getPrecision[]] := { // if x <= 0 // return "Error: Arbitrary logs of negative numbers not yet implemented." origPrec = getPrecision[] try { setPrecision[digits+3] eps = 10.^-(digits+1) sum = 1. * z term = 0 k = 2 do { term = z^k / k^s sum = sum + term k = k + 1 // println[sum] } while abs[term] > eps setPrecision[digits] retval = 1. sum return retval } finally setPrecision[origPrec] } // Binary logarithm (that is, logarithm to base 2.) binaryLog[x, digits = getPrecision[]] := { origPrec = getPrecision[] try { setPrecision[digits+3] x = 1. x y = 0 b = .5 while x < 1 { x = 2 x y = y - 1 } while x >= 2 { x = x / 2 y = y + 1 } setPrecision[15] epsilon = y * 10.^-(digits+3) setPrecision[digits+3] println["Epsilon is $epsilon"] do { x = x^2 if x >= 2 { x = x/2 y = y + b } b = b/2 //println["$y $x $b"] } while b > epsilon setPrecision[digits] return 1. * y } finally setPrecision[origPrec] } /* digits = 1 setPrecision[digits] pi = Pi.getPi[digits] num = pi/4 collapseIntervals[false] inum = new interval[num, num, num] println[arbitrarySin[num,digits]] println[arbitrarySin[-num,digits]] println[sin[num]] println[sin[inum]] println[] println[arbitraryCos[num,digits]] //println[arbitraryCosAround2Pi[num,digits]] println[arbitraryCos[-num,digits]] println[cos[num]] println[cos[inum]] println[] println[arbitraryTan[num,digits]] println[arbitraryTan[-num,digits]] println[tan[num]] println[tan[-num]] println[tan[inum]] setPrecision[2] g = new graphics ps = new polyline pc = new polyline for x=-20 pi to 20 pi step .1 { ps.addPoint[x,-arbitrarySin[x]] pc.addPoint[x,-arbitraryCos[x]] } g.add[ps] g.color[0,0,1] g.add[pc] g.show[] */