Springer Undergraduate Mathematics Series
Advisory Board M.A.J. Chaplain University of Dundee K. Erdmann University of...

Author:
Roozbeh Hazrat

This content was uploaded by our users and we assume good faith they have the permission to share this book. If you own the copyright to this book and it is wrongfully on our website, we offer a simple DMCA procedure to remove your content from our site. Start by pressing the button below!

Springer Undergraduate Mathematics Series

Advisory Board M.A.J. Chaplain University of Dundee K. Erdmann University of Oxford A. MacIntyre Queen Mary, University of London E. S¨uli University of Oxford J.F. Toland University of Bath

For other titles published in this series, go to www.springer.com/series/3423

Roozbeh Hazrat

Mathematica®: A Problem-Centered Approach

Roozbeh Hazrat Department of Pure Mathematics Queen’s University Belfast BT7 1NN United Kingdom [email protected]

ISSN 1615-2085 ISBN 978-1-84996-250-6 e-ISBN 978-1-84996-251-3 DOI 10.1007/978-1-84996-251-3 Springer London Dordrecht Heidelberg New York British Library Cataloguing in Publication Data A catalogue record for this book is available from the British Library Library of Congress Control Number: 2010929694 Mathematics Subject Classification (2000): 68-01, 68N15 © Springer-Verlag London Limited 2010 Apart from any fair dealing for the purposes of research or private study, or criticism or review, as permitted under the Copyright, Designs and Patents Act 1988, this publication may only be reproduced, stored or transmitted, in any form or by any means, with the prior permission in writing of the publishers, or in the case of reprographic reproduction in accordance with the terms of licenses issued by the Copyright Licensing Agency. Enquiries concerning reproduction outside those terms should be sent to the publishers. The use of registered names, trademarks, etc., in this publication does not imply, even in the absence of a specific statement, that such names are exempt from the relevant laws and regulations and therefore free for general use. The publisher makes no representation, express or implied, with regard to the accuracy of the information contained in this book and cannot accept any legal responsibility or liability for any errors or omissions that may be made. Mathematica and the Mathematica logo are registered trademarks of Wolfram Research, Inc (“WRI” – www.wolfram.com) and are used herein with WRI’s permission. WRI did not participate in the creation of this work beyond the inclusion of the accompanying software, and it offers no endorsement beyond the inclusion of the accompanying software Cover design: deblik Printed on acid-free paper Springer is part of Springer Science+Business Media (www.springer.com)

Preface

Teaching the mechanical performance of routine mathematical operations and nothing else is well under the level of the cookbook because kitchen recipes do leave something to the imagination and judgment of the cook but the mathematical recipes do not. G. P´ olya

This book grew out of a course I gave at Queen’s University Belfast during the period of 2004 to 2009. Although there are many books already written about how to use Wolfram Mathematica® , I noticed they fall into two categories: either they provide an explanation about the commands, in the style of: enter the command, push the button and see the result; or they study some problems and write several-paragraph codes in Mathematica. The books in the ﬁrst category did not inspire me (or my imagination) and the second category were too diﬃcult to understand and not suitable for learning (or teaching) Mathematica’s abilities to do programming and solve problems. I could not ﬁnd a book that I could follow to teach this module. In class one cannot go on forever showing students just how commands in Mathematica work; on the other hand it would be very diﬃcult to follow the codes if one writes a program having more than ﬁve lines (especially as Mathematica’s style of programming provides a condensed code). Thus this book. This book promotes Mathematica’s style of programming. I tried to show when we adopt this approach, how naturally one can solve (nice) problems with (Mathematica) style. Here is an example: Does the Euler formula n2 + n + 41 produce prime numbers for n = 1 to 39?

Or in another problem we tried to show how one can eﬀectively use pattern

viii

Preface

matching to check that for two matrices A and B, (ABA−1 )5 = AB 5 A−1 . One only needs to introduce the fact that AA−1 = 1 and then Mathematica will check the problem by cancelling the inverse elements instead of direct calculation. Although the meaning of the code above may not be clear yet, the reader will observe as we proceed how the codes start making sense, as if this is the most natural way to approach the problems. (People who approach problems with a procedural style of programming (such as C++) will experience that this style replaces their way of thinking!) We have tried to let the reader learn from the codes and avoid long and exhausting explanations, as the codes will speak for themselves. Also we have tried to show that in Mathematica (as in the real world) there are many ways to approach a problem and solve it. We have tried to inspire the imagination! Someone once rightly said that the Mathematica programming language is rather like a “Swiss army knife” containing a vast array of features. Mathematica provides us with powerful mathematical functions. Along with this, one can freely mix diﬀerent styles of programming, functional, list-based and procedural, to achieve a lot. This m´elange of programming styles is what we promote in this note. Thus this book could be considered for a course in Mathematica, or for self study. It mainly concentrates on programming and problem solving in Mathematica. I have mostly chosen problems having something to do with numbers as they do not need any particular background. Many of these problems were taken from or inspired by those collected in [3]. I would like to thank Ilan Vardi for answering my emails and Brian McMaster and Judith Millar for polishing the English of this note. Naoko Morita encouraged me to make my notes into this book. I thank her for this and for always smiling and having a little Geschichte zu erz¨ ahlen. Roozbeh Hazrat [email protected] Belfast, October 2009

How to use this book

Each chapter of the book starts with a description of a new topic and some basic examples. We will then demonstrate the use of new commands in several problems and their solutions. We have tried to let the reader learn from the codes and avoided long and exhausting explanations, as the codes will speak for themselves. There are three diﬀerent categories of problems, shown by diﬀerent frames:

Problem 0.1 These problems are the essential parts of the text where new commands are introduced to solve the problem and where it is demonstrated how the commands are used in Wolfram Mathematica® . These problems should not be skipped. =⇒ Solution.

Problem 0.2 These problems further demonstrate how one can use commands already introduced to tackle diﬀerent situations. The readers are encouraged to try out these problems by themselves ﬁrst and then look at the solution. =⇒ Solution.

x

How to use this book

Problem 0.3 These are more challenging problems that could be skipped in the ﬁrst reading. =⇒ Solution.

Most commands in Mathematica have several extensions. When we introduce a command, this is just the tip of the iceberg. In TIPS, we give further indications about the commands, or new commands to explore. Once one has enough competence, then it is quite easy to learn further abilities of a command by using the Mathematica Help and the examples available there.1

1

Included with this book is a free 30 day trial of the Mathematica software. To access your free download, simply go to http://www.wolfram.com/books/resources and enter license number L3280-8445. You will be guided to download and install the latest version of Mathematica.

The Mathematica philosophy

In the beginning is the expression. Wolfram Mathematica® transforms the expression dictated by the rules to another expression. And that’s how a new idea comes into the world! The rules that will be used frequently can be given a name (we call them functions) r[x_]:=1+x^2 r[arrow] 1+arrow^2 r[5] 26

And the transformation could take place immediately or with a delay {x,x}/.x->Random[] {0.0474307, 0.0474307} {x,x}/.x:>Random[] {0.368461, 0.588353}

The most powerful transformations are those which look for a certain pattern in an expression and morph that to a new pattern. (a + b)^c /. (x_ + y_)^z_ -> (x^z + y^z) a^c + b^c

And one of the most important expressions that we will work with is a list. As the name implies this is just a collection of elements (collection of other expressions). We then apply transformations to each element of the list: x^ {1, 5, 7} {x, x^5, x^7}

The Mathematica philosophy

xii

Any expression is considered as having a “head” followed by several arguments, head[arg1,arg2,...]. And one of the transformations which provide a capable tool is to replace the head of an expression by another head! Plus @@ {a,b,c} a+b+c Power @@ (x+y) x^y

Putting these concepts together creates a powerful way to solve problems. In the chapters of this book, we decipher these concepts.

Contents

1.

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Mathematica as a calculator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Numbers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Algebraic computations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Trigonometric computations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Equalities =, :=, == . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 Dynamic variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1 1 5 9 11 12 13 15

2.

Deﬁning functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.1 Formulas as functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.2 Anonymous functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

3.

Lists . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Functions producing lists . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Listable functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Selecting from a list . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.

Changing heads! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

5.

A bit of logic and set theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Being logical . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Handling sets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Decision making, If and Which . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6.

Sums and products . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 6.1 Sum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

25 28 32 34

54 54 58 61

xiv

Contents

6.2 Product . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 7.

Loops and repetitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1 Do, For a While . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Nested loops . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3 Nest, NestList and more . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.4 Fold and FoldList . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.5 Inner and Outer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

72 72 81 84 90 93

8.

Substitution, Mathematica rules . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

9.

Pattern matching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

10. Functions with multiple deﬁnitions . . . . . . . . . . . . . . . . . . . . . . . . . 112 10.1 Functions with local variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 10.2 Functions with conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 11. Recursive functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 12. Linear algebra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 12.1 Vectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 12.2 Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 13. Graphics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 13.1 Two-dimensional graphs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 13.2 Three-dimensional graphs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 14. Calculus and equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158 14.1 Solving equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158 14.2 Calculus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163 15. Solutions to the Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169 Further reading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185 Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 186

1 Introduction

This chapter introduces basic capabilities of Mathematica, which include simple arithmetic, handling algebraic and trigonometric expressions and assigning values to variables. We will also look at dynamic objects, allowing us to see changes in the variables as they happen.

In this chapter we give a quick introduction to the very basic things one can perform with Wolfram Mathematica® . We let the reader learn from reading the codes and avoid long and exhausting explanations, as the codes will speak for themselves.

1.1 Mathematica as a calculator Mathematica can be used as a powerful calculator with the basic arithmetic operations; +, − for addition and subtraction, ∗, / for multiplication and division and ˆ for powers. 10^9 - 987654321 12345679 2682440^4 + 15365639^4 + 18796760^4 180630077292169281088848499041 20615673^4 180630077292169281088848499041

The last two calculations show that 26824404 + 153656394 + 187967604 = 206156734 ,

R. Hazrat, Mathematica ®: A Problem-Centered Approach, Springer Undergraduate Mathematics Series, DOI 10.1007/978-1-84996-251-3 1, © Springer-Verlag London Limited 2010

2

1. Introduction

disproving a conjecture by Euler that three fourth powers can never sum to a fourth power.1 Mathematica can handle large calculations: 2^9941-1 346088282490851215242960395767413316722628668900238547790489283445006220809834 114464364375544153707533664486747635050186414707093323739706083766904042292657 896479937097603584695523190454849100503041498098185402835071596835622329419680 597622813345447397208492609048551927706260549117935903890607959811638387214329 942787636330953774381948448664711249676857988881722120330008214696844649561469 971941269212843362064633138595375772004624420290646813260875582574884704893842 439892702368849786430630930044229396033700105465953863020090730439444822025590 974067005973305707995078329631309387398850801984162586351945229130425629366798 595874957210311737477964188950607019417175060019371524300323636319342657985162 360474512090898647074307803622983070381934454864937566479918042587755749738339 033157350828910293923593527586171850199425548346718610745487724398807296062449 119400666801128238240958164582617618617466040348020564668231437182554927847793 809917495802552633233265364577438941508489539699028185300578708762293298033382 857354192282590221696026655322108347896020516865460114667379813060562474800550 717182503337375022673073441785129507385943306843408026982289639865627325971753 720872956490728302897497713583308679515087108592167432185229188116706374484964 985490944305412774440794079895398574694527721321665808857543604774088429133272 929486968974961416149197398454328358943244736013876096437505146992150326837445 270717186840918321709483693962800611845937461435890688111902531018735953191561 073191960711505984880700270887058427496052030631941911669221061761576093672419 481606259890321279847480810753243826320939137964446657006013912783603230022674 342951943256072806612601193787194051514975551875492521342643946459638539649133 096977765333294018221580031828892780723686021289827103066181151189641318936578 454002968600124203913769646701839835949541124845655973124607377987770920717067 108245037074572201550158995917662449577680068024829766739203929954101642247764 456712221498036579277084129255555428170455724308463899881299605192273139872912 009020608820607337620758922994736664058974270358117868798756943150786544200556 034696253093996539559323104664300391464658054529650140400194238975526755347682 486246319514314931881709059725887801118502811905590736777711874328140886786742 863021082751492584771012964518336519797173751709005056736459646963553313698192 960002673895832892991267383457269803259989559975011766642010428885460856994464 428341952329487874884105957501974387863531192042108558046924605825338329677719 469114599019213249849688100211899682849413315731640563047254808689218234425381 995903838524127868408334796114199701017929783556536507553291382986542462253468 272075036067407459569581273837487178259185274731649705820951813129055192427102 805730231455547936284990105092960558497123779789849218399970374158976741548307 086291454847245367245726224501314799926816843104644494390222505048592508347618 947888895525278984009881962000148685756402331365091456281271913548582750839078 91469979019426224883789463551

If a number of the form 2n − 1 happens to be prime, it is called a Mersenne prime. Recall that a prime number is a number which is divisible only by 1 and itself. It is easy to see 22 − 1 and 23 − 1 and 25 − 1 are Mersenne primes. The list continues. In 1963, Gillies found that the above number, 29941 − 1, is a Mersenne prime. With my laptop it takes about 3 seconds for Mathematica to check that this is a prime number.2 PrimeQ[2^9941-1] True

Back to easier calculations: 24/17 24/17

24 Mathematica always tries to give a precise value, thus gives back instead 17 of attempting to evaluate the fraction. 1

2

This conjecture remained open for almost 200 years, until Noam Elkies at Harvard came up with the above counterexample in 1988. The largest Mersenne prime found so far is 243,112,609 − 1 which was discovered in August 2008 at UCLA and has 12, 978, 189 digits.

1.1 Mathematica as a calculator

3

Sin[Pi/5] √ 1 1 (5 − 5) 2 2 √ Mathematica gives 12 12 (5 − 5) as the value of sin(π/5), which is the precise equality. This shows Mathematica is not approaching the expressions numerically. In order to get the numerical value, one can use the function N. N[24/17] 1.41176 N[24/17, 20] 1.4117647058823529412 ?N N[expr] gives the numerical value of expr. N[expr, n] attempts to give a result with n-digit precision.

As the above line shows, in order to get a quick description of a command, use ?Command. If you need more explanation use ??Command. If you need even more, use Mathematica’s Document Center in the Help Menu and search for the command. Mathematica comes with an excellent Help which contains explanations and many nice examples demonstrating how to use each of the functions available in this software. All elementary mathematical functions are available here, Log, Exp, Sqrt, Sin, Cos, Tan, ArcSin,.... Here we evaluate π ( )2 + (0.5 Log[2])2 , 4 Sqrt[(Pi/4)^2 + (0.5 Log[2])^2] 0.858466

Problem 1.1 Use Mathematica to show that tan

3π 2π √ + 4 sin = 11. 11 11

=⇒ Solution. Tan[3 Pi/11] + 4 Sin[2 Pi/11] Cot[(5 Pi)/22] + 4 Sin[(2 Pi)/11]

We didn’t get

√

11 as an answer. We ask Mathematica to try a bit harder.

4

1. Introduction

?Simplify Simplify[expr] performs a sequence of algebraic and other transformations on expr, and returns the simplest form it finds. Simplify[expr,assum] does simplification using assumptions. >> Simplify[Tan[3 Pi/11] + 4 Sin[2 Pi/11]] Cot[(5 Pi)/22] + 4 Sin[(2 Pi)/11] ?FullSimplify FullSimplify[expr] tries a wide range of transformations on expr involving elementary and special functions, and returns the simplest form it finds. FullSimplify[expr,assum] does simplification using assumptions. >> FullSimplify[Tan[3 Pi/11] + 4 Sin[2 Pi/11]] Sqrt[11]

This problem introduces the commands Simplify and FullSimplify. If one is not happy with the result, one can always use Simplify and even FullSimplify to have Mathematica work a bit harder to come up with more simpliﬁcation.

For a complete list of elementary functions have a look at Functional Navigator: Mathematics and Algorithms: Mathematical Functions in Mathematica’s Help.

Exercise 1.1 Show that

√ 3 64 22 + (1/2)2 − 1 = 4.

Problem 1.2 Using Mathematica, explain why 4 + 6/4 ∗ 3ˆ− 2 + 1 =

31 . 6

=⇒ Solution. If you look at Mathematica’s Help, under “Special Ways to Input Expressions”, you will see the following note: “The Mathematica language has a definite grammar which speciﬁes how your input should be converted to internal form.” One aspect of the grammar is that it speciﬁes how pieces of your input should be grouped. The general rule is that if ⊗ has higher precedence than ⊕, then a⊕b⊗c is interpreted as a⊕(b⊗c), and a⊗b⊕c is interpreted as (a⊗b)⊕c. You will then ﬁnd a long table listing which operation has a higher precedence

1.2 Numbers

5

and thus, based on that, you will be able to explain why 4 + 6/4 ∗ 3ˆ− 2 + 1 . amounts to 31 6 However, common sense tells us that instead of creating an ambiguous expression such as 4 + 6/4 ∗ 3ˆ− 2 + 1, one should use parentheses () to group objects together and make the expression more clear. For example, one could write 4 + (6/4) ∗ 3ˆ(−2) + 1, or even better use Mathematica’s Palette (Basic Math Assistance) and type 4+

♣

6 × 3−2 + 1. 4

TIPS

– Comments can be added to the codes using (* comment *). (* the most beautiful theorem in Mathematics *) E^(I Pi) + 1 0

– The symbol % refers to the previous output produced. %% refers to the second previous output and so on. – If in calculations you don’t get what you are expecting, use Simplify or even FullSimplify (see Problem 1.1). – To get the numerical approximation, use N[expr] or alternatively, expr//N (see Problem 2.2 for diﬀerent ways of applying a function to a variable). Use EngineeringForm[expr,n] and ScientificForm[expr,n] to get other forms of numerical approximations to n signiﬁcant digits. – The mathematical constant e, the exponential number, is deﬁned in Mathematica as E, or Exp. To get en use either E^ n or Exp[n].

1.2 Numbers There are several standard ways to start with an integer number and produce new numbers out of it. For example, starting from 4, one can form 4 × 3 × 2 × 1, which is represented by 4!.

6

1. Introduction

4! 24 123! 1214630436702532967576624324188129585545421708848338231532891 8161829235892362167668831156960612640202170735835221294047782 59109157041165147218602951990626164673073390741981495296000 0000000000000000000000000

The fundamental theorem of arithmetic states that one can decompose any number n as a product of powers of primes and this decomposition is unique, i.e., n = pk11 · · · pkt t where pi ’s are prime. Thus 12 = 22 × 31 and 37534 = 2 × 72 × 383. Mathematica can do all of these: FactorInteger[12] {{2, 2}, {3, 1}} FactorInteger[37534] {{2,1},{7,2},{383,1}} 2^1 * 7^2 * 383^1 37534 FactorInteger[6473434456376432] {{2,4},{3239053,1},{124909859,1}} PrimeQ[124909859] True Prime[8] 19

Prime[n] produces the n-th prime number. PrimeQ[n] determines whether n is a prime number. More than 2200 years ago Euclid proved that the set of prime numbers is inﬁnite. His proof is used even today in modern books. However, it is not that long ago that we also learned that there is no simple formula that produces only prime numbers. n In 1640 Fermat conjectured that the formula 22 +1 always produces a prime number. Almost a hundred years later the ﬁrst counterexample was found. PrimeQ[2^(2^1)+1] True PrimeQ[2^(2^2)+1] True PrimeQ[2^(2^3)+1] True

