Lecture Four:

I. BASIC PROGRAMMING So far, we've been using J in an interactive mode, perhaps saving interim results in variables. When we end the session, our work goes away. For example, we might define polynomial multiplication in one line, and use it to generate a table: mp =: +//.@(*/) NB. polynomial multiplication 1 1 mp^:(i.11) 1 NB. (x+1)^n, powers 0 through 10 1 0 0 0 0 0 0 0 0 0 0 Typically, however, we'll want to define a verb as a series of programmed steps, and have these verbs invoke one another. And we won't want to redefine these verbs from scratch every time we use them, meaning we'll want to save our code in persistent files. The interactive execution window is where to test fragments of a longer program, which, once debugged, may be copied to a text window and saved as a runnable file. When the file is run, the operators it defines (verbs, adverbs, nouns, conjunctions) will become available to us in some locale (more on locales another time). II. BERNOULLI NUMBERS For example, suppose we'd like to compute the Bernoulli numbers. Several algorithms are available, including one based on the above expansion of (x + 1)^n. The algorithm for deriving each successive Bernoulli number uses the binomial coefficients in each successive row of Pascal's Triangle, multiplied by all the Bernoulli numbers derived so far, beginning with B[0]=1: 0 = (2 * B[1]) + (1 * B[0]) 0 = (3 * B[2]) + (3 * B[1]) + (1*B[0]) 0 = (4 * B[3]) + (6 * B[2]) + (4*B[1]) + (1*B[0])or... B[0] =: 1 B[1] =: (_1r2) * (1 * B[0]) B[2] =: (_1r3) * (3 * B[1] + 1*B[0]) B[3] =: (_1r4) * (6 * B[2] + 4*B[1] + 1*B[0]) So one strategy is to collect our Bernoulli numbers in a list, and to use (x+1)^n, with n increasing, to generate the coefficients. Each step of the way, we'll use all our Bernoullis so far, with the binomial coefficients, to generate the next Bernoulli. Let's look at a single row of the binomial expansion from the table at the top of this page (ignore any zero padding on the right): 1 9 36 84 126 126 84 36 9 1 What we need to do is multiply all coefficients after _9 by the Bernoulli numbers collected so far. Then add the resulting terms and multiply by 1r9. This will give our next Bernoulli number, which we prefix to a growing list. Here's a program which codifies this algorithm: mp =: +//.@(*/) NB. polynomial multiplication bernoullis =: 3 : 0 coeffs =. x: 1 1 NB. beginning coefficients bns =. 1 NB. first Bernoulli while. y.>0 do. coeffs =. 1 1 mp coeffs NB. (x1)^n b =. ( % 1{coeffs ) * (+/ bns * 2}.coeffs) bns =. b,bns NB. growing list of Bernoullis y. =. <:y. end. bns ) III. ANALYSIS OF A SHORT PROGRAM If you've programmed before, your eye will likely pick out the looping structure, bracketed by the while and end statements. Automatically, J assumes y. is the name of the verb's right argument, which in this case is named bernoullis. Were this a dyadic verb, the left argument would be named x. (again, J simply assumes this). Asking for bernoullis 10 at the command line will initialize y. to 10. After a couple of variable initializations, the loop begins. It contains a decrementby1 statement: y. =. <:y. so that each time through the loop, y. gets smaller, and at zero, finally fails the y.>0 test (part of the while statement) and escapes the loop. Each time through the loop, we move to the next row of coefficients in the binomial expansion of (x1)^n. The 2nd line in the loop is the most complicated, and implements the algorithm specified above. Let's look at its details: coeffs =. 1 9 36 84 126 126 84 36 9 1 b =. ( % 1{coeffs ) * (+/ bns * 2}.coeffs) The first parenthesized expression takes element 1 of coeffs. Given zerobased indexing, this would be the 9 if coeffs were as shown above. This element (9) is then reciprocated (%) and negated (), giving _1r9. The result is in rational format because our coefficients were initialized as x: 1 1 before the loop started, meaning all subsequent polynomial multiplications have continued to use extended precision format. The second parenthesized expression uses dyadic curtail (}.) to take all but the first two elements of coeffs, i.e. in this case: 36 84 126 126 84 36 9 1. These terms are multiplied by the corresponding eight elements in bns, and the results are summed (+/). Both of these expressions are then multiplied together, giving our next Bernoulli number b. This local variable b is then prefixed to bns in the next line, y. is decremented by one (<:), and the loop begins another iteration. When we finally escape the loop, bns is what gets returned by this verb (whatever is on the last line constitutes the return value). bns contains a list of Bernoulli numbers in reverse order, i.e. the last Bernoulli computed begins the list. bernoullis 10 5r66 0 _1r30 0 1r42 0 _1r30 0 1r6 1r2 1 bernoullis 11 0 5r66 0 _1r30 0 1r42 0 _1r30 0 1r6 1r2 1 Note that oddindexed Bernoulli numbers are 0, after the first, which is _1r2. If we want the nth Bernoulli number, we can just take the first element of what's returned, using head ({.): load 'd:\program files\j501a\user\projects\berns.ijs' {. bernoullis 50 NB. 50th Bernoulli number 495057205241079648212477525r66 We could also codify this as another verb: bernoulli =: 3 : '{. bernoullis y.' bernoulli 20 NB. 20th Bernoulli number _174611r330 IV. PARTS OF SPEECH Also in need of explanation in the above programs, is the use of the number 3, or the expression 3 : 0, at the outset of a verb's definition. This 3 is one of four possible values, 1, 2, 3 or 4, and it refers to the part of speech that's being defined. 1 signifies an adverb, 2 a conjunction, 3 a monadic verb with an optional dyadic meaning, and 4 an exclusively dyadic verb. The : 0 after the 3 indicates that the rest of the definition follows on subsequent lines, until a concluding right parenthesis, alone on a line, is reached. In place of this 0, the complete definition enclosed in single quotes, may be provided on a single line. This was the approach taken in our definition of bernoulli. V. BERNOULLI POLYNOMIALS What might we do with Bernoulli numbers? Originally, they were developed, not by Jakob Bernoulli (16541705), but by Johann Faulhaber (15801635) writing in 1631. Johann wanted to find the sum of the first n1 consecutive integers all to the kth power, e.g. 1^4 + 2^4 + 3^4 ... + 100^4. Let's compute the answer using J in simple calculator mode: +/ (i.>:100x)^4 NB. i.e. 1^4 + 2^4 + ... + 100^4 This sum may be expressed in terms of the 5th Bernoulli polynomial, the coefficients of which are comprised of the first five Bernoulli numbers weighted (multiplied) by the binomial coefficients k ! 5 (k=1...5). bcoeffs =: 3 : 0 The polynomial we seek, giving the sum of the first n1 consecutive integers to the kth power, is comprised of the k terms of the k+1 Bernoulli polynomial, multiplied by 1/(k+1). A prefix of 0 ensures that our coefficients align with corresponding powers, i.e. the last term aligns with n^5. ] b5 =: 0, }. bcoeffs 5 And if the question were about consecutive integers to the 10th power, we would use Bernoulli numbers to get the coefficients of the 11th degree polynomial, and so on. Computing the sum of the first 1000 integers to the 10th power would look like this: (1r11 * 0,}.bcoeffs 11) p. 1001 91409924241424243424241924242500 +/ (>:i.1000x)^10 NB. checking the result 91409924241424243424241924242500 Suppose we'd like to accomplish this step of building the summation polynomial in code. For example, if we wanted the polynomial for the sum of k powers of n1 consecutive integers, we could say pk =: k sumpoly and then evaluate pk, already a polynomial. This is a case of defining a nounmodified verb (a kind of adverb) to return a tacit verb. sumpoly =: 1 : '(0,}.(% >: x:m.) * bcoeffs >:m.) & p.' p4 =: 4 sumpoly NB. returns a tacit verb p4 101 NB. apply returned verb to n=100 Although the last expression above makes it looks like sumpoly is just a dyadic verb, it's really not, as 4 sumpoly wouldn't be a legal expression if it were. The purpose of sumpoly is to bind the appropriate coefficients as the left argument to the polynomial verb. Remember that & means "bind" or "curry." The m. is optionally used in place of x. in the definition of sumpoly in order to specify that sumpoly expects a noun for a left argument  an exception will be raised if something different is passed. The arguments m. and n. stand for left and right nouns, while u. and v. stand for left and right verbs. x. and y. stand for left and right arguments of whatever type. VI. TAYLOR EXPANSIONS J knows how to generate the Taylor expansion of many primitive functions. The Taylor expansion of f(x) is a polynomial of increasing powers of x with rational coefficients. Deriving this polynomial involves taking successively higher derivatives of f(x) around f(a). We call it the Maclaurin expansion if a=0. The verb t. returns the coefficients of this polynomial. It's left argument is the function f, and the right argument is the index of the term we want. Using i.n gives us terms 0 through n1. For example: cos =: 2&o. NB. o. is the circular functions verb cf =: cos t. x:i.10 NB. first 10 Taylor coefficients 1 0 _1r2 0 1r24 0 _1r720 0 1r40320 0 cf p. 1p1 NB. use polynomial to get cos(pi) _0.976022 cos 1p1 NB. the discrepancy is significant _1 (cos t. i.20) p. 1p1 NB. 20 terms instead of just 10 _1 cos t. x:20 NB. the 20th coefficient as a rational 1r2432902008176640000 Note the Taylor expansion of e^x and its derivative (recall p.. from Lecture 2): ] e =: (^ t. x:i.10) NB. Taylor expansion of e^x 1 1 1r2 1r6 1r24 1r120 1r720 1r5040 1r40320 1r362880 p.. e NB. given an endless expansion, p.. e equals e 1 1 1r2 1r6 1r24 1r120 1r720 1r5040 1r40320 The graph below shows the cosine curve (red) plotted with the curve generated by the first 10 terms of the Taylor expansion of cosine (green). load 'plot' x =: 0.1*i:50 ] cf =: cos t. x:i.10 1 0 _1r2 0 1r24 0 _1r720 0 1r40320 0 plot x;>(cf p. x);cos x This next graph shows the cosine curve (red) plotted with a Taylor expansion of cosine to 15 terms. Notice the curves coincide over a much broader range (note also the difference in scale between the two graphs). x =: 0.1*i:100 cf2 =: cos t. x:i.15 plot x;>(cf2 p. x);cos x Confirmation that e^(i x) = cos x + i sin x may also be obtained by generating the Taylor expansions of both sides: f =: ^@(0j1&*) NB. f(x) = e^(i x) g =: cos + (0j1&*@sin) NB. g(x) = cos(x) + i sin(x) f t. i.7 NB. Taylor expansion (to 7 terms) of f(x) 1 0j1 _0.5 0j_0.166667 0.0416667 0j0.00833333 _0.00138889 g t. i.7 NB. ditto for g(x) 1 0j1 _0.5 0j_0.166667 0.0416667 0j0.00833333 _0.00138889 The Bernoulli numbers may be generated from the Taylor expansion of the
function Using 1: allows us to think of the verb train ^  1: as a fork, instead of writing <:@^ (i.e. decrement exp(x) by 1). Then the initial dyadic % "hooks" x, to make it the numerator as well. Let's look at the Taylor expansion of this function: (% ^  1:) t. x:i.11 NB. Taylor expansion of x/(e^x1) What the above lines demonstrate is the the coefficients of the Taylor
expansion of This pattern of k! occuring in the denominator of a Taylor expansion is part of its definition, and the verb t: (note colon), produces what's called the "weighted Taylor" expansion, meaning successive terms 0...k are automatically multiplied by k!. So t: gives us the Bernoulli numbers directly. . bernoullis 14 NB. Bernoullis 0 through 14 1 _1r2 1r6 0 _1r30 0 1r42 0 _1r30 0 5r66 0 _691r2730 0 7r6 (% ^  1:) t: x:i.15 NB. ditto 1 _1r2 1r6 0 _1r30 0 1r42 0 _1r30 0 5r66 0 _691r2730 0 7r6 Given we have direct access to the Bernoullis in this weighted Taylor expansion, it might seem that our earlier algorithm based on binomial coefficients is redundant. However, you'll discover that finding the 100th Bernoulli number using the Taylor expansion takes significantly longer than the binomialsbased approach. However, J does come with a module (./system/packages/math/numbers.ijs) which includes a Bernoulli numbers algorithm that's faster yet, based in turn on socalled tangent numbers.
