From 5b95219e094cfd8806ecad8398b28fa451ff45d4 Mon Sep 17 00:00:00 2001 From: John Maddock Date: Wed, 4 Oct 2006 09:57:05 +0000 Subject: [PATCH] Fixed digamma 128-bit error rates with a more accurate value for the root. Updated png's: hopefully they'll all render OK now. Tweeked docs to reflect recent advances. [SVN r3237] --- doc/digamma.qbk | 9 +++++++-- doc/equations/chf.png | Bin 873 -> 922 bytes doc/equations/dist_tutorial1.png | Bin 2218 -> 1883 bytes doc/equations/dist_tutorial2.png | Bin 796 -> 817 bytes doc/equations/dist_tutorial3.png | Bin 1359 -> 1436 bytes doc/equations/dist_tutorial4.png | Bin 609 -> 630 bytes doc/equations/hazard.png | Bin 776 -> 826 bytes doc/equations/students_t_dist.png | Bin 1192 -> 465 bytes doc/erf.qbk | 9 --------- .../boost/math/special_functions/digamma.hpp | 4 ++-- minimax/f.cpp | 19 ++++++++++++++++-- test/test_digamma.cpp | 11 +--------- 12 files changed, 27 insertions(+), 25 deletions(-) diff --git a/doc/digamma.qbk b/doc/digamma.qbk index 1ab430f6e..6d92d6fa9 100644 --- a/doc/digamma.qbk +++ b/doc/digamma.qbk @@ -40,9 +40,14 @@ than the one shown will have __zero_error. [[53] [Win32 Visual C++ 8] [Peak=0.98 Mean=0.36] [Peak=0.99 Mean=0.5] [Peak=0.95 Mean=0.5] [Peak=214 Mean=16] ] [[64] [Linux IA32 / GCC] [Peak=1.4 Mean=0.4] [Peak=1.3 Mean=0.45] [Peak=0.98 Mean=0.35] [Peak=180 Mean=13] ] [[64] [Linux IA64 / GCC] [Peak=0.92 Mean=0.4] [Peak=1.3 Mean=0.45] [Peak=0.98 Mean=0.4] [Peak=180 Mean=13] ] -[[113] [HPUX IA64, aCC A.06.06] [Peak=0.9 Mean=0.4] [Peak=200 Mean=14] [Peak=0.99 Mean=0.4] [Peak=200 Mean=14] ] +[[113] [HPUX IA64, aCC A.06.06] [Peak=0.9 Mean=0.4] [Peak=1.1 Mean=0.5] [Peak=0.99 Mean=0.4] [Peak=64 Mean=6] ] ] +As shown above, error rates for positive arguments are generally very low. +For negative arguments there are an infinite number of irrational roots: +relative errors very close to these can be arbitrarily large, although +absolute error will remain very low. + [h4 Testing] There are two sets of tests: spot values are computed using @@ -119,7 +124,7 @@ rational approximation. Note that since X0 is irrational, we need twice as many digits in X0 as in x in order to avoid cancellation error during the subtraction (this assumes that /x/ is an exact value, if it's not then all bets are off). That means that even when x is the value of the root rounded to the nearest -representable value, the result of digamma(x) ['[*will not be zero*]]. +representable value, the result of digamma(x) ['[*will not be zero]]. [endsect][/section:tgamma The Gamma Function] diff --git a/doc/equations/chf.png b/doc/equations/chf.png index e13dcf00d23916e0e105314af892bfd2ca31713a..2923d3137390378a2e1f693dfbd19acc4e40cdd4 100644 GIT binary patch delta 859 zcmaFKHj90NBnKM<1H;y`<(nodn%6Vyc)B=-RNQ)V*SF}gfs-!ebDPpUZP-b4>&c@Kf!@wcgAbS6CdiZPxg_Hv+O5f`L zR&BlMYD|6!)D! zKH09HF{k$Ztj#s?Z%!o(E5;wc)_BRP>Au;m*&BU&mxo>7*7dRa?1EUc>gc+@maOw1 z-2Eh%3dL5opS!^IX~xv28>N|V)iZ8gw?65y&C_2o@>gTDit;NTU--Iocl0!=1!V~w z+e&`#R}Q#v&wTNZ7k0B9tEaBMwE1HXQ|T2(4Z-y@(!2IPS6 z&4%YicBTF8SsD8_T=)F0o~M7WG1i0>x`|H+H%q$z{kMXMyk_CI|CnG<^F?S8*h#W>8%>A=3Fdlz1c zJTt|pD9PsQynaUBfA3#EHGV$-TY;guJDl=Zr<*;%+l{PE83H8M*xvz>F(|Hg^Vn|v>q-^XLF?M*x9 zQ)@C#`DCtL_l;*|ct(+G`iGgdx0f6=mFk-ipZs^4^y7`+|NQjnpM1=H^^~pGZg6ZV z+m>{HN88SA`BD9?^;fl)pKiJ_F)T}X?VCl>}q=75fmhG6!K*2X#~HXp?}P;z3Tx@l<|)h=r43uPTFrH$w{Nh0sHmk2nN2v7DY2Q#X0$2_c3FstzshT@Dx2 zU*wCRm8LQVHc=gO3 zvo1$?=`LpqGeO`rqydthqxKh&QQ)K=3SuA@mUo*RV?eppDs}J_v_sQ-`n^!zB zev(w{?UPsE9b(F^eVqSITR7RS9Qke)xO#7 zOa2%?+p^K9>1E54Z`*Ed&^fQ#w14*$qm`EDr?cNa^S@$F{@dgHzb^cmT^9beXIkEL zLF<`YoR3x1pFg@=%XMMX>en}?hUsnmaa1};_~NGxXWeFm%N_KSduL^OZS&!^tL+s| ziyM2|-+8rX@5F-wWjk1ESMF?BCVbAr@XzXRa)iXI#n@y9K*lshdW57 z&zm1)d~fTK_VmBA@9}+j#;PLA*8b06Mp4D#9$C3Ji%kSy3eBpWkbPmtu9&y=8@<+E z+O+qx+TxsYh04l4-pI85fjLr-94~9GU$u44f}31ICNp#Erzsuj`L<4N`l{76w={%( zpWT%A{r1Olac;t{3w!HttSu&2 z+jSb(2}Ls59FO>2^!1I2ZqGWOYJJyyP=aHKWK!^E&}dbZmyrA4A^GrvWD4g6kIB=E z3s3O}xJ_fY_3wh9wd&1dHJqDT7%#|+Pb+aS-o%n2c>Gqx_rsBP40k_@s~=$asagCr U>c?b31_lNOPgg&ebxsLQ0QC!er~m)} diff --git a/doc/equations/dist_tutorial1.png b/doc/equations/dist_tutorial1.png index 80d54afe8882e4755290cd423ca77a810e2ff43e..172eaab0b275f585c1c1731b55299109867b75a4 100644 GIT binary patch delta 1879 zcmZ1_c$;s6WIY=L1H*E@ZMh5#3@qu6zK#qG>ra@ocD-a^V2~_vjVKAuPb(=;EJ|f? zOvz75Rq)JBOiv9;O-!jQJeg|4z`*v_)5S5Q;?|qHuk+>@2(&#^JLUWJU;UR&m)5Py z6*Nj{&OP>{a_K6|2}hUzv@tgO@%Z=e-|P%q9yA|dP%wC5ZBTD#SGRZnKV!GLef#SF zuIJNu=q;>}I>XyJMrxJSf);<@2HV3m=l{K8U}{dYIq1-NhBrx{y?F}zhV%*f{E4-9 zFWz6t+{jqgTQb+eH+K7UpSQUOw!WJ1OIk+xWcT3-Ntrp(s||O$ZT})w)bcRvMDSuh z&5epiENr5!%eDtz-FJOzy_LY^&mNx^-6#%V7mx@szHrg}=G8ULA69JGt5P}9Xu;Yg zX<6lO?tL$p$XsnDAkNg>=#q2ANLY%gnbB!_V2{bv%ZdgpZaKUH5|5Y!B=o9Hzt39D z%XGMz(IKJZbmV)^T!To#8F@^~E=-Mz20bzg25*nFUSMEZ7Wu;|=GT_bIrV4k)~!); zl;C1we$@2w!!EWKp~Wj6hCJCMD521NVjiD<&jAJngBhAn3_lA?Fen%>aPSCd*=+jd zoX}A3BoHW|G+oz-FK+ggj{7cNhoY9L8B_$^xRhnKt%`-I`Q7FsZb$E{ily_F&#w)? zxoq|on{Jczjw_EJY+f->`j_HK>xCQY*ZAwl*riYYe5Th$or9+<_UELTb>1m`c8_%v z0>is(u0)7i&fwfTW37OudrZ)CyZM}-WlS&*P(TyhmZhg$mI5ZlI;ItFpk z41bPGzXi!MSABI-c^G)*+nui0hJN>JGWA@!-RGDnz?bxC&OSwhxv{6de?Q%P z$=xc*bo-j6rF$0hy1koT`qS&wt+7%o`};8b-1*}5(i`)_U(fheo4NZ*-2aw3%gbryUwQ>R&nc$JD$ki> z;jOnK^}fIizl9es_1YF$PcJNb;<)qqs^?Q`U+{D7LanFU!nTh(5OCUtRt+-L;QHn|IvO z|J8P&r6jYMJtVK)eO9dZ>k6~-fX{NOsd{&7%$)C^IB3%Ia{ZYY|BU{dr)<1dl(^p) zm@(V>XXot){jYX>I}{|nuSz{Q{l)ZKwaO=@EmXYD*1b+i(%M;ZJi+H~{d%A3QlHq# zMVHqftN&O2)qDNETUyS36=h5Fz1~01eAT8Hv_5vlMd!mA-!JI@eW&j7%0l~`!m7>R zyjyx31%J1)UtZ{(UH0ehd(W4a9ntq!UwfXPYf-anM~dxlmxO0?Pj1zZ`ra$R>9DzA z0(){B@9k~nZA~kL`GjV0GdU$R)b8`Tvg-S%>-XMEni(9ZSK0gYV8zFxf5OWd{3;A= z+RW~We{h*3{4tJC?fX>*f1ZicZXc>Na5yTEsc#cyzxYGR`?doN&o%F{8~>8t%kbr~ zwz|vBsz*gp+xvGgJ$GIb^LW<53)}8pw0hlj%gCnh65GVn)2tKypC6neCHw5-4W+e{ zdJ7V6or$*dv+e$pC6iu%?2I(ux0;1PM>buR2;Q8N`eV)Amou7Aw8^Z^G4>Cdzm!*1 zXO_RmA~~;bOPDXqy=JVwu~EloMkQE> z<8~!}&t!4F?fC(2N}t8wH!JG#ZC`)Jr~S6@j6;t32}jIR+HZ4Dl(Rqmea6vh`QIDM zbF$8r8c*e&nH)3edxTF~=StZ;`D&*}2OkI~UcBT`6qZ|5#JBSP`uehx{!5#ts>Nr^ z6!^z*_{_Me;ncB3UBXy4FMjdfPVPod(}Vyeo?|E63@05Fnr;#{kw>Eawz-6Mm~y4T zkyPdOM+LusJ$$!iTZ$;BzgTMfZD9!;jUz=nY&m!u^!|t?F#NBulV>yFPG(?WVDNPH Kb6Mw<&;$TXU15m; literal 2218 zcmeAS@N?(olHy`uVBq!ia0y~yV3c5BVCdmsW?*2LaKSB}fq{Xuz$3Dlfr0M`2s2LA z=96Y%U|>mi^mSxl*f=lI$MHJ@1A{`cN02WALzOB6LqjtI!_WT=3=J3JMB}ii%1~O3KR0Dk>_fs;X*gYU=9h8X6j!nwnZ#TH4y$IyySKy1IIL zdiwhM1_lO(hK5E)M#jd*CMG7Prlw|QX6EMR78VwkmX=mlR@TUEST?Jv}|Wy}f;Xef|CY6DCZUIC0{nNs}f|o;+pBl&Mpv zPMbDu`t<2DX3Us5bLOmBvu4kpJ!j6GxpU{vn>TO%{P_zOELga3;i5&07B61BWXY1H zOP4NNwru(G4q*RET)ZvFc88#Zj%xN+mAO`A4v-n?bY zmaSX2ZriqP`}XZScI?=>bLXyIyLRv1y=Tv!y?giW+qZB3{{06I95{IJ;Gsi@4j(>z zIC=8qsZ*y;pFVx&%$c)i&z?JX?)>@l7cN}5c=6(;OP4NR zzI^4%m8(~;Ub}Yf`t|EKZrr$e^X9Evw{G9Qedo@dyLa#2yLa#Y{re9dJb3u<;iE^7 z9zTBkt5>gHzkdDZ&6~Gx-@beI?*04sA3l8e`0?YX zPoF-2{`}?3m#<&He*5{pPZa-&kP_u`PiTTj=wML6Yjvxng#g)PSO4>KvT9JDaKU1+#- z`Jt~&JTH=0>|1bs#gT%`JmN#_E@lZF}=8zw0QTdw^2O1%v7aQFU2#d9qX8MnT?4}*XO+Et#b{TaTz;~Sr|Ao z1Sfcw^lILDxI`nY__n$p*CyXXU+>(A%P>~kx;&%sdUwL#`L6d@EVA>hT$*C;+pif7vW_J{OhepoV z_IxP5?M#A8Fjvprj63%?zTp47rR3q3e+{1O$KEUJYfY_sbJ(Phfv0%cW)0>~Y@27g zF~~4TEcw0Y2TPP@ee$}FCl`Mk7<30c+|p-s|MP>{ewp5rV>_4JxVMct_e(s_j%M}A zAHK5q-1M6gV3xc-`?81fB->uYPs`@079_|>upetIM5c}xW_<9v`~8A5=L`mO)$fMy zl38MUnp8i^{_Z()CU%eAqp7;*BYfnSwR}J7Gb>BdRjwJ_~`V^=DMxs{7m~*#_DN7jw(jijkJGkvgn99yW!&dTK===G^cE1Ql2zP-1z2K z{u_&YRrRL`7{8Rf^zY5aMZ5oK{eEAdm8KV$)((77XPik|qJc3bI`r}saNQ#TLWG%uTtjiu%M z{KY{kOFlUH-0_glxO7C`-V$4N4;DW^7`+N z{lAm#-lf`BC-48ea$$quwYGhpQ7#7>Dr469WPUgH?^$HXVz*F>Rlv(ZF3ab5uC)Bd z&v*Z`+$yR2C1-GazJoc#(hsK{j>-$Hdc6Dj-L1Fwesa3EW=7oVC%d<+GN!zkd7P{L z^pE=Y%ub^Fjpya0$~U!hhwfYH=@!qTp!xW$+=Xl*yJe}WJD;DOBzd(i>Rq(+1_qzf zh+Ci1ZhTR6aebE-^(3i(iD;rrGRqsz=jW2XA60mCWb%nRWw*S=tIW(r^!;C(+&IS< zrPuoDap?#%yC#n&@?67%C_PY&$RgmC%ImtIP2P~G7A>^9>uiJAGwjz=heR&5l# zQM%Rn+3DFr^K1EQ-u_D7=D4X~mhSV}!7(Cx5`$CrEnE@e{Y?I5O;nec=Atbg#j8&I zIlR@C!{V)cC) b_~ZW!PSs9HIwG}BpeBN+tDnm{r-UW|PdpYi diff --git a/doc/equations/dist_tutorial2.png b/doc/equations/dist_tutorial2.png index 8b2bf93486089cfb3cacebc6bf95900d609d77ec..3bf4ff04891fad162f7f3790c04468fe07ea3470 100644 GIT binary patch delta 805 zcmbQkwvlauWIY=L1B0Y?%`XN929|V3Uq=Rp^(V|(yIz8XOI#yLg7ec#$`gxH85~pc zlTsBta}(23gHjVyDhp4h+AuIM9rko_45_&FW^Qa>vjI=bbkXUze%pVE$obtMvXt}Z z*5t$HLbr0ShFz{m>UY+-T_|JE@a55TMFtKJhTGfro_kzh!}C7>jfbjiOL^vhZu!Hr zd9H}I@IN~K!Oh|92ZcqVdjtcL4mk5qoc(s;{f@uKCs*We`FK55ao*8~4)f0+Uw^t% zaIKEn?C+xGJ|(Wd_}RbBd8u3XxW+&CfxGu5o_m^xOXIp5q)s)Lu+9AYQ+ZGJ{bu%8 zCzw9CIrtvDaKdg{Wqf@tACp#ziDh!+@9Wm@Cvm(yB(k}A`_?5Biu>%Dg7kFbT&A2_ z5Xuk+7#3aJS! zn6e_M@xrYJ{ZlVD%XcjPU-3A!%I^DruXUR(0uQyM|GjIesBkdvIp>@mYtlD$>TNx( zn-;dZRsLa1LyOghYY(ELrsi3SyVY|VKk#0&F8OMO8VgGsOK$^1tVG1~j{+B*m}}eL zw>QjTkBqXJk)f`%VH3m5-*4+MPc!w({BZIY{NIFCa0x63UPDf z(B3w`>J^J-VuO&H>E_+}4oPp496jw)7Mz@W@haQ5sgoVM_?e3CEZSbTk42uvO=u<5X=vX$A%cmUKs7M+SzC^Add= zzcVl}C?tCX`7$t6sWLD$G&3;#{LjG9@PdIMh=YMqG=hP_YB~dhc>eV$`3VdROp`oa z978H@y@~PUJFLLt6Y%eUYE_4nB1g2rgVdT8yOwNy{jr{{%wq3%b`L(w6$}lz%Rc*m zFHW5l^Xo!Y)1R^}X&#TZEGm^be(Z~lztB_59)C|))y`nH^OKE(Cf85$nvrnCvuKsh z>-YcXFPrdt?uW}JmoI3gxJ(L9T<5Ips)+r3t4U$ zgsy^!XL+F>k;q;D0&Y2(TV;)ZoEt~Gk zBegaD-8Fwvo#T9p69eaEC$b-yaJXgpitC}gp*Aw<&wCa>-P<1$JkR>}xpOyJR@zpc zx$#RfG2KdK+2f{(voe;Kc-)9QpSVyc__f-4q2o6fdad!vC{9k_`$FLOxjPyTKh-7% zEPpafF)^Kc13&ApM`^*P^A+1`^$*rRU(3y1$=ASX%*UYC!!Tn*%rlwnu@82un zda-fBi3GvRh1aefvJ4c8<2wKHg=_gu4}DJk78kGYSwH7G#;Gx!zRDE%f#>{l)p+TE zkM%`ej<2pJfBPV*^ygNk=nX}cW83bVZuK~QI*)nt-)ACWa{X*+ljRoke)^yh*S2Sk z|G#H(=Oh$!pQr9>)O|6H_V+Po~;1FtF-*x;TbZ+J#f~|JMEa-S1W}Xsx6D zCA~4feqVQeAos-!T9bsUHz-EM9&fy%ICbg6o7|>x{YTc8D^+ekr`ziHhO_maa-nqK zyw$0-b@i*O>Lw@5duZXLA|E*Kc*ezB?}e>-UQ47~*KGGYa_mh1EBAfhKmOjk|7~Vs z)wJ|wUyY=NmngsCTu@qnX!=H*Bl}i_@0XISKA-5W&3}I1is^TS{nj7JO{rX0QJ$o) z@MCd?rI*6RQ-_b~WF#uxY{+=3_|4VIW$ULiuRNB?sk?tGzg)9e_wb5OXFj|=B6s9i zjA*jR{+jED&kI;fPd1p@du^ut`WJHZZwR}V%J$z-%ygHu6%w$XB6+ey_Edd@S(3im zsd-M*7dVL@bQ8E%9BH&>+b!wI1{;ILCdxXbtbee$-;aOq@4GYC?pf_pn-qV8>ka46 zyPw12IFgeadK)dn)J`9DSF*k%bnPNz`cuWSf3+8$uUD%`;#c_amgmH0rL!+sXMM^66t2b({* zm33u2%?#6LU-!Ff+4s3Uk*`l5Yv}WRV>;=IXiMF)eVeV=Hmdt&$v)V@R-HVlPvVr4 zBn!h!=e}F+m{3*us_cknyue`80O6PW@JZSukXu+tY5r2!sqa2LxrEq zn}DzvNB9z3^c-XNEOa+5`tqy)(Sb#mR6^!gnMgA$YQ0 z*wgK{UoCbmi`c|0Xx$>Zi)E{W)AJ8)how&+e%Sl*>+}U|U#DN3bXn9QA$b8$$KOX6 zUR%A|BO_paI;4#`N&N2V0*5!Q^^2=YHkL(hG?V^3Ga@ZVvE$R692tS#7Z>yBXsOCR zyI+|l^XAP$@tdXhuFX}^VLx*2SV-ed35BvB))7j)Gmz?Sd&!98&THTDR|iEx1(4^}oNxkCYUXuNIl#5=*nMHaE|EYV-N#5x&FLPu3i8 yd}D5ttg_$svD|e3SN4;cqiYuO{1Z-K*k=@J#d(GAe>npK1B0ilpUXO@geCyQ>bR-^ delta 1351 zcmbQkeV%KAWIYoD1H&Ev{$K_M2F?PH$YKTtz9S&aI8~cZnt_3VCEd~2k%3|3yhIF!Z|2^f^jLw%W%88&|EKTeU+I~e%g($@;ljnsK5F&CcP@D7>Rl|p|DC-p^me!2 zfmu0=xi{VX{Ill&1@X89AFEw89iRRh2XgK<+nLj9XufVXdz)cIy@%8B7hjeIc&~bK zt@i$F`Q5Yien`mfdQ+b!D|IS2)b0R#!u`vSl9xCyJ*r*BMY9X zG2W}Tr%KI8@VV_UZ~2aQ2P(F?hHT>X%nLlQSyv+2X_@5hO^am%SC%Q<_j^&TqO+wb zJ=)-ymfM|uMLScHZe4Tmbq?9&+otX5dkbswvXeR4cIrl%31|L3s9kO_|Eh6!s>%6n z2|dD@eCgZHFMfSZY&T2b9*?utS$y@UGi}d5%{?*8In+&hn!z%&O3q8y3ObpSSFwHC zko>1Mq59soUp-%EY_t6PclEE5{H*LNnvykV*5vgYS6n=r@32c`-p#%s1GnHUX%-=^ zeTEYixJ_F6Bee{8Lo5y0nhkf}?J1RFO*W92eB-Ren#%_yc#fr%8ryLnG)NP?wBxOG zz0FFUuFe{t70gW$4!iHGUMMQx-_bF7_H4s4MVYjekA{~y9oVjxY|GE%x}&c3@TA)> ziMuVm8)N1e3fHD*p4niL|9+}s&IQ-)>sv)%KX$&sA-`gU8pj)Wm8|`&fB*UJoHE~aEv*}R z?9#fNo+K)!r1-Q|XOy{x9Qt6W8&`Yv_m>|MPxtl3ci8B(yh;Diz>vJE%Cb)J(*}V_ zza-8tdVW~(#iA3^Q*FBE@i!jm)cUKFoXhb>NaK6O2ast`_SOq){`t;PuGBjJ!uOMF zey)4(YZ=&TdYSd<>*acufu#n1A(hXiHrxNt{uwy`m-R6wCk3rUgJY?E%UT=yIb|m? zIe8qExWi+Y#>#Ho#puKlDe-yU;zN#7()~>i3;KkrR-E0Y>aH{;KbnWhPK|4Jx+F&# zXWQB9o{NQc3b7^s>u0JzH_Kt~&J2N|H$pZ!eplo9mLCm@*5J@Ke)eJRU$I4wpJtl= zmlkl+lYP9{#5ep)yYq=1>yM~b#cBcT`Yc-!bmn*PllNUu&dtv8;#N zP8?cm#qC;Zmzs*sz8m5hsCYMcO2iwFaJxy@xO^pra@ocD-a^V2~_vjVKAuPb(=;EJ|f? zOvz75Rq)JBOiv9;O-!jQJeg|4z`(@e>EaktaqCS?G#|4f#~UgAtH0~NMBH}Y)Oqt& z*8=TB6TU3=jb?R|m{4~6`(MV!B|-^ICnA-665O}tRxdqoUmy8^Gr)q~g{#TfVS=5) z8;y_5A|1>2eo$VnRi!BKB+i6t|qDV()XnA+7oqQ?7_onry+RV9uflpV-yF_lDdciGw z`JVQR*=}=8-7aier}A~Dzx?7f^QzkXl}-<==6>3larNXCVOuS!o_*p{ng^v6`X;CL z8h8b}ciSH7I<|B6i?gekcBE$9TK+Kc28TVHQJCKM9V<-kUffc}X%@9UwZ7Ev}`Z{n@6|LW3SA{QQYgr3mI(3R2di=ni&{={%2rlc)`FB#KFKQ8o|I|HJyP$JpX!>`~(ID#_OIg zjv*Dd-poGOcUXbP&16Jyi& zbYK45e@B9ylOt>^xU}U`KYBcW$Yp2|d85hZP7%u)@t%%(^{Iz--PcC!zT#n{WZu8s zWNUce(}M-&CfUq;>us2xA8vCj%A2s{WX6gq4qWOV7hIYeIlJ(K1K%ejo3yvbk1xJ# zExE1tLP_9<4cUKhTRf`Do1q=^WIj*u`SgcUm#_8TI`=85@lNu2>-m0ig3c|UJ$hYp zzuWG;)M@kR|A~U^W~F@p-@fAc{8wh&eE&P&Jos|!^$!_s7I&Cs-~Do^H!rd8oqxRX zpvRNddjvo2kem4ZW5L=97ujStU0Ev~HQ)UPyWRQ4ui`T*E3~C&A1^fx$cf1}jZON@ zvcABsYKqO7#m{EtsK(DfZmQQPeXJm_!dXgxq0XN(HWG4iU(Jie9vW^pI(~j#E$h`O zJ+qD<`!$CvKhs>7b4K(dLx&q|X5W3{4mAJVB(=}^oBNVK25j{j#{C8Voo$Y_KL`oD z{r8AzLcr%463b>zS9!RdBWdX;CIhp-|C)3E@H6bud-%jxthbGUfq}u()z4*}Q$iB} D$q5fL diff --git a/doc/equations/hazard.png b/doc/equations/hazard.png index 7a5978611c9374bad55fb93c3298cc054f7da0e7..b65aebd46b1b476e91e0d02d3ab100c8e609c896 100644 GIT binary patch delta 761 zcmeBR+r>6Pl7o$bfkE{i&)JEJ7WHR5T^vIyZoQd%F|SpDr)~L7@3(*TzeS|)HdJhJ zYGe}Krysegc*;wyznL3feEnIokAY#$P62TSiG#hI40GHWu78!TudBP#;W3Ar?6k= zq8g6uc3pjCoqMC{^EV&=&$3O53d@_9iiY%`NDDOpy!oAJ@LkedeF$>@$@$ zDivCp6JP=m#zv3ekU_6biL2>Fw1IRPK$3W#f6=xCu}{vFZ))@!bQt? zcDZdVHvY7`zAwyFu;<#c8^ULfy6&2=)5BFrcxA82qp9!yUzWVhy7ai(jbL%diQh`k z>^!pC%yeh%y0a_i7A)KROv&?2=Lx_0#TU;N_hdXyTyk%v@V7bK-&&&bxYf^fc5HJp zxY`+YdcjWZN!#7u&f)&XVqDlMYx(rX$0u$V&08`QrXzNOUQ8?@ohr9nfOKT;9%JILpXsU Zr+n+Yr&5V885kHCJYD@<);T3K0RRumciR8} delta 711 zcmdnR*1d_Mo`w||bx zAA?%U?>wIJR&2-SHCYpHeVdc&?N*X?dw*BvTF?63&-)75r@vZ!Yeu=$Hs2-Uto|oo zuM+xxqjauC>_#cus<_Rk-u8X@a`&$8w`a*!!8gB37ak0#T=0F>(H^C+94)mAciBuo zwTW-ME7A9(=A~%T3*9Fk=87*aT=$>p8nwb~-`5#adbw9`TYpIDZ`9t8?#Axrmp7+n zuqU)m+;?(vYG1wN?X_w1+}}kD)UB62bUXK+;ktJR73VIkyv^&rN4fmX@9F2>KYO2e z+~T#{v7GH)e_Q>lErst(izQE6_jtpWMWR|#;Tvjo#dhZD#@x{@KDvJHmvWU7b6=J4 z;}a)kUJ1?Rab0ya$GlCD{qf;9g16g^HoI~OcqYYjuXhVs{4FrIr=IikHv3aZn5#rbLly^>nsb;-(I>!oNN8JiSggNozptDta)q9Br7{(-z1*a>GLNsR+Oec zJ)|AM@@(0J@EO;B20L!x;E{?h+I(}uxzAsd+^RnDbRRG=eYoMlZf=Gj!biU-uyr>a zOl>JG<@j*&1%t+s`3%eq ZS8Dri%ce~1VPIfj@O1TaS?83{1OQYtSRDWW diff --git a/doc/equations/students_t_dist.png b/doc/equations/students_t_dist.png index 5d4657d8c4e667f2ea3a082205331f38ad2552f4..3a0a239ae6d2cb6b241c882642729a27075e0270 100644 GIT binary patch delta 450 zcmZ3%d69X7WIY=L149w3eLMpL153K2uOkD)`V;1?T`w6J7$i$vBT9nv(@M${i&7aJ zQ}UBi6+Ckj(^G>|6H_V+Po~;1Ffb;2x;TbZ+PZfx9UOMIa_aBN+>Nl64y;sV*<(R?dbd!vw zZugX}tDc*@FQ1aDeB$6C^lkYc_5HW$yBC!Pml!g6SPU*+_c)!&78KJV6t#GL zP*2R0aJ6~TnPNhcmD=^8K0P#Jd6L3iprF0+bkaGPXkm3mmY zP4v!F4cVR-J-%PGd?jY4ST|2RmR7ZJ>h*}3+vPhN8A=O+(z&A?wVBo~2$x+eF*D&c zW0a$|)#Tpl!^Lm?Gf#N1zm>dtF7u(d^3>TwR+8x zn}J5f-^4{zmK3w^{`;6qnK|LwuluL=R~wp~U$~utnL%rwwA#cwDgPN57#KWV{an^L HB{Ts5Gts^4 delta 1183 zcmcb}yn=IrWIZzj1B2`B^SKNR44efXk;M!Qd`Cc-ajG_-Gy?+zOS+@4BLl<6d5J!b z-x(Mf6p}rHd>I(3R2di=ni&{={%2rlc)`FB#KFKQ8o|I|HJyP$JpX!>`~(IDW`+Qt z5LX5;U}9ooW@ct#VPR!uWn*JwXJ_Z&;Naxstmoq5;^yY&;o;%s<>lk!)z#C})7RHGFfcGQG&C|YGB!3g zsW&k(H8nLeGcz|gx3I9Vw6wIcva+_ewz09XwY9agv$MCicW`iUbaZrba&mTdc5!iW zb#--fb8~lh_wexW^z`)d^78if_VMxY_4W1h^Yi!j4+sbd3=9kk3JMMm4habf4Gj$o z3kweqkBEqfjEszmii(boj){qhjg5_qtB;G1k55QQNK8yjN=iyjPEJWlNli^nOG`^n zPtVB6$jr>l%F4>l&d$ln$<58p%gf8p&o3w_C@d^2Dk>^2E-oo4DJ?B6D=RB6FR!Sm zsI084s;a84uCA%6sjaQ8tE;Q8uWx8*Xl!h3YHDh3ZfDz7Yinz7Z|~^nsPF9T z?CR?3?(XjC>FMq5?d$97@9&>5VZy|T6DLiYG|;B zXU&>5d-m)(bLPyQJ9pl^dGqJbU$9`o!i5VLEn2jA@!};*mMmSmblI|H%a<=-v0}x_ zl`B`RTD5xh>NRWDtX;cy-MV$_*RS8OVZ+9a8#is*w7Gus<}F*cY~8wb+qP}nw{PFE zW5>>&J9q8cwR`vOJ$v@--Me?+zJ2@m?>}(hz`=tD4;?ym`0(K)M~)mldi2<_W5~S3{Qdj)-@kwV|NqxbT4>I|z`#@zF!cV>9uME_Cf(X~p!3MOeal@mO0J!i%3;&}x93v1QRfZc z8G99$wk7xfn#%^@#3)ND!y}DmQHEv?%-Z>iF@8V;hOrdXD61%zf{>4uibt; z?q1GCe}Pg%j;CuJ?<7U;Ia(ODQOVHLEoo2S?F7!G^~px7Di`K>>@`?hE@Yb#_58MV zXP>NIXHQ-J&0x8W>TNye1+=pg?@Jdv+H>U96&e1OmAXMeI;9$m4?fCrJnUQlb3wC^ sh?;!jvAili$EeWV-jn}J@V_@_2>xT+qVchI1E{3%boFyt=akR{0M@l8?f?J) diff --git a/doc/erf.qbk b/doc/erf.qbk index ccf337f81..82b4dda35 100644 --- a/doc/erf.qbk +++ b/doc/erf.qbk @@ -70,15 +70,6 @@ versions of these functions use differing implementations internally, so this gives us reasonably independent test data. Using our test data to test other "known good" implementations also provides an additional sanity check. -One should note that our tests rely on decimal to binary conversion of -floating-point numbers, and assume that the result will be correctly rounded. -In practice, it appears that in a few very rare cases the test data may be incorrect in the last bit: -this depends upon the compiler and standard library used, and means that the -relative errors quoted above have to treated somewhat circumspectly. Using -binary or hexadecimal coded test data would remove this issue (or at least -confirm whether it is actually an issue or not), but would make the test data -unportable, so is not used at present. - [h4 Implementation] All versions of these functions first use the usual reflection formulas diff --git a/include/boost/math/special_functions/digamma.hpp b/include/boost/math/special_functions/digamma.hpp index 0a925890e..a22f80c72 100644 --- a/include/boost/math/special_functions/digamma.hpp +++ b/include/boost/math/special_functions/digamma.hpp @@ -158,8 +158,8 @@ T digamma_imp_1_2(T x, const mpl::int_<0>*) static const T root1 = 1569415565.0 / 1073741824uL; static const T root2 = (381566830.0 / 1073741824uL) / 1073741824uL; static const T root3 = ((111616537.0 / 1073741824uL) / 1073741824uL) / 1073741824uL; - static const T root4 = (((503992063.0 / 1073741824uL) / 1073741824uL) / 1073741824uL) / 1073741824uL; - static const T root5 = 0.75231638452626400509999138382223723380394592598054e-36L; + static const T root4 = (((503992070.0 / 1073741824uL) / 1073741824uL) / 1073741824uL) / 1073741824uL; + static const T root5 = 0.52112228569249997894452490385577338504019838794544e-36L; static const T P[] = { 0.25479851061131551526977464225335883769L, diff --git a/minimax/f.cpp b/minimax/f.cpp index ed59131ce..862fd97c6 100644 --- a/minimax/f.cpp +++ b/minimax/f.cpp @@ -122,6 +122,10 @@ NTL::RR f(const NTL::RR& x, int variant) } case 10: { + // + // Digamma over the interval [1,2], set x-offset to 1 + // and optimise for absolute error over [0,1]. + // int current_precision = NTL::RR::precision(); if(current_precision < 1000) NTL::RR::SetPrecision(1000); @@ -129,11 +133,16 @@ NTL::RR f(const NTL::RR& x, int variant) // This value for the root of digamma is calculated using our // differentiated lanczos approximation. It agrees with Cody // to ~ 25 digits and to Morris to 35 digits. See: - // TOLMS ALGORITHM 708 (Didonato and Morris). + // TOMS ALGORITHM 708 (Didonato and Morris). // and Math. Comp. 27, 123-127 (1973) by Cody, Strecok and Thacher. // //NTL::RR root = boost::lexical_cast("1.4616321449683623412626595423257213234331845807102825466429633351908372838889871"); - + // + // Actually better to calculate the root on the fly, it appears to be more + // accurate: convergence is easier with the 1000-bit value, the approximation + // produced agrees with functions.mathworld.com values to 35 digits even quite + // near the root. + // static boost::math::tools::eps_tolerance tol(1000); static boost::uintmax_t max_iter = 1000; static const NTL::RR root = boost::math::tools::bracket_and_solve_root(&boost::math::digamma, NTL::RR(1.4), NTL::RR(1.5), true, tol, max_iter).first; @@ -142,6 +151,12 @@ NTL::RR f(const NTL::RR& x, int variant) double lim = 1e-65; if(fabs(x2 - root) < lim) { + // + // This is a problem area: + // x2-root suffers cancellation error, so does digamma. + // That gets compounded again when Remez calculates the error + // function. This cludge seems to stop the worst of the problems: + // static const NTL::RR a = boost::math::digamma(root - lim) / -lim; static const NTL::RR b = boost::math::digamma(root + lim) / lim; NTL::RR fract = (x2 - root + lim) / (2*lim); diff --git a/test/test_digamma.cpp b/test/test_digamma.cpp index 3f7dd4578..ecf244de0 100644 --- a/test/test_digamma.cpp +++ b/test/test_digamma.cpp @@ -42,14 +42,6 @@ void expected_results() // Define the max and mean errors expected for // various compilers and platforms. // - add_expected_result( - ".*aCC.*", // compiler - ".*", // stdlib - ".*", // platform - ".*", // test type(s) - ".*Root.*", // test data group - ".*", 220, 40); // test function - add_expected_result( ".*", // compiler ".*", // stdlib @@ -137,12 +129,11 @@ void test_spots(T, const char* t) // Special tolerance (200eps) for when we're very near the root, // and T has more than 64-bits in it's mantissa: // - T tolerance2 = boost::math::tools::digits() > 64 ? boost::math::tools::epsilon() * 20000 : tolerance; BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast(0.125)), static_cast(-8.3884926632958548678027429230863430000514460424495L), tolerance); BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast(0.5)), static_cast(-1.9635100260214234794409763329987555671931596046604L), tolerance); BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast(1)), static_cast(-0.57721566490153286060651209008240243104215933593992L), tolerance); BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast(1.5)), static_cast(0.036489973978576520559023667001244432806840395339566L), tolerance); - BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast(1.5) - static_cast(1)/32), static_cast(0.00686541147073577672813890866512415766586241385896200579891429L), tolerance2); + BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast(1.5) - static_cast(1)/32), static_cast(0.00686541147073577672813890866512415766586241385896200579891429L), tolerance); BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast(2)), static_cast(0.42278433509846713939348790991759756895784066406008L), tolerance); BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast(8)), static_cast(2.0156414779556099965363450527747404261006978069172L), tolerance); BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast(12)), static_cast(2.4426616799758120167383652547949424463027180089374L), tolerance);