1.2 Numbers

7

PrimeQ[2^(2^4)+1] True PrimeQ[2^(2^5)+1] False 2^(2^5)+1 4294967297 FactorInteger[2^(2^5)+1] {{641,1},{6700417,1}} 5

This shows that 22 + 1 is not a prime number. In fact it decomposes into 5 two prime numbers 22 + 1 = 641 × 7600417.

Problem 1.3 What is the probability that a randomly chosen 12-digit number will be a prime? =⇒ Solution. The probability is the number of 12-digit prime numbers over the number of all 12-digit numbers. So we start with ﬁnding how many 12-digit numbers exist: 10^13 - 10^12 9000000000000

Next, we will ﬁnd how many 12-digit prime numbers exist. We will use the following built-in function of Mathematica. ?PrimePi PrimePi[x] gives the number of primes less than or equal to x. >> PrimePi[10^13] 346065536839 PrimePi[10^12] 37607912018 N[(346065536839 - 37607912018)/9000000000000]*100 3.42731

So the probability that we randomly pick a 12-digit prime number is only 3.42 percent.

8

1. Introduction

Exercise 1.2 Check that 123456789098765432111 is prime.

Exercise 1.3 Check that for any given positive integer n ≥ 3, the least common multiple of the numbers 1, 2, · · · , n is greater than 2n−1 . (Hint: see LCM.)

Exercise 1.4 Check that the number 32! ends with 7 zeros.

Problem 1.4 The function Mod[m,n] gives the remainder on division of m by n. Let a = 98 and b = 75. Show that aa bb ends with exactly 98 zeros. =⇒ Solution. Clearly you can directly compute aa bb and count the number of zeros as in Exercise 1.4. But this time you have to count whether the number ends with 98 zeros. Here is a more clever way: Check whether the remainder of this number when divided by 1098 is zero, but its remainder on division by 1099 is not. Mod[98^98 75^75, 10^98] 0 Mod[98^98 75^75, 10^99] 5000000000000000000000000000000000000000000000000000000000000 00000000000000000000000000000000000000

Problem 1.5 Recall that for integers m and n, the binomial

n! n is deﬁned as . m m!(n − m)!

Using Mathematica, check that (m + n)! m+n . = m!n! m =⇒ Solution. Not only is the binomial function available in Mathematica, but also Mathematica can perfectly handle it symbolically as the solution to this problem shows. We will talk more about symbolic computations in Section 1.3.

1.3 Algebraic computations

9

Binomial[m + n, n] == (n + m)!/(n! m!) Binomial[m + n, n] == (m + n)!/(m! n!) FullSimplify[Binomial[m + n, n] == (n + m)!/(n! m!)] True

This is another instance where we need to use FullSimplify to make Mathematica work harder to come up with the result. We will discuss the diﬀerent equalities available in Mathematica in Section 1.6. However, for the time being, note that == is used to compare both sides of equations.

There are several more integer functions available in Mathematica, which can be found in Functional Navigator: Mathematics and Algorithms: Mathematical Functions: Integer Functions. ♣

TIPS

– The command NextPrime[n] gives the next prime larger than n and PrimePi[n] gives the number of primes less than or equal to n (see Problem 1.3). – For integers m and n, one can ﬁnd unique numbers q and r such that r is positive, m = qn + r and r > Divisible[2^2008, b[23]] True

Later, in Problem 3.4, we will determine all the positive integers n between 3 and 50 for which 22008 is divisible by b(n) and will see that in fact there are very few n with this property.

Later, in Chapter 10, we will deﬁne functions with conditions, functions with several deﬁnitions and functions containing several lines of code (a procedure).

