Login [Register]
Don't have an account? Register now to chat, post, use our tools, and much more.
This is an archived, read-only copy of the United-TI subforum , including posts and topic from May 2003 to April 2012. If you would like to discuss any of the topics in this forum, you can visit Cemetech's TI-BASIC subforum. Some of these topics may also be directly-linked to active Cemetech topics. If you are a Cemetech member with a linked United-TI account, you can link United-TI topics here with your current Cemetech topics.

This forum is locked: you cannot post, reply to, or edit topics. TI-Basic => TI-BASIC
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. Razz )

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. Very Happy


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? Surprised 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. Very Happy  [post="69445"]<{POST_SNAPBACK}>[/post]
I wouldn't either, but a 20 foot link cable is just long enough to chop off a couple of parentheses.


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

Posted: 10 Feb 2006 08:04:50 am    Post subject:

Headache is right. Surprised I had to make a truth table to make sense of your patch. Razz Thank you, though. :D

And we shave off one byte more if it becomes


Code:
∟BN(N+1)(N≠1 xor fPart(.5N


Nice. Very Happy Thanks to you two.

thornahawk
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. Smile 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. Razz

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. :0x5:

[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. Cool

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! Smile


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. Smile 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
Display posts from previous:   
Register to Join the Conversation
Have your own thoughts to add to this or any other topic? Want to ask a question, offer a suggestion, share your own programs and projects, upload a file to the file archives, get help with calculator and computer programming, or simply chat with like-minded coders and tech and calculator enthusiasts via the site-wide AJAX SAX widget? Registration for a free Cemetech account only takes a minute.

» Go to Registration page
    »
» View previous topic :: View next topic  
Page 1 of 1 » All times are GMT - 5 Hours

 

Advertisement