f20-5

902 Chapter20. Less-NumericalAlgorithmsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).20.5 Arithmetic Coding Wesawintheprevioussectionthataperfect(entropy-bounded)codingscheme would use Li=−log2pibits to encode character i(in the range 1≤i≤Nch), ifpiis its probability of occurrence. Huffman coding gives a way of rounding the Li’s to close integer values and constructing a code with those lengths. Arithmetic coding[1], which we now discuss, actually does manage to encode characters using noninteger numbers of bits! It also provides a convenient way to output the resultnot as a stream of bits, but as a stream of symbols in any desired radix. This latter property is particularly useful if you want, e.g., to convert data from bytes (radix 256) to printable ASCII characters (radix 94), or to case-independentalphanumericsequences containing only A-Z and 0-9 (radix 36). In arithmetic coding, an input message of any length is represented as a real number Rin the range 0≤R< 1. The longer the message, the more precision requiredof R. This is bestillustratedbyanexample,so letus returnto thefictitious language,Vowellish,oftheprevioussection. RecallthatVowellishhasa5characteralphabet (A, E, I, O, U), with occurrence probabilities 0.12, 0.42, 0.09, 0.30, and 0.07,respectively. Figure20.5.1showshowamessagebeginning“IOU”isencoded: The interval [0,1)is divided into segments corresponding to the 5 alphabetical characters; the length of a segment is the corresponding p i. We see that the first messagecharacter,“I”,narrowstherangeof Rto0.37≤R< 0.46. This intervalis nowsubdividedintofivesubintervals,againwithlengthsproportionaltothe pi’s. The second message character, “O”, narrows the range of Rto0.3763≤R< 0.4033. The“U” characterfurthernarrowsthe rangeto 0.37630 ≤R< 0.37819.Anyvalue ofRin this range can be sent as encoding “IOU”. In particular, the binary fraction .011000001 is in this range, so “IOU” can be sent in 9 bits. (Huffman coding took 10 bits for this example, see §20.4.) Ofcoursethereistheproblemofknowingwhentostopdecoding. Thefraction .011000001 representsnotsimply“IOU,”but“IOU ...,”wheretheellipsesrepresent an infinite string of successor characters. To resolve this ambiguity, arithmetic coding generally assumes the existence of a special Nch+1th character, EOM (end of message), which occurs only once at the end of the input. Since EOMhas a low probability of occurrence, it gets allocated only a very tiny piece of the number line. In the above example, we gave Ras a binary fraction. We could just as well have output it in any other radix, e.g., base 94 or base 36, whatever is convenient for the anticipated storage or communication channel. You might wonder how one deals with the seemingly incredible precision requiredof Rforalongmessage. Theansweristhat Risneveractuallyrepresented all at once. At any give stage we have upper and lower bounds for Rrepresented as a finite number of digits in the output radix. As digits of the upper and lower boundsbecomeidentical,we can left-shift them away andbringin new digits at the low-significance end. The routines below have a parameter NWKfor the number of working digits to keep around. This must be large enough to make the chance of an accidental degeneracy vanishingly small. (The routines signal if a degeneracy ever occurs.) Since the process of discarding old digits and bringingin new ones is performedidentically onencodingand decoding,everythingstays synchronized. 20.5ArithmeticCoding 903Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).1.0 0.9 0.8 0.7 0.60.5 0.4 0.3 0.2 0.1 0.0A E I O UA E I O UA E I O UA E I O U0.46 0.42 0.41 0.370.385 0.3800.4033 0.37630.37819 0.376300.3780 0.37640.44 0.43 0.3900.3950.400 0.37660.37680.37720.37740.3778 0.37760.45 0.40 0.39 0.380.3770 Figure 20.5.1. Arithmetic coding of the message “IOU...”in thefictitious language Vowellish. Successive characters givesuccessively finersubdivisions oftheinitial interval between 0and1. The final value can be output as the digits of a fraction in any desired radix. Note how the subinterval allocated to a character is proportional to its probability of occurrence. Theroutine arcmakconstructsthecumulativefrequencydistributiontableused to partition the interval at each stage. In the principal routine arcode, when an intervalofsize jdifis to bepartitionedin theproportionsofsome ntosome ntot, say,thenwemustcompute (n*jdif)/ntot . Withintegerarithmetic,thenumerator is likely to over flow; and, unfortunately, an expression like jdif/(ntot/n) is not equivalent. In the implementation below, we resort to double precision floating arithmetic for this calculation. Not only is this inef ficient, but different roundoff errorscan(albeitveryrarely)makedifferentmachinesencodedifferently,thoughanyone type of machine will decode exactly what it encoded, since identical roundoff errorsoccurin the two processes. For serious use, one needsto replacethis floating calculation with an integer computation in a double register (not available to the FORTRAN programmer). The internally set variable minint, which is the minimum allowed number of discrete steps between the upper and lower bounds, determines when new low- significancedigits areadded. minintmustbelargeenoughtoprovideresolutionof all the input characters. That is, we must have p i×minint >1for all i. A value of100Nch,o r 1.1/minpi, whicheveris larger,is generallyadequate. However, for safety, the routine below takes minintto be as large as possible, with the product minint*nradd just smaller than over flow. This results in some time inef ficiency, and in a few unnecessary characters being output at the end of a message. You can 904 Chapter20. Less-NumericalAlgorithmsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).decrease minintif you want to live closer to the edge. Afinal safety featurein arcmakis its refusalto believezerovaluesin the table nfreq;a 0is treated as if it were a 1. If this were not done, the occurrence in a message ofa single characterwhose nfreqentryis zerowouldresult in scrambling the entire rest of the message. If you want to live dangerously, with a very slightlymore efficient coding, you can delete the max( ,1) operation. SUBROUTINE arcmak(nfreq,nchh,nradd) INTEGER nchh,nradd,nfreq(nchh),MC,NWK,MAXINTPARAMETER (MC=512,NWK=20,MAXINT=2147483647) Given a table nfreq(1:nchh) of the frequency of occurrence of nchhsymbols, and given a desired output radix nradd, initialize the cumulative frequency table and other variables for arithmetic compression.Parameters: MCis largest anticipated value of nchh;NWKis the number of working digits (see text); MAXINT is a large positive integer that does not overflow. INTEGER j,jdif,minint,nc,nch,nrad,ncum, * ncumfq(MC+2),ilob(NWK),iupb(NWK) COMMON /arccom/ ncumfq,iupb,ilob,nch,nrad,minint,jdif,nc,ncum SAVE /arccom/if(nchh.gt.MC)pause ’MC too small in arcmak’if(nradd.gt.256)pause ’nradd may not exceed 256 in arcmak’ minint=MAXINT/nradd nch=nchhnrad=nradd ncumfq(1)=0 do 11j=2,nch+1 ncumfq(j)=ncumfq(j-1)+max(nfreq(j-1),1) enddo 11 ncumfq(nch+2)=ncumfq(nch+1)+1ncum=ncumfq(nch+2)return END Individualcharactersinamessagearecodedordecodedbytheroutine arcode, which in turn uses the utility arcsum. SUBROUTINE arcode(ich,code,lcode,lcd,isign) INTEGER ich,isign,lcd,lcode,MC,NWK CHARACTER*1 code(lcode) PARAMETER (MC=512,NWK=20) C USES arcsum Compress ( isign =1) or decompress ( isign =−1) the single character ichinto or out of the character array code(1:lcode) , starting with byte code(lcd) and (if necessary) incrementing lcdso that, on return, lcdpoints to the first unused byte in code.N o t e that this routine saves the result of previous calls until a new byte of code is produced, and only then increments lcd. An initializing call with isign=0 is required for each different array code. The routine arcmak must have previously been called to initialize the common block /arccom/ .Ac a l lw i t h ich=nch (as set in arcmak ) has the reserved meaning “end of message.” INTEGER ihi,j,ja,jdif,jh,jl,k,m,minint,nc,nch,nrad,ilob(NWK), * iupb(NWK),ncumfq(MC+2),ncum,JTRY COMMON /arccom/ ncumfq,iupb,ilob,nch,nrad,minint,jdif,nc,ncum SAVE /arccom/ The following statement function is used to calculate (k*j)/m without overflow. Program efficiency can be improved by substituting an assembly language routine that does integermultiply to a double register. JTRY(j,k,m)=int((dble(k)*dble(j))/dble(m)) if (isign.eq.0) then Initialize enough digits of the upper and lower bounds. jdif=nrad-1 do 11j=NWK,1,-1 20.5ArithmeticCoding 905Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).iupb(j)=nrad-1 ilob(j)=0 nc=j if(jdif.gt.minint)return Initialization complete. jdif=(jdif+1)*nrad-1 enddo 11 pause ’NWK too small in arcode’ else if (isign.gt.0) then If encoding, check for valid input character. if(ich.gt.nch.or.ich.lt.0)pause ’bad ich in arcode’ else If decoding, locate the character ichby bisection. ja=ichar(code(lcd))-ilob(nc)do 12j=nc+1,NWK ja=ja*nrad+(ichar(code(j+lcd-nc))-ilob(j)) enddo 12 ich=0 ihi=nch+1 1 if(ihi-ich.gt.1) then m=(ich+ihi)/2if (ja.ge.JTRY(jdif,ncumfq(m+1),ncum)) then ich=m else ihi=m endif goto 1endifif(ich.eq.nch)return Detected end of message. endif Following code is common for encoding and decoding. Convert character ichto a new subrange [ilob,iupb) . jh=JTRY(jdif,ncumfq(ich+2),ncum) jl=JTRY(jdif,ncumfq(ich+1),ncum)jdif=jh-jlcall arcsum(ilob,iupb,jh,NWK,nrad,nc) call arcsum(ilob,ilob,jl,NWK,nrad,nc) How many leading digits to output (if encoding) or skip over? do 13j=nc,NWK if(ich.ne.nch.and.iupb(j).ne.ilob(j))goto 2if(lcd.gt.lcode)pause ’lcode too small in arcode’ if(isign.gt.0) code(lcd)=char(ilob(j)) lcd=lcd+1 enddo 13 return Ran out of message. Did someone forget to encode a terminating ncd? 2 nc=j j=0 How many digits to shift? 3 if (jdif.lt.minint) then j=j+1 jdif=jdif*nrad goto 3 endif if (nc-j.lt.1) pause ’NWK too small in arcode’if(j.ne.0)then Shift them. do 14k=nc,NWK iupb(k-j)=iupb(k) ilob(k-j)=ilob(k) enddo 14 endifnc=nc-jdo 15k=NWK-j+1,NWK iupb(k)=0 ilob(k)=0 enddo 15 endifreturn Normal return. END 906 Chapter20. Less-NumericalAlgorithmsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).SUBROUTINE arcsum(iin,iout,ja,nwk,nrad,nc) INTEGER ja,nc,nrad,nwk,iin(*),iout(*) Used by arcode . Add the integer jato the radix nradmultiple-precision integer iin(nc..nwk) . Return the result in iout(nc..nwk) . INTEGER j,jtmp,karrykarry=0do 11j=nwk,nc+1,-1 jtmp=ja ja=ja/nradiout(j)=iin(j)+(jtmp-ja*nrad)+karryif (iout(j).ge.nrad) then iout(j)=iout(j)-nrad karry=1 else karry=0 endif enddo 11 iout(nc)=iin(nc)+ja+karry return END If radix-changing, rather than compression, is your primary aim (for example to convert an arbitrary file into printable characters) then you are of course free to set all the components of nfreqequal, say, to 1. CITED REFERENCES AND FURTHER READING: Bell,T.C.,Cleary,J.G.,andWitten,I.H.1990, TextCompression (EnglewoodCliffs,NJ:Prentice- Hall). Nelson, M. 1991, The Data Compression Book (Redwood City, CA: M&T Books). Witten, I.H., Neal, R.M., and Cleary, J.G. 1987, Communications of the ACM , vol. 30, pp. 520– 540. [1] 20.6 Arithmetic at Arbitrary Precision Let’s compute the number πto a couple of thousand decimal places. In doing so, we’ll learn some things about multiple precision arithmetic on computers and meet quite an unusual application of the fast Fourier transform (FFT). We ’ll also develop a set of routines that you can use for other calculations at any desired level of arithmetic precision. To start with, we need an analytic algorithm for π. Useful algorithms are quadratically convergent, i.e., they double the number of signi ficant digits at each iteration. Quadratically convergent algorithms for πare based on the AGM (arithmeticgeometric mean) method,which also finds applicationto the calculation of elliptic integrals (cf. §6.11)and in advancedimplementationsof the ADI method for elliptic partial differential equations ( §19.5). Borwein and Borwein [1]treat this subject, which is beyond our scope here. One of their algorithms for πstarts with the initializations X0=√ 2 π0=2+√ 2 Y0=4√ 2(20.6.1 )