• ulrich_y's avatar
    New MPL routine · 9ff6a342
    ulrich_y authored
    This is the algorithm by GiNaC. In theory one
    could extend this to add a caching mechanism such as
    
        complex(kind=prec) :: cache(size(x),MPLMaxQ)
    
        do q=1,j
          cache(:,q) = x**q/q**m
        enddo
        do q=1,MPLMaxQ
          res = t(1)
    
          ! Fortran uses Column-major order, hence cache(:,q) is
          ! faster than cache(q,:).
          cache(:,q+j-1) = x**(q+j-1)/(q+j-1)**m
          t(j) = t(j) + cache(j,q)
          do k=1,j-1
            t(j-k) = t(j-k) + t(j-k+1) * cache(j-k,k+q)
          enddo
    
          if (mod(q,2) .eq. 1) then
            if (abs(t(1)-res).lt.MPLdel) exit
          endif
        enddo
    
    In practice this doesn't really help because any
    time saved with the cache is paid back through the
    allocation and clearing of cache(:,:). Both
    variations work similarly well now. If at some
    point we might need MPLs with many more arguments
    (size(x)), this might change.
    9ff6a342
globals.f90 1.31 KB