Bernoulli numbers: Difference between revisions
Added AppleScript.
(Added AppleScript.) |
|||
Line 220:
B(60) -1215233140483755572040304994079820246041491 / 56786730
</pre>
=={{header|AppleScript}}==
I've only been able to get accurate results up to B(31) using AppleScript's number classes and operators and up to B(57) with the Foundation framework's NSDecimalNumber class in AppleScriptObjC. To handle up to B(60) and beyond, the script here replaces numbers with lists whose items are integers representing the numbers' decimal digits — which of course requires custom math routines.
<lang applescript>on bernoullis(n) -- Return a list of "numerator / denominator" texts representing Bernoulli numbers B(0) to B(n).
set listMathScript to getListMathScript(10) -- Script object providing custom list math routines.
set output to {}
-- Akiyama–Tanigawa algorithm for the "second Bernoulli numbers".
-- List 'a' will contain {numerator, denominator} lists for the fractions.
-- The numerators and denominators will in turn be lists containing integers as their (decimal) digits.
set a to {}
repeat with m from 0 to n
-- Append the structure for 1 / (m + 1) to the end of a.
set {numerator2, denominator2} to {{1}, listMathScript's intToList(m + 1)}
set a's end to result
repeat with j from m to 1 by -1
-- Retrieve the preceding fraction's numerator and denominator.
set {numerator1, denominator1} to a's item j
tell listMathScript
-- Reset {numerator2, denominator2} to (preceding fraction - new fraction) * j as follows:
-- Work out the two fractions' lowest common denominator and adjust the numerators accordingly.
set lcd to its lcm(denominator1, denominator2)
set numerator1 to its multiply(numerator1, its |div|(lcd, denominator1))
set numerator2 to its multiply(numerator2, its |div|(lcd, denominator2))
-- Subtract numerator2 from numerator1 and multiply the result by j.
set numerator2 to its multiply(its subtract(numerator1, numerator2), its intToList(j))
set denominator2 to lcd
end tell
-- Store the resulting fraction in a's slot j. Don't bother to reduce it.
set a's item j to {numerator2, denominator2}
end repeat
-- The fraction just stored in a's first slot is Bernoulli(m).
-- Reduce it and append a text representation to the output.
tell listMathScript
set gcd to its hcf(numerator2, denominator2)
set numerator2 to its |div|(numerator2, gcd)
set denominator2 to its |div|(denominator2, gcd)
set end of output to its listToText(numerator2) & (" / " & its listToText(denominator2))
end tell
end repeat
return output
end bernoullis
on getListMathScript(base)
script
on multiply(a, b) -- Multiply list a by list b.
set aLength to (count a)
set bLength to (count b)
set productLength to aLength + bLength - 1
set product to {}
repeat productLength times
set product's end to 0
end repeat
repeat with bIndex from -1 to -bLength by -1
set bDigit to b's item bIndex
if (bDigit is not 0) then
set carry to 0
set productIndex to bIndex
repeat with aIndex from aLength to 1 by -1
set subresult to bDigit * (a's item aIndex) + carry + (product's item productIndex)
set product's item productIndex to (subresult mod base)
set carry to subresult div base
set productIndex to productIndex - 1
end repeat
if (carry = 0) then
else if (productIndex < -productLength) then
set product's beginning to carry
else
set product's item productIndex to (product's item productIndex) + carry
end if
end if
end repeat
return product
end multiply
on subtract(a, b) -- Subtract list b from list a.
set aLength to (count a)
set bLength to (count b)
copy a to difference
copy b to b
repeat (bLength - aLength) times
set difference's beginning to 0
end repeat
repeat (aLength - bLength) times
set b's beginning to 0
end repeat
set differenceLength to (count difference)
repeat with i from 1 to differenceLength
set aDigit to difference's item i
set bDigit to b's item i
set |b > a| to (bDigit > aDigit)
if ((|b > a|) or (aDigit > bDigit)) then exit repeat
end repeat
if (|b > a|) then set {b, difference} to {difference, b}
set carry to 0
repeat with i from differenceLength to 1 by -1
set subresult to (difference's item i) - (b's item i) + carry + base
set difference's item i to (subresult mod base)
set carry to subresult div base - 1
end repeat
if (carry > 0) then set difference's beginning to carry
if (|b > a|) then invert(difference)
return difference
end subtract
on |div|(a, b) -- List a div list b.
return divide(a, b)'s quotient
end |div|
on |mod|(a, b) -- List a mod list b.
return divide(a, b)'s remainder
end |mod|
on divide(a, b) -- Divide list a by list b. Return a record containing separate lists for the quotient and remainder.
set dividend to trim(a)
set divisor to trim(b)
set dividendLength to (count dividend)
set divisorLength to (count divisor)
if (divisorLength > dividendLength) then return {quotient:{0}, remainder:a}
set dividendNegative to (dividend's beginning < 0)
if (dividendNegative) then invert(dividend)
set divisorNegative to (divisor's beginning < 0)
if (divisorNegative) then invert(divisor)
set quotient to {}
if (divisorLength is 1) then
set segment to {}
else
set segment to items 1 thru (divisorLength - 1) of dividend
end if
repeat with nextSlot from divisorLength to dividendLength
set end of segment to item nextSlot of dividend
set q to 0
repeat
set newSegment to trim(subtract(segment, divisor))
if (newSegment's beginning < 0) then exit repeat
set q to q + 1
set segment to newSegment
end repeat
set end of quotient to q
end repeat
-- The quotient's negative if the input signs are different. Otherwise positive.
if (dividendNegative ≠ divisorNegative) then invert(quotient)
-- The remainder's sign is always that of dividend.
if (dividendNegative) then invert(segment)
return {quotient:quotient, remainder:segment}
end divide
on lcm(a, b) -- Lowest common multiple of list a and list b.
return multiply(b, |div|(a, hcf(a, b)))
end lcm
on hcf(a, b) -- Highest common factor of list a and list b.
set a to trim(a)
set b to trim(b)
repeat until (b = {0})
set x to a
set a to b
set b to trim(|mod|(x, b))
end repeat
if (a's beginning < 0) then invert(a)
return a
end hcf
on invert(lst) -- Invert the sign of all lst's "digits".
repeat with thisDigit in lst
set thisDigit's contents to -thisDigit
end repeat
end invert
on trim(lst) -- Return a copy of the input list with no leading zeros..
repeat with i from 1 to (count lst)
if (lst's item i is not 0) then exit repeat
end repeat
return lst's items i thru end
end trim
on intToList(n) -- Return a list of numbers representing n's digits.
set lst to {n mod base}
set n to n div base
repeat until (n = 0)
set beginning of lst to n mod base as integer
set n to n div base
end repeat
return lst
end intToList
on listToText(lst) -- Return the number represented by the input list as text.
-- This lazily assumes 2 <= base <= 10. :)
set lst to trim(lst)
if (lst's beginning < 0) then
invert(lst)
set lst's beginning to "-"
end if
return join(lst, "")
end listToText
end script
return result
end getListMathScript
on join(lst, delim)
set astid to AppleScript's text item delimiters
set AppleScript's text item delimiters to delim
set txt to lst as text
set AppleScript's text item delimiters to astid
return txt
end join
on task()
set maxN to 60
set output to {""}
set padding to " = "
set bernoulliNumbers to bernoullis(maxN) -- Get all the results in one sweep.
repeat with n from 0 to maxN
set bernie to bernoulliNumbers's item (n + 1)
if (bernie does not start with "0") then
set Bn to "B(" & n & ")"
set output's end to Bn & ¬
text ((count Bn) - 3) thru (50 - (offset of "/" in bernie)) of padding & ¬
bernie
end if
end repeat
return join(output, linefeed)
end task
task()</lang>
{{output}}
<lang applescript>"
B(0) = 1 / 1
B(1) = 1 / 2
B(2) = 1 / 6
B(4) = -1 / 30
B(6) = 1 / 42
B(8) = -1 / 30
B(10) = 5 / 66
B(12) = -691 / 2730
B(14) = 7 / 6
B(16) = -3617 / 510
B(18) = 43867 / 798
B(20) = -174611 / 330
B(22) = 854513 / 138
B(24) = -236364091 / 2730
B(26) = 8553103 / 6
B(28) = -23749461029 / 870
B(30) = 8615841276005 / 14322
B(32) = -7709321041217 / 510
B(34) = 2577687858367 / 6
B(36) = -26315271553053477373 / 1919190
B(38) = 2929993913841559 / 6
B(40) = -261082718496449122051 / 13530
B(42) = 1520097643918070802691 / 1806
B(44) = -27833269579301024235023 / 690
B(46) = 596451111593912163277961 / 282
B(48) = -5609403368997817686249127547 / 46410
B(50) = 495057205241079648212477525 / 66
B(52) = -801165718135489957347924991853 / 1590
B(54) = 29149963634884862421418123812691 / 798
B(56) = -2479392929313226753685415739663229 / 870
B(58) = 84483613348880041862046775994036021 / 354
B(60) = -1215233140483755572040304994079820246041491 / 56786730"</lang>
=={{header|Bracmat}}==
| |||