2.2 Anonymous functions Sometimes we need to “deﬁne a function as we go” and use it on the spot. Mathematica enables us to deﬁne a function without giving it a name (nor any reference to any speciﬁc variables) use it, and then move on! These functions are called anonymous or pure functions. Obviously if we need to use a speciﬁc function frequently, then the best way is to give it a name and deﬁne it as we did in Section 2.1. Here is an anonymous function equivalent to f (x) = x2 + 4: (#^2+4)&

24

2. Deﬁning functions

The expression (#^2+4)& deﬁnes a nameless function. As usual we can plug in data in place of #. The symbol & determines where the deﬁnition of the function is completed. (#^2+4)&[5] 29

Compare the following: r[x_] := x (x + 1) r[4] 20 # (# + 1) &[4] 20 4 // # (# + 1) & 20

Exercise 2.2 What do the following pure functions do? Fibonacci[15]//Mod[#,5]& PrimeQ[#! + 1] [email protected] 18//2^#+#&

Anonymous functions can handle several variables. Here is an example of an anonymous function for f (x, y) = x2 + y 2 . Sqrt[#1^2+#2^2]&[3,4] 5

As you might guess, #1 and #2 refer to the ﬁrst and second variables in the function.

3 Lists

A list is a collection of data. In this chapter we study how to handle lists and have access to the elements in the lists. Functions which generate lists are studied and methods of selecting elements from a list with a speciﬁc property are examined. Get ready for some serious programming!

One can think of a computer program as a function which accepts some (crude) data or information and gives back the data we would like to obtain. Lists are the way Wolfram Mathematica® handles information. Roughly speaking, a list is a collection of objects. The objects could be of any type and pattern. Let us start with an example of a list: {1, -4/7, stuff, 1 - 2 x + x^2}

This looks like a mathematical set. One diﬀerence is that lists respect order: {1, 2} == {2, 1} False

The other diﬀerence is that a list can contain a copy of the same object several times: {1,2,1} == {1,2} False

The natural thing here is to be able to access the elements of a list. Let us deﬁne a list p as follows: p = {x, 1, -8/3, a, b, {c, d, e}, radio} {x, 1, -(8/3), a, b, {c, d, e}, radio}

As this example shows, the list p can have any sort of data, even another list as one of its elements. Here is how we can access the elements of the list:

R. Hazrat, Mathematica ®: A Problem-Centered Approach, Springer Undergraduate Mathematics Series, DOI 10.1007/978-1-84996-251-3 3, © Springer-Verlag London Limited 2010

26

3. Lists

p[[1]] x p[[5]] b p[[-1]] radio p[[{2, 4}]] {1, a} p[[{-2, 5}]] {{c, d, e}, b} p[[-2, {2, 3}]] {d, e}

Examining the above examples reveals that p[[i]] gives the i-th member of the list. There are other commands to access the elements of a list as follows: p = {x, 1, -8/3, a, b, {c, d, e}, radio} {x, 1, -(8/3), a, b, {c, d, e}, radio} First[p] x Last[p] radio Drop[p, 3] {a, b, {c, d, e}, radio} Take[p, 2] {x, 1} Rest[p] {1, -(8/3), a, b, {c, d, e}, radio} Rest[%] {-(8/3), a, b, {c, d, e}, radio} Rest[%] {a, b, {c, d, e}, radio} Most[p] {x, 1, -(8/3), a, b, {c, d, e}} Most[Rest[p]] == Rest[Most[p]] True

Most of these commands are self-explanatory and a close look at the above examples shows what each of them will do (see the Tips on page 5 for %). All these commands and more are listed in the Mathematica Help under Functional Navigator: Core Language: List Manipulation: Elements of Lists.

3. Lists

27

Problem 3.1 Let p={a,b,{c,d},e}. From this list produce the list {a,b,c,d,e}. =⇒ Solution. {p[[1]], p[[2]], p[[3, 1]], p[[3, 2]], p[[4]]} {a, b, c, d, e}

One can also use the command Flatten to get rid of extra { and }. We will discuss this command in Section 7.2, where we will be looking at nested loops which create nested lists. p = {a, b, c, {d, e}, f} {a, b, c, {d, e}, f} Flatten[p] {a, b, c, d, e, f}

One of the secrets of writing codes comfortably is that one should be able to manipulate lists easily. Often in applications, situations like the following arise: – Given {x1 , x2 , · · · , xn } and {y1 , y2 , · · · , yn }, produce {x1 , y1 , x2 , y2 , · · · , xn , yn }. – Given {x1 , x2 , · · · , xn } and {y1 , y2 , · · · , yn }, produce {x1 + y1 , x2 + y2 , · · · , xn + yn }. – Given {x1 , x2 , · · · , xn } and {y1 , y2 , · · · , yn }, produce

{x1 , y1 }, {x1 , y2 }, {x1 , y3 }, · · · , {x1 , yn }

{x2 , y1 }, {x2 , y2 }, · · · , {xn , y1 }, {xn , y2 }, · · · , {xn , yn } .

– Given {x1 , x2 , · · · , xn }, produce {x1 , x1 + x2 , · · · , x1 + x2 + · · · + xn }. – Given {x1 , x2 , · · · , xn }, produce {x1 }, {x2 , . . . , xn } , {x1 , x2 }, {x3 , . . . , xn } . . . {x1 , . . . xn−1 }, {xn } .

28

3. Lists

Mathematica provides us with commands to obtain the above arrangements easily. We will look at these commands in Sections 7.4, 7.5, 12.1 and Problem 9.4. ♣

TIPS

– One can add elements to a list. There are several commands to handle this, including Append, Prepend and Insert. ?Insert Insert[list,elem,n] insert elem at position n in list. If n is negative, the position is counted from the end. Insert[{i, think, i, am}, "therefore", 3] {i, think, therefore, i, am}

See also Example 5.3 for the use of Append. – Commands Sort, Reverse, RotateLeft and RotateRight are available to rearrange the order of a list. – Commands Delete and Drop are available to remove elements from a list.

3.1 Functions producing lists Mathematica provides us with commands of which the output is a list. These commands have a nature of repetition and replace loops in procedural programming (more on this in Chapter 7). Let us look at some of them here before starting to write more serious codes. Range[10] {1, 2, 3, 4, 5, 6, 7, 8, 9, 10} Range[3, 11] {3, 4, 5, 6, 7, 8, 9, 10, 11} Range[2, 17, 4] {2, 6, 10, 14} ?Range Range[imax] generates the list {1,2,...,imax}. Range[imin, imax] generates the list {imin,...,imax}. Range[imin,imax,di] uses step di.>>

One of the most useful commands which produces a list is Table.

3.1 Functions producing lists

29

Table[2 n + 1, {n, 1, 13}] {3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27}

In Table[2n+1,{n,1,13}], n runs from 1 to 13 and each time the function 2n+1 is evaluated. The following example shows how easily we can work symbolically in Mathematica. Table[x^i + y^i, {i, 2, 17, 4}] {x^2 + y^2, x^6 + y^6, x^10 +y^10, x^14 + y^14}

As in the last example of Range, here in Table, i starts from 2 with steps 4 and thus takes the values 2,6,10,14. Here is one more example demonstrating how beautifully Mathematica can handle symbols Table[xi ,{i,1,10}] {x1 , x2 , x3 , x4 , x5 , x6 , x7 , x8 , x9 , x10 }

Problem 3.2 Produce the list of the ﬁrst 30 Fibonacci numbers. =⇒ Solution. All we have to do is to let i run from 1 to 30 and each time plug i into the function Fibonacci[i] which produces the i-th Fibonacci number. This can be done with Table as follows: Table[Fibonacci[i], {i, 1, 30}] {1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025, 121393, 196418, 317811, 514229, 832040}

Problem 3.3 Find all natural numbers n between 1 and 15 for which the polynomial xn + 64 can be written as a product of two nonconstant polynomials with integer coeﬃcients. =⇒ Solution. We have seen that Factor will give a factorization of an expression (if it is possible).

30

3. Lists

Table[Factor[x^i + 64], {i, 1, 10}] {64 + x, 64 + x^2, (4 + x) (16 - 4 x + x^2), (8 - 4 x + x^2) (8 + 4 x + x^2), 64 + x^5, (4 + x^2) (16 - 4 x^2 + x^4), 64 + x^7, (8 - 4 x^2 + x^4) (8 + 4 x^2 + x^4), (4 + x^3) (16 - 4 x^3 + x^6), 64 + x^10}

Although this solves the problem (for 1 ≤ n ≤ 10), the output is not arranged nicely and it is hard to understand the answer. It would be easier if we keep track of i as well. We change the code slightly. For this we need Print ?Print Print[expr] prints expr as output. >>

Here is the improved code: Table[Print[i, "

", Factor[x^i + 64]], {i, 1, 15}];

1

64+x

2

64+x^2

3

(4+x) (16-4 x+x^2)

4

(8-4 x+x^2) (8+4 x+x^2)

5

64+x^5

6

(4+x^2) (16-4 x^2+x^4)

7

64+x^7

8

(8-4 x^2+x^4) (8+4 x^2+x^4)

9

(4+x^3) (16-4 x^3+x^6)

10

64+x^10

11

64+x^11

12 (2-2 x+x^2) (2+2 x+x^2) (4-4 x+2 x^2-2 x^3+x^4) (4+4 x+2 x^2+2 x^3+x^4) 13

64+x^13

14

64+x^14

15

(4+x^5) (16-4 x^5+x^10)

From the list, for n = 3, 4, 6, 8, 9, 12, 15 we have a factorization into two nonconstant polynomials (in fact, for n = 3k and 4k the polynomial can be factorized into two polynomials).

3.1 Functions producing lists

31

Problem 3.4 Determine all the positive integers n between 3 and 50 for which 22008 is divisible by n n n 1+ + + . 1 2 3 =⇒ Solution. Recall that if Mod[m,n] returns zero, then m is divisible by n. We ﬁrst deﬁne a function b(n) = 1 + n1 + n2 + n3 (see Problem 2.3) and then use a Table to check when 22008 is divisible by b(n) for diﬀerent n. b[n_] := Binomial[n, 1] + Binomial[n, 2] + Binomial[n, 3] + 1 Table[Mod[2^2008, b[n]], {n, 3, 50}] {0, 1, 16, 16, 0, 70, 16, 80, 168, 133, 268, 316, 448, 256, 706, 796, 1096, 723, 1092, 1030, 0, 256, 458, 1240, 2704, 2604, 606, 2922, 640, 1664, 2704, 4076, 2824, 1936, 1024, 8336, 256, 7882, 6974, 4192, 4568, 6061, 1076, 6896, 704, 16669, 6856, 16032}

It is clear from the list that only for n = 3, 7 and 23, is 22008 divisible by the above expression (in fact, one can show that these are the only numbers with this property).

Here is a nice example showing the diﬀerence between the two equalities = and :=.

Example 3.5 This example uses BarChart which is a graphic function. The example is based on the discussion in Section 1.6. ?BarChart BarChart[{y1,y2,...}] makes a bar chart with bar lengths y1,y2,... x=RandomInteger[{1,1000}]; BarChart[Table[x,{200}]] x:=RandomInteger[{1,1000}] BarChart[Table[x,{200}]]

In order to understand this example better, get the list generated by Table[x,1000] for each of the deﬁnitions of x individually and compare them.

32

3. Lists

Figure 3.1 Using = as equality

Figure 3.2 Using := as equality

3.2 Listable functions There are times when we would like to apply a function to all the elements of a list. Suppose f is a function and {a,b,c} is a list. We want to be able to “push” the function f inside the list and get {f[a],f[b],f[c]}. Many of Mathematica’s built-in functions have the property that they simply “go inside” a list. This property of a function is called listable. For example Sqrt is a listable function. Sqrt[{a, b, c, d, e}]

√ √ √ √ √ { a, b, c, d, e} All the arithmetic functions are listable: 1+ {a, b, c, d, e} {1 + a, 1 + b, 1 + c, 1 + d, 1 + e} {a, b, c, d, e}^3 {a^3,b^3,c^3,d^3,e^3} 1/{a, b, c, d, e} {1/a,1/b,1/c,1/d,1/e}

We will use the function Sqrt in the following, to show that the product of four consecutive numbers plus one is always a squared number. Table[n (n + 1) (n + 2) (n + 3) + 1, {n, 1, 10}] {25, 121, 361, 841, 1681, 3025, 5041, 7921, 11881, 17161} Sqrt[%] {5, 11, 19, 29, 41, 55, 71, 89, 109, 131}

3.2 Listable functions

33

Not all the functions are listable. If we want to push a function into a list, the command to use is Map. f[{a, b, c, d, e}] f[{a, b, c, d, e}] Map[f, {a, b, c, d, e}] {f[a], f[b], f[c], f[d], f[e]}

The equivalent shorthand to apply a function to a list is /@ as follows: f /@ {a, b, c, d, e} {f[a], f[b], f[c], f[d], f[e]}

Here is an example: Table[(1 + x)^i, {i, 5}] {1 + x, (1 + x)^2, (1 + x)^3, (1 + x)^4, (1 + x)^5} Expand /@ % {1 + x, 1 + 2 x + x^2, 1 + 3 x + 3 x^2 + x^3, 1 + 4 x + 6 x^2 + 4 x^3 + x^4, 1 + 5 x + 10 x^2 + 10 x^3 + 5 x^4 + x^5}

Problem 3.6 The formula n2 + n + 41 has a very interesting property. Observe that this formula produces prime numbers for all n between 0 and 39. =⇒ Solution. First we produce the numbers: Table[n^2 + n + 41, {n, 0, 40}] {41, 43, 47, 53, 61, 71, 83, 97, 113, 131, 151, 173, 197, 223, 251, 281, 313, 347, 383, 421, 461, 503, 547, 593, 641, 691, 743, 797, 853, 911, 971, 1033, 1097, 1163, 1231, 1301, 1373, 1447, 1523, 1601, 1681}

Then we apply PrimeQ to this list. This function is listable. PrimeQ[%] {True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, False}

The list shows that for n between 0 and 39, the above formula produces prime numbers, however, for n = 40, the last number in the list 402 + 40 + 41 is not prime. This formula was found by Euler around 1772. One wonders, if one changes 41 to another number in the formula n2 + n + 41, whether one gets more consecutive prime numbers. We will examine this in Problem 7.2. Of course we could also simply write Table[PrimeQ[n^2 + n + 41], {n, 0, 40}]

34

3. Lists

as well. As we said, there are many ways to approach a problem in Mathematica.

Exercise 3.1 Generate 6 random integers n between 1 and 3095 and show that n6 + 1091 is not prime. (Hint: see RandomInteger for generating random numbers.)

3.3 Selecting from a list So far we have been able to produce a list of data by using functions producing lists. The next step is to be able to choose, from a list, certain data which ﬁt a speciﬁc description. This can be achieved by the command Select as the following problem demonstrates.

Problem 3.7 How many numbers of the form 3n5 + 11, when n varies from 1 to 2000, are prime? =⇒ Solution. First, let us produce the ﬁrst 20 numbers of this form. plist=Table[3 n^5 + 11, {n, 1, 20}] {14, 107, 740, 3083, 9386, 23339, 50432, 98315, 177158, 300011, 483164, 746507, 1113890, 1613483, 2278136, 3145739, 4259582, 5668715, 7428308, 9600011}

The next step, in the spirit of Section 3.2, would be to apply PrimeQ to all the numbers and ﬁnd out which ones are prime. Since this is a listable function it is enough: PrimeQ[plist] {False, True, False, True, False, True, False, False, False, False, False, True, False, True, False, True, False, False, False, False}

Now we are left to count the number of True ones. This is do-able here; 6 of these numbers are prime. However, if we are dealing with a list with 2000 elements, we need to select, from the elements of the list, those with a desired property,1 here being prime. The command Select does just this: 1

Or a desired pattern, more about this later.

3.3 Selecting from a list

35

Select[plist,PrimeQ] {107, 3083, 23339, 746507, 1613483, 3145739}

These are prime numbers in the list plist. The command Length gives the length of a list. If we assemble all the steps in one line we have Length[Select[Table[3 n^5 + 11, {n, 1, 20}], PrimeQ]] 6

Thus to ﬁnd out how many numbers of the form 3n5 + 11 are prime when n runs from 1 to 2000, all we have to do is to change 20 to 2000: Length[Select[Table[3 n^5 + 11, {n, 1, 2000}], PrimeQ]] 97

In a nutshell, Select[list,f], will apply the function f (which returns True or False) to all the elements, say x, of the list and return those elements for which f[x] is true. A function which returns either True or False is called a Boolean function. So far we have seen one Boolean function, that is, PrimeQ. (More on Boolean statements in Chapter 5.)

Problem 3.8 Show that among the ﬁrst 500 Fibonacci numbers, 18 of them are prime. =⇒ Solution. We do this step by step. First we generate the ﬁrst 10 Fibonacci numbers. Then among those we select the ones which are prime. Then using Length we get the size of this list, as follows: fiblist=Table[Fibonacci[i], {i, 1, 10}] {1, 1, 2, 3, 5, 8, 13, 21, 34, 55} Select[fiblist, PrimeQ] {2, 3, 5, 13} Length[Select[Table[Fibonacci[i], {i, 1, 10}], PrimeQ]] 4

Now instead of 10 we generate 500 numbers: Length[Select[Table[Fibonacci[i], {i, 1, 500}], PrimeQ]] 18

36

3. Lists

Exercise 3.2 Show that, among the ﬁrst 450 Fibonacci numbers, the number of odd Fibonacci numbers is twice the number of even ones. (Hint: see OddQ and EvenQ for odd and even numbers.)

Exercise 3.3 Let m be a natural number and A=

(m + 3)3 + 1 . 3m

Find all the integers m less than 500 such that A is an integer. Show that A is always odd. (Hint: see IntegerQ.) Let us look at another example, slightly diﬀerent but of the same nature as Problem 3.7. The following example shows that anonymous functions ﬁt very well with Select.

Problem 3.9 For which 1 ≤ n ≤ 1000 does the formula 2n + 1 produce a prime number? =⇒ Solution. Here is the solution: Select[Range[1000], PrimeQ[2^(#) + 1] &] {1, 2, 4, 8, 16}

Let us take a deep breath and go through this one-liner code slowly. The function PrimeQ[2^(#)+1]& is an anonymous function which gives True if the number 2n + 1 is prime and False otherwise. PrimeQ[2^(#) + 1] &[12] False

Range[1000] creates a list containing the numbers from 1 to 1000. The command Select applies the anonymous function PrimeQ[2^(#)+1]& to each element of this list and in the cases when the result is true, i.e., when 2n + 1 is prime, the element will be selected. Thus {1, 2, 4, 8, 16} are the only numbers that make 2n + 1 a prime number.

Exercise 3.4 Find the number of positive integers 0 < n < 20000 such that 1997 divides n2 + (n + 1)2 . Try the same code for 2009.

3.3 Selecting from a list

37

Exercise 3.5 For integers 2 ≤ n ≤ 200, ﬁnd all n such that n divides (n − 1)! + 1. Show that there are 46 such n.

Problem 3.10 Notice that 122 = 144 and 212 = 441, namely the numbers and their squares are reverses of each other. Find all the numbers up to 10000 with this property. =⇒ Solution. We need to introduce some new built-in functions. IntegerDigits[n] gives a list of the decimal digits of the integer n. We also need Reverse and FromDigits: IntegerDigits[80972] {8, 0, 9, 7, 2} Reverse[%] {2, 7, 9, 0, 8} FromDigits[%] 27908

Thus the above shows we can easily produce the reverse of a number: re[n_] := FromDigits[Reverse[IntegerDigits[n]]] re[12345] 54321 re[6548629] 9268456

Having this function under our belt, the solution to the problem is just one line. Notice that the problem is asking for the numbers n such that re[nˆ2]=re[n]ˆ2 . Select[Range[10000], re[#]^2 == re[#^2] &] {1, 2, 3, 10, 11, 12, 13, 20, 21, 22, 30, 31, 100, 101, 102, 103, 110, 111, 112, 113, 120, 121, 122, 130, 200, 201, 202, 210, 211, 212, 220, 221, 300, 301, 310, 311, 1000, 1001, 1002, 1003, 1010, 1011, 1012, 1013, 1020, 1021, 1022, 1030, 1031, 1100, 1101, 1102, 1103, 1110, 1111, 1112, 1113, 1120, 1121, 1122, 1130, 1200, 1201, 1202, 1210, 1211, 1212, 1220, 1300, 1301, 2000, 2001, 2002, 2010, 2011, 2012, 2020, 2021, 2022, 2100, 2101, 2102, 2110, 2111, 2120, 2121, 2200, 2201, 2202, 2210, 2211, 3000, 3001, 3010, 3011, 3100, 3101, 3110, 3111, 10000}

38

3. Lists

Here is one more example using the command FromDigits. We know that 11 is a prime number. One wonders what is the next prime number consisting only of ones. A wild guess: a number with 23 digits all one? All we have to do is to produce this number then, with PrimeQ, test whether this is prime. Here is one way to generate this number. Table[1, {23}] {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1} FromDigits[%] 11111111111111111111111 PrimeQ[%] True

Here is the code to ﬁnd out which numbers of this kind up to 500 digits are prime. Select[Range[500], PrimeQ[FromDigits[Table[1, {#}]]] &] {2, 19, 23, 317}

Problem 3.11 Show that the number of k between 0 and 1000 for which

1000 is odd is a k

power of 2. =⇒ Solution. First, Range[0,1000] creates a list of numbers from 0 to 1000. Using an anonymous function, each time we plug these numbers into Binomial[1000, #] and check right away if this is an odd number by OddQ[Binomial[1000, #]] &. If this is the case, Select will pick up these numbers. The rest is clear from the code below. Select[Range[0, 1000], OddQ[Binomial[1000, #]] &] {0, 8, 32, 40, 64, 72, 96, 104, 128, 136, 160, 168, 192, 200, 224, 232, 256, 264, 288, 296, 320, 328, 352, 360, 384, 392, 416, 424, 448, 456, 480, 488, 512, 520, 544, 552, 576, 584, 608, 616, 640, 648, 672, 680, 704, 712, 736, 744, 768, 776, 800, 808, 832, 840, 864, 872, 896, 904, 928, 936, 960, 968, 992, 1000} Length[Select[Range[0, 1000], OddQ[Binomial[1000, #]] &]] 64 FactorInteger[%] {{2, 6}}

3.3 Selecting from a list

39

Exercise 3.6 Show that among the ﬁrst 200 primes p, the ones such that the remainder when 19p−1 is divided by p2 is 1 are {3, 7, 13, 43, 137}.

Exercise 3.7 An integer dn dn−1 . . . d1 is called prime-palindromic if dn dn−1 . . . d1 and d1 . . . dn−1 dn are both prime (for example 941). Write a code to ﬁnd all primepalindromic numbers up to 5000. Observe that there are 167 such numbers.

Problem 3.12 Find all positive integers 0 < m < 105 such that the fourth power of the number of positive divisors of m equals m. (Hint: see Divisors.) =⇒ Solution. The function Divisors[n] gives all the positive numbers which divide n. This is a list which includes 1 and n also. We need to get the number of these divisors, thus the function Length. We are looking for each number m such that the fourth power of Length[Divisors[m]] is m. The rest is clear from the code. Select[Range[10^5], Length[Divisors[#]]^4 == # &] {1, 625, 6561}

Exercise 3.8 Show that there is only one positive integer n smaller than 1000 such that n! + (n + 1)! is the square of an integer. The idea of sending a function into a list, i.e., applying a function to each element of a list, seems to be a good one. We have already mentioned that the listable built-in functions are able to go inside a list, like PrimeQ or Prime. Have a look at the Attributes of Prime in the following: ??Prime Prime[n] gives the nth prime number. Attributes[Prime] = {Listable, Protected}

Also recall that if a function is not listable, Map or /@ will push the function into the list as the following problem demonstrates:

40

3. Lists

Problem 3.13 What digit does not appear as the last digit of the ﬁrst 20 Fibonacci numbers? =⇒ Solution. This one-liner code collects all the digits which appears as the last digit: Union[Last /@ (IntegerDigits /@ (Fibonacci /@ Range[20]))] {0, 1, 2, 3, 4, 5, 7, 8, 9}

Thus 6 is the only digit which is not present. Let us understand this code. As the command /@ applies a function to all elements of a list, Fibonacci /@ Range[20] produces the ﬁrst 20 Fibonacci numbers. Fibonacci /@ Range[20] {1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765}

Then IntegerDigits would go inside this list and get the digits of each number. IntegerDigits/@ Fibonacci /@ Range[20] {{1}, {1}, {2}, {3}, {5}, {8}, {1, 3}, {2, 1}, {3, 4}, {5, 5}, {8, 9}, {1, 4, 4}, {2, 3, 3}, {3, 7, 7}, {6, 1, 0}, {9, 8, 7}, {1, 5, 9, 7}, {2, 5, 8, 4}, {4, 1, 8, 1}, {6, 7, 6, 5}}

Then the function Last will get the last digits as required. Last /@ IntegerDigits /@ Fibonacci /@ Range[20] {1, 1, 2, 3, 5, 8, 3, 1, 4, 5, 9, 4, 3, 7, 0, 7, 7, 4, 1, 5}

Union will get rid of any repetitions in the list. (More on this command in Section 5.2.) Since Fibonacci and IntegerDigits are listable functions, one can also write the above code as follows: Union[Last /@ IntegerDigits[Fibonacci[Range[20]]]]

If one does not want to use IntegerDigits then one can use the Mod function to get access to the last digit of a number. ?Mod Mod[m, n] gives the remainder on division of m by n. Mod[264,10] 4 ?Quotient Quotient[m, n] gives the integer quotient of m and n. Quotient[264,10] 26 26*10+4 264

3.3 Selecting from a list

41

Thus another way to write the code is as follows. Note that Mod is also a listable function. Union[Mod[Fibonacci[Range[20]],10]] {0, 1, 2, 3, 4, 5, 7, 8, 9}

In Section 5.2 we will see how to use Mathematica to get directly the digit 6 that we are after, namely to handle sets.

Recall that one can decompose any number n as a product of powers of primes and this decomposition is unique, i.e., n = pk11 · · · pkt t where pi ’s are prime. Let us call a number a square free number if, in its decomposition to primes, all the ki ’s are 1. Namely, no power of primes can divide this number. Thus 15 = 3 × 5 is a square free number but 48 = 24 × 3 is not. Recall that Select[list,f] will apply the function f (which returns True or False) to all the elements, say x, of the list and return those elements for which f[x] is true. There is an option in Select which makes it possible to get only the ﬁrst n elements of list that satisfy f, i.e., f returns True, as follows: Select[list,f,n]. This comes in very handy, as in some problems we are only interested in a certain number of data which satisfy f. Also this can be used in circumstances where we want to test the elements until something goes wrong or some desirable element comes up. The following example demonstrates this.

Problem 3.14 Write a function squareFreeQ[n] that returns True if the number n is a square free number, and False otherwise. =⇒ Solution. Here is the code: squareFreeQ[n_]:= Select[Last/@ FactorInteger[n], # != 1&,1] == {} squareFreeQ[2*5*7] True squareFreeQ[2*5*7*7] False

Let us decipher this code. Recall from Section 1.2 that if n = pk11 · · · pkt t is the decomposition of n into its prime factors then FactorInteger[n]

42

3. Lists

{{p1 , k1 }, {p2 , k2 }, · · · , {pt , kt }} Now all we have to do is to go through this list and see if all ki ’s are one. So the ﬁrst step is to apply Last to each list to discard pi ’s and get ki ’s. Last /@ FactorInteger[n]

{k1 , k2 , · · · , kt } Having this list, we shall go through the list one by one and examine if ki ’s are one. The anonymous function #= 1& does exactly this. So Select[{k1 , k2 , · · · , kt },#= 1&] gives the list of ki ’s which are not one. But in our case, looking for square free primes, it is enough if only one ki is not one. Then the number is not square free. Thus we use an option of Select which goes through the list until it ﬁnds an element such that ki is not one. So we need to modify the code to Select[{k1 , k2 , · · · , kt },#= 1&,1]. We are almost done. All we have to do is to see if this list is empty or not (namely is there any ki not equal to one). And for this we compare Select[{k1 , k2 , · · · , kt },#= 1&,1]=={}. We can solve this problem later with a slightly diﬀerent method (see Problem 4.3).

Exercise 3.9 Find the ﬁrst 5 positive integers n such that n6 + 1091 is prime. Show that all these n are between 3500 and 8500.

Exercise 3.10 A number with n digits is called cyclic if multiplication by 1, 2, 3, · · · , n produces the same digits in a diﬀerent order. Find the only 6-digit cyclic number.

Problem 3.15 Find out how many primes bigger than n and smaller than 2n exist, when n goes from 1 to 30. =⇒ Solution. We deﬁne an anonymous function which ﬁnds all the primes bigger than n and smaller than 2n and then gets the size of this list. Once we are done with this, we apply this function to a list of numbers from 1 to 30. Our anonymous function looks like this: Length[Select[Range[# + 1, 2 # - 1], PrimeQ]] &

3.3 Selecting from a list

43

Here, Range[# + 1, 2 # - 1] produces all the numbers between n and 2n. Then Select ﬁnds out which of them are in fact prime. Then we use the command Length to get the number of elements of this list. One can optimize this a bit, as we don’t need to look at the whole range of n to 2n, as clearly even numbers are not prime so we can ignore them right from the beginning. But we leave it to the reader to do this. All we have to do now is to apply this function with Map or /@ to numbers from 1 to 30. Length[Select[Range[# + 1, 2 # - 1], PrimeQ]] & /@ Range[30] {0, 1, 1, 2, 1, 2, 2, 2, 3, 4, 3, 4, 3, 3, 4, 5, 4, 4, 4, 4, 5, 6, 5, 6, 6, 6, 7, 7, 6, 7}

One can also write a more innocent code, using Table as follows: Table[Length[Select[Range[n + 1, 2n - 1], PrimeQ]], {n, 1, 30}] {0, 1, 1, 2, 1, 2, 2, 2, 3, 4, 3, 4, 3, 3, 4, 5, 4, 4, 4, 4, 5, 6, 5, 6, 6, 6, 7, 7, 6, 7}

This doesn’t seem to be the most eﬃcient way to write this problem as each time we test the same numbers repeatedly to see whether they are prime. We will write another code for this problem in Problem 7.1 using a Do Loop. One can also use the built-in function PrimePi to write a very short solution for this problem. For the usage of PrimePi see Problem 1.3.

A word is called palindromic if it reads the same backwards as forwards, e.g., madam. In the next problem we are going to ﬁnd all the palindromic words in a lot of languages! The approach is similar to Problem 3.10, however, in Problem 3.16 we will be working with strings of characters.

Problem 3.16 Write a function to check whether a word is palindromic. Using DictionaryLookup, ﬁnd all the palindromic words in the English language. Then draw a bar chart showing how many palindromic words there are in diﬀerent languages which are supported by Mathematica. =⇒ Solution. We ﬁrst deﬁne a function to check whether a word is palindromic. Using StringReverse, this is easy: pal[n_] := n == StringReverse[n] pal["test"] False pal["kayak"] True

44

3. Lists

Now we select from the English dictionary available in Mathematica all the words which are palindromic. Select[DictionaryLookup[], pal] {"a", "aha", "aka", "bib", "bob", "boob", "bub", "CFC", "civic", "dad", "deed", "deified", "did", "dud", "DVD", "eke", "ere", "eve", "ewe", "eye", "gag", "gig", "huh", "I", "kayak", "kook", "level", "ma’am", "madam", "mam", "MGM", "minim", "mom", "mum", "nan", "non", "noon", "nun", "oho", "pap", "peep", "pep", "pip", "poop", "pop", "pup", "radar", "redder", "refer", "repaper", "reviver", "rotor", "sagas", "sees", "seres", "sexes", "shahs", "sis", "solos", "SOS", "stats", "stets", "tat", "tenet", "TNT", "toot", "tot", "tut", "wow", "WWW"} Length[Select[DictionaryLookup[], pal]] 70

This shows there are 70 palindromic words in this dictionary. We will do the same with other languages, but ﬁrst let us see what languages are available. DictionaryLookup[All] {"Arabic", "BrazilianPortuguese", "Breton", "BritishEnglish", "Catalan", "Croatian", "Danish", "Dutch", "English", "Esperanto", "Faroese", "Finnish", "French", "Galician", "German", "Hebrew", "Hindi", "Hungarian", "IrishGaelic", "Italian", "Latin", "Polish", "Portuguese", "Russian", "ScottishGaelic", "Spanish", "Swedish"}

We will choose all the palindromic words from all these dictionaries. The following is going to take about a minute. t = Table[{lan, Length[Select[DictionaryLookup[{lan, All}], pal]]}, {lan,DictionaryLookup[All]}] {{"Arabic", 55}, {"BrazilianPortuguese", 81}, {"Breton", 57}, {"BritishEnglish", 94}, {"Catalan", 150}, {"Croatian", 106}, {"Danish", 157}, {"Dutch", 83}, {"English", 70}, {"Esperanto", 14}, {"Faroese", 85}, {"Finnish", 72}, {"French", 41}, {"Galician", 92}, {"German", 20}, {"Hebrew", 159}, {"Hindi", 27}, {"Hungarian", 146}, {"IrishGaelic", 42}, {"Italian", 22}, {"Latin", 21}, {"Polish", 133}, {"Portuguese", 121}, {"Russian", 25}, {"ScottishGaelic", 37}, {"Spanish", 58}, {"Swedish", 60}}

We are now ready to put these on a Bar chart. palnum = Last /@ t {55, 81, 57, 94, 150, 106, 157, 83, 70, 14, 85, 72, 41, 92, 20, 159, 27, 146, 42, 22, 21, 133, 121, 25, 37, 58, 60} BarChart[palnum, BarOrigin -> Left, ChartStyle -> "DarkRainbow", ChartLabels -> DictionaryLookup[All]]

3.3 Selecting from a list

♣

45

TIPS

– In Problem 3.16 we used the command StringReverse, which reverses the order of the characters, so StringReverse["this is great"] "taerg si siht"

Mathematica has several more commands to work with strings, including StringLength, StringTake, StringDrop, StringReplace, ToString, etc.

Problem 3.17 Find all the words in the Mathematica dictionary such that the reverses of the words are in the dictionary as well (for example, loots and stool). =⇒ Solution. There are 496 words in the Mathematica dictionary such that their reverses have a meaning as well. Here is the code: Select[DictionaryLookup[], {StringReverse[#]} == DictionaryLookup[StringReverse[#] ] &]

However, as we don’t want to print all the words, we use Short to get one line of the list.

46

3. Lists

?Short Short[expr] prints as a short form of expr, less than about one line long. Short[expr,n] prints as a form of expr about n lines long. >> Short[ Select[DictionaryLookup[], {StringReverse[#]} == DictionaryLookup[StringReverse[#] ] &]] {a,abut,agar,ah,aha,,yaps,yard,yaw,yaws,yob}

4 Changing heads! Mathematica considers expressions in a uniﬁed manner. Any expression consists of a head and its arguments. For example {a,b,c} is considered as List[a,b,c] with the function List as the head and a,b,c as its arguments. Mathematica enables us to replace the head of an expression with another expression. For example, replacing the head of {a,b,c} with Plus gives Plus[a,b,c] which is a+b+c. This simple idea provides a powerful method to approach solving problems as this chapter demonstrates.

Let us for a moment be a bit abstract. Wolfram Mathematica® has a very consistent way of dealing with expressions. Any expression in Mathematica has the following presentation head[arg1,arg2,...,argn] where head and arg could be expressions themselves. To make this point clear let us use the command FullForm which shows how Mathematica considers an expression. FullForm[a + b + c] Plus[a, b, c] FullForm[a*b*c] Times[a, b, c] FullForm[{a, b, c}] List[a, b, c]

Notice that the expressions a+b+c and {a,b,c} which represent very diﬀerent things have such close presentations. Here Plus is a function and a,b,c are plugged into this function. Plus is the head of the expression a+b+c. One can see from the FullForm that the only diﬀerence between a+b+c and {a,b,c} is their heads! We can get the head of any expression: Head[{a, b, c}] List Head[a + b + c] Plus {a, b, c}[[0]] List R. Hazrat, Mathematica ®: A Problem-Centered Approach, Springer Undergraduate Mathematics Series, DOI 10.1007/978-1-84996-251-3 4, © Springer-Verlag London Limited 2010

48

4. Changing heads!

Mathematica gives us the ability to replace the head of an expression with another head. The consequence of this is simply (head and) mind-blowing! This can be done with the command Apply. Here is the traditional example: Apply[Plus,{a,b,c}] a+b+c

It is not diﬃcult to explain this. The full form of {a,b,c} is List[a,b,c] with the head List. All Mathematica does is to change the head List to Plus, thus we have Plus[a,b,c] which is a+b+c. The shorthand for Apply is @@, as the following example shows: Plus @@ Range[10] 55

This gives the sum of 1 to 10.

Problem 4.1 Deﬁne the following functions: 13 + 23 + · · · + n3 1 1 1 ep(n) = 1 + + + · · · + 1 2! n! p(n) = (1 + x)(1 + x2 ) · · · (1 + xn ). cs(n)

=

=⇒ Solution. Recall that Range[n] gives the list {1, 2, ..., n} and since ˆ is a listable function, Range[n]ˆ 3 produces {13 , 23 , ..., n3 }. Therefore, as in the example above, if we replace the head of the list with Plus, i.e., Plus @@ Range[n]ˆ 3 we will have 13 + 23 + ... + n3 . The other functions use a similar approach. cs[n_]:=Plus @@ Range[n]^3 ep[n_]:=1.+Plus @@ (1/Range[n]!) p[n_]:=Times @@ (1+x^Range[n])

Let us look at the second example also. Again, Range[n] produces a list {1, 2, 3, · · · , n}. Note that the factorial function ! is listable, thus Range[n]! would produce {1!, 2!, 3!, · · · , n!}. Recall that all the arithmetic operations are 1 1 1 1 also listable, including / . Thus 1/Range[n]! produces { 1! , 2! , 3! , · · · , n! }. We 1 1 1 1 are almost there, all we have to do is to replace the head of { 1! , 2! , 3! , · · · , n! } 1 1 1 which is a List with Plus and as a result we get 1! + 2! + · · · + n! . All these functions are classical examples of using Sum and Product which are available in Mathematica. We will meet these commands in Chapter 6.

4. Changing heads!

49

As the expressions get more complicated, it is quite diﬃcult to analyze the FullForm of an expression (if this is necessary at all). The command TreeForm is an excellent facility to shed more light on how the functions are composed to get the expression. Here is an example: FullForm[Sqrt[2 + x]/5 + x/y] Plus[Times[Rational[1,5], Power[Plus[2,x], Rational[1,2]]], Times[x, Power[y, -1]]]

To understand this, start from the bottom of the tree and make your way to the top and you will get the expression √ 2+x x + . 5 y

Problem 4.2 Show that the only n less than 1000 such that 3n + 4n + · · · + (n + 2)n = (n + 3)n are the numbers 2 and 3. =⇒ Solution. We ﬁrst need to write a code to evaluate 3n + 4n + · · · + (n + 2)n . This is very similar to the functions in Problem 4.1 and it looks like this Apply[Plus, Range[3, n+2]^n]

50

4. Changing heads!

Or, if using the shorthand, the equivalent code is Plus @@ Range[3, n + 2]ˆ n. We then select from the range of 1 to 1000, those in which the result of the above is the same as (n+3)ˆ n. Note that in order to plug these into Select we need to make these two into one anonymous function: Select[Range[1000],Apply[Plus, Range[3, # + 2]^#] == (# + 3)^# &] {2, 3}

Exercise 4.1 A number is called a Harshad number1 if it is divisible by the sum of its digits (e.g., 12 is Harshad as it is divisible by 1+2=3). Find all 2-digit Harshad numbers. How many 5-digit Harshad numbers are there? Let’s look at Problem 3.14 again.

Problem 4.3 Write a function squareFreeQ[n] that returns True if the number n is a square free number, and False otherwise. =⇒ Solution. Here is the whole code: squareFree1Q[n_] := Times @@ Last /@ FactorInteger[n] == 1 squareFree1Q /@ {12, 13, 14, 25, 26} {False, True, True, False, True}

It might take a while to understand how this code works. If n = pk11 · · · pkt t , then FactorInteger[n] will produce {{p1 , k1 }, {p2 , k2 }, · · · , {pt , kt }}. We are after numbers such that all the ki are 1 in the decomposition. Thus we can get all the ki , multiply them and if we get anything other than 1, then this would be a non-square free number. Thus the ﬁrst step is to apply Last to the list Last /@ FactorInteger[n] to get {k1 , k2 , · · · , kt }. Then all we have to do is to multiply them all together, and here comes the Times @@ to change the head of {k1 , k2 , · · · , kt } from List to Times.

1

Harshad means “giving joy” in Sanskrit, named and deﬁned by the Indian mathematician D. Kaprekar.

4. Changing heads!

51

Problem 4.4 Find all the numbers up to one million which have the following property: if n = d1 d2 · · · dk then n = d1 ! + d2 ! + · · · + dk ! (e.g. 145 = 1! + 4! + 5!). =⇒ Solution. Select[Range[1000000], Plus @@ Factorial /@ IntegerDigits[#] == # &] {1, 2, 145, 40585}

The code consists of an anonymous function which, for any n, checks whether it has the desired property of the problem. Then, by using Select, we check the list of all the numbers from 1 to one million, Range[1000000]. Our anonymous function is Plus @@ Factorial /@ IntegerDigits[#] == # &. Let us look at the left-hand side of ==. The built-in function IntegerDigits[#] applied to n = d1 d2 · · · dk produces the list of digits of n, namely {d1 , d2 , · · · , dk }. Next, applying Factorial /@ to this list, we get {d1 !, d2 !, · · · , dk !}. Now all we need is to get the sum of elements of this list, and this is possible by changing the head from List to Plus by Plus @@. Once this is done, we compare the left-hand side of == with the right-hand side which is the original number #.

Problem 4.5 A number is perfect if it is equal to the sum of its proper divisors, e.g., 6 = 1 + 2 + 3 but 18 = 1 + 2 + 3 + 6 + 9. Write a program to ﬁnd all the perfect numbers up to 10000. (Hint, have a look at the command Divisors.) =⇒ Solution. Here is a step-by-step approach to the solution. Divisors[6] {1,2,3,6} Most[Divisors[6]] {1,2,3} Apply[Plus,Most[Divisors[6]]] 6 Select[Range[10000], # == Apply[Plus, Most[Divisors[#]]] &] {6, 28, 496, 8128}

The numbers 6, 28 and 496 were already known as being perfect numbers 2000 years before Christ. A glance at the list shows that all the perfect num-

52

4. Changing heads!

bers we have found are even. It is still unknown whether there is an odd perfect number. Probably this is the oldest unsolved question in mathematics!

Problem 4.6 Among the ﬁrst 100000 numbers, what is the largest number n which is divisible √ by all positive integers ≤ n? =⇒ Solution. Select[Range[100000],(Mod[#,LCM @@ Range[Floor[Sqrt[#]]]]==0)&] {1, 2, 3, 4, 6, 8, 12, 24}

Exercise 4.2 Decipher what the following codes do: g[n_] := Times @@ Apply[Plus , Inner[List, x^Range[n], 1/x^Range[n], List], 1] t[n_] := Times @@ Apply[Plus, Thread[List[x^Range[n], 1/x^Range[n]]], 1]

Exercise 4.3 Find all the numbers up to 100 which have the following property: if n = pk11 · · · pkt t is the prime decomposition of n then n = k1 × p1 + k2 × p 2 + · · · + kt × p t .

Exercise 4.4 Consider all numbers of the form 3n2 + n + 1, where 0 < n < 10000. How small can the sum of the digits of such a number be? (Hint, see Min.)

Exercise 4.5 Decipher what the following codes do: Power @@ (x + y) Plus @@@ (x^y + y^z)

4. Changing heads!

53

Exercise 4.6 Show that the sum of all the divisors of the number 608655567023837898967037173424316962265783077 3351885970528324860512791691264 is a perfect number. (This number is the only known sublime number besides 12. A sublime number is a positive integer which has a perfect number of positive divisors (including itself), and whose positive divisors add up to another perfect number. See Problem 4.5 for perfect numbers and MathWorld [7] for sublime numbers.)

Exercise 4.7 A weird number is a number such that the sum of the proper divisors (divisors including 1 but not itself) of the number is greater than the number, but no subset of those divisors sums to the number itself. Find all the weird numbers up to 10000.

5 A bit of logic and set theory

This chapter focuses on statements which take values of True or False, i.e., Boolean statements. We can combine statements of this type to make decisions based on their values. We introduce the If statement available in Mathematica.

5.1 Being logical In mathematical logic, statements can have a value of True, False or undeﬁned.1 These are called Boolean expressions. This helps us to “make a decision” and write programs based on the value of a statement (I am thinking of the classical If-Then statement; If something is True, Then do this, Otherwise do that). We have seen == which compares the left-hand side with the righthand side. Studying the following examples carefully will tell us how Wolfram Mathematica® approaches logical statements: 3^2+4^2==5^2 True 3^2+4^2>5^2 False 9Sqrt[10!] < 10Sqrt[9!] False (x-1)(x+1)==x^2-1 (x-1)(x+1)==x^2-1 1

We don’t want to go into detail here mainly because I don’t know the detail!

R. Hazrat, Mathematica ®: A Problem-Centered Approach, Springer Undergraduate Mathematics Series, DOI 10.1007/978-1-84996-251-3 5, © Springer-Verlag London Limited 2010

5.1 Being logical

55

Simplify[%] True x==5 x==5 {1,2}=={2,1} False {a,b}=={b,a} {a,b}=={b,a}

As one notices, Mathematica echoes back the expressions that she can’t evaluate (e.g., x==5). Among them {a,b}=={b,a}, although one expects to get False as lists respect order. This is because Mathematica does not know about the values of a and b, and in case a and b are the same then {a,b}=={b,a} is True, and False otherwise. If you want Mathematica to judge from the face value, then use ===, x==5 False {a,b}=={b,a} False ?=== ihs===rhs yields True if the expression lhs is identical to rhs and yields False otherwise. >>

One can combine logical statements with usual operations And, Or, Not, or the equivalent &&, ||, ! as the following examples show: 2 > 3 && 3 > 2 False And[2 > 3, 3 > 2] False 1 < 2 < 3 True 2 > 3 || 2 < 3 True Or[2 > 3, 3 > 2] True 3^2+4^2>= 5^2 True

In general, for two statements A and B, the statement A&&B is false if one of A or B is false and A || B is true if one of them is true. In order to produce all possible combinations of true and false cases, we use the command Outer

56

5. A bit of logic and set theory

as the following example shows (we will look at the command Outer in more detail in Section 7.5). Outer[f, {a, b}, {x, y}] {{f[a, x], f[a, y]}, {f[b, x], f[b, y]}}

Thus, if in the above we replace f by And or Or we will get all the possible combinations of True and False. Outer[And, {True, False}, {True, False}] {{True, False},{False,False}} Outer[Or, {True, False}, {True, False}] {{True, True}, {True,False}}

In Mathematica, for a variable, one can specify certain domains. This means that the variable takes its values from a speciﬁc type of data. The domains available are Algebraics, Booleans, Complexes, Integers, Primes, Rationals and Reals. One of the fundamental theorems in number theory is to show that π is not a rational number, i.e., is not of the form m/n, where m and n are integers. Look at the following examples: Pi ∈ Rationals False Sqrt[7] ∈ Integers False

Problem 5.1

√ √ √ Show that 1 + 3 + 5 + 7 is an algebraic number (i.e., it is a solution of a polynomial equation with integer coeﬃcients). =⇒ Solution. Plus @@ Sqrt[Range[1, 7, 2]] ∈ Algebraics True √ √ √ The last example shows that 1 + 3 + 5 + 7 is an algebraic number (i.e. it is a solution of a polynomial equation with integer coeﬃcients). One can use membership (∈) to write neat solutions to several problems.

Problem 5.2 Does the formula (n!)2 + 1 give prime numbers for n = 1 to 6?

5.1 Being logical

57

=⇒ Solution. (#!^2 + 1) & /@ Range[6] ∈ Primes False Here we ﬁrst apply the anonymous function (#!^2 + 1) which is the formula (n!)2 + 1 to the list containing 1 to 6. Then we ask Mathematica if the elements of this list belong to the domain Primes. The answer is False. The following code shows that the above formula does not produce a prime number for n = 6: PrimeQ /@ ((#!^2 + 1) & /@ Range[1, 6]) {True, True, True, True, True, False}

One should be careful that Mathematica cannot(yet) perform miracles. For √ √ 3 3 example, one can actually prove that 2 + 5 + 2 − 5 is an integer, but (2 + 5^(1/2))^(1/3) + (2 - 5^(1/2))^(1/3) ∈ Integers False FullSimplify[(2 + 5^(1/2))^(1/3) + (2 - 5^(1/2))^(1/3) ∈ Integers] False Mathematica provides the logical quantiﬁers ∀, ∃ and ⇒ with ForAll, Exist and Implies commands. But these seem not to be that powerful yet. For example, one cannot prove Fermat’s little theorem which says 2p−1 ≡ 1(mod p) where p > 2 is a prime number with them! ForAll[p, p ∈ Primes, Mod[2^(p - 1), p] == 1] Or even an easy fact that the product of four consecutive numbers plus one is a squared number. Implies[n ∈ Integers && n > 0, Sqrt[n(n + 1)(n + 2)(n + 3) + 1] ∈ Integers] In both cases Mathematica gives back the same expression, indicating she cannot decide on them.

58

5. A bit of logic and set theory

5.2 Handling sets Now it has been agreed among mathematicians that any mathematics starts by considering sets, i.e., collections of objects.2 As we mentioned, the diﬀerence between mathematical sets and lists in Mathematica is that lists respect order and repetition, which is to say one can have several copies of one object in a list (see Chapter 3). Sets are not sensitive about repeated objects, e.g., the set {a, b} is the same as the set {a, b, b, a}. There is no concept of sets in Mathematica and if necessary one considers a list as a set. If one wants to get rid of duplications in a list, one can use Union[{a,b,b,a}] {a,b}

Considering two sets, the natural operations between them are union and intersection. Mathematica provides Union to collect all elements from diﬀerent lists in one list (after removing all the duplications) and Intersection for collecting common elements (again discarding repeated elements). The following examples show how these commands work. u= {1, 2, 3, 4, 5, 2, 4, 7, 4}; a = {1, 4, 7, 3}; b = {5, 4, 3, 2}; Union[u] {1, 2, 3, 4, 5, 7} a ∪ b {1, 2, 3, 4, 5, 7} a ∩ b {3,4} Complement[u, a] {2, 5} ?Complement Complement[eall, e1, e2, ... ] gives the elements in eall which are not in any of the ei. 2

To be precise, one ﬁrst considers classes.

5.2 Handling sets

59

Complement[u, a ∩ b] == Complement[u, a] ∪ Complement[u, b] True The ﬁrst example shows Union[list] will get rid of repetition in a list. The command Complement[u,a] will give the elements of u which are not in a. From the example one can see that a ∩ b is acceptable in Mathematica and is a shorthand for Intersection[a,b]. In the last example we checked a theorem of set theory that (A ∩ B)c = Ac ∪ B c where c stands for complement.

Example 5.3 The following trick will be used later (in Chapter 7, inside a loop) to collect data. A={} A=A {x} {x} A=A {y} {x,y} A=A {z} {x,y,z} This is similar to the traditional trick sum=sum+i. Each time sum=sum+i is performed, i will be added to sum and this result will be the new value of sum. There are other ways to add an element to a list. Append[{a,b,c},d] {a,b,c,d} A={};A=Append[A,x] {x} A=Append[A,y] {x,y} A=Append[A,z] {x,y,z}

Also note that the command AppendTo[s, elem] is equivalent to s = Append[s, elem].

Problem 5.4 How many positive integers n are there such that n is a divisor of at least one of the numbers 1040 and 2030 ?

60

5. A bit of logic and set theory

=⇒ Solution. Recall that Divisors[n] will produce a list of all the positive integers which divide n. So all we need to do is to get the collection of divisors of 1040 and 2030 . This can be done with Union. Note that we are not going to print out all the divisors as this is a long list. Length[Union[Divisors[10^40], Divisors[20^30]]] 2301

Problem 5.5 Find out how many common words there are between the English language with French, Dutch and German respectively. =⇒ Solution. Recall that DictionaryLookup contains all the words in 26 languages (see Problem 3.16). All we have to do is to intersect the list of words in the English language with the list of French words and so on as follows: Length[Intersection[DictionaryLookup[{"English", All}], DictionaryLookup[{"French", All}]]] 6897 Length[Intersection[DictionaryLookup[{"English", All}], DictionaryLookup[{"Dutch", All}]]] 6166 Length[ Intersection[DictionaryLookup[{"English", All}], DictionaryLookup[{"German", All}]]] 1286

♣

TIPS

– The command Tally[list] tallies the elements in list, listing all distinct elements together with their multiplicities. – The command Join[list1,list2] concatenates lists or other expressions that share the same head.

5.3 Decision making, If and Which

61

5.3 Decision making, If and Which The statement If[stat,this,that] where stat is a Boolean expression, i.e., has the value of True or False, will execute this if the stat value is True and that otherwise. That means, in either case, one of the statements this or that will be performed (but not both). So this gives an ability to make a decision about which part one wants to perform. Here is an example: If[12^13+13>13^12+12, Print["12^13+13>13^12+12"], Print["12^13+1313^12+12"], 12^13 + 13 < 13^12 + 12, Print["12^13+13

Advisory Board M.A.J. Chaplain University of Dundee K. Erdmann University of Oxford A. MacIntyre Queen Mary, University of London E. S¨uli University of Oxford J.F. Toland University of Bath

For other titles published in this series, go to www.springer.com/series/3423

Roozbeh Hazrat

Mathematica®: A Problem-Centered Approach

Roozbeh Hazrat Department of Pure Mathematics Queen’s University Belfast BT7 1NN United Kingdom [email protected]

ISSN 1615-2085 ISBN 978-1-84996-250-6 e-ISBN 978-1-84996-251-3 DOI 10.1007/978-1-84996-251-3 Springer London Dordrecht Heidelberg New York British Library Cataloguing in Publication Data A catalogue record for this book is available from the British Library Library of Congress Control Number: 2010929694 Mathematics Subject Classification (2000): 68-01, 68N15 © Springer-Verlag London Limited 2010 Apart from any fair dealing for the purposes of research or private study, or criticism or review, as permitted under the Copyright, Designs and Patents Act 1988, this publication may only be reproduced, stored or transmitted, in any form or by any means, with the prior permission in writing of the publishers, or in the case of reprographic reproduction in accordance with the terms of licenses issued by the Copyright Licensing Agency. Enquiries concerning reproduction outside those terms should be sent to the publishers. The use of registered names, trademarks, etc., in this publication does not imply, even in the absence of a specific statement, that such names are exempt from the relevant laws and regulations and therefore free for general use. The publisher makes no representation, express or implied, with regard to the accuracy of the information contained in this book and cannot accept any legal responsibility or liability for any errors or omissions that may be made. Mathematica and the Mathematica logo are registered trademarks of Wolfram Research, Inc (“WRI” – www.wolfram.com) and are used herein with WRI’s permission. WRI did not participate in the creation of this work beyond the inclusion of the accompanying software, and it offers no endorsement beyond the inclusion of the accompanying software Cover design: deblik Printed on acid-free paper Springer is part of Springer Science+Business Media (www.springer.com)

Preface

Teaching the mechanical performance of routine mathematical operations and nothing else is well under the level of the cookbook because kitchen recipes do leave something to the imagination and judgment of the cook but the mathematical recipes do not. G. P´ olya

This book grew out of a course I gave at Queen’s University Belfast during the period of 2004 to 2009. Although there are many books already written about how to use Wolfram Mathematica® , I noticed they fall into two categories: either they provide an explanation about the commands, in the style of: enter the command, push the button and see the result; or they study some problems and write several-paragraph codes in Mathematica. The books in the ﬁrst category did not inspire me (or my imagination) and the second category were too diﬃcult to understand and not suitable for learning (or teaching) Mathematica’s abilities to do programming and solve problems. I could not ﬁnd a book that I could follow to teach this module. In class one cannot go on forever showing students just how commands in Mathematica work; on the other hand it would be very diﬃcult to follow the codes if one writes a program having more than ﬁve lines (especially as Mathematica’s style of programming provides a condensed code). Thus this book. This book promotes Mathematica’s style of programming. I tried to show when we adopt this approach, how naturally one can solve (nice) problems with (Mathematica) style. Here is an example: Does the Euler formula n2 + n + 41 produce prime numbers for n = 1 to 39?

Or in another problem we tried to show how one can eﬀectively use pattern

viii

Preface

matching to check that for two matrices A and B, (ABA−1 )5 = AB 5 A−1 . One only needs to introduce the fact that AA−1 = 1 and then Mathematica will check the problem by cancelling the inverse elements instead of direct calculation. Although the meaning of the code above may not be clear yet, the reader will observe as we proceed how the codes start making sense, as if this is the most natural way to approach the problems. (People who approach problems with a procedural style of programming (such as C++) will experience that this style replaces their way of thinking!) We have tried to let the reader learn from the codes and avoid long and exhausting explanations, as the codes will speak for themselves. Also we have tried to show that in Mathematica (as in the real world) there are many ways to approach a problem and solve it. We have tried to inspire the imagination! Someone once rightly said that the Mathematica programming language is rather like a “Swiss army knife” containing a vast array of features. Mathematica provides us with powerful mathematical functions. Along with this, one can freely mix diﬀerent styles of programming, functional, list-based and procedural, to achieve a lot. This m´elange of programming styles is what we promote in this note. Thus this book could be considered for a course in Mathematica, or for self study. It mainly concentrates on programming and problem solving in Mathematica. I have mostly chosen problems having something to do with numbers as they do not need any particular background. Many of these problems were taken from or inspired by those collected in [3]. I would like to thank Ilan Vardi for answering my emails and Brian McMaster and Judith Millar for polishing the English of this note. Naoko Morita encouraged me to make my notes into this book. I thank her for this and for always smiling and having a little Geschichte zu erz¨ ahlen. Roozbeh Hazrat [email protected] Belfast, October 2009

How to use this book

Each chapter of the book starts with a description of a new topic and some basic examples. We will then demonstrate the use of new commands in several problems and their solutions. We have tried to let the reader learn from the codes and avoided long and exhausting explanations, as the codes will speak for themselves. There are three diﬀerent categories of problems, shown by diﬀerent frames:

Problem 0.1 These problems are the essential parts of the text where new commands are introduced to solve the problem and where it is demonstrated how the commands are used in Wolfram Mathematica® . These problems should not be skipped. =⇒ Solution.

Problem 0.2 These problems further demonstrate how one can use commands already introduced to tackle diﬀerent situations. The readers are encouraged to try out these problems by themselves ﬁrst and then look at the solution. =⇒ Solution.

x

How to use this book

Problem 0.3 These are more challenging problems that could be skipped in the ﬁrst reading. =⇒ Solution.

Most commands in Mathematica have several extensions. When we introduce a command, this is just the tip of the iceberg. In TIPS, we give further indications about the commands, or new commands to explore. Once one has enough competence, then it is quite easy to learn further abilities of a command by using the Mathematica Help and the examples available there.1

1

Included with this book is a free 30 day trial of the Mathematica software. To access your free download, simply go to http://www.wolfram.com/books/resources and enter license number L3280-8445. You will be guided to download and install the latest version of Mathematica.

The Mathematica philosophy

In the beginning is the expression. Wolfram Mathematica® transforms the expression dictated by the rules to another expression. And that’s how a new idea comes into the world! The rules that will be used frequently can be given a name (we call them functions) r[x_]:=1+x^2 r[arrow] 1+arrow^2 r[5] 26

And the transformation could take place immediately or with a delay {x,x}/.x->Random[] {0.0474307, 0.0474307} {x,x}/.x:>Random[] {0.368461, 0.588353}

The most powerful transformations are those which look for a certain pattern in an expression and morph that to a new pattern. (a + b)^c /. (x_ + y_)^z_ -> (x^z + y^z) a^c + b^c

And one of the most important expressions that we will work with is a list. As the name implies this is just a collection of elements (collection of other expressions). We then apply transformations to each element of the list: x^ {1, 5, 7} {x, x^5, x^7}

The Mathematica philosophy

xii

Any expression is considered as having a “head” followed by several arguments, head[arg1,arg2,...]. And one of the transformations which provide a capable tool is to replace the head of an expression by another head! Plus @@ {a,b,c} a+b+c Power @@ (x+y) x^y

Putting these concepts together creates a powerful way to solve problems. In the chapters of this book, we decipher these concepts.

Contents

1.

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Mathematica as a calculator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Numbers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Algebraic computations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Trigonometric computations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Equalities =, :=, == . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 Dynamic variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1 1 5 9 11 12 13 15

2.

Deﬁning functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.1 Formulas as functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.2 Anonymous functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

3.

Lists . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Functions producing lists . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Listable functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Selecting from a list . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.

Changing heads! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

5.

A bit of logic and set theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Being logical . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Handling sets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Decision making, If and Which . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6.

Sums and products . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 6.1 Sum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

25 28 32 34

54 54 58 61

xiv

Contents

6.2 Product . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 7.

Loops and repetitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1 Do, For a While . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Nested loops . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3 Nest, NestList and more . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.4 Fold and FoldList . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.5 Inner and Outer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

72 72 81 84 90 93

8.

Substitution, Mathematica rules . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

9.

Pattern matching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

10. Functions with multiple deﬁnitions . . . . . . . . . . . . . . . . . . . . . . . . . 112 10.1 Functions with local variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 10.2 Functions with conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 11. Recursive functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 12. Linear algebra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 12.1 Vectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 12.2 Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 13. Graphics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 13.1 Two-dimensional graphs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 13.2 Three-dimensional graphs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 14. Calculus and equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158 14.1 Solving equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158 14.2 Calculus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163 15. Solutions to the Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169 Further reading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185 Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 186

1 Introduction

This chapter introduces basic capabilities of Mathematica, which include simple arithmetic, handling algebraic and trigonometric expressions and assigning values to variables. We will also look at dynamic objects, allowing us to see changes in the variables as they happen.

In this chapter we give a quick introduction to the very basic things one can perform with Wolfram Mathematica® . We let the reader learn from reading the codes and avoid long and exhausting explanations, as the codes will speak for themselves.

1.1 Mathematica as a calculator Mathematica can be used as a powerful calculator with the basic arithmetic operations; +, − for addition and subtraction, ∗, / for multiplication and division and ˆ for powers. 10^9 - 987654321 12345679 2682440^4 + 15365639^4 + 18796760^4 180630077292169281088848499041 20615673^4 180630077292169281088848499041

The last two calculations show that 26824404 + 153656394 + 187967604 = 206156734 ,

R. Hazrat, Mathematica ®: A Problem-Centered Approach, Springer Undergraduate Mathematics Series, DOI 10.1007/978-1-84996-251-3 1, © Springer-Verlag London Limited 2010

2

1. Introduction

disproving a conjecture by Euler that three fourth powers can never sum to a fourth power.1 Mathematica can handle large calculations: 2^9941-1 346088282490851215242960395767413316722628668900238547790489283445006220809834 114464364375544153707533664486747635050186414707093323739706083766904042292657 896479937097603584695523190454849100503041498098185402835071596835622329419680 597622813345447397208492609048551927706260549117935903890607959811638387214329 942787636330953774381948448664711249676857988881722120330008214696844649561469 971941269212843362064633138595375772004624420290646813260875582574884704893842 439892702368849786430630930044229396033700105465953863020090730439444822025590 974067005973305707995078329631309387398850801984162586351945229130425629366798 595874957210311737477964188950607019417175060019371524300323636319342657985162 360474512090898647074307803622983070381934454864937566479918042587755749738339 033157350828910293923593527586171850199425548346718610745487724398807296062449 119400666801128238240958164582617618617466040348020564668231437182554927847793 809917495802552633233265364577438941508489539699028185300578708762293298033382 857354192282590221696026655322108347896020516865460114667379813060562474800550 717182503337375022673073441785129507385943306843408026982289639865627325971753 720872956490728302897497713583308679515087108592167432185229188116706374484964 985490944305412774440794079895398574694527721321665808857543604774088429133272 929486968974961416149197398454328358943244736013876096437505146992150326837445 270717186840918321709483693962800611845937461435890688111902531018735953191561 073191960711505984880700270887058427496052030631941911669221061761576093672419 481606259890321279847480810753243826320939137964446657006013912783603230022674 342951943256072806612601193787194051514975551875492521342643946459638539649133 096977765333294018221580031828892780723686021289827103066181151189641318936578 454002968600124203913769646701839835949541124845655973124607377987770920717067 108245037074572201550158995917662449577680068024829766739203929954101642247764 456712221498036579277084129255555428170455724308463899881299605192273139872912 009020608820607337620758922994736664058974270358117868798756943150786544200556 034696253093996539559323104664300391464658054529650140400194238975526755347682 486246319514314931881709059725887801118502811905590736777711874328140886786742 863021082751492584771012964518336519797173751709005056736459646963553313698192 960002673895832892991267383457269803259989559975011766642010428885460856994464 428341952329487874884105957501974387863531192042108558046924605825338329677719 469114599019213249849688100211899682849413315731640563047254808689218234425381 995903838524127868408334796114199701017929783556536507553291382986542462253468 272075036067407459569581273837487178259185274731649705820951813129055192427102 805730231455547936284990105092960558497123779789849218399970374158976741548307 086291454847245367245726224501314799926816843104644494390222505048592508347618 947888895525278984009881962000148685756402331365091456281271913548582750839078 91469979019426224883789463551

If a number of the form 2n − 1 happens to be prime, it is called a Mersenne prime. Recall that a prime number is a number which is divisible only by 1 and itself. It is easy to see 22 − 1 and 23 − 1 and 25 − 1 are Mersenne primes. The list continues. In 1963, Gillies found that the above number, 29941 − 1, is a Mersenne prime. With my laptop it takes about 3 seconds for Mathematica to check that this is a prime number.2 PrimeQ[2^9941-1] True

Back to easier calculations: 24/17 24/17

24 Mathematica always tries to give a precise value, thus gives back instead 17 of attempting to evaluate the fraction. 1

2

This conjecture remained open for almost 200 years, until Noam Elkies at Harvard came up with the above counterexample in 1988. The largest Mersenne prime found so far is 243,112,609 − 1 which was discovered in August 2008 at UCLA and has 12, 978, 189 digits.

1.1 Mathematica as a calculator

3

Sin[Pi/5] √ 1 1 (5 − 5) 2 2 √ Mathematica gives 12 12 (5 − 5) as the value of sin(π/5), which is the precise equality. This shows Mathematica is not approaching the expressions numerically. In order to get the numerical value, one can use the function N. N[24/17] 1.41176 N[24/17, 20] 1.4117647058823529412 ?N N[expr] gives the numerical value of expr. N[expr, n] attempts to give a result with n-digit precision.

As the above line shows, in order to get a quick description of a command, use ?Command. If you need more explanation use ??Command. If you need even more, use Mathematica’s Document Center in the Help Menu and search for the command. Mathematica comes with an excellent Help which contains explanations and many nice examples demonstrating how to use each of the functions available in this software. All elementary mathematical functions are available here, Log, Exp, Sqrt, Sin, Cos, Tan, ArcSin,.... Here we evaluate π ( )2 + (0.5 Log[2])2 , 4 Sqrt[(Pi/4)^2 + (0.5 Log[2])^2] 0.858466

Problem 1.1 Use Mathematica to show that tan

3π 2π √ + 4 sin = 11. 11 11

=⇒ Solution. Tan[3 Pi/11] + 4 Sin[2 Pi/11] Cot[(5 Pi)/22] + 4 Sin[(2 Pi)/11]

We didn’t get

√

11 as an answer. We ask Mathematica to try a bit harder.

4

1. Introduction

?Simplify Simplify[expr] performs a sequence of algebraic and other transformations on expr, and returns the simplest form it finds. Simplify[expr,assum] does simplification using assumptions. >> Simplify[Tan[3 Pi/11] + 4 Sin[2 Pi/11]] Cot[(5 Pi)/22] + 4 Sin[(2 Pi)/11] ?FullSimplify FullSimplify[expr] tries a wide range of transformations on expr involving elementary and special functions, and returns the simplest form it finds. FullSimplify[expr,assum] does simplification using assumptions. >> FullSimplify[Tan[3 Pi/11] + 4 Sin[2 Pi/11]] Sqrt[11]

This problem introduces the commands Simplify and FullSimplify. If one is not happy with the result, one can always use Simplify and even FullSimplify to have Mathematica work a bit harder to come up with more simpliﬁcation.

For a complete list of elementary functions have a look at Functional Navigator: Mathematics and Algorithms: Mathematical Functions in Mathematica’s Help.

Exercise 1.1 Show that

√ 3 64 22 + (1/2)2 − 1 = 4.

Problem 1.2 Using Mathematica, explain why 4 + 6/4 ∗ 3ˆ− 2 + 1 =

31 . 6

=⇒ Solution. If you look at Mathematica’s Help, under “Special Ways to Input Expressions”, you will see the following note: “The Mathematica language has a definite grammar which speciﬁes how your input should be converted to internal form.” One aspect of the grammar is that it speciﬁes how pieces of your input should be grouped. The general rule is that if ⊗ has higher precedence than ⊕, then a⊕b⊗c is interpreted as a⊕(b⊗c), and a⊗b⊕c is interpreted as (a⊗b)⊕c. You will then ﬁnd a long table listing which operation has a higher precedence

1.2 Numbers

5

and thus, based on that, you will be able to explain why 4 + 6/4 ∗ 3ˆ− 2 + 1 . amounts to 31 6 However, common sense tells us that instead of creating an ambiguous expression such as 4 + 6/4 ∗ 3ˆ− 2 + 1, one should use parentheses () to group objects together and make the expression more clear. For example, one could write 4 + (6/4) ∗ 3ˆ(−2) + 1, or even better use Mathematica’s Palette (Basic Math Assistance) and type 4+

♣

6 × 3−2 + 1. 4

TIPS

– Comments can be added to the codes using (* comment *). (* the most beautiful theorem in Mathematics *) E^(I Pi) + 1 0

– The symbol % refers to the previous output produced. %% refers to the second previous output and so on. – If in calculations you don’t get what you are expecting, use Simplify or even FullSimplify (see Problem 1.1). – To get the numerical approximation, use N[expr] or alternatively, expr//N (see Problem 2.2 for diﬀerent ways of applying a function to a variable). Use EngineeringForm[expr,n] and ScientificForm[expr,n] to get other forms of numerical approximations to n signiﬁcant digits. – The mathematical constant e, the exponential number, is deﬁned in Mathematica as E, or Exp. To get en use either E^ n or Exp[n].

1.2 Numbers There are several standard ways to start with an integer number and produce new numbers out of it. For example, starting from 4, one can form 4 × 3 × 2 × 1, which is represented by 4!.

6

1. Introduction

4! 24 123! 1214630436702532967576624324188129585545421708848338231532891 8161829235892362167668831156960612640202170735835221294047782 59109157041165147218602951990626164673073390741981495296000 0000000000000000000000000

The fundamental theorem of arithmetic states that one can decompose any number n as a product of powers of primes and this decomposition is unique, i.e., n = pk11 · · · pkt t where pi ’s are prime. Thus 12 = 22 × 31 and 37534 = 2 × 72 × 383. Mathematica can do all of these: FactorInteger[12] {{2, 2}, {3, 1}} FactorInteger[37534] {{2,1},{7,2},{383,1}} 2^1 * 7^2 * 383^1 37534 FactorInteger[6473434456376432] {{2,4},{3239053,1},{124909859,1}} PrimeQ[124909859] True Prime[8] 19

Prime[n] produces the n-th prime number. PrimeQ[n] determines whether n is a prime number. More than 2200 years ago Euclid proved that the set of prime numbers is inﬁnite. His proof is used even today in modern books. However, it is not that long ago that we also learned that there is no simple formula that produces only prime numbers. n In 1640 Fermat conjectured that the formula 22 +1 always produces a prime number. Almost a hundred years later the ﬁrst counterexample was found. PrimeQ[2^(2^1)+1] True PrimeQ[2^(2^2)+1] True PrimeQ[2^(2^3)+1] True

1.2 Numbers

7

PrimeQ[2^(2^4)+1] True PrimeQ[2^(2^5)+1] False 2^(2^5)+1 4294967297 FactorInteger[2^(2^5)+1] {{641,1},{6700417,1}} 5

This shows that 22 + 1 is not a prime number. In fact it decomposes into 5 two prime numbers 22 + 1 = 641 × 7600417.

Problem 1.3 What is the probability that a randomly chosen 12-digit number will be a prime? =⇒ Solution. The probability is the number of 12-digit prime numbers over the number of all 12-digit numbers. So we start with ﬁnding how many 12-digit numbers exist: 10^13 - 10^12 9000000000000

Next, we will ﬁnd how many 12-digit prime numbers exist. We will use the following built-in function of Mathematica. ?PrimePi PrimePi[x] gives the number of primes less than or equal to x. >> PrimePi[10^13] 346065536839 PrimePi[10^12] 37607912018 N[(346065536839 - 37607912018)/9000000000000]*100 3.42731

So the probability that we randomly pick a 12-digit prime number is only 3.42 percent.

8

1. Introduction

Exercise 1.2 Check that 123456789098765432111 is prime.

Exercise 1.3 Check that for any given positive integer n ≥ 3, the least common multiple of the numbers 1, 2, · · · , n is greater than 2n−1 . (Hint: see LCM.)

Exercise 1.4 Check that the number 32! ends with 7 zeros.

Problem 1.4 The function Mod[m,n] gives the remainder on division of m by n. Let a = 98 and b = 75. Show that aa bb ends with exactly 98 zeros. =⇒ Solution. Clearly you can directly compute aa bb and count the number of zeros as in Exercise 1.4. But this time you have to count whether the number ends with 98 zeros. Here is a more clever way: Check whether the remainder of this number when divided by 1098 is zero, but its remainder on division by 1099 is not. Mod[98^98 75^75, 10^98] 0 Mod[98^98 75^75, 10^99] 5000000000000000000000000000000000000000000000000000000000000 00000000000000000000000000000000000000

Problem 1.5 Recall that for integers m and n, the binomial

n! n is deﬁned as . m m!(n − m)!

Using Mathematica, check that (m + n)! m+n . = m!n! m =⇒ Solution. Not only is the binomial function available in Mathematica, but also Mathematica can perfectly handle it symbolically as the solution to this problem shows. We will talk more about symbolic computations in Section 1.3.

1.3 Algebraic computations

9

Binomial[m + n, n] == (n + m)!/(n! m!) Binomial[m + n, n] == (m + n)!/(m! n!) FullSimplify[Binomial[m + n, n] == (n + m)!/(n! m!)] True

This is another instance where we need to use FullSimplify to make Mathematica work harder to come up with the result. We will discuss the diﬀerent equalities available in Mathematica in Section 1.6. However, for the time being, note that == is used to compare both sides of equations.

There are several more integer functions available in Mathematica, which can be found in Functional Navigator: Mathematics and Algorithms: Mathematical Functions: Integer Functions. ♣

TIPS

– The command NextPrime[n] gives the next prime larger than n and PrimePi[n] gives the number of primes less than or equal to n (see Problem 1.3). – For integers m and n, one can ﬁnd unique numbers q and r such that r is positive, m = qn + r and r > Divisible[2^2008, b[23]] True

Later, in Problem 3.4, we will determine all the positive integers n between 3 and 50 for which 22008 is divisible by b(n) and will see that in fact there are very few n with this property.

Later, in Chapter 10, we will deﬁne functions with conditions, functions with several deﬁnitions and functions containing several lines of code (a procedure).

2.2 Anonymous functions Sometimes we need to “deﬁne a function as we go” and use it on the spot. Mathematica enables us to deﬁne a function without giving it a name (nor any reference to any speciﬁc variables) use it, and then move on! These functions are called anonymous or pure functions. Obviously if we need to use a speciﬁc function frequently, then the best way is to give it a name and deﬁne it as we did in Section 2.1. Here is an anonymous function equivalent to f (x) = x2 + 4: (#^2+4)&

24

2. Deﬁning functions

The expression (#^2+4)& deﬁnes a nameless function. As usual we can plug in data in place of #. The symbol & determines where the deﬁnition of the function is completed. (#^2+4)&[5] 29

Compare the following: r[x_] := x (x + 1) r[4] 20 # (# + 1) &[4] 20 4 // # (# + 1) & 20

Exercise 2.2 What do the following pure functions do? Fibonacci[15]//Mod[#,5]& PrimeQ[#! + 1] [email protected] 18//2^#+#&

Anonymous functions can handle several variables. Here is an example of an anonymous function for f (x, y) = x2 + y 2 . Sqrt[#1^2+#2^2]&[3,4] 5

As you might guess, #1 and #2 refer to the ﬁrst and second variables in the function.

3 Lists

A list is a collection of data. In this chapter we study how to handle lists and have access to the elements in the lists. Functions which generate lists are studied and methods of selecting elements from a list with a speciﬁc property are examined. Get ready for some serious programming!

One can think of a computer program as a function which accepts some (crude) data or information and gives back the data we would like to obtain. Lists are the way Wolfram Mathematica® handles information. Roughly speaking, a list is a collection of objects. The objects could be of any type and pattern. Let us start with an example of a list: {1, -4/7, stuff, 1 - 2 x + x^2}

This looks like a mathematical set. One diﬀerence is that lists respect order: {1, 2} == {2, 1} False

The other diﬀerence is that a list can contain a copy of the same object several times: {1,2,1} == {1,2} False

The natural thing here is to be able to access the elements of a list. Let us deﬁne a list p as follows: p = {x, 1, -8/3, a, b, {c, d, e}, radio} {x, 1, -(8/3), a, b, {c, d, e}, radio}

As this example shows, the list p can have any sort of data, even another list as one of its elements. Here is how we can access the elements of the list:

R. Hazrat, Mathematica ®: A Problem-Centered Approach, Springer Undergraduate Mathematics Series, DOI 10.1007/978-1-84996-251-3 3, © Springer-Verlag London Limited 2010

26

3. Lists

p[[1]] x p[[5]] b p[[-1]] radio p[[{2, 4}]] {1, a} p[[{-2, 5}]] {{c, d, e}, b} p[[-2, {2, 3}]] {d, e}

Examining the above examples reveals that p[[i]] gives the i-th member of the list. There are other commands to access the elements of a list as follows: p = {x, 1, -8/3, a, b, {c, d, e}, radio} {x, 1, -(8/3), a, b, {c, d, e}, radio} First[p] x Last[p] radio Drop[p, 3] {a, b, {c, d, e}, radio} Take[p, 2] {x, 1} Rest[p] {1, -(8/3), a, b, {c, d, e}, radio} Rest[%] {-(8/3), a, b, {c, d, e}, radio} Rest[%] {a, b, {c, d, e}, radio} Most[p] {x, 1, -(8/3), a, b, {c, d, e}} Most[Rest[p]] == Rest[Most[p]] True

Most of these commands are self-explanatory and a close look at the above examples shows what each of them will do (see the Tips on page 5 for %). All these commands and more are listed in the Mathematica Help under Functional Navigator: Core Language: List Manipulation: Elements of Lists.

3. Lists

27

Problem 3.1 Let p={a,b,{c,d},e}. From this list produce the list {a,b,c,d,e}. =⇒ Solution. {p[[1]], p[[2]], p[[3, 1]], p[[3, 2]], p[[4]]} {a, b, c, d, e}

One can also use the command Flatten to get rid of extra { and }. We will discuss this command in Section 7.2, where we will be looking at nested loops which create nested lists. p = {a, b, c, {d, e}, f} {a, b, c, {d, e}, f} Flatten[p] {a, b, c, d, e, f}

One of the secrets of writing codes comfortably is that one should be able to manipulate lists easily. Often in applications, situations like the following arise: – Given {x1 , x2 , · · · , xn } and {y1 , y2 , · · · , yn }, produce {x1 , y1 , x2 , y2 , · · · , xn , yn }. – Given {x1 , x2 , · · · , xn } and {y1 , y2 , · · · , yn }, produce {x1 + y1 , x2 + y2 , · · · , xn + yn }. – Given {x1 , x2 , · · · , xn } and {y1 , y2 , · · · , yn }, produce

{x1 , y1 }, {x1 , y2 }, {x1 , y3 }, · · · , {x1 , yn }

{x2 , y1 }, {x2 , y2 }, · · · , {xn , y1 }, {xn , y2 }, · · · , {xn , yn } .

– Given {x1 , x2 , · · · , xn }, produce {x1 , x1 + x2 , · · · , x1 + x2 + · · · + xn }. – Given {x1 , x2 , · · · , xn }, produce {x1 }, {x2 , . . . , xn } , {x1 , x2 }, {x3 , . . . , xn } . . . {x1 , . . . xn−1 }, {xn } .

28

3. Lists

Mathematica provides us with commands to obtain the above arrangements easily. We will look at these commands in Sections 7.4, 7.5, 12.1 and Problem 9.4. ♣

TIPS

– One can add elements to a list. There are several commands to handle this, including Append, Prepend and Insert. ?Insert Insert[list,elem,n] insert elem at position n in list. If n is negative, the position is counted from the end. Insert[{i, think, i, am}, "therefore", 3] {i, think, therefore, i, am}

See also Example 5.3 for the use of Append. – Commands Sort, Reverse, RotateLeft and RotateRight are available to rearrange the order of a list. – Commands Delete and Drop are available to remove elements from a list.

3.1 Functions producing lists Mathematica provides us with commands of which the output is a list. These commands have a nature of repetition and replace loops in procedural programming (more on this in Chapter 7). Let us look at some of them here before starting to write more serious codes. Range[10] {1, 2, 3, 4, 5, 6, 7, 8, 9, 10} Range[3, 11] {3, 4, 5, 6, 7, 8, 9, 10, 11} Range[2, 17, 4] {2, 6, 10, 14} ?Range Range[imax] generates the list {1,2,...,imax}. Range[imin, imax] generates the list {imin,...,imax}. Range[imin,imax,di] uses step di.>>

One of the most useful commands which produces a list is Table.

3.1 Functions producing lists

29

Table[2 n + 1, {n, 1, 13}] {3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27}

In Table[2n+1,{n,1,13}], n runs from 1 to 13 and each time the function 2n+1 is evaluated. The following example shows how easily we can work symbolically in Mathematica. Table[x^i + y^i, {i, 2, 17, 4}] {x^2 + y^2, x^6 + y^6, x^10 +y^10, x^14 + y^14}

As in the last example of Range, here in Table, i starts from 2 with steps 4 and thus takes the values 2,6,10,14. Here is one more example demonstrating how beautifully Mathematica can handle symbols Table[xi ,{i,1,10}] {x1 , x2 , x3 , x4 , x5 , x6 , x7 , x8 , x9 , x10 }

Problem 3.2 Produce the list of the ﬁrst 30 Fibonacci numbers. =⇒ Solution. All we have to do is to let i run from 1 to 30 and each time plug i into the function Fibonacci[i] which produces the i-th Fibonacci number. This can be done with Table as follows: Table[Fibonacci[i], {i, 1, 30}] {1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025, 121393, 196418, 317811, 514229, 832040}

Problem 3.3 Find all natural numbers n between 1 and 15 for which the polynomial xn + 64 can be written as a product of two nonconstant polynomials with integer coeﬃcients. =⇒ Solution. We have seen that Factor will give a factorization of an expression (if it is possible).

30

3. Lists

Table[Factor[x^i + 64], {i, 1, 10}] {64 + x, 64 + x^2, (4 + x) (16 - 4 x + x^2), (8 - 4 x + x^2) (8 + 4 x + x^2), 64 + x^5, (4 + x^2) (16 - 4 x^2 + x^4), 64 + x^7, (8 - 4 x^2 + x^4) (8 + 4 x^2 + x^4), (4 + x^3) (16 - 4 x^3 + x^6), 64 + x^10}

Although this solves the problem (for 1 ≤ n ≤ 10), the output is not arranged nicely and it is hard to understand the answer. It would be easier if we keep track of i as well. We change the code slightly. For this we need Print ?Print Print[expr] prints expr as output. >>

Here is the improved code: Table[Print[i, "

", Factor[x^i + 64]], {i, 1, 15}];

1

64+x

2

64+x^2

3

(4+x) (16-4 x+x^2)

4

(8-4 x+x^2) (8+4 x+x^2)

5

64+x^5

6

(4+x^2) (16-4 x^2+x^4)

7

64+x^7

8

(8-4 x^2+x^4) (8+4 x^2+x^4)

9

(4+x^3) (16-4 x^3+x^6)

10

64+x^10

11

64+x^11

12 (2-2 x+x^2) (2+2 x+x^2) (4-4 x+2 x^2-2 x^3+x^4) (4+4 x+2 x^2+2 x^3+x^4) 13

64+x^13

14

64+x^14

15

(4+x^5) (16-4 x^5+x^10)

From the list, for n = 3, 4, 6, 8, 9, 12, 15 we have a factorization into two nonconstant polynomials (in fact, for n = 3k and 4k the polynomial can be factorized into two polynomials).

3.1 Functions producing lists

31

Problem 3.4 Determine all the positive integers n between 3 and 50 for which 22008 is divisible by n n n 1+ + + . 1 2 3 =⇒ Solution. Recall that if Mod[m,n] returns zero, then m is divisible by n. We ﬁrst deﬁne a function b(n) = 1 + n1 + n2 + n3 (see Problem 2.3) and then use a Table to check when 22008 is divisible by b(n) for diﬀerent n. b[n_] := Binomial[n, 1] + Binomial[n, 2] + Binomial[n, 3] + 1 Table[Mod[2^2008, b[n]], {n, 3, 50}] {0, 1, 16, 16, 0, 70, 16, 80, 168, 133, 268, 316, 448, 256, 706, 796, 1096, 723, 1092, 1030, 0, 256, 458, 1240, 2704, 2604, 606, 2922, 640, 1664, 2704, 4076, 2824, 1936, 1024, 8336, 256, 7882, 6974, 4192, 4568, 6061, 1076, 6896, 704, 16669, 6856, 16032}

It is clear from the list that only for n = 3, 7 and 23, is 22008 divisible by the above expression (in fact, one can show that these are the only numbers with this property).

Here is a nice example showing the diﬀerence between the two equalities = and :=.

Example 3.5 This example uses BarChart which is a graphic function. The example is based on the discussion in Section 1.6. ?BarChart BarChart[{y1,y2,...}] makes a bar chart with bar lengths y1,y2,... x=RandomInteger[{1,1000}]; BarChart[Table[x,{200}]] x:=RandomInteger[{1,1000}] BarChart[Table[x,{200}]]

In order to understand this example better, get the list generated by Table[x,1000] for each of the deﬁnitions of x individually and compare them.

32

3. Lists

Figure 3.1 Using = as equality

Figure 3.2 Using := as equality

3.2 Listable functions There are times when we would like to apply a function to all the elements of a list. Suppose f is a function and {a,b,c} is a list. We want to be able to “push” the function f inside the list and get {f[a],f[b],f[c]}. Many of Mathematica’s built-in functions have the property that they simply “go inside” a list. This property of a function is called listable. For example Sqrt is a listable function. Sqrt[{a, b, c, d, e}]

√ √ √ √ √ { a, b, c, d, e} All the arithmetic functions are listable: 1+ {a, b, c, d, e} {1 + a, 1 + b, 1 + c, 1 + d, 1 + e} {a, b, c, d, e}^3 {a^3,b^3,c^3,d^3,e^3} 1/{a, b, c, d, e} {1/a,1/b,1/c,1/d,1/e}

We will use the function Sqrt in the following, to show that the product of four consecutive numbers plus one is always a squared number. Table[n (n + 1) (n + 2) (n + 3) + 1, {n, 1, 10}] {25, 121, 361, 841, 1681, 3025, 5041, 7921, 11881, 17161} Sqrt[%] {5, 11, 19, 29, 41, 55, 71, 89, 109, 131}

3.2 Listable functions

33

Not all the functions are listable. If we want to push a function into a list, the command to use is Map. f[{a, b, c, d, e}] f[{a, b, c, d, e}] Map[f, {a, b, c, d, e}] {f[a], f[b], f[c], f[d], f[e]}

The equivalent shorthand to apply a function to a list is /@ as follows: f /@ {a, b, c, d, e} {f[a], f[b], f[c], f[d], f[e]}

Here is an example: Table[(1 + x)^i, {i, 5}] {1 + x, (1 + x)^2, (1 + x)^3, (1 + x)^4, (1 + x)^5} Expand /@ % {1 + x, 1 + 2 x + x^2, 1 + 3 x + 3 x^2 + x^3, 1 + 4 x + 6 x^2 + 4 x^3 + x^4, 1 + 5 x + 10 x^2 + 10 x^3 + 5 x^4 + x^5}

Problem 3.6 The formula n2 + n + 41 has a very interesting property. Observe that this formula produces prime numbers for all n between 0 and 39. =⇒ Solution. First we produce the numbers: Table[n^2 + n + 41, {n, 0, 40}] {41, 43, 47, 53, 61, 71, 83, 97, 113, 131, 151, 173, 197, 223, 251, 281, 313, 347, 383, 421, 461, 503, 547, 593, 641, 691, 743, 797, 853, 911, 971, 1033, 1097, 1163, 1231, 1301, 1373, 1447, 1523, 1601, 1681}

Then we apply PrimeQ to this list. This function is listable. PrimeQ[%] {True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, False}

The list shows that for n between 0 and 39, the above formula produces prime numbers, however, for n = 40, the last number in the list 402 + 40 + 41 is not prime. This formula was found by Euler around 1772. One wonders, if one changes 41 to another number in the formula n2 + n + 41, whether one gets more consecutive prime numbers. We will examine this in Problem 7.2. Of course we could also simply write Table[PrimeQ[n^2 + n + 41], {n, 0, 40}]

34

3. Lists

as well. As we said, there are many ways to approach a problem in Mathematica.

Exercise 3.1 Generate 6 random integers n between 1 and 3095 and show that n6 + 1091 is not prime. (Hint: see RandomInteger for generating random numbers.)

3.3 Selecting from a list So far we have been able to produce a list of data by using functions producing lists. The next step is to be able to choose, from a list, certain data which ﬁt a speciﬁc description. This can be achieved by the command Select as the following problem demonstrates.

Problem 3.7 How many numbers of the form 3n5 + 11, when n varies from 1 to 2000, are prime? =⇒ Solution. First, let us produce the ﬁrst 20 numbers of this form. plist=Table[3 n^5 + 11, {n, 1, 20}] {14, 107, 740, 3083, 9386, 23339, 50432, 98315, 177158, 300011, 483164, 746507, 1113890, 1613483, 2278136, 3145739, 4259582, 5668715, 7428308, 9600011}

The next step, in the spirit of Section 3.2, would be to apply PrimeQ to all the numbers and ﬁnd out which ones are prime. Since this is a listable function it is enough: PrimeQ[plist] {False, True, False, True, False, True, False, False, False, False, False, True, False, True, False, True, False, False, False, False}

Now we are left to count the number of True ones. This is do-able here; 6 of these numbers are prime. However, if we are dealing with a list with 2000 elements, we need to select, from the elements of the list, those with a desired property,1 here being prime. The command Select does just this: 1

Or a desired pattern, more about this later.

3.3 Selecting from a list

35

Select[plist,PrimeQ] {107, 3083, 23339, 746507, 1613483, 3145739}

These are prime numbers in the list plist. The command Length gives the length of a list. If we assemble all the steps in one line we have Length[Select[Table[3 n^5 + 11, {n, 1, 20}], PrimeQ]] 6

Thus to ﬁnd out how many numbers of the form 3n5 + 11 are prime when n runs from 1 to 2000, all we have to do is to change 20 to 2000: Length[Select[Table[3 n^5 + 11, {n, 1, 2000}], PrimeQ]] 97

In a nutshell, Select[list,f], will apply the function f (which returns True or False) to all the elements, say x, of the list and return those elements for which f[x] is true. A function which returns either True or False is called a Boolean function. So far we have seen one Boolean function, that is, PrimeQ. (More on Boolean statements in Chapter 5.)

Problem 3.8 Show that among the ﬁrst 500 Fibonacci numbers, 18 of them are prime. =⇒ Solution. We do this step by step. First we generate the ﬁrst 10 Fibonacci numbers. Then among those we select the ones which are prime. Then using Length we get the size of this list, as follows: fiblist=Table[Fibonacci[i], {i, 1, 10}] {1, 1, 2, 3, 5, 8, 13, 21, 34, 55} Select[fiblist, PrimeQ] {2, 3, 5, 13} Length[Select[Table[Fibonacci[i], {i, 1, 10}], PrimeQ]] 4

Now instead of 10 we generate 500 numbers: Length[Select[Table[Fibonacci[i], {i, 1, 500}], PrimeQ]] 18

36

3. Lists

Exercise 3.2 Show that, among the ﬁrst 450 Fibonacci numbers, the number of odd Fibonacci numbers is twice the number of even ones. (Hint: see OddQ and EvenQ for odd and even numbers.)

Exercise 3.3 Let m be a natural number and A=

(m + 3)3 + 1 . 3m

Find all the integers m less than 500 such that A is an integer. Show that A is always odd. (Hint: see IntegerQ.) Let us look at another example, slightly diﬀerent but of the same nature as Problem 3.7. The following example shows that anonymous functions ﬁt very well with Select.

Problem 3.9 For which 1 ≤ n ≤ 1000 does the formula 2n + 1 produce a prime number? =⇒ Solution. Here is the solution: Select[Range[1000], PrimeQ[2^(#) + 1] &] {1, 2, 4, 8, 16}

Let us take a deep breath and go through this one-liner code slowly. The function PrimeQ[2^(#)+1]& is an anonymous function which gives True if the number 2n + 1 is prime and False otherwise. PrimeQ[2^(#) + 1] &[12] False

Range[1000] creates a list containing the numbers from 1 to 1000. The command Select applies the anonymous function PrimeQ[2^(#)+1]& to each element of this list and in the cases when the result is true, i.e., when 2n + 1 is prime, the element will be selected. Thus {1, 2, 4, 8, 16} are the only numbers that make 2n + 1 a prime number.

Exercise 3.4 Find the number of positive integers 0 < n < 20000 such that 1997 divides n2 + (n + 1)2 . Try the same code for 2009.

3.3 Selecting from a list

37

Exercise 3.5 For integers 2 ≤ n ≤ 200, ﬁnd all n such that n divides (n − 1)! + 1. Show that there are 46 such n.

Problem 3.10 Notice that 122 = 144 and 212 = 441, namely the numbers and their squares are reverses of each other. Find all the numbers up to 10000 with this property. =⇒ Solution. We need to introduce some new built-in functions. IntegerDigits[n] gives a list of the decimal digits of the integer n. We also need Reverse and FromDigits: IntegerDigits[80972] {8, 0, 9, 7, 2} Reverse[%] {2, 7, 9, 0, 8} FromDigits[%] 27908

Thus the above shows we can easily produce the reverse of a number: re[n_] := FromDigits[Reverse[IntegerDigits[n]]] re[12345] 54321 re[6548629] 9268456

Having this function under our belt, the solution to the problem is just one line. Notice that the problem is asking for the numbers n such that re[nˆ2]=re[n]ˆ2 . Select[Range[10000], re[#]^2 == re[#^2] &] {1, 2, 3, 10, 11, 12, 13, 20, 21, 22, 30, 31, 100, 101, 102, 103, 110, 111, 112, 113, 120, 121, 122, 130, 200, 201, 202, 210, 211, 212, 220, 221, 300, 301, 310, 311, 1000, 1001, 1002, 1003, 1010, 1011, 1012, 1013, 1020, 1021, 1022, 1030, 1031, 1100, 1101, 1102, 1103, 1110, 1111, 1112, 1113, 1120, 1121, 1122, 1130, 1200, 1201, 1202, 1210, 1211, 1212, 1220, 1300, 1301, 2000, 2001, 2002, 2010, 2011, 2012, 2020, 2021, 2022, 2100, 2101, 2102, 2110, 2111, 2120, 2121, 2200, 2201, 2202, 2210, 2211, 3000, 3001, 3010, 3011, 3100, 3101, 3110, 3111, 10000}

38

3. Lists

Here is one more example using the command FromDigits. We know that 11 is a prime number. One wonders what is the next prime number consisting only of ones. A wild guess: a number with 23 digits all one? All we have to do is to produce this number then, with PrimeQ, test whether this is prime. Here is one way to generate this number. Table[1, {23}] {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1} FromDigits[%] 11111111111111111111111 PrimeQ[%] True

Here is the code to ﬁnd out which numbers of this kind up to 500 digits are prime. Select[Range[500], PrimeQ[FromDigits[Table[1, {#}]]] &] {2, 19, 23, 317}

Problem 3.11 Show that the number of k between 0 and 1000 for which

1000 is odd is a k

power of 2. =⇒ Solution. First, Range[0,1000] creates a list of numbers from 0 to 1000. Using an anonymous function, each time we plug these numbers into Binomial[1000, #] and check right away if this is an odd number by OddQ[Binomial[1000, #]] &. If this is the case, Select will pick up these numbers. The rest is clear from the code below. Select[Range[0, 1000], OddQ[Binomial[1000, #]] &] {0, 8, 32, 40, 64, 72, 96, 104, 128, 136, 160, 168, 192, 200, 224, 232, 256, 264, 288, 296, 320, 328, 352, 360, 384, 392, 416, 424, 448, 456, 480, 488, 512, 520, 544, 552, 576, 584, 608, 616, 640, 648, 672, 680, 704, 712, 736, 744, 768, 776, 800, 808, 832, 840, 864, 872, 896, 904, 928, 936, 960, 968, 992, 1000} Length[Select[Range[0, 1000], OddQ[Binomial[1000, #]] &]] 64 FactorInteger[%] {{2, 6}}

3.3 Selecting from a list

39

Exercise 3.6 Show that among the ﬁrst 200 primes p, the ones such that the remainder when 19p−1 is divided by p2 is 1 are {3, 7, 13, 43, 137}.

Exercise 3.7 An integer dn dn−1 . . . d1 is called prime-palindromic if dn dn−1 . . . d1 and d1 . . . dn−1 dn are both prime (for example 941). Write a code to ﬁnd all primepalindromic numbers up to 5000. Observe that there are 167 such numbers.

Problem 3.12 Find all positive integers 0 < m < 105 such that the fourth power of the number of positive divisors of m equals m. (Hint: see Divisors.) =⇒ Solution. The function Divisors[n] gives all the positive numbers which divide n. This is a list which includes 1 and n also. We need to get the number of these divisors, thus the function Length. We are looking for each number m such that the fourth power of Length[Divisors[m]] is m. The rest is clear from the code. Select[Range[10^5], Length[Divisors[#]]^4 == # &] {1, 625, 6561}

Exercise 3.8 Show that there is only one positive integer n smaller than 1000 such that n! + (n + 1)! is the square of an integer. The idea of sending a function into a list, i.e., applying a function to each element of a list, seems to be a good one. We have already mentioned that the listable built-in functions are able to go inside a list, like PrimeQ or Prime. Have a look at the Attributes of Prime in the following: ??Prime Prime[n] gives the nth prime number. Attributes[Prime] = {Listable, Protected}

Also recall that if a function is not listable, Map or /@ will push the function into the list as the following problem demonstrates:

40

3. Lists

Problem 3.13 What digit does not appear as the last digit of the ﬁrst 20 Fibonacci numbers? =⇒ Solution. This one-liner code collects all the digits which appears as the last digit: Union[Last /@ (IntegerDigits /@ (Fibonacci /@ Range[20]))] {0, 1, 2, 3, 4, 5, 7, 8, 9}

Thus 6 is the only digit which is not present. Let us understand this code. As the command /@ applies a function to all elements of a list, Fibonacci /@ Range[20] produces the ﬁrst 20 Fibonacci numbers. Fibonacci /@ Range[20] {1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765}

Then IntegerDigits would go inside this list and get the digits of each number. IntegerDigits/@ Fibonacci /@ Range[20] {{1}, {1}, {2}, {3}, {5}, {8}, {1, 3}, {2, 1}, {3, 4}, {5, 5}, {8, 9}, {1, 4, 4}, {2, 3, 3}, {3, 7, 7}, {6, 1, 0}, {9, 8, 7}, {1, 5, 9, 7}, {2, 5, 8, 4}, {4, 1, 8, 1}, {6, 7, 6, 5}}

Then the function Last will get the last digits as required. Last /@ IntegerDigits /@ Fibonacci /@ Range[20] {1, 1, 2, 3, 5, 8, 3, 1, 4, 5, 9, 4, 3, 7, 0, 7, 7, 4, 1, 5}

Union will get rid of any repetitions in the list. (More on this command in Section 5.2.) Since Fibonacci and IntegerDigits are listable functions, one can also write the above code as follows: Union[Last /@ IntegerDigits[Fibonacci[Range[20]]]]

If one does not want to use IntegerDigits then one can use the Mod function to get access to the last digit of a number. ?Mod Mod[m, n] gives the remainder on division of m by n. Mod[264,10] 4 ?Quotient Quotient[m, n] gives the integer quotient of m and n. Quotient[264,10] 26 26*10+4 264

3.3 Selecting from a list

41

Thus another way to write the code is as follows. Note that Mod is also a listable function. Union[Mod[Fibonacci[Range[20]],10]] {0, 1, 2, 3, 4, 5, 7, 8, 9}

In Section 5.2 we will see how to use Mathematica to get directly the digit 6 that we are after, namely to handle sets.

Recall that one can decompose any number n as a product of powers of primes and this decomposition is unique, i.e., n = pk11 · · · pkt t where pi ’s are prime. Let us call a number a square free number if, in its decomposition to primes, all the ki ’s are 1. Namely, no power of primes can divide this number. Thus 15 = 3 × 5 is a square free number but 48 = 24 × 3 is not. Recall that Select[list,f] will apply the function f (which returns True or False) to all the elements, say x, of the list and return those elements for which f[x] is true. There is an option in Select which makes it possible to get only the ﬁrst n elements of list that satisfy f, i.e., f returns True, as follows: Select[list,f,n]. This comes in very handy, as in some problems we are only interested in a certain number of data which satisfy f. Also this can be used in circumstances where we want to test the elements until something goes wrong or some desirable element comes up. The following example demonstrates this.

Problem 3.14 Write a function squareFreeQ[n] that returns True if the number n is a square free number, and False otherwise. =⇒ Solution. Here is the code: squareFreeQ[n_]:= Select[Last/@ FactorInteger[n], # != 1&,1] == {} squareFreeQ[2*5*7] True squareFreeQ[2*5*7*7] False

Let us decipher this code. Recall from Section 1.2 that if n = pk11 · · · pkt t is the decomposition of n into its prime factors then FactorInteger[n]

42

3. Lists

{{p1 , k1 }, {p2 , k2 }, · · · , {pt , kt }} Now all we have to do is to go through this list and see if all ki ’s are one. So the ﬁrst step is to apply Last to each list to discard pi ’s and get ki ’s. Last /@ FactorInteger[n]

{k1 , k2 , · · · , kt } Having this list, we shall go through the list one by one and examine if ki ’s are one. The anonymous function #= 1& does exactly this. So Select[{k1 , k2 , · · · , kt },#= 1&] gives the list of ki ’s which are not one. But in our case, looking for square free primes, it is enough if only one ki is not one. Then the number is not square free. Thus we use an option of Select which goes through the list until it ﬁnds an element such that ki is not one. So we need to modify the code to Select[{k1 , k2 , · · · , kt },#= 1&,1]. We are almost done. All we have to do is to see if this list is empty or not (namely is there any ki not equal to one). And for this we compare Select[{k1 , k2 , · · · , kt },#= 1&,1]=={}. We can solve this problem later with a slightly diﬀerent method (see Problem 4.3).

Exercise 3.9 Find the ﬁrst 5 positive integers n such that n6 + 1091 is prime. Show that all these n are between 3500 and 8500.

Exercise 3.10 A number with n digits is called cyclic if multiplication by 1, 2, 3, · · · , n produces the same digits in a diﬀerent order. Find the only 6-digit cyclic number.

Problem 3.15 Find out how many primes bigger than n and smaller than 2n exist, when n goes from 1 to 30. =⇒ Solution. We deﬁne an anonymous function which ﬁnds all the primes bigger than n and smaller than 2n and then gets the size of this list. Once we are done with this, we apply this function to a list of numbers from 1 to 30. Our anonymous function looks like this: Length[Select[Range[# + 1, 2 # - 1], PrimeQ]] &

3.3 Selecting from a list

43

Here, Range[# + 1, 2 # - 1] produces all the numbers between n and 2n. Then Select ﬁnds out which of them are in fact prime. Then we use the command Length to get the number of elements of this list. One can optimize this a bit, as we don’t need to look at the whole range of n to 2n, as clearly even numbers are not prime so we can ignore them right from the beginning. But we leave it to the reader to do this. All we have to do now is to apply this function with Map or /@ to numbers from 1 to 30. Length[Select[Range[# + 1, 2 # - 1], PrimeQ]] & /@ Range[30] {0, 1, 1, 2, 1, 2, 2, 2, 3, 4, 3, 4, 3, 3, 4, 5, 4, 4, 4, 4, 5, 6, 5, 6, 6, 6, 7, 7, 6, 7}

One can also write a more innocent code, using Table as follows: Table[Length[Select[Range[n + 1, 2n - 1], PrimeQ]], {n, 1, 30}] {0, 1, 1, 2, 1, 2, 2, 2, 3, 4, 3, 4, 3, 3, 4, 5, 4, 4, 4, 4, 5, 6, 5, 6, 6, 6, 7, 7, 6, 7}

This doesn’t seem to be the most eﬃcient way to write this problem as each time we test the same numbers repeatedly to see whether they are prime. We will write another code for this problem in Problem 7.1 using a Do Loop. One can also use the built-in function PrimePi to write a very short solution for this problem. For the usage of PrimePi see Problem 1.3.

A word is called palindromic if it reads the same backwards as forwards, e.g., madam. In the next problem we are going to ﬁnd all the palindromic words in a lot of languages! The approach is similar to Problem 3.10, however, in Problem 3.16 we will be working with strings of characters.

Problem 3.16 Write a function to check whether a word is palindromic. Using DictionaryLookup, ﬁnd all the palindromic words in the English language. Then draw a bar chart showing how many palindromic words there are in diﬀerent languages which are supported by Mathematica. =⇒ Solution. We ﬁrst deﬁne a function to check whether a word is palindromic. Using StringReverse, this is easy: pal[n_] := n == StringReverse[n] pal["test"] False pal["kayak"] True

44

3. Lists

Now we select from the English dictionary available in Mathematica all the words which are palindromic. Select[DictionaryLookup[], pal] {"a", "aha", "aka", "bib", "bob", "boob", "bub", "CFC", "civic", "dad", "deed", "deified", "did", "dud", "DVD", "eke", "ere", "eve", "ewe", "eye", "gag", "gig", "huh", "I", "kayak", "kook", "level", "ma’am", "madam", "mam", "MGM", "minim", "mom", "mum", "nan", "non", "noon", "nun", "oho", "pap", "peep", "pep", "pip", "poop", "pop", "pup", "radar", "redder", "refer", "repaper", "reviver", "rotor", "sagas", "sees", "seres", "sexes", "shahs", "sis", "solos", "SOS", "stats", "stets", "tat", "tenet", "TNT", "toot", "tot", "tut", "wow", "WWW"} Length[Select[DictionaryLookup[], pal]] 70

This shows there are 70 palindromic words in this dictionary. We will do the same with other languages, but ﬁrst let us see what languages are available. DictionaryLookup[All] {"Arabic", "BrazilianPortuguese", "Breton", "BritishEnglish", "Catalan", "Croatian", "Danish", "Dutch", "English", "Esperanto", "Faroese", "Finnish", "French", "Galician", "German", "Hebrew", "Hindi", "Hungarian", "IrishGaelic", "Italian", "Latin", "Polish", "Portuguese", "Russian", "ScottishGaelic", "Spanish", "Swedish"}

We will choose all the palindromic words from all these dictionaries. The following is going to take about a minute. t = Table[{lan, Length[Select[DictionaryLookup[{lan, All}], pal]]}, {lan,DictionaryLookup[All]}] {{"Arabic", 55}, {"BrazilianPortuguese", 81}, {"Breton", 57}, {"BritishEnglish", 94}, {"Catalan", 150}, {"Croatian", 106}, {"Danish", 157}, {"Dutch", 83}, {"English", 70}, {"Esperanto", 14}, {"Faroese", 85}, {"Finnish", 72}, {"French", 41}, {"Galician", 92}, {"German", 20}, {"Hebrew", 159}, {"Hindi", 27}, {"Hungarian", 146}, {"IrishGaelic", 42}, {"Italian", 22}, {"Latin", 21}, {"Polish", 133}, {"Portuguese", 121}, {"Russian", 25}, {"ScottishGaelic", 37}, {"Spanish", 58}, {"Swedish", 60}}

We are now ready to put these on a Bar chart. palnum = Last /@ t {55, 81, 57, 94, 150, 106, 157, 83, 70, 14, 85, 72, 41, 92, 20, 159, 27, 146, 42, 22, 21, 133, 121, 25, 37, 58, 60} BarChart[palnum, BarOrigin -> Left, ChartStyle -> "DarkRainbow", ChartLabels -> DictionaryLookup[All]]

3.3 Selecting from a list

♣

45

TIPS

– In Problem 3.16 we used the command StringReverse, which reverses the order of the characters, so StringReverse["this is great"] "taerg si siht"

Mathematica has several more commands to work with strings, including StringLength, StringTake, StringDrop, StringReplace, ToString, etc.

Problem 3.17 Find all the words in the Mathematica dictionary such that the reverses of the words are in the dictionary as well (for example, loots and stool). =⇒ Solution. There are 496 words in the Mathematica dictionary such that their reverses have a meaning as well. Here is the code: Select[DictionaryLookup[], {StringReverse[#]} == DictionaryLookup[StringReverse[#] ] &]

However, as we don’t want to print all the words, we use Short to get one line of the list.

46

3. Lists

?Short Short[expr] prints as a short form of expr, less than about one line long. Short[expr,n] prints as a form of expr about n lines long. >> Short[ Select[DictionaryLookup[], {StringReverse[#]} == DictionaryLookup[StringReverse[#] ] &]] {a,abut,agar,ah,aha,,yaps,yard,yaw,yaws,yob}

4 Changing heads! Mathematica considers expressions in a uniﬁed manner. Any expression consists of a head and its arguments. For example {a,b,c} is considered as List[a,b,c] with the function List as the head and a,b,c as its arguments. Mathematica enables us to replace the head of an expression with another expression. For example, replacing the head of {a,b,c} with Plus gives Plus[a,b,c] which is a+b+c. This simple idea provides a powerful method to approach solving problems as this chapter demonstrates.

Let us for a moment be a bit abstract. Wolfram Mathematica® has a very consistent way of dealing with expressions. Any expression in Mathematica has the following presentation head[arg1,arg2,...,argn] where head and arg could be expressions themselves. To make this point clear let us use the command FullForm which shows how Mathematica considers an expression. FullForm[a + b + c] Plus[a, b, c] FullForm[a*b*c] Times[a, b, c] FullForm[{a, b, c}] List[a, b, c]

Notice that the expressions a+b+c and {a,b,c} which represent very diﬀerent things have such close presentations. Here Plus is a function and a,b,c are plugged into this function. Plus is the head of the expression a+b+c. One can see from the FullForm that the only diﬀerence between a+b+c and {a,b,c} is their heads! We can get the head of any expression: Head[{a, b, c}] List Head[a + b + c] Plus {a, b, c}[[0]] List R. Hazrat, Mathematica ®: A Problem-Centered Approach, Springer Undergraduate Mathematics Series, DOI 10.1007/978-1-84996-251-3 4, © Springer-Verlag London Limited 2010

48

4. Changing heads!

Mathematica gives us the ability to replace the head of an expression with another head. The consequence of this is simply (head and) mind-blowing! This can be done with the command Apply. Here is the traditional example: Apply[Plus,{a,b,c}] a+b+c

It is not diﬃcult to explain this. The full form of {a,b,c} is List[a,b,c] with the head List. All Mathematica does is to change the head List to Plus, thus we have Plus[a,b,c] which is a+b+c. The shorthand for Apply is @@, as the following example shows: Plus @@ Range[10] 55

This gives the sum of 1 to 10.

Problem 4.1 Deﬁne the following functions: 13 + 23 + · · · + n3 1 1 1 ep(n) = 1 + + + · · · + 1 2! n! p(n) = (1 + x)(1 + x2 ) · · · (1 + xn ). cs(n)

=

=⇒ Solution. Recall that Range[n] gives the list {1, 2, ..., n} and since ˆ is a listable function, Range[n]ˆ 3 produces {13 , 23 , ..., n3 }. Therefore, as in the example above, if we replace the head of the list with Plus, i.e., Plus @@ Range[n]ˆ 3 we will have 13 + 23 + ... + n3 . The other functions use a similar approach. cs[n_]:=Plus @@ Range[n]^3 ep[n_]:=1.+Plus @@ (1/Range[n]!) p[n_]:=Times @@ (1+x^Range[n])

Let us look at the second example also. Again, Range[n] produces a list {1, 2, 3, · · · , n}. Note that the factorial function ! is listable, thus Range[n]! would produce {1!, 2!, 3!, · · · , n!}. Recall that all the arithmetic operations are 1 1 1 1 also listable, including / . Thus 1/Range[n]! produces { 1! , 2! , 3! , · · · , n! }. We 1 1 1 1 are almost there, all we have to do is to replace the head of { 1! , 2! , 3! , · · · , n! } 1 1 1 which is a List with Plus and as a result we get 1! + 2! + · · · + n! . All these functions are classical examples of using Sum and Product which are available in Mathematica. We will meet these commands in Chapter 6.

4. Changing heads!

49

As the expressions get more complicated, it is quite diﬃcult to analyze the FullForm of an expression (if this is necessary at all). The command TreeForm is an excellent facility to shed more light on how the functions are composed to get the expression. Here is an example: FullForm[Sqrt[2 + x]/5 + x/y] Plus[Times[Rational[1,5], Power[Plus[2,x], Rational[1,2]]], Times[x, Power[y, -1]]]

To understand this, start from the bottom of the tree and make your way to the top and you will get the expression √ 2+x x + . 5 y

Problem 4.2 Show that the only n less than 1000 such that 3n + 4n + · · · + (n + 2)n = (n + 3)n are the numbers 2 and 3. =⇒ Solution. We ﬁrst need to write a code to evaluate 3n + 4n + · · · + (n + 2)n . This is very similar to the functions in Problem 4.1 and it looks like this Apply[Plus, Range[3, n+2]^n]

50

4. Changing heads!

Or, if using the shorthand, the equivalent code is Plus @@ Range[3, n + 2]ˆ n. We then select from the range of 1 to 1000, those in which the result of the above is the same as (n+3)ˆ n. Note that in order to plug these into Select we need to make these two into one anonymous function: Select[Range[1000],Apply[Plus, Range[3, # + 2]^#] == (# + 3)^# &] {2, 3}

Exercise 4.1 A number is called a Harshad number1 if it is divisible by the sum of its digits (e.g., 12 is Harshad as it is divisible by 1+2=3). Find all 2-digit Harshad numbers. How many 5-digit Harshad numbers are there? Let’s look at Problem 3.14 again.

Problem 4.3 Write a function squareFreeQ[n] that returns True if the number n is a square free number, and False otherwise. =⇒ Solution. Here is the whole code: squareFree1Q[n_] := Times @@ Last /@ FactorInteger[n] == 1 squareFree1Q /@ {12, 13, 14, 25, 26} {False, True, True, False, True}

It might take a while to understand how this code works. If n = pk11 · · · pkt t , then FactorInteger[n] will produce {{p1 , k1 }, {p2 , k2 }, · · · , {pt , kt }}. We are after numbers such that all the ki are 1 in the decomposition. Thus we can get all the ki , multiply them and if we get anything other than 1, then this would be a non-square free number. Thus the ﬁrst step is to apply Last to the list Last /@ FactorInteger[n] to get {k1 , k2 , · · · , kt }. Then all we have to do is to multiply them all together, and here comes the Times @@ to change the head of {k1 , k2 , · · · , kt } from List to Times.

1

Harshad means “giving joy” in Sanskrit, named and deﬁned by the Indian mathematician D. Kaprekar.

4. Changing heads!

51

Problem 4.4 Find all the numbers up to one million which have the following property: if n = d1 d2 · · · dk then n = d1 ! + d2 ! + · · · + dk ! (e.g. 145 = 1! + 4! + 5!). =⇒ Solution. Select[Range[1000000], Plus @@ Factorial /@ IntegerDigits[#] == # &] {1, 2, 145, 40585}

The code consists of an anonymous function which, for any n, checks whether it has the desired property of the problem. Then, by using Select, we check the list of all the numbers from 1 to one million, Range[1000000]. Our anonymous function is Plus @@ Factorial /@ IntegerDigits[#] == # &. Let us look at the left-hand side of ==. The built-in function IntegerDigits[#] applied to n = d1 d2 · · · dk produces the list of digits of n, namely {d1 , d2 , · · · , dk }. Next, applying Factorial /@ to this list, we get {d1 !, d2 !, · · · , dk !}. Now all we need is to get the sum of elements of this list, and this is possible by changing the head from List to Plus by Plus @@. Once this is done, we compare the left-hand side of == with the right-hand side which is the original number #.

Problem 4.5 A number is perfect if it is equal to the sum of its proper divisors, e.g., 6 = 1 + 2 + 3 but 18 = 1 + 2 + 3 + 6 + 9. Write a program to ﬁnd all the perfect numbers up to 10000. (Hint, have a look at the command Divisors.) =⇒ Solution. Here is a step-by-step approach to the solution. Divisors[6] {1,2,3,6} Most[Divisors[6]] {1,2,3} Apply[Plus,Most[Divisors[6]]] 6 Select[Range[10000], # == Apply[Plus, Most[Divisors[#]]] &] {6, 28, 496, 8128}

The numbers 6, 28 and 496 were already known as being perfect numbers 2000 years before Christ. A glance at the list shows that all the perfect num-

52

4. Changing heads!

bers we have found are even. It is still unknown whether there is an odd perfect number. Probably this is the oldest unsolved question in mathematics!

Problem 4.6 Among the ﬁrst 100000 numbers, what is the largest number n which is divisible √ by all positive integers ≤ n? =⇒ Solution. Select[Range[100000],(Mod[#,LCM @@ Range[Floor[Sqrt[#]]]]==0)&] {1, 2, 3, 4, 6, 8, 12, 24}

Exercise 4.2 Decipher what the following codes do: g[n_] := Times @@ Apply[Plus , Inner[List, x^Range[n], 1/x^Range[n], List], 1] t[n_] := Times @@ Apply[Plus, Thread[List[x^Range[n], 1/x^Range[n]]], 1]

Exercise 4.3 Find all the numbers up to 100 which have the following property: if n = pk11 · · · pkt t is the prime decomposition of n then n = k1 × p1 + k2 × p 2 + · · · + kt × p t .

Exercise 4.4 Consider all numbers of the form 3n2 + n + 1, where 0 < n < 10000. How small can the sum of the digits of such a number be? (Hint, see Min.)

Exercise 4.5 Decipher what the following codes do: Power @@ (x + y) Plus @@@ (x^y + y^z)

4. Changing heads!

53

Exercise 4.6 Show that the sum of all the divisors of the number 608655567023837898967037173424316962265783077 3351885970528324860512791691264 is a perfect number. (This number is the only known sublime number besides 12. A sublime number is a positive integer which has a perfect number of positive divisors (including itself), and whose positive divisors add up to another perfect number. See Problem 4.5 for perfect numbers and MathWorld [7] for sublime numbers.)

Exercise 4.7 A weird number is a number such that the sum of the proper divisors (divisors including 1 but not itself) of the number is greater than the number, but no subset of those divisors sums to the number itself. Find all the weird numbers up to 10000.

5 A bit of logic and set theory

This chapter focuses on statements which take values of True or False, i.e., Boolean statements. We can combine statements of this type to make decisions based on their values. We introduce the If statement available in Mathematica.

5.1 Being logical In mathematical logic, statements can have a value of True, False or undeﬁned.1 These are called Boolean expressions. This helps us to “make a decision” and write programs based on the value of a statement (I am thinking of the classical If-Then statement; If something is True, Then do this, Otherwise do that). We have seen == which compares the left-hand side with the righthand side. Studying the following examples carefully will tell us how Wolfram Mathematica® approaches logical statements: 3^2+4^2==5^2 True 3^2+4^2>5^2 False 9Sqrt[10!] < 10Sqrt[9!] False (x-1)(x+1)==x^2-1 (x-1)(x+1)==x^2-1 1

We don’t want to go into detail here mainly because I don’t know the detail!

R. Hazrat, Mathematica ®: A Problem-Centered Approach, Springer Undergraduate Mathematics Series, DOI 10.1007/978-1-84996-251-3 5, © Springer-Verlag London Limited 2010

5.1 Being logical

55

Simplify[%] True x==5 x==5 {1,2}=={2,1} False {a,b}=={b,a} {a,b}=={b,a}

As one notices, Mathematica echoes back the expressions that she can’t evaluate (e.g., x==5). Among them {a,b}=={b,a}, although one expects to get False as lists respect order. This is because Mathematica does not know about the values of a and b, and in case a and b are the same then {a,b}=={b,a} is True, and False otherwise. If you want Mathematica to judge from the face value, then use ===, x==5 False {a,b}=={b,a} False ?=== ihs===rhs yields True if the expression lhs is identical to rhs and yields False otherwise. >>

One can combine logical statements with usual operations And, Or, Not, or the equivalent &&, ||, ! as the following examples show: 2 > 3 && 3 > 2 False And[2 > 3, 3 > 2] False 1 < 2 < 3 True 2 > 3 || 2 < 3 True Or[2 > 3, 3 > 2] True 3^2+4^2>= 5^2 True

In general, for two statements A and B, the statement A&&B is false if one of A or B is false and A || B is true if one of them is true. In order to produce all possible combinations of true and false cases, we use the command Outer

56

5. A bit of logic and set theory

as the following example shows (we will look at the command Outer in more detail in Section 7.5). Outer[f, {a, b}, {x, y}] {{f[a, x], f[a, y]}, {f[b, x], f[b, y]}}

Thus, if in the above we replace f by And or Or we will get all the possible combinations of True and False. Outer[And, {True, False}, {True, False}] {{True, False},{False,False}} Outer[Or, {True, False}, {True, False}] {{True, True}, {True,False}}

In Mathematica, for a variable, one can specify certain domains. This means that the variable takes its values from a speciﬁc type of data. The domains available are Algebraics, Booleans, Complexes, Integers, Primes, Rationals and Reals. One of the fundamental theorems in number theory is to show that π is not a rational number, i.e., is not of the form m/n, where m and n are integers. Look at the following examples: Pi ∈ Rationals False Sqrt[7] ∈ Integers False

Problem 5.1

√ √ √ Show that 1 + 3 + 5 + 7 is an algebraic number (i.e., it is a solution of a polynomial equation with integer coeﬃcients). =⇒ Solution. Plus @@ Sqrt[Range[1, 7, 2]] ∈ Algebraics True √ √ √ The last example shows that 1 + 3 + 5 + 7 is an algebraic number (i.e. it is a solution of a polynomial equation with integer coeﬃcients). One can use membership (∈) to write neat solutions to several problems.

Problem 5.2 Does the formula (n!)2 + 1 give prime numbers for n = 1 to 6?

5.1 Being logical

57

=⇒ Solution. (#!^2 + 1) & /@ Range[6] ∈ Primes False Here we ﬁrst apply the anonymous function (#!^2 + 1) which is the formula (n!)2 + 1 to the list containing 1 to 6. Then we ask Mathematica if the elements of this list belong to the domain Primes. The answer is False. The following code shows that the above formula does not produce a prime number for n = 6: PrimeQ /@ ((#!^2 + 1) & /@ Range[1, 6]) {True, True, True, True, True, False}

One should be careful that Mathematica cannot(yet) perform miracles. For √ √ 3 3 example, one can actually prove that 2 + 5 + 2 − 5 is an integer, but (2 + 5^(1/2))^(1/3) + (2 - 5^(1/2))^(1/3) ∈ Integers False FullSimplify[(2 + 5^(1/2))^(1/3) + (2 - 5^(1/2))^(1/3) ∈ Integers] False Mathematica provides the logical quantiﬁers ∀, ∃ and ⇒ with ForAll, Exist and Implies commands. But these seem not to be that powerful yet. For example, one cannot prove Fermat’s little theorem which says 2p−1 ≡ 1(mod p) where p > 2 is a prime number with them! ForAll[p, p ∈ Primes, Mod[2^(p - 1), p] == 1] Or even an easy fact that the product of four consecutive numbers plus one is a squared number. Implies[n ∈ Integers && n > 0, Sqrt[n(n + 1)(n + 2)(n + 3) + 1] ∈ Integers] In both cases Mathematica gives back the same expression, indicating she cannot decide on them.

58

5. A bit of logic and set theory

5.2 Handling sets Now it has been agreed among mathematicians that any mathematics starts by considering sets, i.e., collections of objects.2 As we mentioned, the diﬀerence between mathematical sets and lists in Mathematica is that lists respect order and repetition, which is to say one can have several copies of one object in a list (see Chapter 3). Sets are not sensitive about repeated objects, e.g., the set {a, b} is the same as the set {a, b, b, a}. There is no concept of sets in Mathematica and if necessary one considers a list as a set. If one wants to get rid of duplications in a list, one can use Union[{a,b,b,a}] {a,b}

Considering two sets, the natural operations between them are union and intersection. Mathematica provides Union to collect all elements from diﬀerent lists in one list (after removing all the duplications) and Intersection for collecting common elements (again discarding repeated elements). The following examples show how these commands work. u= {1, 2, 3, 4, 5, 2, 4, 7, 4}; a = {1, 4, 7, 3}; b = {5, 4, 3, 2}; Union[u] {1, 2, 3, 4, 5, 7} a ∪ b {1, 2, 3, 4, 5, 7} a ∩ b {3,4} Complement[u, a] {2, 5} ?Complement Complement[eall, e1, e2, ... ] gives the elements in eall which are not in any of the ei. 2

To be precise, one ﬁrst considers classes.

5.2 Handling sets

59

Complement[u, a ∩ b] == Complement[u, a] ∪ Complement[u, b] True The ﬁrst example shows Union[list] will get rid of repetition in a list. The command Complement[u,a] will give the elements of u which are not in a. From the example one can see that a ∩ b is acceptable in Mathematica and is a shorthand for Intersection[a,b]. In the last example we checked a theorem of set theory that (A ∩ B)c = Ac ∪ B c where c stands for complement.

Example 5.3 The following trick will be used later (in Chapter 7, inside a loop) to collect data. A={} A=A {x} {x} A=A {y} {x,y} A=A {z} {x,y,z} This is similar to the traditional trick sum=sum+i. Each time sum=sum+i is performed, i will be added to sum and this result will be the new value of sum. There are other ways to add an element to a list. Append[{a,b,c},d] {a,b,c,d} A={};A=Append[A,x] {x} A=Append[A,y] {x,y} A=Append[A,z] {x,y,z}

Also note that the command AppendTo[s, elem] is equivalent to s = Append[s, elem].

Problem 5.4 How many positive integers n are there such that n is a divisor of at least one of the numbers 1040 and 2030 ?

60

5. A bit of logic and set theory

=⇒ Solution. Recall that Divisors[n] will produce a list of all the positive integers which divide n. So all we need to do is to get the collection of divisors of 1040 and 2030 . This can be done with Union. Note that we are not going to print out all the divisors as this is a long list. Length[Union[Divisors[10^40], Divisors[20^30]]] 2301

Problem 5.5 Find out how many common words there are between the English language with French, Dutch and German respectively. =⇒ Solution. Recall that DictionaryLookup contains all the words in 26 languages (see Problem 3.16). All we have to do is to intersect the list of words in the English language with the list of French words and so on as follows: Length[Intersection[DictionaryLookup[{"English", All}], DictionaryLookup[{"French", All}]]] 6897 Length[Intersection[DictionaryLookup[{"English", All}], DictionaryLookup[{"Dutch", All}]]] 6166 Length[ Intersection[DictionaryLookup[{"English", All}], DictionaryLookup[{"German", All}]]] 1286

♣

TIPS

– The command Tally[list] tallies the elements in list, listing all distinct elements together with their multiplicities. – The command Join[list1,list2] concatenates lists or other expressions that share the same head.

5.3 Decision making, If and Which

61

5.3 Decision making, If and Which The statement If[stat,this,that] where stat is a Boolean expression, i.e., has the value of True or False, will execute this if the stat value is True and that otherwise. That means, in either case, one of the statements this or that will be performed (but not both). So this gives an ability to make a decision about which part one wants to perform. Here is an example: If[12^13+13>13^12+12, Print["12^13+13>13^12+12"], Print["12^13+1313^12+12"], 12^13 + 13 < 13^12 + 12, Print["12^13+13

Our partners will collect data and use cookies for ad personalization and measurement. Learn how we and our ad partner Google, collect and use data. Agree & close