| Author |
Message |
|
thornahawk
μολών λαβέ
Active Member
Joined: 27 Mar 2005
Posts: 569
|
 Posted: 08 Feb 2006 08:39:46 pm Post subject: |
|
|
More info here.
The following TI-BASIC program uses a recurrence with caching for computing the Bernoulli numbers:
Code:
PROGRAM:BERN
Prompt N
If fPart(N)
Stop
If dim(∟BN)=1
{1,-.5,6ֿ¹}→BN
For(M,dim(∟BN),N)
.5–sum(augment({M+1}ֿ¹,seq((M nCr K)∟BN(K+1)/(1+M–K),K,2,M–1)))→∟BN(M+1)
End
max(not(fPart(.5N)),1=N)∟BN(N+1)
Usage:
When running the program for the first time, initialize with
{1}→BN
(dirty trick, I know.  )
before running BERN. You are prompted for the index, and the program will return the corresponding Bernoulli number. Do not alter list ∟BN in any way.
thornahawk |
|
| Back to top |
|
|
DarkerLine
ceci n'est pas une |
Super Elite (Last Title)
Joined: 04 Nov 2003
Posts: 8328
|
| Posted: 08 Feb 2006 09:22:04 pm Post subject: |
|
|
Some non-calculation-altering optimizations and suggestions:
Use the sequential variable n in place of N, and take out the fPart() check. If the user tries to enter any negative or fractional value in that case, the program will exit with an ERR:DOMAIN error.
It looks like you're using the size of the list as a check to see if it's been initialized. Use SetUpEditor instead in the following way:
SetUpEditor BN
If not(dim(∟BN
{1,-.5,6ֿ¹→BN
This will create the list if it hasn't been created; on the TI-83 Plus, it will also prevent crashes in case the list is archived.
As for that complicated-looking equation, I wouldn't go near that with a ten-foot link cable. 
Last edited by Guest on 08 Feb 2006 09:22:25 pm; edited 1 time in total |
|
| Back to top |
|
|
thornahawk
μολών λαβέ
Active Member
Joined: 27 Mar 2005
Posts: 569
|
| Posted: 09 Feb 2006 02:01:07 am Post subject: |
|
|
I'll incorporate your suggestions after I finish this shift I'm currently in. :)
"...that complicated-looking equation..."
You think?  That was a straightforward application of formula 30 in the link I gave above. ;)
(edit: For those too lazy to click
)
thornahawk
Last edited by Guest on 09 Feb 2006 02:03:55 am; edited 1 time in total |
|
| Back to top |
|
|
alexrudd
pm me if you read this
Bandwidth Hog
Joined: 06 Oct 2004
Posts: 2335
|
| Posted: 09 Feb 2006 07:20:56 pm Post subject: |
|
|
DarkerLine wrote: As for that complicated-looking equation, I wouldn't go near that with a ten-foot link cable. [post="69445"]<{POST_SNAPBACK}>[/post]
Code:
(fPart(.5N) xor 1!=N)∟BN(N+1
is two bytes smaller than
Code:
max(not(fPart(.5N)),1=N)∟BN(N+1
Whew! That's my headache for the day. |
|
| Back to top |
|
|
thornahawk
μολών λαβέ
Active Member
Joined: 27 Mar 2005
Posts: 569
|
|
| Back to top |
|
|
Weregoose
Authentic INTJ
Super Elite (Last Title)
Joined: 25 Nov 2004
Posts: 3976
|
| Posted: 14 Jul 2006 12:30:55 am Post subject: |
|
|
After glancing between my calculator and a couple math-related websites for about a month now, I've finally made a program which will calculate Bernoulli numbers:
If Ans<0 or fPart(Ans
Stop
.5Ans→X
If 2Ans≤9 or fPart(Ans
Then
{0,‾.5,1,6‾¹,‾30‾¹,42‾¹,‾30‾¹
Ans(1+(X=.5)+int(X+2)not(fPart(X
Else
(‾1)^(Ans-1)2(2Ans)!/(2[font="times new roman"]π)^(2Ans)sum(seq(X^(‾2Ans),X,1,16
End
Ans►Frac
Alright, so the only real math calculations it does is for even, integer inputs greater than nine.† While this is the closest I've been, I'm not really satisfied. Results are only given in the form of p/q for a maximum input of 22. But, oh well—it's incredibly fast considering the other things I've tried. Theoretically, this program puts the calculator one step away from the Riemann zeta function, Euler polynomials, and other stuff most people have never heard about. If you can comment on this sort of thing, great. I only wish I knew how Charles Babbage had put it together.
†Odd numbers—excepting one—do return zero in "real life", so odd-number inputs remain accurate.
–Goose
Last edited by Guest on 11 Jul 2010 06:42:56 pm; edited 1 time in total |
|
| Back to top |
|
|
thornahawk
μολών λαβέ
Active Member
Joined: 27 Mar 2005
Posts: 569
|
| Posted: 14 Jul 2006 06:54:12 am Post subject: |
|
|
If rationalizing the output is needed, it can be run through DE2FR. Although, it might not always return the expected answer.
*looks at Goose's snippet* I see, a table lookup and Riemann zeta.  I'll post the Riemann zeta program I have a little bit later...
thornahawk
Last edited by Guest on 14 Jul 2006 08:17:00 am; edited 1 time in total |
|
| Back to top |
|
|
Weregoose
Authentic INTJ
Super Elite (Last Title)
Joined: 25 Nov 2004
Posts: 3976
|
| Posted: 14 Jul 2006 07:13:53 am Post subject: |
|
|
I do keep that fraction program very handy. I use it as often as I'm dealing with what I believe to be rational values. Sometimes otherwise.
Last edited by Guest on 21 Jul 2006 07:52:49 pm; edited 1 time in total |
|
| Back to top |
|
|
thornahawk
μολών λαβέ
Active Member
Joined: 27 Mar 2005
Posts: 569
|
| Posted: 14 Jul 2006 08:35:40 am Post subject: |
|
|
About implementing Riemann ζ:
http://www.cecm.sfu.ca/personal/pborwein/PAPERS/P155.pdf
The Riemann ζ function program:
Code:
PROGRAM:ZETA
17→N:1→S
For(J,0,2N–1)
S/(J+1)^X→∟S(J+1)
0→∟E(J+1)
-S→S
If J≥N
∟E(J)+N nCr (J–N)→∟E(J+1)
End
(∟E/2^N–1)∟S→S
DelVar ∟E
SortA(∟S)
sum(∟S)→Y
DelVar ∟S
Y/(2^(1–X)–1)→Y
which works for any complex argument X.
N is an adjustable parameter that controls the program's accuracy. The default value in the program typically gives results accurate up to the last digit.
If you wonder what all the complication in the above code was for, this is to mitigate round-off error in the algorithm for certain complex values of the argument.
thornahawk
Last edited by Guest on 05 Jun 2007 03:32:43 am; edited 1 time in total |
|
| Back to top |
|
|
Weregoose
Authentic INTJ
Super Elite (Last Title)
Joined: 25 Nov 2004
Posts: 3976
|
| Posted: 15 Jul 2006 11:07:48 pm Post subject: |
|
|
How about some Euler numbers: int(.5+4^(Ans+1)(2Ans)!/(Ans³/910-.014Ans²+.215Ans+.8+[font="times new roman"]π^(2Ans+1
This was crafted from a lower bound and the contorted result of a cubic regression.
The error for each term (as far as the calculator's precision goes) is shown below:
[font="courier new"]Input______Error
0__________0
1__________0
…__________…
10_________‾210
11_________‾11000
12_________0
13_________0
…__________…
Therefore, suffixing the expression accordingly...
int(.5+4^(Ans+1)(2Ans)!/(Ans³/910-.014Ans²+.215Ans+.8+[font="times new roman"]π^(2Ans+1))+(Ans>9)(210+10790(Ans-10
...should remove all error.
The highest input it will allow is 29, because the calculator doesn't like the value 4^31*60! (ERR:OVERFLOW).
[EDIT]
It seems to flunk on most inputs of int(x)+1/2. These are the expected terms. Sorry. :lol:
[EDIT×2]
Although much slower (especially with higher terms)...
{1,1
For(X,1,77
augment(Ans,{.5sum(seq(X nCr KAns(K+1)Ans(X-K+1),K,0,X
Disp Ans(X
End
...this recursion gives all the numbers (up to X=77) flawlessly.
It should be easy to speed this up remarkably. I'll have that next edit.
[EDIT×3]
Okay, maybe I won't. 
Last edited by Guest on 11 Jul 2010 06:44:27 pm; edited 1 time in total |
|
| Back to top |
|
|
thornahawk
μολών λαβέ
Active Member
Joined: 27 Mar 2005
Posts: 569
|
| Posted: 16 Jul 2006 04:34:14 am Post subject: |
|
|
Ah, Euler numbers...
In the spirit of the Bernoulli number program I posted:
Code:
PROGRAM:EULR
Prompt N
If fPart(N)
Stop
If dim(∟EN)=1
{1,-1}→EN
For(M,dim(∟EN),N/2)
-1–sum(seq(((2M) nCr (2K))∟EN(K+1),K,1,M–1))→∟EN(M+1)
End
fPart(N/2)→M
(M=0)∟EN(N/2–M+1)
which is also initialized with the snippet
{1}→EN;
otherwise, the program can be modified (according to DarkerLine's suggestion) to use SetUpEditor and the sequence variable n. ;)
thornahawk |
|
| Back to top |
|
|
Weregoose
Authentic INTJ
Super Elite (Last Title)
Joined: 25 Nov 2004
Posts: 3976
|
| Posted: 23 Jan 2009 03:07:04 pm Post subject: |
|
|
Ans→X
{1:For(A,2,X,2
If not(fPart(X/A
augment(Ans,A{1=max(gcd(A+1,seq(A,A,1,A
End
-int(-prod(Ans+1){2X!/(2[font="times new roman"]π)^X,1
This is a work in progress. It combines the von Staudt–Clausen theorem with an asymptotic series for B2n to compute both the numerator and denominator of Bn directly. Non-even inputs which should return zero are incorrect—easily fixable—and the numerators of B24 and Bn≥28 are off, but not entirely; choosing a different series might improve the accuracy for more numbers. And from what I gather, the denominators are 100% perfect, although seeing many of them requires separating prod(Ans+1 from the series and displaying it in the line above, since the numerators tend to overflow so easily (B500 has 743 digits in the numerator, but just seven in the divisor!).
Last edited by Guest on 11 Jul 2010 06:38:45 pm; edited 1 time in total |
|
| Back to top |
|
|
thornahawk
μολών λαβέ
Active Member
Joined: 27 Mar 2005
Posts: 569
|
| Posted: 26 Jan 2009 12:29:36 pm Post subject: |
|
|
I had an experimental routine that also made use of the von Staudt–Clausen theorem, only that I used a precomputed list of primes. It was essentially the algorithm suggested here and used by the program YACAS (scroll down to section 6.7). It happened to be substantially longer than the routine in the first post, so I never bothered to post it...
thornahawk |
|
| Back to top |
|
|
Weregoose
Authentic INTJ
Super Elite (Last Title)
Joined: 25 Nov 2004
Posts: 3976
|
| Posted: 26 Jan 2009 01:40:44 pm Post subject: |
|
|
Thanks for your reply... Oops, looks like I left an extra line of code up there. *deletes*
Many of my programs would stretch a great distance if I allowed them to!
Last edited by Guest on 26 Jan 2009 01:44:23 pm; edited 1 time in total |
|
| Back to top |
|
|
thornahawk
μολών λαβέ
Active Member
Joined: 27 Mar 2005
Posts: 569
|
| Posted: 06 Feb 2009 07:20:06 am Post subject: |
|
|
At the moment, I have nothing particularly useful to add to this discussion, so I will instead share this old paper by Knuth and Buckholtz: [attachment=2592:bernnum.pdf]
thornahawk
(edit: I have replaced the external link with an optimized version)
Last edited by Guest on 07 Feb 2009 11:05:12 pm; edited 1 time in total |
|
| Back to top |
|
|
thornahawk
μολών λαβέ
Active Member
Joined: 27 Mar 2005
Posts: 569
|
| Posted: 18 Feb 2010 11:16:30 am Post subject: |
|
|
I had encountered the so-called Akiyama-Tanigawa algorithm for computing the Bernoulli numbers around the time before I posted the first message in this thread, but only realized the *proper* way of implementing it on a TI calculator recently (due to somewhat unrelated work).
Naively, the AT algorithm works this way:
Code:
PROGRAM:BERNAT
Prompt N
{0}→AT
For(K,1,N+1)
Kֿ¹→∟AT(K)
For(J,K-1,1,-1)
J(∟AT(J)-∟AT(J+1))→∟AT(J)
End
Pause ∟AT►Frac \\ if you want to see the whole AT array, remove this line otherwise
End
∟AT(1)
(note that in the AT array, B1 has the value 1/2; P. Luschny gives some convincing arguments as to why this is a more reasonable definition)
However, in this form, the AT algorithm is ill-suited for computing with inexact non-integers (try the program above and you'll see what I mean).
A solution would be to use two arrays, one for the numerator and one for the denominator. Thus:
Code:
PROGRAM:BERNAT
Prompt N
{0}→A1:{0}→A2
For(K,1,N+1)
1→∟A1(K)
K→∟A2(K)
For(J,K-1,1,-1)
lcm(∟A2(J),∟A2(J+1))→V
J((V/∟A2(J))∟A1(J)-(V/∟A2(J+1))∟A1(J+1))→U
gcd(U,V)→W
U/W→∟A1(J)
V/W→∟A2(J)
End
Pause ∟A1/∟A2►Frac \\ if you want to see the whole AT array, remove this line otherwise
End
∟A1(1)/∟A2(1)
This works much better, at least up until you hit the limits of the built-in gcd(/lcm( functions. The way to circumvent this problem is to explicitly implement the Euclidean algorithm to replace the corresponding calls to gcd( and lcm(, but around that point where the built-ins start failing, one might want to consider using the zeta series instead for evaluating the Bernoulli numbers.
thornahawk |
|
| Back to top |
|
|
Weregoose
Authentic INTJ
Super Elite (Last Title)
Joined: 25 Nov 2004
Posts: 3976
|
| Posted: 23 Feb 2010 10:15:49 pm Post subject: |
|
|
Why is this form
[attachment=3097:CodeCogsEqn.png]
never written?
It seems a lot more lucid to me – almost dead-easy. What I mean is, I sure wish I'd seen it written that way when I first encountered Bernoulli numbers a long time ago! :biggrin:
Last edited by Guest on 23 Feb 2010 11:08:59 pm; edited 1 time in total |
|
| Back to top |
|
|
thornahawk
μολών λαβέ
Active Member
Joined: 27 Mar 2005
Posts: 569
|
| Posted: 24 Feb 2010 03:58:05 am Post subject: |
|
|
That's actually derivable from the identity used in the very first program of this thread. As you may have seen, that first algorithm I used can be quite unstable.
On the other hand, it could be profitable to rearrange that identity such that only the even-order Bernoulli numbers need to be cached (to be compensated for by explicitly substituting the value for B1)
thornahawk |
|
| Back to top |
|
|
TheStorm
Calc Guru
Joined: 17 Apr 2007
Posts: 1233
|
| Posted: 24 Feb 2010 10:30:04 am Post subject: |
|
|
*shudders after reading the name Bernoulli and then goes back to just barely passing his Diff Eq. class....
I am still amazed at how you two manage to put such abstract math concepts into TI-Basic when I barely understand them to begin with. |
|
| Back to top |
|
|
kinkoa
Member
Joined: 28 Jul 2009
Posts: 103
|
| Posted: 25 Feb 2010 11:22:05 am Post subject: |
|
|
| hmm this is cool its something i couldnt do |
|
| Back to top |
|
|
|
TI-Basic => TI-BASIC
