!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12 program black ! implicit real(kind=kind(1.d0)) (a-h,o-z) ! implicit real(kind=selected_real_kind(16)) (a-h,o-z) ! Works. implicit real(kind=selected_real_kind(32)) (a-h,o-z) ! Works. ! implicit real(kind=SELECTED_REAL_KIND(64)) (a-h,o-z) ! Don't Work ! See https://software.intel.com/en-us/forums/intel-fortran-compiler-for-linux-and-mac-os-x/topic/269257 ! implicit real(kind=16 (a-h,o-z) ! Doesn't work. ! implicit selected_real_kind(2*precision(1.0_dp)) (a-h,o-z) ! Doesn't work. ! implicit real(kind=32) (a-h,o-z) ! Does not work. ! implicit real_kind_precision24 (a-h,o-z) ! Does not work. ! implicit real(kind=kind(dp)) (a-h,o-z) ! Accepts but what does it do? ! implicit real(kind=kind(qp)) (a-h,o-z) ! Accepts but what does it do? ! integer, parameter :: nprecision=kind(1.d0) ! integer, parameter :: np = selected_real_kind(16) ! Works. integer, parameter :: np = selected_real_kind(32) ! Works. ! integer, parameter :: np = selected_real_kind(64) ! Don't Works. ! integer, parameter :: nprecision=kind(32) ! integer, parameter :: nprecision=kind(dp) ! integer, parameter :: nprecision=kind(qp) integer, parameter :: nfa=16 integer, parameter :: nprime=8 integer, parameter :: npz=20 real (kind=np) :: fa(0:nfa) integer :: iiprime(nprime) real (kind=np) :: pz(npz) data pz/1.e-6_np,1.e-5_np, & & 0.0001_np, 0.001_np, 0.01_np, 0.03_np, 1.1_np, & & 1.15_np, & & 1.2_np, & & 1.3_np, 1.4_np, 1.5_np, 1.6_np, & & 1.7_np, 1.8_np, 1.9_np, 2.0_np, 3.0_np, & & 4.0_np, 5.0_np/ data iiprime/1,2,3,5, 7,11,13,17/ ! From https://en.wikipedia.org/wiki/Factorial data fa/1.0_np, 1.0_np, 2.0_np, 6.0_np, 24.0_np, & & 120.0_np, 720.0_np, 5040.0_np, 40320.0_np, & & 362880.0_np, 3628800.0_np, 39916800.0_np, 479001600.0_np, & & 6227020800_np, 87178291200_np, 1307674368000_np, & & 20922789888000_np/ ! print* print*,'Blackbody radiation distribution and integrals' print*,'Related to Debye function integrals' print*,'which are effectively, Riemann Zeta functions' print*,'See Arfken p. 332' print*,'Max of integrands exponent p=1.01 to 5' write(*,"(1x,'iteration',13x,'iterate',6x,'Newton-Raphson')") ! https://en.wiktionary.org/wiki/iterate the noun con2=1.5936242600400403_np ! Numerical solution for p=2 ! iexact=0 if(iexact .eq. 1) then pde=1.0e-6_np pe=1.0_np+pde ! 23456789a123456789b123456789c123456789d123456789e123456789f123456789g12 xe1=1.99999933333377777745181382492999647E-0006_np ! what np=32 gives ! 1.99999933333377777745181382501483924E-0006 ! After 1 000 000 iterations converged ! 1.99999933333377777745181382499397862E-0006 ! After converged ! 1.99999933333377777745183802478120206E-06 ! one less that total ! 1.99999933333377777745183802479528595E-0006 ! converged to 29 digits ! xe1=1.9933774e-2_np ! Yes, they do iterate xe2=xe1 xe3=xe1 ii=1 000 000 do i=1,ii xe1=pe*(1.0_np-exp(-xe1)) xe2=pe*(exp(xe2)-1.0_np)/exp(xe2) sum=0.0_np do j=nfa,1,-1 sum=sum+(-xe3)**j/fa(j) end do xe3=-pe*sum if(i .gt. ii-11) write(*,910) i,xe1,xe2,xe3 end do 910 format(i7,3e45.36) print*,'pe,xe1,xe2,xe3' print*,pe,xe1,xe2,xe3 stop end if ! print*,'The expansion terms' a1=2.0_np ! determined analytically 1st a2=-2.0_np/3.0_np ! determined analytically 2nd a3=4.0_np/9.0_np ! determined numerically 3rd pd=1.e-5 to relerr 1e-5 a4=-1.0_np/3.0_np ! Approximately right, fiducial value ! a4=-0.325963952847781078141194663752917956_np ! Very small delz fit but not right a5=1.0_np/3.0_np ! Fiducial value at z=1.1 ! a5=0.0_np ! a5 can't be fit unless, a4 were known exactly, without more cleverness. ! a5=1.0_np/4.0_np ! Consistent with a4=-44/135 ?, but it doesn't stablize xcoef6 ! a5=1.0_np ! determined numerically 5th this looks right at delz=.01 to 0.5 % ! a5=1.0_np/3.0_np ! determined numerically 5th this looks right at delz=.1 a6=-2.0_np/9.0_np ! Fiducial value at z=1.1. ! a6=0.0_np ! p = z xlow-x2 ! 1.1 2.43943495082172799567213935641107047E-0009 ! 1.2 1.00027505387241362096291602996416415E-0005 ! 1.3 8.70623203195004906480843155165127304E-0005 ! 1.4 3.31905192729496031072133295792039345E-0004 ! 1.5 7.82534201282920938385984169344795317E-0004 ! 1.6 1.18189225328027400815451880144153590E-0003 ! 1.7 6.07272745384688824701656951475041720E-0004 ! 1.8 -3.04593994617942485727758941976983294E-0003 ! 1.9 -1.37308194737776213527942068124790343E-0002 ! 2.0 -3.80687044844845367674863203196046635E-0002 probably on the brink of divergence ! 3.0 -6.82143937212207889340319133029448507 clearly diverging a7=2.d0/81.d0 ! Fiducial value at z=1.1 . I don't trust. Don't use. ! a7=0.0_np ! Can't determine it a8=0.0_np a9=0.0_np b1=1.0_np ! Determined by Bell polynomials and confirmed. See p. 19 b2=1.0_np b3=1.5_np b4=8.0_np/3.0_np b5=125.0_np/24.0_np b6=54.0_np/5.0_np ! Not confirmed isatz=2 ! Which interpolation formula to use. ! https://en.wikipedia.org/wiki/Ansatz icorrect=0 ! 0 for none, 1 for interpolation, 2 for Newton-Raphson itest=0 ! Test for the low series coefficients if(isatz .eq. 1) then ! Crude interpolation formula acon=2.0_np else if(isatz .eq. 2) then ! Gem formula acon=2.0_np/(1-exp(-1.0_np)) bcon=exp(2.0_np)*(1-exp(-1.0_np))/3.0_np bcon2=(exp(-1.0_np) & & -bcon*exp(-2.0_np))/(1.0_np-exp(-1.0_np)) bcon3=bcon2 & & +(-0.5_np*exp(-1.0_np)+bcon*2.0_np*exp(-2.0_np)) & & /(1.0_np-exp(-1.0_np)) acon2=(1.0_np/6.0_np)*acon**2 acon3=0.5_np*acon*(-1.0_np-bcon2) ccon=-exp(3.0_np)*(1.0_np-exp(-1.0_np)) & & *(a3-bcon3-acon2-bcon2-acon3) print*,'For gem formula' print*,'acon,bcon,ccon,b3' print*,acon,bcon,ccon,b3 x01=0.1937475579949906_np pd=0.1_np p=1.0_np+pd ccon2=-( (x01/(p*(1.0_np-exp(-acon*pd))))-1.0_np+exp(-p) & & +bcon*pd*exp(-2.0_np*p) ) & & /(pd**2*exp(-3.0_np*p)) print*,'ccon,ccon2' print*,ccon,ccon2 ! stop else if(isatz .eq. 3) then ! Series interpolation formula print*,'isatz=',isatz xtrans=1.85_np ! seems good place to sharp transition else if(isatz .eq. 4) then ! Series-smooth interpolation formula pcon=4.0_np ! smooth transition power efold=.9_np/log(2.0_np)**(1.0_np/pcon) ! smooth transition point else if(isatz .eq. 5) then ! The newer brilliant flop formula print*,'isatz = ',isatz acon=2.0_np/(1.0_np-exp(-1.0_np)) bcon=1.5569247568238684_np ccon=1.5442655706398485_np else if(isatz .eq. 9) then ! The new brilliant flop formula print*,'isatz = ',isatz else if(isatz .eq. 10) then ! The brilliant flop formula print*,'isatz = ',isatz end if do i=1,npz ! Run through the chosen p=z's. print* if(pz(i) .le. 0.11_np) then pd=pz(i) p=1.0_np+pd else p=pz(i) pd=p-1.0_np end if write(*,'(a,i4,a,e10.3,a,e10.3)') & & 'Iteration cycle ',i,' for z=p=',p,' delz=pd=',pd yexp=exp(-p) y=p*yexp xlow2=pd*(a1+pd*(a2)) xlow3=pd*(a1+pd*(a2+pd*(a3))) xlow4=pd*(a1+pd*(a2+pd*(a3+pd*(a4)))) xlow5=pd*(a1+pd*(a2+pd*(a3+pd*(a4+pd*(a5))))) xlow6=pd*(a1+pd*(a2+pd*(a3+pd*(a4+pd*(a5+pd*a6))) )) xlow8=pd*(a1+pd*(a2+pd*(a3+pd*(a4+pd*(a5+pd*a6 & & +pd*(a7+pd*a8))) ))) ! High precision form. xlow=xlow6 ! xhigh=p-(y+y**2+b3*y**3+b4*y**4+b5*y**5+b6*y**6) xhigh=p-y*(b1+y*(b2+y*(b3+y*(b4+y*(b5+y*b6))) )) ! High precision form ! if(isatz .eq. 1) then ! Crude interpolation formula x0=p*(1.0_np-exp(-acon*pd)) ! prec max = 4, double prec max =7, relerr max= 9.0014819399899054E-002 at 1.7 good ! prec max = 3, double prec max =6, relerr max= 4.5495036393012151E-003 at 1.6 good with NR correction else if(isatz .eq. 2) then ! Gem formula x0=p*(1.0_np-exp(-p)-bcon*pd*exp(-2.0_np*p) & ! prec max = 3, double prec max =6, relerr max= 6.5047149472154484E-003 at 1.7 good ! & -b3*pd**2*exp(-3.0_np*p) & ! prec max = 3, double prec max =6, relerr max= 7.3877224167747660E-004 at 1.8 much better ! prec max = 2, double prec max =5, relerr max= 3.4485709062780802E-007 at 1.7 excellent with NR correction & -ccon*pd**2*exp(-3.0_np*p) & ! prec max = 3, double prec max =6, relerr max= 5.7797367485568822E-004 at 1.8 slightly significantly better ! & -0.1_np*pd**2*exp(-3.0_np*p) & ! prec max = 3, double prec max =6, relerr max= -4.8571899523800627E-004 at 3.0 marginal improve with fit ! & -b4*pd**3*exp(-4.0_np*p) & ! prec max = 3, double prec max =6, relerr max= 6.8422439795701086E-004 at 1.4 marginally better ! relerr = 5.8364429522015239E-004 at 1.7 probably negligible. ! & -b4*pd**2*p*exp(-4.0_np*p) & ! prec max = 3, double prec max =6, relerr max= -2.9347868449807234E-003 at 1.5 definitely worse ! & -b4*pd**3*exp(-4.0_np*p)-b5*pd**4*exp(-5.0_np*p) & ! prec max = 3, double prec max =6, relerr max= -9.3720229008774933E-004 at 1.6 marginally worse ! & ) & & *(1.0_np-exp(-acon*pd)) else if(isatz .eq. 3) then ! Series interpolation formula if(p .le. xtrans) then ! Numerically clear best just to use a sharp transition for accuracy. x0=xlow else x0=xhigh end if ! prec max = 3, double prec max =6, relerr max= 6.9558650904086279E-003 at 1.9 good ! So the maximum is worse by a factor 10 than the magic formula. ! Maybe the magic formula is just semi-accidentally a good fit. else if(isatz .eq. 4) then ! Series smooth-transition interpolation formula w=exp(-(pd/efold)**pcon) x0=xlow*w+xhigh*(1.0_np-w) else if(isatz .eq. 5) then ! Newer briiliant flop formula. isatz = 5 y=pd*exp(-p) wlow=exp(-acon*pd) w=1.0_np-exp(-acon*pd) ! x0=1.0_np-y*(b1+y*(b2+y*(b3+y*(b4+y*(b5+y*b6))) ))/pd ! High precision form x0=1.0_np-y*(b1+y*(bcon+y*(ccon+y*(b4+y*(b5+y*b6))) ))/pd ! High precision form hybrid x0=p*x0*w ! prec max = 4, double prec max =9, relerr max= 1.8141849762850569E-002 at 1.4 not competitive ! prec max = 4, double prec max =6, relerr max= -1.1862074441575882E-003 at 1.6 hybrid, worse ! x0=x0*w+xlow*wlow ! prec max = 4, double prec max =6, relerr max= 1.4647559416957010E-002 at 1.6 marginal about the same else if(isatz .eq. 9) then ! New briiliant flop formula. isatz = 9 pd2=2.0_np*pd ! Never bad, but up to 10 % error. y2=pd2*exp(-pd2) ! Exactly right only at z=p=2 xpd2=pd2-pd*y2*(b1+y2*(b2+y2*(b3+y2*(b4+y2*(b5+y2*b6))) )) ! The gem formula was better even at x0=p*(1.0_np-exp(-xpd2)) ! z=p=2 where this is supposed to be ! print*,'pd2,y2,xpd2,x0' ! the exact large series value. ! print*,'pd2,y2,xpd2,x0,1.0_np' ! Well the large series isn't that good at 2. ! print*,pd2,y2,xpd2,x0,1.0_np ! stop else if(isatz .eq. 10) then ! Briiliant flop formula. if(p .lt. 2.0_np) then ! Spectacular error. w=pd**2 else w=1.0_np end if wlow=max(1.0_np-w,0.0_np) ! Should never be numerically negative. c0=p*w ! The low series coeficient is 0 and the high is z=1+delz=p. c1=-a1*wlow+b1*w ! The minus are because the it's natural to put a minus in front of the c2=-a2*wlow+b2*w ! whole nonzero order set of terms. c3=-a3*wlow+b3*w c3=-a4*wlow+b4*w c4=-a5*wlow+b5*w c6=-a6*wlow+b6*w s=pd*wlow+y*w ! The interpolation expansion parameter. x0=c0-s*(c1+s*(c2+s*(c3+s*(c4+s*(c5+s*c6)) ))) xlow=pd*(a1+pd*(a2+pd*(a3+pd*(a4+pd*(a5+pd*a6 & +pd*a7))) )) ! High precision form. ! a7=0 except for tests xhigh=p-y*(b1+y*(b2+y*(b3+y*(b4+y*(b5+y*b6))) )) ! High precision form. print*,'pd,y,s' print*,pd,y,s end if if(icorrect .eq. 0) then xcor=0.0_np else if(icorrect .eq. 1) then xcor=p*(1.0_np-exp(-x0)) ! Interation formula correction x0=xcor else if(icorrect .eq. 2) then ! xcor=p*(1.0_np-2.0_np*exp(-x0))/(1.0_np-p*exp(-x0)) ! Newton-Raphson correction, numerically rotten xcor=x0-( (x0-p*(1.0_np-exp(-x0)))/(1.0_np-p*exp(-x0)) ) ! Newton-Raphson equation x0=xcor end if print*,'xlow,xhigh,(xlow-xhigh)/xhigh,xcor,x0' print*,xlow,xhigh,(xlow-xhigh)/xhigh,xcor,x0 ! x1=x0 ! x1 for the iteration equation x2=x0 ! x2 for the Newton-Raphson jcon2a=0 ! Iterations to single precision convergence jcon2b=0 ! Iterations to machine precision convergence dcon1=1.d6 ! Not used at present, but for variable symmetry. dcon2=1.d6 ! Tests to see if Newton-Raphson has reached jitter stage ! ! implying machine precision convergence. j=0 write(*,990) j,p,x1,x2 ! Initial values x1old=x1 x2old=x2 kcon=10 ! Maximum number of iterations. If more needed, something is wrong. do k=1,kcon j=j+1 x1=p*(1.0_np-exp(-x1old)) ! Iteration equation x2=x2old-( (x2old-p*(1.0_np-exp(-x2old))) & & /(1.0_np-p*exp(-x2old)) ) ! Newton-Raphson equation ! x2=p*(1.0_np-2.0_np*exp(-x2old))/(1.0_np-p*exp(-x2old)) ! NR equation compact, but numerically rotten write(*,990) j,p,x1,x2 if(abs(x2-x2old) .le. 1.0e-8_np*x2 & & .and. jcon2a .eq. 0) jcon2a=j ! Single precision test if(abs(x2-x2old) .ge. dcon2 .or. jcon2b .gt. 0) then ! Machine-precsion/jitter tests if(jcon2b .eq. 0) jcon2b=j if(k .ge. 4) exit ! Always do at least 4 iterations end if dcon1=abs(x1-x1old) dcon2=abs(x2-x2old) ! For machine-precision/jitter test x1old=x1 x2old=x2 end do ! write(*,990) j,p,x1,x2 print*,'NR iteration, relerr in fit', & ' relerr in it-function=',kcon print*,jcon2a,jcon2b,(x0-x2)/x2,abs(x1-x2)/x2 if(isatz .ge. 2) then print*,'xlow-x2,xhigh-x2,x0-x2,x2' print*,xlow-x2,xhigh-x2,x0-x2,x2 end if if(itest .eq. 1 .and. & & ( & abs(pd-1.0e-6_np) .lt. 1.0e-10_np .or. & precision 32 looses marbles & abs(pd-1.0e-5_np) .lt. 1.0e-10_np .or. & precision 32 looses marbles & abs(pd-0.0001_np) .lt. 1.0e-10_np .or. & precision 32 looses marbles ! & abs(pd-0.003_np) .lt. 0.0001_np .or. & precision 32 looses marbles & abs(pd-0.01_np) .lt. 0.001_np .or. & ! & abs(pd-0.03_np) .lt. 0.001_np .or. & & abs(pd-0.1_np) .lt. 0.001_np .or. & & abs(pd-0.13_np) .lt. 0.001_np .or. & & abs(pd-0.15_np) .lt. 0.001_np .or. & & abs(pd-0.2_np) .lt. 0.001_np .or. & & abs(pd-0.3_np) .lt. 0.001_np .or. & & abs(pd-0.5_np) .lt. 0.001_np .or. & & abs(pd-0.8_np) .lt. 0.001_np .or. & & abs(pd-1.0_np) .lt. 0.001_np & & ) ) then if(abs(pd-1.0e-6_np) .lt. 1.0e-10_np) then x2_in=1.99999933333377777745183802479528595E-0006_np print*,'Replace yes/no?' print*,x2,x2_in,x2-x2_in ! x2=x2_in end if print*,'For Testing for the low series coefficients.' xcoef2=(x2-xlow)/pd**2 + a2 xcoef3=(x2-xlow2)/pd**3 xcoef4=(x2-xlow3)/pd**4 xcoef5=(x2-xlow4)/pd**5 xcoef6=(x2-xlow5)/pd**6 xcoef7=(x2-xlow)/pd**7 + a7 xcoef8=(x2-xlow)/pd**8 + a8 xcoef9=(x2-xlow)/pd**9 + a9 print*,'xcoef2,a2,(a2-xcoef2)/xcoef2' print*,xcoef2,a2,(a2-xcoef2)/xcoef2 print*,'xcoef3,a3,(a3-xcoef3)/xcoef3' print*,xcoef3,a3,(a3-xcoef3)/xcoef3 print*,'xcoef4,a4,(a4-xcoef4)/xcoef4' ! print*,xcoef4,a4,(a4-xcoef4)/xcoef4 print* print*,'xcoef4*9.0_np,xcoef4*18.0_np' print*,xcoef4*9.0_np,xcoef4*18.0_np print*,'xcoef4*27.0_np,xcoef4*81.0_np' print*,xcoef4*27.0_np,xcoef4*81.0_np print*,'xcoef4*11.0_np,xcoef4*14.0_np' print*,xcoef4*11.0_np,xcoef4*14.0_np print*,'xcoef4*21.0_np,xcoef4*24.0_np' print*,xcoef4*21.0_np,xcoef4*24.0_np xc=abs(xcoef4) xsign=sign(1.0_np,xcoef4) in=0 if(in .eq. 1) then do k=1,2000 xd=real(k) xn=xd*xc xdiff=xn-real(int(xn+0.5_np)) relerr=xsign*(xdiff/xd)/pd itell=0 iprime=nprime do kk=2,iprime ! Must have one prime factor in this range jtell=mod(k,iiprime(kk)) if(jtell .eq. 0) itell=itell+1 end do if(abs(relerr) .lt. 20.0_np & & .and. itell .gt. 0) & print*,itell,k,xn,xdiff,relerr end do stop end if ! print* print*,'xcoef5,a5,(a5-xcoef5)/xcoef5' print*,xcoef5,a5,(a5-xcoef5)/xcoef5 xc=abs(xcoef5) xsign=sign(1.0_np,xcoef5) in=0 if(in .eq. 1) then do k=1,10 000 xd=real(k) xn=xd*xc xdiff=xn-real(int(xn+0.5_np)) relerr=xsign*(xdiff/xd)/pd if(abs(relerr) .lt. 20.0_np) print*,k,xn,xdiff,relerr end do stop end if ! print*,'xcoef6,a6,(a6-xcoef6)/xcoef6' print*,xcoef6,a6,(a6-xcoef6)/xcoef6 ! print* print*,'Remainder ',(x2-xlow)/pd**6 print* print*,'xcoef7,a7,(a7-xcoef7)/xcoef7' print*,xcoef7,a7,(a7-xcoef7)/xcoef7 print*,'xcoef7*9.0_np,xcoef7*18.0_np' print*,xcoef7*9.0_np,xcoef7*18.0_np print*,'xcoef7*27.0_np,xcoef7*81.0_np' print*,xcoef7*27.0_np,xcoef7*81.0_np ! print* print*,'xcoef8,a8,(a8-xcoef8)/xcoef8' print*,xcoef8,a8,(a8-xcoef8)/xcoef8 print*,'xcoef8*9.0_np,xcoef8*18.0_np' print*,xcoef8*9.0_np,xcoef8*18.0_np print*,'xcoef8*27.0_np,xcoef8*81.0_np' print*,xcoef8*27.0_np,xcoef8*81.0_np ! print* print*,'xcoef9,a9,(a9-xcoef9)/xcoef9' print*,xcoef9,a9,(a9-xcoef9)/xcoef9 print*,'xcoef9*9.0_np,xcoef9*18.0_np' print*,xcoef9*9.0_np,xcoef9*18.0_np print*,'xcoef9*27.0_np,xcoef9*81.0_np' print*,xcoef9*27.0_np,xcoef9*81.0_np end if end do 990 format(i10,f10.5,2f20.16) ! ! http://physics.nist.gov/cuu/Constants/Table/allascii.txt boltev=8.6173303e-5_np ! 8 digits boltzmann's constant NIST ryd=13.6056930090_np ! 11 digits rydberg NIST ! ryd=13.605692530_np ! http://en.wikipedia.org/wiki/Rydberg_constant t=ryd/(boltev*x) print*,'Comparison of Liddle reionization temperature and mine.' print*,5700.0_np,t ! end