Chapter link: https://sarabander.github.io/sicp/html/1_002e2.xhtml
(* n (factorial (- n 1)))
is
a linear recursive process.(fact-iter (* n count) (- count 1))
is linear recursive.Again, a footnote is critical here:
These statements mask a great deal of oversimplification. For instance, if we count process steps as “machine operations” we are making the assumption that the number of machine operations needed to perform, say, a multiplication is independent of the size of the numbers to be multiplied, which is false if the numbers are sufficiently large. Similar remarks hold for the estimates of space. Like the design and description of a process, the analysis of a process can be carried out at various levels of abstraction.
So, you know, we aren't aiming for perfectly formal proofs that algorithms are of a certain order of growth, rather we want to come up with heuristic arguments of our own and understand the heuristic arguments of others, using models of computation.
fast-expt
is such a cool algorithm
(define (even? n)
(= (remainder n 2) 0))
(define (fast-expt b n)
(cond ((= n 0)
1)
((even? n)
(square (fast-expt b (/ n 2))))
(else
(* b (fast-expt b (- n 1))))))
Okay... in C this can be implemented easily as well as a recursive algorithm...
int fastexpt(int b, int n){
if(n==0)
return 1;
if(n%2==0){
int k=fastexpt(b,n/2);
return k*k;
}
return b*fastexpt(b,n-1);
}
The tricky thing is exercise 1.16 turning this into an iterative process.
Euclidean algorithm, also known as "a banger".
Lamé’s Theorem: If Euclid’s Algorithm requires $k$ steps to compute the $\textrm{GCD}$ of some pair, then the smaller number in the pair must be greater than or equal to the $k^{th}$ Fibonacci number.
[43] This theorem was proved in 1845 by Gabriel Lamé, a French mathematician and engineer known chiefly for his contributions to mathematical physics. To prove the theorem, we consider pairs ${(a_k, b_k)}$, where ${a_k \ge b_k}$, for which Euclid's Algorithm terminates in $k$ steps. The proof is based on the claim that, if ${(a_{k+1}, b_{k+1})} \to {(a_k, b_k)} \to {(a_{k-1}, b_{k-1})}$ are three successive pairs in the reduction process, then we must have $b_{k+1} \ge b_k + b_{k-1}$. To verify the claim, consider that a reduction step is defined by applying the transformation ${a_{k-1} = b_k}$, ${b_{k-1} =}$ remainder of $a_k$ divided by $b_k$. The second equation means that $a_k = {qb_k} + {b_{k-1}}$ for some positive integer $q$. And since $q$ must be at least 1 we have $a_k = {qb_k} + b_{k-1} \ge b_k + b_{k-1}$. But in the previous reduction step we have $b_{k+1} = a_k$. Therefore, $b_{k+1} = a_k \ge b_k + b_{k-1}$. This verifies the claim. Now we can prove the theorem by induction on $k$, the number of steps that the algorithm requires to terminate. The result is true for ${k = 1}$, since this merely requires that $b$ be at least as large as ${\text{Fib}(1) = 1}$. Now, assume that the result is true for all integers less than or equal to $k$ and establish the result for ${k + 1}$. Let ${(a_{k+1}, b_{k+1})} \to {(a_k, b_k)} \to {(a_{k-1}, b_{k-1})}$ be successive pairs in the reduction process. By our induction hypotheses, we have $b_{k-1} \ge {\text{Fib}(k - 1)}$ and $b_k \ge {\text{Fib}(k)}$. Thus, applying the claim we just proved together with the definition of the Fibonacci numbers gives $b_{k+1} \ge b_k + b_{k-1} \ge {\text{Fib}(k)} + {\text{Fib}(k-1)} = {\text{Fib}(k+1)}$, which completes the proof of Lamé's Theorem.
[44] If d is a divisor of n , then so is n/d . But d and n/d cannot both be greater than n⎯⎯√ .
Fermat's Little Theorem. If $n$ is a prime number and $a$ is any positive integer less than $n,$ then $a^n\equiv a\textrm{ mod }n.$
Other algorithms worth implementing might be writing the digits of the decimal form of $1/a$ in base $b$. IIRC period of repetition is equal to $b^n \textrm{ mod }a.$ For example, $10^n\textrm{ mod }7$ takes on six values, so $7^{-1}=0.\overline{142857}$ has period six. Compare that to $10^n\textrm{ mod }11$ which only ever takes on values $1,10$. So $11^{-1}=0.\overline{09}$ has period two.
Another algorithm is the extended Euclidean algorithm, which actually finds the integers such that $ax+by=\textrm{gcd}(a,b).$ I've needed this a few times in Advent of Code problems.
It's tempting to mention continued fractions here, but they get their time in the limelight in the next chapter.
Definition of $O(f(n))$: We write that $R(n)$ is $O(f(n))$ as $n\to \infty$ if $$\exists_{c,N}:\forall_{n\gt N}\;\; R(n)\leq c \cdot f(n)$$ $R(n)=O(f(n))$ is a statement that the growth rate of $f(n)$ times a constant is an upper bound on $R(n)$ for sufficiently large $n.$
Definition of $\Omega(f(n))$: We write that $R(n)$ is $\Omega(f(n))$ as $n\to\infty$ if $$\exists_{c,N}:\forall_{n\gt N}\;\; 0 \leq c \cdot f(n) \leq R(n)$$ $R(n)=\Omega(f(n))$ is a statement that the growth rate of $f(n)$ times a constant is a lower bound on $R(n)$ for sufficiently large $n.$
Once we have those two definitions, something is $\Theta(f(n))$ if it is both $O(f(n))$ and $\Omega(f(n)).$
Definition of $f(n)\gg g(n)$: $f(n)\gg g(n)$ as $n\to a$ if $\lim_{n\to a} \frac{g(n)}{f(n)}=0.$ This definition gives: $$n!\gg 2^n \gg n^3\gg n^2\gg n\log n\gg n\gg \log n \gg 1$$
Definition of $f(n)\sim g(n)$: $f(n)\sim g(n)$ as $n\to a$ if $\lim_{n\to a} \frac{g(n)}{f(n)}=1.$
I have to write this because it's a physicist's bread and butter! For example, Stirling's formula is an asymptotic series, and in quantum field theories Feynman diagrams are used to calculate terms in an asymptotic series.
An asymptotic series is defined as follows. Note that this is an abuse of notation, and that we never actually carry out the infinite sum.
Definition of $\sum_{i=1}^\infty f_i(n) \sim g(n)$: We write $\sum_{i=1}^\infty f_i(n) \sim g(n)$ as $n\to a$ if:
For example, Stirling's formula can be written using this abuse of notation as: $$\log(n!)\sim n\log(n)-n+\frac{1}{2}\log(2\pi n) + \sum_{k=2}^\infty \frac{(-1)^k B_k}{k(k-1)n^{k-1}}$$ The coefficients $B_k$ grow so rapidly that for any fixed $n$, the infinite sum doesn't converge! In fact we really mean to have an asymptotic relation. Another way of writing this is that for any fixed $M$ and as $n\to\infty:$ $$\log(n!)=n\log(n)-n+\frac{1}{2}\log(2\pi n) + \sum_{k=2}^M \frac{(-1)^k B_k}{k(k-1)n^{k-1}} + O\left(\frac{1}{n^M}\right)$$
There's another simpler example of asymptotic series that I really like. Say that you're stuck on a desert island and you have to compute $f(0.1)$, where $f$ is a transcendental function defined by the integral $$f(a)=\int_0^\infty \frac{e^{-a x}}{1+x} \mathrm{d}x.$$ Well, let's use a series expansion and totally ignore issues of convergence and real analysis:
$$\begin{align*} \int_0^\infty \frac{e^{-a x}}{1+x} \mathrm{d}x &= \int_0^\infty e^{-x/a}\sum_{n=0}^\infty (-1)^n x^n \mathrm{d}x \\ &= \sum_{n=0}^\infty \int_0^\infty (-1)^n e^{-x/a} x^n \mathrm{d}x \\ &= \sum_{n=0}^\infty (-1)^n n!a^{n+1} \tag{integrate by parts}\\ \end{align*}$$
This infinite sum has a radius of convergence of zero! However if we plug in $a=0.1$ and keep the first few terms, we get $$f(0.1)\approx 0.1 - 0.01 + 0.002 - 0.0006 + 0.00024=0.09164$$ compared to the true value of $f(0.1)=0.0915633\ldots$
Really what we've proven is that $$\int_0^\infty \frac{e^{-a x}}{1+x} \mathrm{d}x \sim \sum_{n=0}^\infty (-1)^n n!a^{n+1} \tag{as $a\to 0$}$$
TLDR this is the simplest example of coming up with an asymptotic series (as $a\to 0$) for a crazy integral.
I learned this from Bender and Orszag's "Advanced Mathematical Methods for Scientists and Engineers: Asymptotic Methods and Perturbation Theory". The title is a mouthful, and it comes across a lot better in video format. Carl Bender writes a lot of these equations with a smile on his face, because $\sum_{n-1}^\infty (-1)^n n!$ is an absurd sum and is kind of a joke! To a physicist or applied mathematician, what we mean by the infinite sum is really $f(1)$ where $f$ is the transcendental function defined above.
No official meeting today, this was just a problem session.
We spent a lot of time on exercise 1.14.
Read section 1.3
I'll try to do all the practice problems. No worries if someone doesn't do em.
Each of the following two procedures defines a method for adding two positive integers in terms of the procedures inc, which increments its argument by 1, and dec, which decrements its argument by 1.
(define (+ a b)
(if (= a 0)
b
(inc (+ (dec a) b))))
(define (+ a b)
(if (= a 0)
b
(+ (dec a) (inc b))))
Using the substitution model, illustrate the process generated by each procedure in evaluating (+ 4 5)
. Are these processes iterative or recursive?
Using the first definition:
; (define (+ a b)
; (if (= a 0)
; b
; (inc (+ (dec a) b))))
(+ 4 5)
(inc (+ 3 5))
(inc (inc (+ 2 5)))
(inc (inc (inc (+ 1 5))))
(inc (inc (inc (inc (+ 0 5)))))
(inc (inc (inc (inc 5))))
(inc (inc (inc 6)))
(inc (inc 7))
(inc 8)
9
; (define (+ a b)
; (if (= a 0)
; b
; (+ (dec a) (inc b))))
(+ 4 5)
(+ 3 6)
(+ 2 7)
(+ 1 8)
(+ 0 9)
9
The following procedure computes a mathematical function called Ackermann’s function.
(define (A x y)
(cond ((= y 0) 0)
((= x 0) (* 2 y))
((= y 1) 2)
(else (A (- x 1)
(A x (- y 1))))))
What are the values of the following expressions?
(A 1 10)
(A 2 4)
(A 3 3)
Consider the following procedures, where A is the procedure defined above:
(define (f n) (A 0 n))
(define (g n) (A 1 n))
(define (h n) (A 2 n))
(define (k n) (* 5 n n))
Give concise mathematical definitions for the functions computed by the procedures f, g, and h for positive integer values of n. For example, (k n)
computes $5n^2$.
#lang sicp
(define (A x y)
(cond ((= y 0) 0)
((= x 0) (* 2 y))
((= y 1) 2)
(else (A (- x 1)
(A x (- y 1))))))
(A 1 10)
(A 2 4)
(A 3 3)
(A 2 5)
1024 65536 65536 2003529930406846464979072351560255750447825475569751419265016973710894059556311453089506130880933348101038234342907263181822949382118812668869506364761547029165041871916351587966347219442930927982084309104855990570159318959639524863372367203002916969592156108764948889254090805911457037675208500206671563702366126359747144807111774815880914135742720967190151836282560618091458852699826141425030123391108273603843767876449043205960379124490905707560314035076162562476031863793126484703743782954975613770981604614413308692118102485959152380195331030292162800160568670105651646750568038741529463842244845292537361442533614373729088303794601274724958414864915930647252015155693922628180691650796381064132275307267143998158508811292628901134237782705567421080070065283963322155077831214288551675554073345107213112427399562982719769150054883905223804357045848197956393157853510018992000024141963706813559840464039472194016069517690156119726982337890017641517190051133466306898140219383481435426387306539552969691388024158161859561100640362119796101859534802787167200122604642492385111393400464351623867567078745259464670903886547743483217897012764455529409092021959585751622973333576159552394885297579954028471943529913543763705986928913757153740001986394332464890052543106629669165243419174691389632476560289415199775477703138064781342309596190960654591300890188887588084733625956065444888501447335706058817090162108499714529568344061979690565469813631162053579369791403236328496233046421066136200220175787851857409162050489711781820400187282939943446186224328009837323764931814789848119452713007440220765680910376203999203492023906626264491909167985461515778839060397720759279378852241294301017458086862263369284725851403039615558564330385450688652213114813638408384778263790459607186876728509763471271988890680478243230394718650525660978150729861141430305816927924971409161059417185352275887504477592218301158780701975535722241400019548102005661773589781499532325208589753463547007786690406429016763808161740550405117670093673202804549339027992491867306539931640720492238474815280619166900933805732120816350707634351669869625020969023162859350071874190579161241536897514808261904847946571736601005892476655445840838334790544144817684255327207315586349347605137419779525190365032198020108764738368682531025183377533908861426184800374008082238104076468878471647552945326947661700424461063311238021134588694532200116564076327023074292426051582811070387018345324567635625951430032037432740780879056283663406965030844225855967039271869461158513793386475699748568670079823960604393478850861649260304945061743412365828352144806726676841807083754862211408236579802961200027441324438432402331257403545019352428776430880232850855886089962774458164680857875115807014743763867976955049991643998284357290415378143438847303484261903388841494031366139854257635577105335580206622185577060082551288893332226436281984838613239570676191409638533832374343758830859233722284644287996245605476932428998432652677378373173288063210753211238680604674708428051166488709084770291208161104912555598322366244868556651402684641209694982590565519216188104341226838996283071654868525536914850299539675503954938371853405900096187489473992880432496373165753803673586710175783994818471798498246948060532081996066183434012476096639519778021441199752546704080608499344178256285092726523709898651539462193004607364507926212975917698293892367015170992091531567814439791248475706237804600009918293321306880570046591458387208088016887445835557926258465124763087148566313528934166117490617526671492672176128330845273936469244582892571388877839056300482483799839692029222215486145902373478222682521639957440801727144146179559226175083889020074169926238300282286249284182671243405751424188569994272331606998712986882771820617214453142574944015066139463169197629181506579745526236191224848063890033669074365989226349564114665503062965960199720636202603521917776740668777463549375318899587866282125469797102065747232721372918144666659421872003474508942830911535189271114287108376159222380276605327823351661555149369375778466670145717971901227117812780450240026384758788339396817962950690798817121690686929538248529830023476068454114178139110648560236549754227497231007615131870024053910510913817843721791422528587432098524957878034683703337818421444017138688124249984418618129271198533315382567321870421530631197748535214670955334626336610864667332292409879849256691109516143618601548909740241913509623043612196128165950518666022030715613684732364660868905014263913906515063908199378852318365059897299125404479443425166774299659811849233151555272883274028352688442408752811283289980625912673699546247341543333500147231430612750390307397135252069338173843322950701049061867539433130784798015655130384758155685236218010419650255596181934986315913233036096461905990236112681196023441843363334594927631946101716652913823717182394299216272538461776065694542297877071383198817036964588689811863210976900355735884624464835706291453052757101278872027965364479724025405448132748391794128826423835171949197209797145936887537198729130831738033911016128547415377377715951728084111627597186384924222802373441925469991983672192131287035585307966942713416391033882754318613643490100943197409047331014476299861725424423355612237435715825933382804986243892498222780715951762757847109475119033482241412025182688713728193104253478196128440176479531505057110722974314569915223451643121848657575786528197564843508958384722923534559464521215831657751471298708225909292655638836651120681943836904116252668710044560243704200663709001941185557160472044643696932850060046928140507119069261393993902735534545567470314903886022024639948260501762431969305640666366626090207048887438898907498152865444381862917382901051820869936382661868303915273264581286782806601337500096593364625146091723180312930347877421234679118454791311109897794648216922505629399956793483801699157439700537542134485874586856047286751065423341893839099110586465595113646061055156838541217459801807133163612573079611168343863767667307354583494789788316330129240800836356825939157113130978030516441716682518346573675934198084958947940983292500086389778563494693212473426103062713745077286156922596628573857905533240641849018451328284632709269753830867308409142247659474439973348130810986399417379789657010687026734161967196591599588537834822988270125605842365589539690306474965584147981310997157542043256395776070485100881578291408250777738559790129129407309462785944505859412273194812753225152324801503466519048228961406646890305102510916237770448486230229488966711380555607956620732449373374027836767300203011615227008921843515652121379215748206859356920790214502277133099987729459596952817044582181956080965811702798062669891205061560742325686842271306295009864421853470810407128917646906550836129916694778023822502789667843489199409657361704586786242554006942516693979292624714524945408858422726153755260071904336329196375777502176005195800693847635789586878489536872122898557806826518192703632099480155874455575175312736471421295536494084385586615208012115079075068553344489258693283859653013272046970694571546959353658571788894862333292465202735853188533370948455403336565356988172582528918056635488363743793348411845580168331827676834646291995605513470039147876808640322629616641560667508153710646723108461964247537490553744805318226002710216400980584497526023035640038083472053149941172965736785066421400842696497103241919182121213206939769143923368374709228267738708132236680086924703491586840991153098315412063566123187504305467536983230827966457417620806593177265685841681837966106144963432544111706941700222657817358351259821080769101961052229263879745049019254311900620561906577452416191913187533984049343976823310298465893318373015809592522829206820862230332585280119266496314441316442773003237792274712330696417149945532261035475145631290668854345426869788447742981777493710117614651624183616680254815296335308490849943006763654806102940094693750609845588558043970485914449584445079978497045583550685408745163316464118083123079704389849190506587586425810738422420591191941674182490452700288263983057950057341711487031187142834184499153456702915280104485145176055306971441761368582384102787659324662689978418319620312262421177391477208004883578333569204533935953254564897028558589735505751235129536540502842081022785248776603574246366673148680279486052445782673626230852978265057114624846595914210278122788941448163994973881884622768244851622051817076722169863265701654316919742651230041757329904473537672536845792754365412826553581858046840069367718605020070547247548400805530424951854495267247261347318174742180078574693465447136036975884118029408039616746946288540679172138601225419503819704538417268006398820656328792839582708510919958839448297775647152026132871089526163417707151642899487953564854553553148754978134009964854498635824847690590033116961303766127923464323129706628411307427046202032013368350385425360313636763575212604707425311209233402837482949453104727418969287275572027615272268283376741393425652653283068469997597097750005560889932685025049212884068274139881631540456490350775871680074055685724021758685439053228133770707415830756269628316955687424060527726485853050611356384851965918968649596335568216975437621430778665934730450164822432964891270709898076676625671517269062058815549666382573829274182082278960684488222983394816670984039024283514306813767253460126007269262969468672750794346190439996618979611928750519442356402644303271737341591281496056168353988188569484045342311424613559925272330064881627466723523751234311893442118885085079358163848994487544756331689213869675574302737953785262542329024881047181939037220666894702204258836895840939998453560948869946833852579675161882159410981624918741813364726965123980677561947912557957446471427868624053750576104204267149366084980238274680575982591331006919941904651906531171908926077949119217946407355129633864523035673345588033313197080365457184791550432654899559705862888286866606618021882248602144999973122164138170653480175510438406624412822803616648904257377640956326482825258407669045608439490325290526337532316509087681336614242398309530806549661879381949120033919489494065132398816642080088395554942237096734840072642705701165089075196155370186264797456381187856175457113400473810762763014953309735174180655479112660938034311378532532883533352024934365979129341284854970946826329075830193072665337782559314331110963848053940859283988907796210479847919686876539987477095912788727475874439806779824968278272200926449944559380414608770641941810440758269805688038949654616587983904660587645341810289907194293021774519976104495043196841503455514044820928933378657363052830619990077748726922998608279053171691876578860908941817057993404890218441559791092676862796597583952483926734883634745651687016166240642424241228961118010615682342539392180052483454723779219911228595914191877491793823340010078128326506710281781396029120914720100947878752551263372884222353869490067927664511634758101193875319657242121476038284774774571704578610417385747911301908583877890152334343013005282797038580359815182929600305682612091950943737325454171056383887047528950563961029843641360935641632589408137981511693338619797339821670761004607980096016024823096943043806956620123213650140549586250615282588033022908385812478469315720323233601899469437647726721879376826431828382603564520699468630216048874528424363593558622333506235945002890558581611275341783750455936126130852640828051213873177490200249552738734585956405160830583053770732533971552620444705429573538361113677523169972740292941674204423248113875075631319078272188864053374694213842169928862940479635305150560788126366206497231257579019598873041195626227343728900516561111094111745277965482790471250581999077498063821559376885546498822938985408291325129076478386322494781016753491693489288104203015610283386143827378160946341335383578340765314321417150655877547820252454780657301342277470616744241968952613164274104695474621483756288299771804186785084546965619150908695874251184435837306590951460980451247409411373899927822492983367796011015387096129749705566301637307202750734759922943792393824427421186158236161317886392553095117188421298508307238259729144142251579403883011359083331651858234967221259621812507058113759495525022747274674369887131926670769299199084467161228738858457584622726573330753735572823951616964175198675012681745429323738294143824814377139861906716657572945807804820559511881687188075212971832636442155336787751274766940790117057509819575084563565217389544179875074523854455200133572033332379895074393905312918212255259833790909463630202185353848854825062897715616963860712382771725621313460549401770413581731931763370136332252819127547191443450920711848838366818174263342949611870091503049165339464763717766439120798347494627397822171502090670190302469762151278521956142070806461631373236517853976292092025500288962012970141379640038055734949269073535145961208674796547733692958773628635660143767964038430796864138563447801328261284589184898528048048844180821639423974014362903481665458114454366460032490618763039502356402044530748210241366895196644221339200757479128683805175150634662569391937740283512075666260829890491877287833852178522792045771846965855278790447562192663992008409302075673925363735628390829817577902153202106409617373283598494066652141198183810884515459772895164572131897797907491941013148368544639616904607030107596818933741217575988165127000761262789169510406315857637534787420070222051070891257612361658026806815858499852631465878086616800733264676830206391697203064894405628195406190685242003053463156621891327309069687353181641094514288036605995220248248886711554429104721929134248346438705368508648749099178812670565665387191049721820042371492740164460943459845392536706132210616533085662021188968234005752675486101476993688738209584552211571923479686888160853631615862880150395949418529489227074410828207169303387818084936204018255222271010985653444817207470756019245915599431072949578197878590578940052540122867517142511184356437184053563024181225473266093302710397968091064939272722683035410467632591355279683837705019855234621222858410557119921731717969804339317707750755627056047831779844447637560254637033369247114220815519973691371975163241302748712199863404548248524570118553342675264715978310731245663429805221455494156252724028915333354349341217862037007260315279870771872491234494477147909520734761385425485311552773301030342476835865496093722324007154518129732692081058424090557725645803681462234493189708138897143299831347617799679712453782310703739151473878692119187566700319321281896803322696594459286210607438827416919465162267632540665070881071030394178860564893769816734159025925194611823642945652669372203155504700213598846292758012527715422016629954863130324912311029627923723899766416803497141226527931907636326136814145516376656559839788489381733082668779901962886932296597379951931621187215455287394170243669885593888793316744533363119541518404088283815193421234122820030950313341050704760159987985472529190665222479319715440331794836837373220821885773341623856441380700541913530245943913502554531886454796252260251762928374330465102361057583514550739443339610216229675461415781127197001738611494279501411253280621254775810512972088465263158094806633687670147310733540717710876615935856814098212967730759197382973441445256688770855324570888958320993823432102718224114763732791357568615421252849657903335093152776925505845644010552192644505312073756287744998163646332835816140330175813967359427327690448920361880386754955751806890058532927201493923500525845146706982628548257883267398735220457228239290207144822219885587102896991935873074277815159757620764023951243860202032596596250212578349957710085626386118233813318509014686577064010676278617583772772895892746039403930337271873850536912957126715066896688493880885142943609962012966759079225082275313812849851526902931700263136328942095797577959327635531162066753488651317323872438748063513314512644889967589828812925480076425186586490241111127301357197181381602583178506932244007998656635371544088454866393181708395735780799059730839094881804060935959190907473960904410150516321749681412100765719177483767355751000733616922386537429079457803200042337452807566153042929014495780629634138383551783599764708851349004856973697965238695845994595592090709058956891451141412684505462117945026611750166928260250950770778211950432617383223562437601776799362796099368975191394965033358507155418436456852616674243688920371037495328425927131610537834980740739158633817967658425258036737206469351248652238481341663808061505704829059890696451936440018597120425723007316410009916987524260377362177763430621616744884930810929901009517974541564251204822086714586849255132444266777127863728211331536224301091824391243380214046242223349153559516890816288487989988273630445372432174280215755777967021666317047969728172483392841015642274507271779269399929740308072770395013581545142494049026536105825409373114653104943382484379718606937214444600826798002471229489405761853892203425608302697052876621377373594394224114707074072902725461307358541745691419446487624357682397065703184168467540733466346293673983620004041400714054277632480132742202685393698869787607009590048684650626771363070979821006557285101306601010780633743344773073478653881742681230743766066643312775356466578603715192922768440458273283243808212841218776132042460464900801054731426749260826922155637405486241717031027919996942645620955619816454547662045022411449404749349832206807191352767986747813458203859570413466177937228534940031631599544093684089572533438702986717829770373332806801764639502090023941931499115009105276821119510999063166150311585582835582607179410052528583611369961303442790173811787412061288182062023263849861515656451230047792967563618345768105043341769543067538041113928553792529241347339481050532025708728186307291158911335942014761872664291564036371927602306283840650425441742335464549987055318726887926424102147363698625463747159744354943443899730051742525110877357886390946812096673428152585919924857640488055071329814299359911463239919113959926752576359007446572810191805841807342227734721397723218231771716916400108826112549093361186780575722391018186168549108500885272274374212086524852372456248697662245384819298671129452945515497030585919307198497105414181636968976131126744027009648667545934567059936995464500558921628047976365686133316563907395703272034389175415267500915011198856872708848195531676931681272892143031376818016445477367518353497857924276463354162433601125960252109501612264110346083465648235597934274056868849224458745493776752120324703803035491157544831295275891939893680876327685438769557694881422844311998595700727521393176837831770339130423060958999137314684569010422095161967070506420256733873446115655276175992727151877660010238944760539789516945708802728736225121076224091810066700883474737605156285533943565843756271241244457651663064085939507947550920463932245202535463634444791755661725962187199279186575490857852950012840229035061514937310107009446151011613712423761426722541732055959202782129325725947146417224977321316381845326555279604270541871496236585252458648933254145062642337885651464670604298564781968461593663288954299780722542264790400616019751975007460545150060291806638271497016110987951336633771378434416194053121445291855180136575558667615019373029691932076120009255065081583275508499340768797252369987023567931026804136745718956641431852679054717169962990363015545645090044802789055701968328313630718997699153166679208958768572290600915472919636381673596673959975710326015571920237348580521128117458610065152598883843114511894880552129145775699146577530041384717124577965048175856395072895337539755822087777506072339445587895905719156736
For the first function,
(A 0 n)
(cond ((= n 0) 0)
((= 0 0) (* 2 n))
((= n 1) 2)
(else (A (- 0 1)
(A x (- n 1)))))
; simplify assuming n >= 0
(* 2 n)
The second function,
(A 1 n)
(cond ((= n 0) 0)
((= 1 0) (* 2 n))
((= n 1) 2)
(else (A (- 1 1)
(A 1 (- n 1)))))
; (A 1 0) is 0
; (A 1 1) is 2
; (A 1 n) for n>=2 is
; =(* 2 (A 1 (- n 1)))
; so (A 1 2) is (* 2 2)=4
; so (A 1 3) is (* 2 (* 2 2))=8
; so (A 1 n) is 2^n, simple proof by induction.
The third function,
(A 2 n)
(cond ((= n 0) 0)
((= 2 0) (* 2 n))
((= n 1) 2)
(else (A (- 2 1)
(A 2 (- n 1)))))
; (A 2 0) is 0
; (A 2 1) is 2
; (A 2 n) for n>=2 is
; =(A 1 (A 2 (- n 1)))
; =(pow 2 (A 2 (- n 1)))
; eg. (A 2 2) is 2^2
; eg. (A 2 3) is 2^(2^2)
; eg. (A 2 4) is 2^(2^(2^2))
So This can be described with Knuth's up arrow notation, (A 2 n)
is $2\uparrow \uparrow n$.
A function $f$ is defined by the rule that $f(n)=n$ if $n\lt 3$ and $f(n)=f(n-1)+2f(n-2)+3f(n-3)$ if $n\geq 3.$ Write a procedure that computes $f$ by means of a recursive process. Write a procedure that computes $f$ by means of an iterative process.
#lang sicp
(define (f-rec n)
(if (< n 3)
n
(+ (f-rec (- n 1))
(* 2 (f-rec (- n 2)))
(* 3 (f-rec (- n 3))))))
(f-rec 1)
(f-rec 2)
(f-rec 3)
(f-rec 4)
(f-rec 5)
(f-rec 6)
1 2 4 11 25 59
#lang sicp
(define (f-it n)
(f-it-helper 2 1 0 n))
(define (f-it-helper a b c count)
(if (= count 0)
c
(f-it-helper
(+ a (* 2 b) (* 3 c))
a
b
(- count 1))))
(f-it 1)
(f-it 2)
(f-it 3)
(f-it 4)
(f-it 5)
(f-it 6)
1 2 4 11 25 59
Write a recursive procedure that computes elements of Pascal's triangle.
The implementation of the procedure to compute the elements is just the first
function choose
, the rest is just to make it print nicely!
#lang sicp
(define (choose n m)
(if (or (< m 1) (> m (- n 1)))
1
(+ (choose (- n 1) m)
(choose (- n 1) (- m 1)))))
(define (print-n char n)
(display char)
(if (> n 1) (print-n char (- n 1))))
(define (print-binom-seq n m)
(display (choose n m))
(display " ")
(if (< m n) (print-binom-seq n (+ m 1))))
(define (print-pascal-line n middle)
(print-n " " (- middle n))
(print-binom-seq n 0)
(newline))
(define (print-pascal-triangle n)
(define (print-pascal-loop nprime middle)
(print-pascal-line nprime middle)
(if (< nprime n) (print-pascal-loop (+ nprime 1) middle)))
(print-pascal-loop 0 (+ n 1)))
(print-pascal-triangle 6)
1 1 1 1 2 1 1 3 3 1 1 4 6 4 1 1 5 10 10 5 1 1 6 15 20 15 6 1
Prove that $\textrm{Fib}(n)$ is the closest integer to $\varphi^n/\sqrt{5}$ where $\varphi=(1+\sqrt{5})/2.$ Hint: Let $\psi=(1-\sqrt{5})/2.$ Use induction and the definition of the Fibonacci numbers to prove that $\textrm{Fib}(n)=(\varphi^n-\psi^n)/\sqrt{5}.$
First we prove that the formula gives the Fibonacci numbers:
Base Cases. Let's work out the first few base cases by hand to get a feel for it. $$f(0)=0,$$ $$f(1)=\frac{1+\sqrt{5}-1+\sqrt{5}}{2\sqrt{5}}=1,$$ $$f(2)=\frac{1+2\sqrt{5}+5-(1-2\sqrt{5}+5)}{4\sqrt{5}}=\frac{4\sqrt{5}}{4\sqrt{5}}=1,$$
Inductive Step. For the inductive step, we know the specific values of $\varphi$ and $\psi$ are important, so we'll write out the numerators and denominators explicitly.
$$f(n)+f(n-1)=\frac{(1+\sqrt{5})^n-(1-\sqrt{5})^n}{2^n\sqrt{5}}+\frac{(1+\sqrt{5})^{n-1}-(1-\sqrt{5})^{n-1}}{2^{n-1}\sqrt{5}}$$
The first term can be expanded as $$\frac{(1+\sqrt{5})(1+\sqrt{5})^{n-1}-(1-\sqrt{5})(1-\sqrt{5})^{n-1}}{2^n\sqrt{5}}$$ The second as $$\frac{2(1+\sqrt{5})^{n-1}-2(1-\sqrt{5})^{n-1}}{2^{n}\sqrt{5}}$$ Note also that $(1+\sqrt{5})^2=6+2\sqrt{5}=2(3+\sqrt{5})$ and $(1-\sqrt{5})^2=6-2\sqrt{5}=2(3-\sqrt{5}).$ Combining like terms gives:
Next, we know that $\varphi^n/\sqrt{5}$ is always within a radius of $|\psi^n/\sqrt{5}|$ of an integer. But $|\psi^n/\sqrt{5}|$ is the monotonic decreasing sequence ${0.27\ldots, 0.17\ldots, 0.10\ldots,\ldots}$ so that already at $n=1$ the nearest integer to $\varphi^n/\sqrt{5}$ is equal to $f(n)$.
Draw the tree illustrating the process generated by the count-change
procedure of 1.2.2 in making change for 11 cents. What are the orders of growth of the space and number of steps used by this process as the amount to be changed increases?
For the orders of growth, the memory growth is the simplest. The amount of space needed for the algorithm is the same order as the depth of the tree (to see this, I imagine traversing the tree from left to right, adding $0+0+\ldots$. At each point in the computation we have a cut somewhere down the middle of the tree, and this cut can only cut through at most $\Theta(\textrm{depth})$ etges).
For the running time growth, I count the number of nodes in the graph. Denote this as $T(a,n)$ where $a$ is the amount we're adding up to and $n\in{1,2,3,4,5}$ is the number of types of coins we care about. $T(a,1)$ satisfies the recursion relation $T(a,1)=T(a-1,1)+2,$ so that $T(a,1)=2a+1.$
$T(a,2)$ satisfies another recursion relation $T(a,2)=T(a-5,2)+T(a,1)+1,$ so that $T(a,2)=T(a-5,2)+2a+2.$
Since this increases by $2a+2$ in strides of $5$, it is $\Theta(a^2).$
Similarly, $T(a,3)=T(a-10,3)+T(a,2)+1,$ which increases by an amount $\Theta(a^2)$ in strides of ten, so is $\Theta(a^3).$
In general, $T(a,n)=\Theta(a^n)$ and so our algorithm is of time complexity $\Theta(a^5).$
The sine of an angle (specified in radians) can be computed by making use of the approximation $\sin x\approx x$ if x is sufficiently small, and the trigonometric identity $$ \sin x = 3 \sin \frac{x}{3}-4\sin^3\frac{x}{3}$$ to reduce the size of the argument of $\sin.$ (For purposes of this exercise an angle is considered “sufficiently small” if its magnitude is not greater than 0.1 radians.) These ideas are incorporated in the following procedures:
(define (cube x) (* x x x))
(define (p x) (- (* 3 x) (* 4 (cube x))))
(define (sine angle)
(if (not (> (abs angle) 0.1))
angle
(p (sine (/ angle 3.0)))))
1.How many times is the procedure p
applied when (sine 12.15)
is evaluated?
2.What is the order of growth in space and number of steps (as a function of $a$
) used by the process generated by the sine procedure when (sine a)
is evaluated?
#lang sicp
(define (cube x) (* x x x))
(define (p x)
(display "p call")
(newline)
(- (* 3 x) (* 4 (cube x))))
(define (sine angle)
(if (not (> (abs angle) 0.1))
angle
(p (sine (/ angle 3.0)))))
(sine 12.5)
p call p call p call p call p call -0.060813768577286265
Design a procedure that evolves an iterative exponentiation process that uses successive squaring and uses a logarithmic number of steps, as does fast-expt
. (Hint: Using the observation that $(b^{n/2})^2=(b^2)^{n/2},$
keep, along with the exponent $n$
and the base $b,$
an additional state variable $a,$
and define the state transformation in such a way that the product $ab^n$
is unchanged from state to state. At the beginning of the process $a$
is taken to be $1,$ and the answer is given by the value of $a$
at the end of the process. In general, the technique of defining an invariant quantity that remains unchanged from state to state is a powerful way to think about the design of iterative algorithms.)
This one always confuses the hell out of me! The key is that in the iteration function of three arguments, we change all three arguments.
At each step, we either do:
$a b^n = a (b^2)^{n/2}$, meaning we leave $a$ unchanged, we toss out $b$ and only consider $b^2$ from now on, and we consider $n/2$ from now on.
$a b^n = (a b) b^{n-1}$, meaning we change $a$ to $(ab)$, we leave $b$ unchanged, and we only consider $n-1$ from now on.
#lang sicp
(define (powfast b n)
(powfast-iter 1 b n))
(define (powfast-iter a b n)
(cond ((= n 0) a)
((= (remainder n 2) 0) (powfast-iter a (* b b) (/ n 2)))
(else (powfast-iter (* a b) b (- n 1)))))
(powfast 3 1)
(powfast 3 2)
(powfast 3 3)
(powfast 3 4)
(powfast 3 5)
(powfast 3 6)
(powfast 3 7)
(powfast 3 8)
3 9 27 81 243 729 2187 6561
The exponentiation algorithms in this section are based on performing exponentiation by means of repeated multiplication. In a similar way, one can perform integer multiplication by means of repeated addition. The following multiplication procedure (in which it is assumed that our language can only add, not multiply) is analogous to the expt procedure:
(define (* a b)
(if (= b 0)
0
(+ a (* a (- b 1)))))
This algorithm takes a number of steps that is linear in b
. Now suppose we include, together with addition, operations double
, which doubles an integer, and halve
, which divides an (even) integer by 2. Using these, design a multiplication procedure analogous to fast-expt
that uses a logarithmic number of steps.
The $O(n)$ solution relies on rewriting $ab=a+a(b-1).$ So we probably want to do something like...
$ab=(2a)(b/2)$ (if b is even)
$ab=a+a(b-1)$ (if b is odd)
Here we see this computes 232*168
in ten steps.
#lang sicp
(define (double arg) (* arg 2))
(define (halve arg) (/ arg 2))
(define (times a b)
(display "evaluating times")
(newline)
(cond ((= b 0) 0)
((= (remainder b 2) 0) (times (double a) (halve b)))
(else (+ a (times a (- b 1))))))
(times 232 168)
evaluating times evaluating times evaluating times evaluating times evaluating times evaluating times evaluating times evaluating times evaluating times evaluating times evaluating times 38976
And checking that zero is handled correctly:
#lang sicp
(define (double arg) (* arg 2))
(define (halve arg) (/ arg 2))
(define (times a b)
(times-inner 0 a b))
(define (times-inner prod a b)
(cond ((= b 0) prod)
((= (remainder b 2) 0) (times-inner prod (double a) (halve b)))
(else (times-inner (+ a prod) a (- b 1)))))
(times 32 1)
(times 1 32)
(times 32 0)
(times 0 32)
32 32 0 0
Using the results of Exercise 1.16 and Exercise 1.17, devise a procedure that generates an iterative process for multiplying two integers in terms of adding, doubling, and halving and uses a logarithmic number of steps.
We'll work with an invariant $\textrm{prod}+ab.$
If $b$ is even, we take $\textrm{prod}+(2a)(b/2).$ That is, $(\textrm{prod},a,b)\mapsto (\textrm{prod},2a,b/2).$
If $b$ is odd, we take $\textrm{prod}+a+a(b-1).$ That is, $(\textrm{prod},a,b)\mapsto (\textrm{prod}+a,a,b-1).$
#lang sicp
(define (double arg) (* arg 2))
(define (halve arg) (/ arg 2))
(define (times a b)
(times-inner 0 a b))
(define (times-inner prod a b)
(display "evaluating times")
(newline)
(cond ((= b 0) prod)
((= (remainder b 2) 0) (times-inner prod (double a) (halve b)))
(else (times-inner (+ a prod) a (- b 1)))))
(times 232 168)
evaluating times evaluating times evaluating times evaluating times evaluating times evaluating times evaluating times evaluating times evaluating times evaluating times evaluating times 38976
And let's just check that zero is handled correctly:
#lang sicp
(define (double arg) (* arg 2))
(define (halve arg) (/ arg 2))
(define (times a b)
(times-inner 0 a b))
(define (times-inner prod a b)
(cond ((= b 0) prod)
((= (remainder b 2) 0) (times-inner prod (double a) (halve b)))
(else (times-inner (+ a prod) a (- b 1)))))
(times 32 1)
(times 1 32)
(times 32 0)
(times 0 32)
32 32 0 0
There is a clever algorithm for computing the Fibonacci numbers in a logarithmic number of steps. Recall the transformation of the state variables $a$ and $b$
in the fib-iter process of 1.2.2: $a\leftarrow a+b$ and $b\leftarrow a.$
Call this transformation $T$
, and observe that applying $T$
over and over again $n$
times, starting with $1$ and $0,$ produces the pair $\textrm{Fib}(n+1)$ and $\textrm{Fib}(n).$
In other words, the Fibonacci numbers are produced by applying $T^n,$ the $n^{\textrm{th}}$
power of the transformation $T$, starting with the pair $(1, 0).$
Now consider $T$ to be the special case of $p=0$ and $q=1$ in a family of transformations $T_{pq}$ where $T_{pq}$ transforms the pair $(a,b)$ according to
$a\leftarrow bq+aq+ap$ and $b\leftarrow bp+aq.$ Show that if we apply such a transformation
$T_{pq}$ twice, the effect is the same as using a single transformation $T_{p'q'}$ of the same form,
and compute $p'$ and $q'$ in terms of $p$ and $q.$ This gives us an explicit way to
square these transformations, and thus we can compute $T^n$ using successive squaring, as in the
fast-expt
procedure. Put this all together to complete the following procedure, which runs in
a logarithmic number of steps:
(define (fib n)
(fib-iter 1 0 0 1 n))
(define (fib-iter a b p q count)
(cond ((= count 0)
b)
((even? count)
(fib-iter a
b
⟨??⟩ ;compute p'
⟨??⟩ ;compute q'
(/ count 2)))
(else
(fib-iter (+ (* b q)
(* a q)
(* a p))
(+ (* b p)
(* a q))
p
q
(- count 1)))))
The rules are better written in matrix form: $$T_{pq}=\begin{bmatrix} p+q & q \\ q & p \end{bmatrix}$$ $$T_{pq}\begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} bq+aq+ap\\ bp+aq \end{bmatrix}$$
We have: $$T_{pq}^2=\begin{bmatrix} p^2+2pq+2q^2 & 2pq+q^2 \\ 2pq+q^2 & p^2+q^2 \end{bmatrix}=T_{(p^2+q^2)(2pq+q^2)}$$
#lang sicp
(define (fib n)
(fib-iter 1 0 0 1 n))
(define (fib-iter a b p q count)
(cond ((= count 0)
b)
((even? count)
(fib-iter a
b
(+ (* p p) (* q q))
(+ (* 2 p q) (* q q))
(/ count 2)))
(else
(fib-iter (+ (* b q)
(* a q)
(* a p))
(+ (* b p)
(* a q))
p
q
(- count 1)))))
(fib 0)
(fib 1)
(fib 2)
(fib 3)
(fib 4)
(fib 5)
(fib 20)
0 1 1 2 3 5 6765
The process that a procedure generates is of course dependent on the rules used by the interpreter. As an example, consider the iterative gcd
procedure given above. Suppose we were to interpret this procedure using normal-order evaluation, as discussed in 1.1.5. (The normal-order-evaluation rule for if is described in Exercise 1.5.) Using the substitution method (for normal order), illustrate the process generated in evaluating (gcd 206 40)
and indicate the remainder operations that are actually performed. How many remainder operations are actually performed in the normal-order evaluation of (gcd 206 40)
? In the applicative-order evaluation?
Let's consider normal order first, since that's more difficult.
I can think of a really silly way for normal order to work, but probably no interpreter would actually be stupid enough to do it like this. We could take "fully expand then reduce" to its extreme, and leave a whole bunch of unevaluated mods (written using "%" to make things less verbose!)
(gcd 206 40)
(if (= 40 0) 206 (gcd 40 (% 206 40)))
(gcd 40 (% 206 40))
(if (= (% 206 40) 0) 40 (gcd (% 206 40) (% 40 (% 206 40)))
(if (= 6 0) 40 (gcd (% 206 40) (% 40 (% 206 40))) ; first eval
This would be extremely stupid, I mean, finding (% 206 40)
everwhere wouldn't even be memoization. The function definition tells us where to place b
. So instead, we shouldn't have this exponentially growing last argument, the evaluation will look like
(gcd 206 40)
(if (= 40 0) 206 (gcd 40 (% 206 40)))
(gcd 40 (% 206 40))
(if (= (% 206 40) 0) 40 (gcd (% 206 40) (% 40 (% 206 40))))
(if (= 6 0) 40 (gcd 6 (% 40 6))) ; still counts as 1 eval
(gcd 6 (% 40 6))
(if (= (% 40 6) 0) 6 (gcd (% 40 6) (% 6 (% 40 6))))
(if (= 4 0) 6 (gcd 4 (% 6 4))) ; 1 eval
(gcd 4 (% 6 4))
(if (= (% 6 4) 0) 4 (gcd (% 6 4) (% 4 (% 6 4))))
(if (= 2 0) 4 (gcd 2 (% 4 2))) ; 1 eval
(if (= 0 0) 2 (gcd 0 (% 2 0))) ; 1 eval
So we need 4 evals.
Okay, so that is precisely the wrong answer! We do the "extremely stupid" method.
Firt off, applicative order uses four evaluations:
#lang sicp
(define (mod a b)
(display "yep")
(newline)
(remainder a b))
(define (gcd a b)
(if (= b 0)
a
(gcd b (mod a b))))
(gcd 206 40)
yep yep yep yep 2
Next up, let's churn through applicative order and pray I don't make a typo.
(gcd 206 40)
(if (= 40 0) 206 (gcd 40 (% 206 40)))
(gcd 40 (% 206 40))
(if (= (% 206 40) 0) 40 (gcd (% 206 40) (% 40 (% 206 40)))
; first eval:
(if (= 6 0) 40 (gcd (% 206 40) (% 40 (% 206 40)))
(gcd (% 206 40) (% 40 (% 206 40)))
; Next two evals:
(if (= (% 40 (% 206 40)) 0) (% 206 40) (gcd (% 40 (% 206 40)) (% (% 206 40) (% 40 (% 206 40)))))
; Next four evals (if statement evaluates to 2)
(if (= (% (% 206 40) (% 40 (% 206 40))) 0) (% 40 (% 206 40)) (gcd (% (% 206 40) (% 40 (% 206 40))) (% (% 40 (% 206 40)) (% (% 206 40) (% 40 (% 206 40))))))
; Next seven evals (if statement evaluates to 0.)
(if (= (% (% 40 (% 206 40)) (% (% 206 40) (% 40 (% 206 40)))) 0) (% (% 206 40) (% 40 (% 206 40))) (gcd b (% a b)))
; Next 4 evals
(% (% 206 40) (% 40 (% 206 40)))
2
So we have 1+2+4+7+4=18 evals. I also checked this with Mathematica which can have its behavior forced in various ways. This prints out 18:
counter=0;
mod[a_,b_]:=(counter++;Mod[a,b]);
gcd[a_,b_]:=If[b===0,a,gcd[Unevaluated[b],Unevaluated[mod[a,b]]]];
result=gcd[206,40];
counter
Use the smallest-divisor procedure to find the smallest divisor of each of the following numbers: 199, 1999, 19999.
#lang sicp
(define (smallest-divisor n)
(find-divisor n 2))
(define (find-divisor n test-divisor)
(cond ((> (* test-divisor test-divisor) n)
n)
((divides? test-divisor n)
test-divisor)
(else (find-divisor
n
(+ test-divisor 1)))))
(define (divides? a b)
(= (remainder b a) 0))
(smallest-divisor 199)
(smallest-divisor 1999)
(smallest-divisor 19999)
199 1999 7
Most Lisp implementations include a primitive called runtime
that returns an integer that specifies the amount of time the system has been running (measured, for example, in microseconds). The following timed-prime-test
procedure, when called with an integer $n,$
prints $n$
and checks to see if $n$
is prime. If $n$
is prime, the procedure prints three asterisks followed by the amount of time used in performing the test.
(define (timed-prime-test n)
(newline)
(display n)
(start-prime-test n (runtime)))
(define (start-prime-test n start-time)
(if (prime? n)
(report-prime (- (runtime)
start-time))))
(define (report-prime elapsed-time)
(display " *** ")
(display elapsed-time))
Using this procedure, write a procedure search-for-primes
that checks the
primality of consecutive odd integers in a specified range. Use your procedure
to find the three smallest primes larger than 1000; larger than 10,000; larger
than 100,000; larger than 1,000,000. Note the time needed to test each prime.
Since the testing algorithm has order of growth of
$\Theta\left(\sqrt{n}\right),$ you should expect that testing for primes around
10,000 should take about $\sqrt{10}$ times as long as testing for primes around
1000. Do your timing data bear this out? How well do the data for 100,000 and
1,000,000 support the $\Theta\left(\sqrt{n}\right)$ prediction? Is your result
compatible with the notion that programs on your machine run in time
proportional to the number of steps required for the computation?
I had to go to $10^8$ (22ms), $10^9$ (70ms), and $10^{10}$ (220ms) to get times above 1ms!
$22\sqrt{10}\approx 73.8,$ so we're right on the money.
#lang sicp
(define (smallest-divisor n)
(find-divisor n 2))
(define (find-divisor n test-divisor)
(cond ((> (* test-divisor test-divisor) n)
n)
((divides? test-divisor n)
test-divisor)
(else (find-divisor
n
(+ test-divisor 1)))))
(define (divides? a b)
(= (remainder b a) 0))
(define (prime? n)
(= n (smallest-divisor n)))
(define (timed-prime-test n start-time)
(if (prime? n)
(begin
(display n)
(display " : ")
(display (- (runtime) start-time))
(display "ms")
(newline)
#t)
#f))
(define (primes-larger-than n k-primes)
(if (not (= k-primes 0))
(if (timed-prime-test n (runtime))
(primes-larger-than (+ n 1) (- k-primes 1))
(primes-larger-than (+ n 1) k-primes))))
(primes-larger-than 100000000 3)
(primes-larger-than 1000000000 3)
(primes-larger-than 10000000000 3)
100000007 : 26ms 100000037 : 23ms 100000039 : 22ms 1000000007 : 70ms 1000000009 : 71ms 1000000021 : 70ms 10000000019 : 220ms 10000000033 : 222ms 10000000061 : 237ms
The smallest-divisor
procedure shown at the start of this section does lots of
needless testing: After it checks to see if the number is divisible by 2 there
is no point in checking to see if it is divisible by any larger even numbers.
This suggests that the values used for test-divisor
should not be 2, 3, 4, 5,
6, …, but rather 2, 3, 5, 7, 9, …. To implement this change, define a procedure
next that returns 3 if its input is equal to 2 and otherwise returns its input
plus 2. Modify the smallest-divisor
procedure to use (next test-divisor)
instead of (+ test-divisor 1)
. With timed-prime-test
incorporating this
modified version of smallest-divisor
, run the test for each of the 12 primes
found in Exercise 1.22. Since this modification halves the number of test
steps, you should expect it to run about twice as fast. Is this expectation
confirmed? If not, what is the observed ratio of the speeds of the two
algorithms, and how do you explain the fact that it is different from 2?
#lang sicp
(define (smallest-divisor n)
(find-divisor n 2))
(define (find-divisor n test-divisor)
(cond ((> (* test-divisor test-divisor) n)
n)
((divides? test-divisor n)
test-divisor)
(else (find-divisor
n
(nextdiv test-divisor)))))
(define (divides? a b)
(= (remainder b a) 0))
(define (prime? n)
(= n (smallest-divisor n)))
(define (nextdiv test-divisor)
(if (= test-divisor 2)
3
(+ test-divisor 2)))
(define (timed-prime-test n start-time)
(if (prime? n)
(begin
(display n)
(display " : ")
(display (- (runtime) start-time))
(display "ms")
(newline)
#t)
#f))
(define (primes-larger-than n k-primes)
(if (not (= k-primes 0))
(if (timed-prime-test n (runtime))
(primes-larger-than (+ n 1) (- k-primes 1))
(primes-larger-than (+ n 1) k-primes))))
(primes-larger-than 100000000 3)
(primes-larger-than 1000000000 3)
(primes-larger-than 10000000000 3)
100000007 : 15ms 100000037 : 12ms 100000039 : 12ms 1000000007 : 36ms 1000000009 : 36ms 1000000021 : 35ms 10000000019 : 111ms 10000000033 : 114ms 10000000061 : 111ms
It's very closely within a factor of two, yep.
Modify the timed-prime-test
procedure of Exercise 1.22 to use
fast-prime?
(the Fermat method), and test each of the 12 primes you found in
that exercise. Since the Fermat test has $\Theta(\log n)$ growth, how would you expect
the time to test primes near 1,000,000 to compare with the time needed to test
primes near 1000? Do your data bear this out? Can you explain any discrepancy
you find?
#lang sicp
(define (square n) (* n n))
(define (expmod base exp m)
(cond ((= exp 0) 1)
((even? exp)
(remainder
(square (expmod base (/ exp 2) m))
m))
(else
(remainder
(* base (expmod base (- exp 1) m))
m))))
(define (fermat-test n)
(define (try-it a)
(= (expmod a n n) a))
(try-it (+ 1 (random (- n 1)))))
(define (fast-prime? n times)
(cond ((= times 0) true)
((fermat-test n)
(fast-prime? n (- times 1)))
(else false)))
(define (timed-prime-test n start-time)
(if (fast-prime? n 600)
(begin
(display n)
(display " : ")
(display (- (runtime) start-time))
(display "ms")
(newline)
#t)
#f))
(define (primes-larger-than n k-primes)
(if (not (= k-primes 0))
(if (timed-prime-test n (runtime))
(primes-larger-than (+ n 1) (- k-primes 1))
(primes-larger-than (+ n 1) k-primes))))
(primes-larger-than 10 3)
(primes-larger-than 100 3)
(primes-larger-than 1000 3)
(primes-larger-than 10000 3)
(primes-larger-than 100000 3)
(primes-larger-than 1000000 3)
(primes-larger-than 10000000 3)
(primes-larger-than 100000000 3)
11 : 49ms 13 : 46ms 17 : 59ms 101 : 64ms 103 : 66ms 107 : 67ms 1009 : 99ms 1013 : 99ms 1019 : 103ms 10007 : 112ms 10009 : 113ms 10037 : 114ms 100003 : 132ms 100019 : 133ms 100043 : 131ms 1000003 : 151ms 1000033 : 158ms 1000037 : 155ms 10000019 : 176ms 10000079 : 184ms 10000103 : 183ms 100000007 : 208ms 100000037 : 209ms 100000039 : 215ms
The graph of timings looks logarithmic as expected. Note that if we make $n$ too large we'll start to run into big number arithmetic. To get measurable running times, I changed the number of checks to $600$ instead of something more reasonable like $10$ or $20.$
Alyssa P. Hacker complains that we went to a lot of extra work in writing expmod. After all, she says, since we already know how to compute exponentials, we could have simply written
(define (expmod base exp m)
(remainder (fast-expt base exp) m))
Is she correct? Would this procedure serve as well for our fast prime tester? Explain.
We'd run into large number arithmetic, and "slow" is an understatement! The number of digits would grow ruining the time and space complexity of the algorithm.
We have $a^n$ where $a$ is of order $n$. If $n$ is $k$ digits long, then (fast-expt a n)
is of order $kn$ digits long, which is $O(n\log(n))$ memory growth.
Similarly, each multiplication will be more expensive and so we expect this to ruin the time complexity of the algorithm.
Louis Reasoner is having great difficulty doing Exercise 1.24. His
fast-prime?
test seems to run more slowly than his prime?
test. Louis calls
his friend Eva Lu Ator over to help. When they examine Louis’s code, they find
that he has rewritten the expmod procedure to use an explicit multiplication,
rather than calling square:
(define (expmod base exp m)
(cond ((= exp 0) 1)
((even? exp)
(remainder
(* (expmod base (/ exp 2) m)
(expmod base (/ exp 2) m))
m))
(else
(remainder
(* base
(expmod base (- exp 1) m))
m))))
“I don’t see what difference that could make,” says Louis. “I do.” says Eva. “By writing the procedure like that, you have transformed the $\Theta(\log n)$ process into a $\Theta( n)$ process.” Explain.
This is another emphasis of applicative order evaluation. When exp
is even, the recursive expmod
step is calculated twice. If those are even, then it's twice again, and so on. For example if exp
is $2^k$ for some $k,$ then this is a binary tree of depth $k$ with $O(2^k)$ function calls. So for this case, the algorithm is linear $O(n)$ with $n=2^k.$
Demonstrate that the Carmichael numbers listed in Footnote 47 really do fool the Fermat test. That is, write a procedure that takes an integer $n$ and tests whether $a^n$ is congruent to $a$ modulo $n$ for every $a\lt n,$ and try your procedure on the given Carmichael numbers.
#lang sicp
(define (square n) (* n n))
(define (expmod base exp m)
(cond ((= exp 0) 1)
((even? exp)
(remainder
(square (expmod base (/ exp 2) m))
m))
(else
(remainder
(* base (expmod base (- exp 1) m))
m))))
(define (test-all-iter n test bool)
(if (= n test)
bool
(test-all-iter n
(+ test 1)
(and bool (= (expmod test n n) test)))))
(define (carm-test n)
(display "Fermat test of ")
(display n)
(display " : ")
(display (test-all-iter n 2 #t))
(newline))
(display "Testing some non-carmichael non-primes") (newline)
(carm-test 231)
(carm-test 1419)
(carm-test 2001)
(carm-test 8631)
(newline) (display "Testing some carmichael non-primes") (newline)
(carm-test 561)
(carm-test 1105)
(carm-test 1729)
(carm-test 2465)
(carm-test 2821)
(carm-test 6601)
Testing some non-carmichael non-primes Fermat test of 231 : #f Fermat test of 1419 : #f Fermat test of 2001 : #f Fermat test of 8631 : #f Testing some carmichael non-primes Fermat test of 561 : #t Fermat test of 1105 : #t Fermat test of 1729 : #t Fermat test of 2465 : #t Fermat test of 2821 : #t Fermat test of 6601 : #t
One variant of the Fermat test that cannot be fooled is called the Miller-Rabin test (Miller 1976; Rabin 1980). This starts from an alternate form of Fermat’s Little Theorem, which states that if $n$ is a prime number and $a$ is any positive integer less than $n,$ then $a^{n-1}\equiv 1\textrm{ mod }n.$
To test the primality of a number $n$ by the Miller-Rabin test, we
pick a random number $a\lt n$ and find $a^{n-1}\textrm{ mod }n$ using
the expmod
procedure. However, whenever we perform the squaring step in expmod
,
we check to see if we have discovered a “nontrivial square root of 1 modulo $n$
,” that is, a number not equal to 1 or $n−1$ whose square is equal to 1 modulo $n.$
It is possible to prove that if such a nontrivial square root of 1 exists,
then $n$ is not prime. It is also possible to prove that if $n$ is an odd number
that is not prime, then, for at least half the numbers $a\lt n$ , computing $a^{n−1}$
in this way will reveal a nontrivial square root of 1 modulo $n.$ (This is why
the Miller-Rabin test cannot be fooled.) Modify the expmod procedure to signal
if it discovers a nontrivial square root of 1, and use this to implement the
Miller-Rabin test with a procedure analogous to fermat-test
. Check your
procedure by testing various known primes and non-primes. Hint: One convenient
way to make expmod
signal is to have it return 0.
#lang sicp
(define (square n) (* n n))
; In the context of our prime test, base is `a` (the randomly chosen test)
; exp is n (The potential prime), and m is also n.
; Consider the piece of code:
; (remainder
; (square (expmod base (/ exp 2) m))
; m))
; Normally I'd want to do something like
; result = (expmod base (/ exp 2) m)
; if(result*result%m == 1 && result !=1 && result != m-1)
; return 0;
; else
; return result*result%m
(define (expmod-fancy base exp m)
; Kinda silly to evaluate remainder square result twice, but whatever.
(define (even-function result)
(if
(and (= (remainder (square result) m) 1)
(not (= result 1))
(not (= result (- m 1))))
0
; if result is zero, it'll be passed through
(remainder (square result) m)))
(cond ((= exp 0) 1)
((even? exp)
(even-function (expmod-fancy base (/ exp 2) m)))
(else
; here, if expmod-fancy returns 0, it's passed through also.
(remainder
(* base (expmod-fancy base (- exp 1) m))
m))))
(define (fermat-test n)
(define (try-it a)
(= (expmod-fancy a (- n 1) n) 1))
(try-it (+ 1 (random (- n 1)))))
(define (miller-rabin-prime? n times)
(cond ((= times 0) true)
((fermat-test n)
(miller-rabin-prime? n (- times 1)))
(else false)))
(define (mrtest n)
(display "Miller-Rabin test of ")
(display n)
(display " : ")
(display (miller-rabin-prime? n 1))
(newline))
(display "Testing some non-carmichael non-primes") (newline)
(mrtest 231)
(mrtest 1419)
(mrtest 2001)
(mrtest 8631)
(newline) (display "Testing some carmichael non-primes") (newline)
(mrtest 561)
(mrtest 1105)
(mrtest 1729)
(mrtest 2465)
(mrtest 2821)
(mrtest 6601)
(newline) (display "Testing some primes") (newline)
(mrtest 1297)
(mrtest 1579)
(mrtest 2671)
(mrtest 6199)
(mrtest 7013)
Testing some non-carmichael non-primes Miller-Rabin test of 231 : #f Miller-Rabin test of 1419 : #f Miller-Rabin test of 2001 : #f Miller-Rabin test of 8631 : #f Testing some carmichael non-primes Miller-Rabin test of 561 : #f Miller-Rabin test of 1105 : #f Miller-Rabin test of 1729 : #f Miller-Rabin test of 2465 : #f Miller-Rabin test of 2821 : #f Miller-Rabin test of 6601 : #f Testing some primes Miller-Rabin test of 1297 : #t Miller-Rabin test of 1579 : #t Miller-Rabin test of 2671 : #t Miller-Rabin test of 6199 : #t Miller-Rabin test of 7013 : #t