Monnet St Etienne and LIPENS Lyon JEANMICHEL MULLER CNRS LIPENS Lyon and ARNAUD TISSERAND INRIA LIPENS Lyon Polynomial approximations are almost always used when implementing functions on a computing system In most cases the polynomial that best app ID: 57321
Download Pdf The PPT/PDF document "Computing MachineEfcient Polynomial Appr..." is the property of its rightful owner. Permission is granted to download and print the materials on this web site for personal, non-commercial use only, and to display it on your personal computer provided you do not modify the materials and that you retain all copyright notices contained in the materials. By downloading content from our website, you accept the terms of this agreement.
ComputingMachine-EfÞcientPolynomialApproximationsNICOLASBRISEBARREUniversit«eJ.Monnet,St-EtienneandLIP-E.N.S.LyonJEAN-MICHELMULLERCNRS,LIP-ENSLyonARNAUDTISSERANDINRIA,LIP-ENSLyon Polynomialapproximationsarealmostalwaysusedwhenimplementingfunctionsonacomputing AuthorsÕaddresses:N.Brisebarre,LArAl,UniversiteJ.Monnet,23,rueduDrP.Michelon,F-42023EtienneCedex,FranceandLIP/Arenaire(CNRS-ENSLyon-INRIA-UCBL),ENSLyon,46eedÕItalie,F-69364LyonCedex07France;email:;J.-M.Muller,LIP/Ar ComputingMachine-EfÞcientPolynomialApproximations1.INTRODUCTIONThebasicßoating-pointoperationsthatareimplementedinhardwareonmod-ernprocessorsareaddition/subtraction,multiplication,andsometimesdivisionand/orfusedmultiply-add,thatis,anexpressionoftheform,computedwithoneÞnalroundingonly.Moreover,divisionisfrequentlymuchslowerthanadditionormultiplication(seeTableI),andsometimes,forinstance,ontheIta-nium,itisnothardwiredatall.Insuchcases,itisimplementedasasequenceoffusedmultiplyandaddoperations,usinganiterativealgorithm.Therefore,whentryingtoimplementagivenregularenoughfunction,itseemsreasonabletotrytoavoiddivisionsandtouseadditions,subtractions,multiplicationsandpossiblyfusedmultiply-adds.SincetheonlyfunctionsofonevariablethatonecanimplementusingaÞnitenumberoftheseoperationsandcomparisonsarepiecewisepolynomials,anaturalchoiceistofocusonpiecewisepolynomialapproximationsto.Indeed,mostrecentsoftware-orientedele-mentaryfunctionalgorithmsusepolynomialapproximations[Markstein2000;Muller1997;StoryandTang1999;Corneaetal.2002].Twokindsofpolynomialapproximationsareused:theapproximationsthatminimizetheaverageerror,calledleastsquaresapproximations,andtheap-proximationsthatminimizetheworst-caseerror,calledleastmaximumap-,orminimaxapproximations.Inbothcases,wewanttominimizeadistance,whereisapolynomialofagivendegree.Forleastsquaresapproximations,thatdistanceis:2,[isacontinuousweightfunctionthatcanbeusedtoselectpartsofofa,b]wherewewanttheapproximationtobemoreaccurate.Forminimaxapproximations,thedistanceis:is:a,b]=supaxb|p(x)f(x)|.Onecouldalsoconsiderdistancessuchasrel,[ Theleastsquaresapproximationsarecomputedbyaprojectionmethodus-ingorthogonalpolynomials.MinimaxapproximationsarecomputedusinganalgorithmcreditedtoRemez[Remes1934;Hartetal.1968].SeeMarkstein[2000]andMuller[1997]forrecentpresentationsofelementaryfunctionInthisarticle,weareconcernedwithminimaxapproximationsusingdis-dis-a,b].Andyet,ourgeneralmethodinSection3.2alsoappliestodistancerel,[.OurapproximationswillbeusedinÞnite-precisionarithmetic.Hence,thecomputedpolynomialcoefÞcientsareusuallyrounded:beaÞxedÞnitesequenceofnaturalintegers,thecoefÞcientoftheminimaxapproximation+···+ACMTransactionsonMathematicalSoftware,Vol.32,No.2,June2006. N.Brisebarreetal.TableI.Latencies(innumberofcycles)ofDoublePrecisionFloating-PointAddition,Multiplication,andDivisiononSomeRecentProcessors Processor FPadd FPmult. FPdiv. PentiumIV 5 7 38 PowerPC750 3 4 31 UltraSPARCIII 4 4 24 Alpha21264 4 4 15 AthlonK6-III 3 3 20 isroundedto,forinstance,thenearestmultipleof2.Bydoingthis,weobtainaslightlydifferentpolynomialapproximationö.Butwehavenoguaranteethatisthebestminimaxapproximationtoamongthepolynomialswhosedegree-coefÞcientisamultipleof2.TheaimofthisarticleistopresentawayofÞndingthisbesttruncatedapproximation.Wehavetwogoalsinmind:Ñratherlowprecision(e.g.,around24bits),hardware-oriented,forspeciÞc-purposeimplementations.Insuchcases,tominimizemultipliersizes(whichincreasesspeedandsavessiliconarea),thevaluesof,for1,shouldbeverysmall.Thedegreesofthepolynomialapproximationsarelow.Typi-calrecentexamplesaregiveninWeietal.[2001]andPineiroetal.[2001].Roughlyspeaking,whatmattershereisreducingthecost(intermsofdelayandarea)withoutmakingtheaccuracyunacceptable;Ñsingle-precisionordouble-precision,software-oriented,forimplementationoncurrentgeneralpurposemicroprocessors.Usingtable-drivenmethods,suchastheonessuggestedbyTang[1989,1990,1991,1992],thedegreeofthepolynomialapproximationscanbemaderatherlow.Roughlyspeaking,whatmattersinthiscaseistogetveryhighaccuracywithoutmakingthecost(intermsofdelayandmemory)unacceptable.OnecouldobjectthattheÕsarenotnecessarilyknownapriori.Forinstance,ifonewishesacoefÞcienttobeexactlyrepresentableindoubleprecisionarith-metic,oneneedstoknowtheorderofmagnitudeofthatcoefÞcienttoknowwhatvalueofcorrespondstothatwish.Andyet,inpractice,goodapproxi-mationsofthesamedegreetoagivenfunctionhavecoefÞcientsthatareveryclose(theapproachgiveninSection3.1showsthis)sothatusingourapproach,withpossiblytwodifferentvaluesofifthedegree-coefÞcientoftheminimaxapproximationisveryclosetoapowerof,sufÞces.Itisimportanttonoticethatourpolynomialapproximationswillbecom-putedonceonlyandwillbeusedveryfrequently(indeed,severalbilliontimesforanelementaryfunctionprograminawidelydistributedlibrary).Hence,ifitremainsreasonablyfast,thespeedofanalgorithmthatcomputesadequateapproximationsisnotextremelyimportant.However,inthepracticalcaseswehavestudiedsofar,ourmethodwillveryquicklygivearesult.Inthisarticle,weprovideageneralandefÞcientmethodforÞndingthebesttruncatedapproximation(s)(itisnotnecessarilyunique).ItconsistsinbuildingapolytopetowhichthenumeratorsofthecoefÞcientsofthis(these)besttruncatedapproximation(s)belong,suchthatcontainsanumberACMTransactionsonMathematicalSoftware,Vol.32,No.2,June2006. ComputingMachine-EfÞcientPolynomialApproximationsassmallaspossibleofpointsof.Onceitisachieved,wedoanexhaustivesearchbycomputingthenorms 2m0+a1 +···+ a,b]with(Themethodpresentedhereisveryßexiblesinceitalsoapplieswhenweimposesupplementaryconstraintsonthetruncatedpolynomials,and/orwhenrel,[isconsidered.Forexample,thesearchcanberestrictedtooddorevenpolynomialsor,moregenerally,topolynomialswithsomeÞxedcoefÞcients.Thisisfrequentlyuseful:onemightforinstance,wishthecom-putedvalueofexp()tobeexactlyoneif0(hence,requiringthedegree-0coefÞcientofthepolynomialapproximationtobe1).Ofcourse,onewouldliketotakeintoaccounttheroundofferrorthatoc-cursduringpolynomialevaluation:gettingthepolynomial,withconstraintsonthesizeofthecoefÞcients,thatminimizesthetotal(approximationplusroundoff)errorwouldbeextremelyuseful.Althoughwearecurrentlywork-ingonthatproblem,wedonotyethaveasolution.First,itisveryalgorithm-and-architecturedependent(forinstance,somearchitectureshaveanextendedinternalprecision).Second,sincethelargestroundofferrorandthelargestap-proximationerrorareextremelyunlikelytobeattainedatexactlythesamepoints,thetotalerrorisdifÞculttopredictaccurately.Andyet,hereareafewobservationsthatleadustobelievethat,inmanypracticalcases,ourapproachwillgiveuspolynomialsthatwillbeveryclosetotheseidealapproximations.Pleasenotethattheseobservationsaremerelyintuitivefeelings,andthatonecanalwaysbuildcasesforwhichtheidealap-proximationsdifferfromtheoneswecompute.(1)Goodapproximationsofthesamedegreetoagivenfunctionhaveco-efÞcientsthatareverycloseinpractice.Indeed,theapproachgiveninSection3.1showsthis.(2)WhenevaluatingtwopolynomialswhosecoefÞcientsareverycloseonvari-ablesthatbelongtothesameinputinterval,thelargestroundofferrorswillbeveryclose,too.(3)Inallpracticalcases,theapproximationerroroscillatesslowly,whereastheroundofferrorvariesveryquicklysothat,iftheinputintervalisrea-sonablysmall,anerrorveryclosetothemaximumerrorisreachednearanypoint.ThisisillustratedinFigures1and2:wehavedeÞnedasthepolynomialobtainedbyroundingtothenearestdoubleprecisionnumbereachcoefÞcientofthedegree-5minimaxapproximationtoin[0,1Figure1showsthedifference(approximationerror),andFigure2showsthedifferencebetweenthecomputedvalueandtheexactvalueof),assumingHornerÕsschemeisused,indoubleprecisionarithmetic. Sofar,wehavecomputedthesenormsusingthefunctionofMaple.OurresearchgroupisworkingonaCimplementationthatwillusemultipleprecisionintervalarithmetictogetcertiÞedupperandlowerboundsontheinÞnitenormofaregularenoughfunction.ACMTransactionsonMathematicalSoftware,Vol.32,No.2,June2006. N.Brisebarreetal. Fig.1.Approximationerror.Wehaveplottedthedifference Fig.2.Roundofferror.Wehaveplottedthedifferencebetweentheexactandcomputedvaluesof).ThecomputationsareperformedusingHornerÕsscheme,indoubleprecision,withoutusingalargerinternalformat.Theseobservationstendtoindicatethat,forallcandidatepolynomials,theroundofferrorswillbeveryclose,andthetotalerrorwillbeclosetothesumoftheapproximationandroundofferrors.Hence,thebestpolynomialwhencon-sideringtheapproximationerroronly,willbeveryclosetothebestpolynomialwhenconsideringapproximationandroundofferrors.Ofcourse,theseobservationsarenotproofs;theyarejustresultfromsomeexperiments,andwearefarfrombeingabletosolvethegeneralproblemofACMTransactionsonMathematicalSoftware,Vol.32,No.2,June2006. ComputingMachine-EfÞcientPolynomialApproximationsÞndingthebestpolynomial,withsizeconstraintsoncoefÞcients,whenconsid-eringapproximationandroundofferrors.Hopefully,theseremarkswillonedayhelptobuildamoregeneralmethod.Theoutlineofthearticleisthefollowing.WegiveanaccountofChebyshevpolynomialsandsomeoftheirpropertiesinSection2.InSection3,weÞrstprovideamethodbasedonChebyshevpolynomialsthatpartiallyanswerstheproblem,andthenwegiveageneralandefÞcientmethodbasedonpolytopesthatÞndsabesttruncatedapproximationofafunctionoveracompactinter-val[].DespitethefactthatitislessefÞcientandgeneralthanthepolytopemethod,wepresentthemethodbasedonChebyshevpolynomialsbecausethisapproachseemsinterestinginitself,issimple,andgivesresultsthatareeasytouseand,moreover,mightbeusefulinotherproblems.WeendSection3witharemarkillustratingtheßexibilityofourmethod.WeÞnishwithsomeexamplesinSection4.Wecompletethearticlewiththreeappendices.IntheÞrstone,wecollecttheproofsofthestatementsgiveninSection2.Inthesecondone,weprovealemmausedinSection3.2thatimplies,inparticular,theexistenceofabesttruncatedpolynomialapproximation.Inthelastone,wegiveaworkedexampleofthemethodspresentedhere.Toendthisintroduction,letusmentionthataCimplementationofourmethodisinprocessandalsothatthemethodappliestosomesignalprocess-ingproblems,namely,ÞndingtherationallinearcombinationofcosineswithconstraintsonthesizeinbitsoftherationalcoefÞcientsinordertoimplement(insoftwareorhardware)digitalFIRÞlters.Thiswillbethepurposeofafuturearticle.Asweonlydealwiththesupremumnorm,whereverthereisnoambiguity,wewillwriteinsteadof,whereisanyrealset.2.SOMEREMINDERSONCHEBYSHEVPOLYNOMIALSChebyshevPolynomialsTheChebyshevpolynomialscanbedeÞnedeitherbytherecurrencerelationorby)for)forApresentationofChebyshevpolynomialscanbefoundinBorweinandelyi[1995]andespeciallyinRivlin[1990].Thesepolynomialsplayacen-tralroleinapproximationtheory.ThefollowingpropertyiseasilyderivedfromDeÞnition1.ACMTransactionsonMathematicalSoftware,Vol.32,No.2,June2006. N.Brisebarreetal.Forn,wehave 2 n/2 k=0(1)k(nk1)! Hence,ThasdegreenanditsleadingcoefÞcientis.Ithasnrealroots,allstrictlybetweenWerecallthatapolynomialisapolynomialwhoseleadingcoefÞcientis1.ThefollowingstatementisawellknownandremarkablepropertyofChebyshevPolynomials.2(MOLYNOMIALSOFLetaThemonicdegree-npolynomialhavingthesmallestsmallesta,b]normis 22n1Tn2xba Inthefollowing,wewillmakeuseofthepolynomialsWehave(seeFoxandParker[1972,Chapter3],e.g.)),henceallthecoefÞcientsofarenonzerointegers.Now,westatetwopropositionsthatgeneralizeProperty2whendealingwithintervalsoftheform[0,]and[Leta,deÞne+···+ Letkbeaninteger,n,thepolynomial kT nx hasthesmallestsmallesta]normamongthepolynomialsofdegreeatmostnwithadegree-kcoefÞcientequalto.ThatnormisRemark1.Moreover,when0or1,wecanshowthatthispolynomialistheonlyonehavingthisproperty.Wedonotgivetheproofofthisuniquenesspropertyinthisarticlesinceweonlyneedtheexistenceresultinthesequel.Leta,deÞne+···+ Letkandnbeintegers,ÑIfkandnarebothevenorodd,thepolynomial k,nTnx ACMTransactionsonMathematicalSoftware,Vol.32,No.2,June2006. ComputingMachine-EfÞcientPolynomialApproximationshasthesmallestsmallesta,a]normamongthepolynomialsofdegreeatmostnwithadegree-kcoefÞcientequalto.ThatnormisÑElse,thepolynomial k,n1Tn1x hasthesmallestsmallesta,a]normamongthepolynomialsofdegreeatmostnwithadegree-kcoefÞcientequalto.Thatnormis3.GETTINGTHETRUNCATEDPOLYNOMIALTHATISCLOSESTTOAFUNCTIONONACOMPACTINTERVALbetworealnumbers,letbeafunctiondeÞnedon[]and1integers.DeÞneDeÞnem0,m1,...,mn]nasthesetofthepolynomialsofdegreelessthanorequaltowhosedegree-coefÞcientisamultipleof2forallbetween0and(wewillcallthesepolynomialstruncatedpolynomials),thatistosay, 2m0+a1 +···+ betheminimaxapproximationtoon[].DeÞneöasthepolynomialwhosedegree-coefÞcientisobtainedbyroundingthedegree-coefÞcientoftothenearestmultipleof2(withanarbitrarychoiceincaseofatie)forisanelementofofm0,m1,...,mn]n.AlsodeÞneandööa,b]andööa,b].WeassumethatöWestateourproblemasfollows.Let,wearelookingforatruncatedtruncatedm0,m1,...,mn]nsuchthatthata,b]=minqP[m0,m1,...,mn]nfq[a,b]andfp[a,b]K.(1)Lemma2inAppendix2impliesthatthenumberoftruncatedpolynomialssatisfying(1)isÞnite.,thisproblemhasasolutionsinceösatisÞes1.Itshouldbenotedthat,inthatcase,isnotnecessarilyequaltoöWecanput,forexample,3.1APartialApproachThroughChebyshevPolynomialsThetermpartialreferstothefactthattheintervalswecandealwithinthissectionmustbeoftheform[0,]or[],where0.Thisrestrictioncomesfromthefollowingtwoproblems.(1)Wedonothaveinthegeneralcase[]aresultanalogoustoPropositions1and2andsimpletostate.ACMTransactionsonMathematicalSoftware,Vol.32,No.2,June2006. N.Brisebarreetal.(2)Fromagivenpolynomialandagiveninterval,asimplechangeofvariableallowsustoreducetheinitialapproximationproblemtoanotherapproxi-mationproblemwithanintervaloftheform[0,]or[].Unfortunately,thischangeofvariables(thatdoesnotkeepthesizeofthecoefÞcientsin-variant)leadstoasystemofinequalitiesoftheformwefacedinSection3.2.Then,wehavetoperformadditionaloperationsinordertoproducecandi-datepolynomials.Indoingso,welosethemaininterestÑthesimplicityÑoftheapproachthroughChebyshevpolynomialsinthecases[0,]and[Wewillonlydealwithintervals[0,]where0sincethemethodpresentedinthefollowingeasilyadaptstointervals[]whereInthissection,wecomputeboundssuchthat,ifthecoefÞcientsofapolyno-polyno-m0,m1,...,mn]narenotwithinthesebounds,thenthena]�+K.Knowingtheseboundswillmakeanexhaustivesearchingofpossible.Todothis,considerapolynomialwhosedegree-coefÞcientis.FromProposition1,wehave isthenonzerodegree-coefÞcientof).Now,ifisatadistancegreaterthan,itcannotbebea]qp[0,a]pf[0,a]�K.Therefore,ifthereexists,suchthatthata]�+Kandtherefore.Hence,thedegree-necessarilyliesintheinterval[].Thuswehave ci2mipi2mi(pi+(+K)|i|) ,(2)since2isaninteger.Notethat,as00a],Condition(1)impliesinthatis,since2isaninteger, c02m0p0m0(f(0)+K) Wereplacewithmax()andwithmin(For,wehave1possiblevaluesfortheinteger2.Thismeansthatwehave1)candidatepolynomials.Ifthisamountissmallenough,wesearchforbycomputingthenormsnormsa],qrunningamongthepossiblepolynomials.Otherwise,weneedanadditionalsteptode-creasethenumberofcandidates.Hence,wenowgiveamethodforthispurpose.ACMTransactionsonMathematicalSoftware,Vol.32,No.2,June2006. ComputingMachine-EfÞcientPolynomialApproximationsItallowsustoreducethenumberofcandidatepolynomialsdramaticallyandappliesmoregenerallytointervalsoftheform[]whereareanyrealnumbers.3.2AGeneralandEfÞcientApproachThroughPolytopesFromnowon,wewilldealwithintervalsoftheform[]whereanyrealnumbers.WerecallthefollowingdeÞnitionsfromSchrijver[2003].1.LetAsubsetiscalledapolyhedronifthereexistsanwithrealcoefÞcientsandavector(forsome0)suchthatAsubsetiscalledapolytopeifisaboundedpolyhedron.Apolyhedron(orapolytope)iscalledrationalifitisdeterminedbyarational,respectively,systemoflinearinequalities.1inequalitiesgivenby2deÞnearationalpolytopeofwhichthenumeratorsof(i.e.the2)belongto.Theideaistobuildapolytope,stillcontainingthe2,suchthatisthesmallestpossible,whichmeansthatthenumberofcandidatepolynomialsisthesmallestpossible,inordertoreduceasmuchaspossibletheÞnalstepofcomputationofsupremumnorms.Oncewegetthispolytope,weneedanefÞcientwayofproducingthesecandidates,thatis,anefÞcientwayofscanningtheintegerpoints(i.e.,topointswithintegercoordinates)oftherationalpolytopewebuilt.Severalalgorithmsallowustoachievethis.TheonegiveninAncourtandIrigoin[1991]usestheFourier-Motzkinpairwiseelimination,theonegiveninFeautrier[1988]andCollardetal.[1995]isaparameterizedversionoftheDualSimplexmethodandtheonegiveninLeVergeetal.[1994]isbasedonthedualrepresentationofpolyhedrausedinPolylib[ThePolylibTeam2004].ThelasttwoalgorithmsallowustoproduceinanoptimizedwaytheloopsinourÞnalprogramofexhaustivesearch.Notethatthesealgorithmshave,atworst,thecomplexityofintegerlinearprogramming.Now,letusgivethedetailsofthemethod.Firstwenoticethatthepreviousapproachhandlestheunknownsratelywhichseemsunnatural.Hence,thebasicaimofthemethodistoconstructapolytopedeÞnedbyinequalitiesthattakeintoaccountinamoresatisfyingwaythedependencerelationsbetweentheunknowns.Thispolytopeshouldcon-tainfewerpointsofthantheonebuiltfromtheChebyshevpolynomialsÕapproach.TheexamplesofSection4indicatethatthisseemstobethecase(andtheimprovementscanbedramatic). Afterthesubmissionofthisarticle,wereadthatthisideahasalreadybeenproposedinHabsiegerandSalvy[1997],apaperdealingwithnumber-theoreticalissues,buttheauthorsdidnothaveanyefÞcientmethodforscanningtheintegerpointsofthepolytope.Byoptimized,wemeanthatallpointsofthepolytopearescannedonlyonce.Infact,thealgorithmgiveninFeautrier[1988]andCollardetal.[1995]has,inthesituationwearefacing,thecomplexityofrationallinearprogramming.ACMTransactionsonMathematicalSoftware,Vol.32,No.2,June2006. N.Brisebarreetal.Condition(1)meansforallalla,b].Theideaistoconsiderinequalities(3)foracertainnumber(chosenbytheuser)ofvaluesofofa,b].Aswewanttobuildrationalpolytopes,thesevaluesmustberationalnumbers.Weproposethefollowingchoiceofvalues.Letarationalapproximationtogreaterthanorequalto,letbearationalapproximationtolessthanorequaltoandsuchthat,letbeanonzeronaturalintegerchosenbytheuser.WeshowinLemma2inAppendix2thatwemusthaveinordertoensurethatwegetapolytope(whichimpliestheÞnitenessofthenumberofthebesttruncatedapproximation(s)).Weconsidertherationalvalues .Thenagain,sincethepolytopehastoberational,wecompute,for,tworationalthatarerationalapproximationsto,respectively,suchthat.TherationalpolytopesearchedisthereforedeÞnedbytheinequalitiesissmallenough(thiscanbeestimatedthankstoPolylib),westartourexhaustivesearchbycomputingthenorms 2m0+a1 2m1x+n a,b],(5)with(.ThissetcanbescannedefÞcientlythankstooneofthealgorithmsgiveninAncourtandIrigoin[1991],Feautrier[1988]andCollardetal.[1995]orLeVergeetal.[1994]thatwehavepreviouslyquoted.Orelse,weincreasethevalueoftheparameterinordertoconstructanotherrationalpolytopethatcontainsfewerelementsof.Wemustpointoutthat,assumingthatthenewparameterisgreaterthan,doesnotnecessarilyleadtoanewpolytopewithfewerelementsofinside(itiseasytoshowsuchcounterexamples),butitisreasonabletoexpectthat,generally,apolytopebuiltwithagreaterparametershouldcontainfewerelementsofinsidesinceitisdeÞnedfromalargernumberofinequalities(4)or,inotherwords,asourdiscretizationmethodisdoneusingalargernumberofrationalpoints,whichshouldallowustorestrictthenumberofpossiblecandidatessincetherearemoreconditionstosatisfy.Itisindeedthecaseifwechooseanypositiveintegermultipleofasnewparameter.Inthiscase,thenewpolytope,associatedwiththeparameter,with,isasubsetofisbuiltfromthesetofrationalpoints whichisasubsetof fromwhichisbuilt. Choosingequally-spacedrationalvaluesseemsanaturalchoicewhendealingwithregularenoughfunctions.Andyet,inafewcases,wegetabetterresultwithveryirregularly-spacedpoints.Thisissomethingweplantoinvestigateinthenearfuture.ACMTransactionsonMathematicalSoftware,Vol.32,No.2,June2006. ComputingMachine-EfÞcientPolynomialApproximationsRemark2.Aswesaidintheintroduction,thismethodisveryßexible.Wegivesomeexamplestoillustrateit.ÑIfwerestrictoursearchtooddtruncatedpolynomials(inthiscase,weput1andwemusthave),itsufÞcestoreplaceinequalities(4)tocreateapolytopewhosepointswithintegercoordinateswescan.ÑIfwerestrictoursearchtotruncatedpolynomialssomeofwhosecoefÞcientshaveÞxedvalues,(e.g.,ifweassumethatthetruncatedpolynomialssoughthaveconstanttermequalto1)itsufÞcestoreplaceinequalities(4)withtocreateapolytopewhosepointswescanwithintegercoordinates(wemusthaveÑOurmethodalsoappliestothesearchforthebesttruncatedpolynomialwithrespecttotherelativeerrordistancerel,[deÞnedintheintroduction.Inthiscase,wecanstatetheproblemasfollows.Let0,wesearchforatruncatedpolynomialpolynomialm0,m1,...,mn]nsuchthatrel,[[m0,m1,...,mn]nfqrel,[rel,[ItsufÞcestoconsidertheinequalitiesforatleast1distinctrationalvaluesofofa,b]tomakeapolytopetowhichweapplythemethodpresented.4.EXAMPLESWeimplementedinMaplethemethodgiveninSection3.1,andwestarteddevelopingaClibrarythatimplementsthemethoddescribedinSection3.2.Forproducingtheresultspresentedinthissection,weimplementedtheapproachthroughpolytopesinaCprogramofourown,basedonthepolyhedrallibraryPolylib[ThePolylibTeam2004].Thisallowsustotreatalotofexamples,butitistooroughlydonetotackleexampleswithapproximationsofdegree15to20,forinstance.OurgoalistoimplementthemethodpresentedhereinaClibrarywhichwouldusePolylibandPIPFeautrier[2003]thatimplementstheparametricintegerlinearprogrammingsolvergiveninFeautrier[1988]andCollardetal.[1995]forscanningtheintegerpointsofthepolytope.ThechoiceofPIPinsteadofthealgorithmgiveninLeVergeetal.[1994]isduetoACMTransactionsonMathematicalSoftware,Vol.32,No.2,June2006. N.Brisebarreetal.TableII.Examples f [a,b] m ö K 1 cos 0, 4 [12,10,6,4] 1.135...104 6.939...104 ö 2 exp 0,1 2 [15,14,12,10] 2.622...105 3.963...105 ö 3 exp 0,log [56,45,33,23] 1.184...10 .362...10 ö 4 arctan(1+x) 0,1 4 [24,21,18,17,16] 2.381...108 3.774...108 ö 5 exp log(2) log(2) [25,17,9] 8.270...10 .310...109 ö 6 exp log(2) log(2) [28,19,9] 8.270...10 .310...109 ö 7 log(3/4+x) log(2) 4,1 [12,9,7,5] .371...104 7.731...104 ö 8 log( 2/2+x) log(2) 1 2 2,2 2 2 [12,9,7,5] .371...104 9.347...104 ö TableIII.CorrespondingResults Chebyshev Polytope(d) T1 T2 GainofAccuracyinBits 1 330 1(4) 0.62 0.26 1.5 2 84357 9(20) 0.51 2.18 0.375 3 9346920 15(20) 0.99 1.77 0.22 4 192346275 1(20) 0.15 0.55 0.08 5 1 0(4) 0.05 0 0 6 4 1(4) 0.03 0.10 0.41 7 38016 2(15) 0.13 0.69 0.06 8 12(20) 0.59 5.16 0.26 thefactthatourpolytopemay,insomecases,havealargeamountofvertices(duetotheneedforalargeamountofpoints,i.e.,constraints,whenbuildingthepolytope)whichcangeneratememoryproblemsifweuseLeVergeetal.[1994].WealsoneedtheMPFRmultiple-precisionlibrary[TheSpacesProject2004]sincewefaceprecisionproblemswhenformingthepolytope,andtheGMPlibrary[Granlund2002]sincethepolytopeisdeÞnedbyamatrixandavectorwhichmayhaveverylargeintegercoefÞcients.TheGMPlibraryisalsousedinsidesomePolylibcomputations.WeputsomeexamplesinTableII,andwegroupthecorrespondingresultsinTableIII.InthelastcolumnofTableII,thenotationmeansthatwechoseavalueslightlysmallerthanö,namelyöroundeddowntofourfractionaldigits.IntheÞrstcolumnofTableIII,onecanÞndthenumberofcandidatesgivenbyChebyshevÕsapproach.Inthesecond,wegivethenumberofcandidatesgivenbytheapproachthroughpolytopesandthecorrespondingparameterdenotesthetimeinsecondsforproducingthecandidatepolynomialswiththepolytopemethod.denotesthetimeinsecondsforcomputingthenorms(5)which,forthemoment(cf.footnote1),isdonewithMaple9inwhichthevalueofwasÞxedequalto20.Inthelastcolumn,wegivethegainofaccuracyinbitsthatwegetifweusethebesttruncatedpolynomialinsteadofpolynomialöAllthecomputationsweredoneona2.53GHzPentium4computerrun-ningDebianLinuxwith256MBRAM.ThecompilerusedwasACMTransactionsonMathematicalSoftware,Vol.32,No.2,June2006. ComputingMachine-EfÞcientPolynomialApproximations Fig.3.ThisÞgureshowshowthenumberofintegerpointsofthepolytopedropswhendecreases,inthecaseofexample2.Thetermallcandidatesmeanstheintegerpointsofthepolytope.Bygoodcandidate,wemeanthepointsthatfulÞlltherequirement(1).4.1ChoiceoftheExamplesTheseexamplescorrespondtocommoncasesthatoccurinelementaryandspecialfunctionapproximation[Muller1997].Examples1,2,and4areofimmediateinterest.Example3isofinterestforimplementingtheexponen-tialfunctionindoubleprecisionarithmetic,usingatable-drivenmethodsuchasTangÕsmethod[Tang1991,1992].Examples5and6alsocorrespondtoatable-drivenimplementationoftheexponentialfunctioninsingleprecision.Examples7and8aimatproducingverycheapapproximationstologarithmsin[12,1],forspecialpurposeapplications.Figure3showshowthenumberofintegerpointsofthepolytopedropswhendecreases,inthecaseofexample2.APPENDIX1.PROOFSOFPROPOSITIONS1AND2ProvingthosepropositionsÞrstrequiresthefollowingtworesults.TheÞrstoneiswellknown.Thereareexactlynvaluesxsuchthatwhichsatisfysatisfy1,1]Thatis,themaximumabsolutevalueofTisattainedatthexÕs,andthesignofTalternatesatthesepoints..Theseextremepointsaresimplythepointscos( ACMTransactionsonMathematicalSoftware,Vol.32,No.2,June2006. N.Brisebarreetal.beanincreasingsequenceofnonnegativeintegers+···++x],theneitherPorPhasatmostnzerosin.Byinductionon.For0,itisstraightforward.Nowweassumethatthepropertyistrueuntiltherank.Let+···++x]with00.Assumethathasatleastzerosin(0,).Thenhasatleast2zerosin(0,Thus,thenonzeropolynomial+···+has,fromRolleÕsTheorem,atleast1zerosin(0,)whichcontradictstheinductionhypothesis. ROOFOF1.FromProperty3,thereexist01suchthatthata]=1k(1)ni.Letq(x)= nj=0,j =kcjxjR[x]satisfysatisfya]|1k|.Thenthepolynomial))hastheformandisnotidenticallyzero.Asitchangessignbetweenanytwoconsecutiveextremaof,ithasatleastzerosin(0,).Hence,fromLemma1,itmustvanishidenticallywhichisthecontradictiondesired.ROOFOF2.Weassumethatisevensinceastraightforwardadaptationoftheproofinthiscasegivestheproofofthecasewhereisodd.First,wesupposeeven.Let+···++···+suchthatthata,a]|1/k,n|.Then,forall |k,n|Q(xP(x)+P(x) 21 thatis,sincearebotheven,forall +···++···+ FromProperty3andtheinequalityinequalitya,a]|1/k,n|=Tn(·/a)[a,a],thechangessignbetweentwoconsecutiveextremaof).Then,ithasatleastdistinctzerosinina,a]and,moreprecisely,in[,0)]since0isanextremaof)asiseven.Hence,thepolynomial hasatleast2distinctzerosin(0, ]:itisidenticallyzerofromLemma1.ACMTransactionsonMathematicalSoftware,Vol.32,No.2,June2006. ComputingMachine-EfÞcientPolynomialApproximationsWehavejustprovedthat k,nTnx a +P(x)P(x) Werecallthatforallalla,a].Therefore,wehave,bysubstituting1and1to |k,n|1 k,nTn1 aP(1)P(1) 21 |k,n|1 k,nTn1 aand1 |k,n|1 k,nTn1 aP(1)P(1) 21 |k,n|1 k,nTn1 iseven,weknowthat1.Hence,weobtain whichisthecontradictiondesired.Now,wesupposeodd.Let+···++···+suchthatthata,a]|1/k,n1|.Then,forall |k,n1|Q(xP(x)P(x) 21 thatis,forall +···++···+ FromProperty3andtheinequalityinequalitya,a]|1/k,n1|=Tn1(·/a)[a,a],thepolynomialchangessignbetweentwoconsecutiveextremaof).Then,ithasatleast1distinctzerosin[]and,inparticular,atleast2distinctzerosin[,0)Hence,thepolynomial x)/ hasatleast1distinctzerosin(0, ]:itisidenticallyzerofromLemma1.Wehavejustobtainedthat k,n1Tn1x a+P(x)+P(x) Hereagain,werecallthatforallalla,a].Thus,itfollows,bysubstituting1and1to |k,n1|1 k,n1Tn1 aP(1)+P(1) 21 |k,n1|1 k,n1Tn11 aand1 |k,n1|1 k,n1Tn11 aP(1)+P(1) 21 |k,n1|1 k,n1Tn11 ACMTransactionsonMathematicalSoftware,Vol.32,No.2,June2006. N.Brisebarreetal.iseven,wehave1.Hence,weobtain whichisthecontradictiondesired.APPENDIX2Wenowprovethefollowingstatement.Letd,letx(resp.)suchthatxifij,letthesetdeÞnedbyresp.foriIfdn,thenisapolytope(resp.rationalpolytope).Ifdn,theneitheremptyorisunbounded..Thesetisobviouslyapolyhedron(respectivelyrationalpolyhe-dron).So,ifwewanttoprovethatisapolytope(respectivelyrationalpoly-tope),itsufÞcestoshowthatitisbounded.First,weassume.Then,iscontainedinthesetdeÞnedbyThelinearapplicationisanisomorphismforthematrixisanonsingularVandermondematrixsince.Thelinearappli-,deÞnedonaÞnitedimensional-vectorspace,iscontinuousandthesetsetli,ui]isacompactof.Thus,thesetwhichisequalto)isacompactof,whichimpliesthatisnecessarilybounded.ACMTransactionsonMathematicalSoftware,Vol.32,No.2,June2006. ComputingMachine-EfÞcientPolynomialApproximationsNow,weassumenonempty.Then,wenoticethatisdeÞnedbytheinequalities,thematrixisnonsingular.Hence,inparticular,thepolyhedroncontainsthefamilyofwhichprovesthatcannotbebounded. APPENDIX3.AWORKEDEXAMPLELetusgivethedetailsoftheÞrstexampleofTableII.Wefocusonthedegree-3approximationtothecosinefunctionin[0,First,wedetermine(hereusingthepackageofMaplewithequalto10)thedegree-3minimaxpolynomialassociatedtocoson[0,4].Wee /4]=0.0001135879209.Thismeansthatsuchanapprox-imationisnotgoodenoughforsingle-precisionimplementationofthecosinefunction.Itcanbeofinterestforsomespecialpurposeimplementations.Wehave 16x317 32x2+5 1andöö /4]=0.0006939707.Letuschoose2.Then,theapproachthatusesChebyshevpolynomialsÑ3possiblevaluesbetween40954096and40974096forthedegree-0term,Ñ22possiblevaluesbetween512and151024forthedegree-1term,Ñ5possiblevaluesbetween16and2forthedegree-2term,Ñ1possiblevalueequalto116forthedegree-3term.ACMTransactionsonMathematicalSoftware,Vol.32,No.2,June2006. N.Brisebarreetal.Hence,theapproachthatusesChebyshevpolynomialsprovides330candi-datepolynomials.Thisisareasonablysmallnumberthetimenecessaryfortheexhaustivecomputationofthenorms5issmall.Thereisnoneedtousetheapproachbasedonpolytopes;weuseitanywayinthefollowingtoshowhowitworks.First,wedeÞneasubinterval[]of[0,4]withrationalbounds.Let0,wechoosearationalapproximation4suchthat4anditsnumeratoranddenominatoraresmallinordertospeedupthecomputations.ThebiggertheintegersoccuringinthedatadeÞningthepolytopeare,theslowerthecomputationsareperformed.Hence,let9,whichisaconvergentofthecontinuedfractionexpansionofWechooseavalueofequalto4,thatis,thepolytopeisbuiltfromÞverationalpointswhichare0, 4·7 9,2 4·7 9,3 4·7 9,7 .Let,3.Therefore,wewanttheinequalitieshereaftertobesatisÞed: 36K3 j=0aj 2mj7i 36jcos7i Westillassumethat2.TheseinequalitiesdeÞneapolytope,butitisnotarationalone.Hence,wecomputerationalapproximationsofthelowerandupperboundsin7:for,4,wecomputesuchthat ,cos( .Hereagain,ourtargetistopreventtheincreaseofthesizeoftheinvolvedintegers.Therefore,wecaneithercomputeconvergentsoftheorwecanimposethat,foreachhavethesamedenominatorasthe 2mj(7i ,3.Wechoosethesecondpossibility,hence,forall,4,weput 36K,/Diandui=-Dicos7i istheleastcommonmultipleofthedenominatorsofthe 2mj(7i ,3.WeobtainarationalpolytopedeÞnedby100010007295671764137272911347056109762763588137272922682822487808 Around1.5secondona2.53GHzPentium4computerrunningDebianLinuxwith256MBRAM.ACMTransactionsonMathematicalSoftware,Vol.32,No.2,June2006. ComputingMachine-EfÞcientPolynomialApproximationsTableIV.NumberofPoints,andComputationalTimes(inCandinMaple)forVariousChoicesofinExample1 Time[s] d K T1 T2 #Polynomials 9 0010 2.60 1.68 7 9 0010 1.59 0.98 4 9 5010 0.25 0.26 1 9 9310 2.55 1.40 6 18 9310 2.84 1.40 6 36 9310 3.37 1.40 6 72 9310 7.39 1.40 6 :(4095,6,34,1)4410(gain15bits) Then,ourcurrentCimplementationreturnsonlyonecandidate(4095,6,34,1).WecheckthatitsatisÞescondition(1).Wethereforedirectly 16x317 32x2+3 512x+4095 /4]=0.0002441406250.Inthisexample,thedistancebetweenisapproximately035timesthedistancebetweenandö.Usingourapproachsavesaround5bitsofaccuracy.TableIVgivessomeÞgures(numberofpointsofthepolytope,delayofcom-putation)dependingonthechoicesofinthisexample.Hereagain,denotesthetimeinsecondsforproducingthecandidatepolynomials,anddenotesthetimeinsecondsforcomputingthenorms5.ACKNOWLEDGMENTSWewouldliketothankthereferees,whogreatlyhelpedtoimprovethequalityofthemanuscript.WewouldalsoliketothankNicolasBoullisandSergeTorres,whogreatlyhelpedcomputingtheexamples,andPaulFeautrierandCBastoul,whoseexpertiseinpolyhedriccomputationshasbeenhelpful.,C.,F.1991.Scanningpolyhedrawithdoloops.InProceedingsofthe3rdACMSIGPLANSymposiumonPrinciplesandPracticeofParallelProgramming(PPoPPÕ91)ACMPress,NewYork,NY,39Ð50.,P.ELYI,T.1995.PolynomialsandPolynomialsInequalities.GraduateTextsinMathematics,161.Springer-Verlag,NewYork,NY.,J.-F.,FEAUTRIER,P.,,T.1995.ConstructionofdoloopsfromsystemsofafÞneconstraints.Parall.Process.Lett.5,421Ð436.,M.,H,J.,,P.T.P.2002.ScientiÞcComputingonItanium-BasedSystemsIntelPress.EAUTRIER,P.1988.Parametricintegerprogramming.RAIRORech.Op«er.22,3,243Ð268.EAUTRIER,P.2003.PIP/Piplib,aparametricintegerlinearprogrammingsolver.,L.,I.B.1972.ChebyshevPolynomialsinNumericalAnalysis.OxfordMathe-maticalHandbooks.OxfordUniversityPress,London,UK.ACMTransactionsonMathematicalSoftware,Vol.32,No.2,June2006. N.Brisebarreetal.,T.2002.GMP,theGNUmultipleprecisionarithmeticlibrary,version4.1.2.,L.ALVY,B.1997.OnintegerChebyshevpolynomials.Math.Computat.66,218,,J.F.,C,E.W.,LAWSON,C.L.,MAEHLY,H.J.,M,C.K.,R,J.R.,T,H.G.,,C.1968.ComputerApproximations.JohnWiley&Sons,NewYork,NY.,H.,VAN,V.,,D.K.1994.Loopnestsynthesisusingthepolyhedrallibrary.Tech.Rep.INRIAResearchReportRR-2288,(May)INRIA.,P.2000.IA-64andElementaryFunctions:SpeedandPrecision.Hewlett-PackardProfessionalBooks.PrenticeHall.,J.-M.1997.ElementaryFunctions,AlgorithmsandImplementation.Birkhauser,Boston,MA.,J.,B,J.,,J.-M.2001.Faithfulpoweringcomputationusingtablelook-upandafusedaccumulationtree.InProceedingsofthe15thIEEESymposiumonComputerArithmetic(Arith-15).BurgessandCiminieraEds.IEEEComputerSocietyPress,LosAlamitos,CA,40Ð58.,E.1934.SurunproceconvergentdÕapproximationssuccessivespourdeterminerlesomesdÕapproximation.C.R.Acad.Sci.Paris198,2063Ð2065.,T.J.1990.ChebyshevPolynomials.FromApproximationTheorytoAlgebra2ndEd.PureandAppliedMathematics.JohnWiley&Sons,NewYork,NY.,A.2003.Combinatorialoptimization.PolyhedraandefÞciency.Vol.A.AlgorithmsandCombinatorics,24.Springer-Verlag,Berlin,Germany.TORY,S.,P.T.P.1999.NewalgorithmsforimprovedtranscendentalfunctionsonIA-64.InProceedingsofthe14thIEEESymposiumonComputerArithmetic.KorenandKornerup,Eds.IEEEComputerSocietyPress,LosAlamitos,CA,4Ð11.,P.T.P.1989.Table-drivenimplementationoftheexponentialfunctioninIEEEßoating-pointarithmetic.ACMTrans.Math.Soft.15,2(June),144Ð157.,P.T.P.1990.Table-drivenimplementationofthelogarithmfunctioninIEEEßoating-pointarithmetic.ACMTrans.Math.Soft.16,4(Dec.),378Ð400.,P.T.P.1991.Tablelookupalgorithmsforelementaryfunctionsandtheirerroranalysis.Proceedingsofthe10thIEEESymposiumonComputerArithmetic.P.KornerupandD.W.Matula,Eds.IEEEComputerSocietyPress,LosAlamitos,CA,232Ð236.,P.T.P.1992.Table-drivenimplementationoftheexpm1functioninIEEEßoating-pointarithmetic.ACMTrans.Math.Soft.18,2(June),211Ð222.OLYLIB.2004.Polylib,alibraryofpolyhedralfunctions,version5.20.0.u-strasbg.fr/polylib/PACESPROJECT.2004.MPFR,themultipleprecisionßoatingpointreliablelibrary,version,B.,C,J.,,J.2001.High-performancearchitecturesforelementaryfunctiongeneration.InProceedingsofthe15thIEEESymposiumonComputerArithmetic(Arith-15)BurgessandCiminieraEds.IEEEComputerSocietyPress,LosAlamitos,CA,136Ð146.ReceivedApril2004;revisedMay2005;acceptedSeptember2005ACMTransactionsonMathematicalSoftware,Vol.32,No.2,June2006.