Movatterモバイル変換

{-# LANGUAGE DeriveDataTypeable #-}{- |Module      :  Fractal.RUFF.Mandelbrot.AddressCopyright   :  (c) Claude Heiland-Allen 2010-2011License     :  BSD3Maintainer  :  claudiusmaximus@goto10.orgStability   :  unstablePortability :  portableExternal angles give rise to kneading sequences under the angle doublingmap.  Internal addresses encode kneading sequences in human-readable form,when extended to angled internal addresses they distinguish hyperboliccomponents in a concise and meaningful way.The algorithms are mostly based on Dierk Schleicher's paper/Internal Addresses Of The Mandelbrot Set And Galois Groups Of Polynomials (version of February 5, 2008)/<http://arxiv.org/abs/math/9411238v2>.-}moduleFractal.RUFF.Mandelbrot.Address(Angle,double,wrap,prettyAngle,prettyAngles,Knead(..),kneadChar,Kneading(..),prettyKneading,kneading,period,unwrap,associated,upper,lower,InternalAddress(..),prettyInternalAddress,internalAddress,internalFromList,internalToList,AngledInternalAddress(..),prettyAngledInternalAddress,angledInternalAddress,angledFromList,angledToList,externalAngles,stripAngles,splitAddress,joinAddress,addressPeriod,parseAngle,parseAngles,parseKnead,parseKneading,parseInternalAddress,parseAngledInternalAddress)whereimportData.Data(Data())importData.Typeable(Typeable())importControl.Monad(guard)importControl.Monad.Identity(Identity())importData.Char(digitToInt)importData.Bits(testBit)importData.List(genericDrop,genericIndex,genericLength,genericReplicate,genericSplitAt,genericTake,foldl')importData.Maybe(isJust,listToMaybe)importData.Ratio((%),numerator,denominator)importText.Parsec(ParsecT(),choice,digit,eof,many,many1,runP,sepEndBy,string,try)-- | Angle as a fraction of a turn, usually in [0, 1).typeAngle=Rational-- | Convert to human readable form.prettyAngle::Angle->StringprettyAnglea=show(numeratora)++" / "++show(denominatora)-- | Convert to human readable form.prettyAngles::[Angle]->StringprettyAngles[]=""prettyAngles[a]=show(numeratora)++"/"++show(denominatora)prettyAngles(a:as)=show(numeratora)++"/"++show(denominatora)++" "++prettyAnglesas-- | Wrap an angle into [0, 1).wrap::Angle->Anglewrapa|f<0=1+f|otherwise=fwhere(_,f)=properFractiona::(Integer,Angle)-- | Angle doubling map.double::Angle->Angledoublea=wrap(2*a)-- | Binary representation of a (pre-)periodic angle.typeBinAngle=([Bool],[Bool])-- | Convert an angle from binary representation.unbinary::BinAngle->Angleunbinary(pre,per)|n==0=bitspre%(2^m)|otherwise=(bitspre%(2^m))+(bitsper%(2^m*(2^n-1)))wherem=lengthpren=lengthper-- | Convert a list of bits to an integer.bits::[Bool]->Integerbits=foldl'(\ab->2*a+ifbthen1else0)0-- | Convert an angle to binary representation.binary::Angle->BinAnglebinarya|a==0=([],[])|even(denominatora)=let(pre,per)=binary(doublea)b=a>=1/2in(b:pre,per)|otherwise=let(t,p)=head.dropWhile((1/=).denominator.fst).map(\q->(a*(2^q-1),q))$[(1::Int)..]s=numeratortn=fromIntegralpper=[s`testBit`i|i<-[n-1,n-2..0]]in([],per)-- | Tuning transformation for binary represented periodic angles.--   Probably only valid for angle pairs presenting ray pairs.btune::BinAngle->(BinAngle,BinAngle)->BinAnglebtune(tpre,tper)(([],per0),([],per1))=(concatMapftpre,concatMapftper)wherefFalse=per0fTrue=per1btune__=error"btune: can't handle pre-periods"-- | Tuning transformation for angles.tune::Angle->(Angle,Angle)->Angletunet(t0,t1)=unbinary$btune(binaryt)(binaryt0,binaryt1)-- | Elements of kneading sequences.dataKnead=Zero|One|Starderiving(Read,Show,Eq,Ord,Enum,Bounded,Data,Typeable)-- | Knead character representation.kneadChar::Knead->CharkneadCharZero='0'kneadCharOne='1'kneadCharStar='*'-- | Kneading sequences.  Note that the 'Aperiodic' case has an infinite list.dataKneading=Aperiodic[Knead]|PrePeriodic[Knead][Knead]|StarPeriodic[Knead]|Periodic[Knead]deriving(Read,Show,Eq,Ord,Data,Typeable)-- | Kneading sequence as a string.  The 'Aperiodic' case is truncated arbitrarily.prettyKneading::Kneading->StringprettyKneading(Aperiodicks)=mapkneadChar(take17ks)++"..."prettyKneading(PrePeriodicusvs)=mapkneadCharus++"("++mapkneadCharvs++")"prettyKneading(StarPeriodicvs)="("++mapkneadCharvs++")"prettyKneading(Periodicvs)="("++mapkneadCharvs++")"-- | The kneading sequence for an external angle.kneading::Angle->Kneadingkneadinga0'|a0==0=StarPeriodic[Star]|otherwise=fstkneadswherea0=wrapa0'lo=a0/2hi=(a0+1)/2kneads=kneading'1(doublea0)ks=(a0,One):sndkneadskneading'::Integer->Angle->(Kneading,[(Angle,Knead)])kneading'na|isJusti=caseiofJust0->caselastqsofStar->(StarPeriodicqs,[])_->(Periodicqs,[])Justj->let(p,q)=genericSplitAtjqsin(PrePeriodicpq,[])_->error"Fractal.Mandelbrot.Address.kneading (isJust -> Nothing?)"|a==lo=((a,Star):)`mapP`k|a==hi=((a,Star):)`mapP`k|lo<a&&a<hi=((a,One):)`mapP`k|hi<a||a<lo=((a,Zero):)`mapP`k|otherwise=error"Fractal.Mandelbrot.Address.kneading (unmatched?)"wherek=kneading'(n+1)(doublea)ps=genericTakenksqs=mapsndpsi=fmapfst.listToMaybe.filter((a==).fst.snd).zip[(0::Integer)..]$psmapPf~(x,y)=(x,fy)-- | The period of a kneading sequence, or 'Nothing' when it isn't periodic.period::Kneading->MaybeIntegerperiod(StarPeriodick)=Just(genericLengthk)period(Periodick)=Just(genericLengthk)period_=Nothingrho::Kneading->Integer->Integerrhovr|r>=1&&fmap(r`mod`)(periodv)/=Just0=((1+r)+).genericLength.takeWhileid.zipWith(==)vs.genericDropr$vs|otherwise=rhov(r+1)wherevs=unwrapv-- | Unwrap a kneading sequence to an infinite list.unwrap::Kneading->[Knead]unwrap(Aperiodicvs)=vsunwrap(PrePeriodicusvs)=us++cyclevsunwrap(StarPeriodicvs)=cyclevsunwrap(Periodicvs)=cyclevsorbit::(a->a)->a->[a]orbit=iterate-- | Internal addresses are a non-empty sequence of strictly increasing--   integers beginning with '1'.dataInternalAddress=InternalAddress[Integer]deriving(Read,Show,Eq,Ord,Data,Typeable)-- | Internal address as a string.prettyInternalAddress::InternalAddress->StringprettyInternalAddress(InternalAddress[])=error"Fractal.Mandelbrot.Address.InternalAddress.pretty"prettyInternalAddress(InternalAddress[x])=showxprettyInternalAddress(InternalAddress(x:ys))=showx++" "++prettyInternalAddress(InternalAddressys)-- | Construct a valid 'InternalAddress', checking the precondition.internalFromList::[Integer]->MaybeInternalAddressinternalFromListx0s@(1:_)=InternalAddress`fmap`fromList'0x0swherefromList'n[x]|x>n=Just[x]fromList'n(x:xs)|x>n=(x:)`fmap`fromList'xxsfromList'__=NothinginternalFromList_=Nothing-- | Extract the sequence of integers.internalToList::InternalAddress->[Integer]internalToList(InternalAddressxs)=xs-- | Construct an 'InternalAddress' from a kneading sequence.internalAddress::Kneading->MaybeInternalAddressinternalAddress(StarPeriodic[Star])=Just(InternalAddress[1])internalAddress(StarPeriodicv@(One:_))=Just.InternalAddress.address'per(genericLengthv)$vinternalAddress(Periodicv@(One:_))=Just.InternalAddress.address'per(genericLengthv)$vinternalAddressk@(Aperiodic(One:_))=Just.InternalAddress.address'inf.unwrap$kinternalAddress_=Nothingaddress'inf::[Knead]->[Integer]address'infv=address'vaddress'per::Integer->[Knead]->[Integer]address'perpv=takeWhile(<=p)$address'vaddress'::[Knead]->[Integer]address'v=address''1[One]whereaddress''skvk=sk:address''sk'vk'wheresk'=(1+).genericLength.takeWhileid.zipWith(==)v.cycle$vkvk'=genericTakesk'(cyclev)-- | A star-periodic kneading sequence's upper and lower associated--   kneading sequences.associated::Kneading->Maybe(Kneading,Kneading)associated(StarPeriodick)=Just(Periodica,Periodicabar)wheren=genericLengthkdivisors=[m|m<-[1..n],n`mod`m==0]abar=head.filter(and.zipWith(==)a'.cycle).map(`genericTake`a')$divisors(a,a')=if((n`elem`).internalToList)`fmap`internalAddress(Periodica1)==JustTruethen(a1,a2)else(a2,a1)a1=map(\s->casesofStar->Zero;t->t)ka2=map(\s->casesofStar->One;t->t)kassociated_=Nothing-- | The upper associated kneading sequence.upper::Kneading->MaybeKneadingupper=fmapfst.associated-- | The lower associated kneading sequence.lower::Kneading->MaybeKneadinglower=fmapsnd.associated-- | Angled internal addresses have angles between each integer in an--   internal address.dataAngledInternalAddress=UnangledInteger|AngledIntegerAngleAngledInternalAddressderiving(Read,Show,Eq,Ord,Data,Typeable)-- | Angled internal address as a string.prettyAngledInternalAddress::AngledInternalAddress->StringprettyAngledInternalAddress(Unangledn)=shownprettyAngledInternalAddress(Anglednra)|r/=1/2=shown++" "++show(numeratorr)++"/"++show(denominatorr)++" "++prettyAngledInternalAddressa|otherwise=shown++" "++prettyAngledInternalAddressa-- | Builds a valid 'AngledInternalAddress' from a list, checking the--   precondition that only the last 'Maybe Angle' should be 'Nothing',--   and the 'Integer' must be strictly increasing.angledFromList::[(Integer,MaybeAngle)]->MaybeAngledInternalAddressangledFromList=fromList'0wherefromList'x[(n,Nothing)]|n>x=Just(Unangledn)fromList'x((n,Justr):xs)|n>x&&0<r&&r<1=Anglednr`fmap`fromList'nxsfromList'__=NothingunsafeAngledFromList::[(Integer,MaybeAngle)]->AngledInternalAddressunsafeAngledFromList=fromList'0wherefromList'x[(n,Nothing)]|n>x=UnanglednfromList'x((n,Justr):xs)|n>x&&0<r&&r<1=Anglednr(fromList'nxs)fromList'__=error"Fractal.Mandelbrot.Address.unsafeAngledFromList"-- | Convert an 'AngledInternalAddress' to a list.angledToList::AngledInternalAddress->[(Integer,MaybeAngle)]angledToList(Unangledn)=[(n,Nothing)]angledToList(Anglednra)=(n,Justr):angledToListadenominators::InternalAddress->Kneading->[Integer]denominatorsav=denominators'(internalToLista)wheredenominators'(s0:ss@(s1:_))=letrr=rs0s1in(((s1-rr)`div`s0)+ifs0`elem`takeWhile(<=s0)(orbitprr)then1else2):denominators'ssdenominators'_=[]rss'=cases'`mod`sof0->st->tp=rhovnumerators::Angle->InternalAddress->[Integer]->[Integer]numeratorsraqs=zipWithnum(internalToLista)qswherenumsq=genericLength.filter(<=r).map(genericIndexrs)$[0..q-2]wherers=iterate(foldr(.)id.genericReplicates$double)r-- | The angled internal address corresponding to an external angle.angledInternalAddress::Angle->MaybeAngledInternalAddressangledInternalAddressr0=doletr=wrapr0k=kneadingri<-internalAddresskletd=denominatorsikn=numeratorsridreturn.unsafeAngledFromList.zip(internalToListi).(++[Nothing]).mapJust.zipWith(%)n$d-- | Split an angled internal address at the last island.splitAddress::AngledInternalAddress->(AngledInternalAddress,[Angle])splitAddressa=let(ps0,rs0)=unzip$angledToListaps1=reverseps0rs1=reverse(Nothing:initrs0)prs1=zipps1rs1f((p,Justr):qrs@((q,_):_))acc|p==denominatorr*q=fqrs(r:acc)fprsacc=gprsaccgprsacc=let(ps2,rs2)=unzipprsps3=reverseps2rs3=reverse(Nothing:initrs2)prs3=zipps3rs3aa=unsafeAngledFromListprs3in(aa,acc)infprs1[]-- | The inverse of 'splitAddress'.joinAddress::AngledInternalAddress->[Angle]->AngledInternalAddressjoinAddress(Unangledp)[]=UnangledpjoinAddress(Unangledp)(r:rs)=Angledpr(joinAddress(Unangled$p*denominatorr)rs)joinAddress(Angledpra)rs=Angledpr(joinAddressars)-- | The period of an angled internal address.addressPeriod::AngledInternalAddress->IntegeraddressPeriod(Unangledp)=paddressPeriod(Angled__a)=addressPerioda-- | Discard angle information from an internal address.stripAngles::AngledInternalAddress->InternalAddressstripAngles=InternalAddress.mapfst.angledToList-- | The pair of external angles whose rays land at the root of the--   hyperbolic component described by the angled internal address.externalAngles::AngledInternalAddress->Maybe(Rational,Rational)externalAngles=externalAngles'1(0,1)externalAngles'::Integer->(Rational,Rational)->AngledInternalAddress->Maybe(Rational,Rational)externalAngles'p0lohia0@(Unangledp)|p0/=p=casewakeeslohipof[lh]->externalAngles'plha0_->Nothing|otherwise=JustlohiexternalAngles'p0lohia0@(Angledpra)|p0/=p=casewakeeslohipof[lh]->externalAngles'plha0_->Nothing|otherwise=do{-      let num = numerator r          den = denominator r          q = p * den          ws = wakees lohi q          nums = [ num' | num' <- [ 1.. den - 1 ], let r' = num' % den, denominator r' == den ]          nws, nnums :: Integer          nws = genericLength ws          nnums = genericLength nums      guard (nws == nnums)      i <- genericElemIndex num nums      lh <- safeGenericIndex ws (i :: Integer)      externalAngles' q lh a-}letnum=numeratorrden=denominatorrws=wakees(0,1)dennums=[num'|num'<-[1..den-1],letr'=num'%den,denominatorr'==den]nws,nnums::Integernws=genericLengthwsnnums=genericLengthnumsguard(nws==nnums)i<-genericElemIndexnumnums(l,h)<-safeGenericIndexws(i::Integer)externalAngles'(p*den)(ifp>1then(tunellohi,tunehlohi)else(l,h))awakees::(Rational,Rational)->Integer->[(Rational,Rational)]wakees(lo,hi)q=letgaps(l,h)n|n==0=[(l,h)]--        | h - l < 1 % (2 ^ n - 1) = [(l, h)]|n>0=letgs=gaps(l,h)(n-1)cs=candidatesngsinaccumulatecsgs|otherwise=error"Fractal.Mandelbrot.Address.gaps !(n >= 0)"candidatesngs=letden=2^n-1in[r|(l,h)<-gs,num<-[ceiling(l*fromIntegerden)..floor(h*fromIntegerden)],letr=num%den,l<r,r<h,period(kneadingr)==Justn]accumulate[]ws=wsaccumulate(l:h:lhs)ws=let(ls,ms@((ml,_):_))=break(l`inside`)ws(_s,(_,rh):rs)=break(h`inside`)msinls++[(ml,l)]++accumulatelhs((h,rh):rs)accumulate__=error"Fractal.Mandelbrot.Address.gaps !even"insidex(l,h)=l<x&&x<hinchunk2.candidatesq.gaps(lo,hi)$(q-1)chunk2::[t]->[(t,t)]chunk2[]=[]chunk2(x:y:zs)=(x,y):chunk2zschunk2_=error"Fractal.Mandelbrot.Address.chunk2 !even"genericElemIndex::(Eqa,Integralb)=>a->[a]->MaybebgenericElemIndex_[]=NothinggenericElemIndexe(f:fs)|e==f=Just0|otherwise=(1+)`fmap`genericElemIndexefssafeGenericIndex::Integralb=>[a]->b->MaybeasafeGenericIndex[]_=NothingsafeGenericIndex(x:xs)i|i<0=Nothing|i>0=safeGenericIndexxs(i-1)|otherwise=Justx-- | Parse an angle.parseAngle::String->MaybeAngleparseAngles=caserunPpFraction()""sofLeft_->NothingRightf->Just(unFractionf)-- | Parse a list of angles.parseAngles::String->Maybe[Angle]parseAngless=caserunP(manypFraction)()""sofLeft_->NothingRightfs->Just(mapunFractionfs)-- | Parse a kneading element.parseKnead::String->MaybeKneadparseKneads=caserunPpKnead()""sofLeft_->NothingRightk->Justk-- | Parse a non-aperiodic kneading sequence.parseKneading::String->MaybeKneadingparseKneadings=caserunPpKneading()""sofLeft_->NothingRightks->Justks-- | Parse an internal address.parseInternalAddress::String->MaybeInternalAddressparseInternalAddresss=caserunP(manypNumber)()""sofLeft_->NothingRightns->internalFromList(mapunNumberns)-- | Parse an angled internal address, accepting some unambiguous--   abbreviations.parseAngledInternalAddress::String->MaybeAngledInternalAddressparseAngledInternalAddresss=caserunPparser()""sofLeft_->NothingRighta->JustadataToken=NumberInteger|FractionIntegerIntegerunFraction::Token->AngleunFraction(Fractiontb)=t%bunFraction_=error"Fractal.Mandelbrot.Address.unFraction"unNumber::Token->IntegerunNumber(Numbern)=nunNumber_=error"Fractal.Mandelbrot.Address.unNumber"typeParset=ParsecTString()Identitytparser::ParseAngledInternalAddressparser=dots<-pTokensaccum1tswhereaccump[]=return$Unangledpaccum_[Numbern]=return$Unanglednaccum_(Numbern:ts@(Number_:_))=doa<-accumntsreturn$Angledn(1%2)aaccum_(Numbern:Fractiontb:ts)=doa<-accum(n*b)tsreturn$Angledn(t%b)aaccump(Fractiontb:ts)=doa<-accum(p*b)tsreturn$Angledp(t%b)apTokens::Parse[Token]pTokens=do_<-pOptionalSpacets<-pToken`sepEndBy`pSpaceeofreturntspToken::ParseTokenpToken=choice[trypFraction,pNumber]pFraction::ParseTokenpFraction=doNumbertop<-pNumber_<-pOptionalSpace_<-string"/"_<-pOptionalSpaceNumberbottom<-pNumberguard$top<bottomreturn$FractiontopbottompNumber::ParseTokenpNumber=don<-foldl(\xy->10*x+y)0`fmap`map(toInteger.digitToInt)`fmap`many1digitguard$0<nreturn$NumbernpSpace::Parse[String]pSpace=many1(string" ")pOptionalSpace::Parse[String]pOptionalSpace=many(string" ")pKnead::ParseKneadpKnead=choice[string"0">>returnZero,string"1">>returnOne,string"*">>returnStar]pKneading::ParseKneadingpKneading=dopre<-manypKnead_<-string"("per<-many1pKnead_<-string")"return$case(nullpre,lastper)of(False,_)->PrePeriodicpreper(True,Star)->StarPeriodicper_->Periodicper