Code covered by the BSD License  

Highlights from
HPF - a big decimal class

HPF - a big decimal class

by

 

04 May 2012 (Updated )

High precision floating point arithmetic, a new class written in MATLAB

updateHPFconstants(maxdigits,constantname)
function updateHPFconstants(maxdigits,constantname)
% updateHPFconstants: generates new values for the HPF stored constants: e, pi, and the natural log of 10
% usage: updateHPFconstants(maxdigits,constantname)
%
% When I wrote the HPF tools, I generated those important constants to
% what seemed like a sufficient number of digits for any reasonable
% computation. However, Moore's law being what it is, perhaps 100,000
% digits or so will be too small next year. This tool allows matlab to
% update those constants programmatically, then saving them out for future
% use.
%
% Recognize that updateHPFconstants will be a cpu hog, taking a serious
% amount of time for large numbers of digits. This is NOT a tool that you
% will use on any daily basis, since those constants are already stored
% internally to plenty of digits for any sane, rational person. In fact,
% expect the running time for this code to be in the hours, or even days
% to update these constants beyond a million digits or so.
%
% arguments: (input)
%  maxdigits - scalar integer - defines a new maximum number of stored
%      digits for the three constants usd internally, e, pi, and the
%      natural log of 10.
%
%      if maxdigits is less than the currently stored number of digits
%      for those constants, then nothing is done.
%
%  constantname - string containing the name of the constant to be updated
%      in value. If the new value has more digits in it for this constant
%      than are currently stored by HPF, then an update is done to the
%      file on disk for that constant.
%
%      constant name may be a string containing one constant name, or a
%      cell array of strings.
%
%      Currently, the list of constants in HPF is {'pi', 'e', 'ln10'}.

% get the path to _special_numbers_.mat
snpath = which('_special_numbers_.mat');

% bring in the current values in _special_numbers_.mat
load(snpath)

% 50 spare digits should be entirely adequate
sparedigits = 50;
NDig = [maxdigits,sparedigits];

if (nargin < 2) || isempty(constantname)
  constantname = {'pi' 'e', 'ln10'};
else
  constantname = cellstr(constantname);
  if any(~ismember(constantname,{'pi', 'e', 'ln10'}))
    error('Legal constants that may be updated are: ''pi'', ''e'', ''ln10''')
  end
end

% do we need to update e?
if ismember('e',constantname) && (maxdigits > numel(edigits))
  % we need to update e. do so by computing exp(2^-20), then raise that to
  % the power 1048576. I'll keep plenty of spare digits on the side, so there
  % is no fear of a loss of precision here.
  z = hpf('0.00000095367431640625',NDig); % 1/1048576 = 2^-20
  
  % How many terms in the Taylor series do we need to compute?
  % Since the exponential series has as its last term
  % z^m/factorial(m), we can quit when the natural log of that
  % expression is less than the log of our exponential goal.
  % Since log(exp(z)) is never less than -0.5, we have a simple
  % cutoff, based on the number of significant digits in our
  % hpf number.
  %
  %  -0.5 + log(10^-NDig) = log(z^m/m!) = m log(z) - log(m!)
  %
  % the logs here are natural logs of course. Using Stirling's
  % approximation for m!, we get
  %
  %  -0.5 - NDig*log(10) = m*log(z) - [1/2*log(2*pi) + 1/2*log(m) + m*log(m) - m]
  %
  % Combining terms, this reduces to
  %
  %  1/2*(log(2*pi)-1) - NDig*log(10) = m*(log(z) + 1) - log(m)*(m + 1/2)
  %
  % Bump up NDig by 1, just to ensure convergence for the full digits.
  % fun is a decreasing function. The zero crossing defines the number
  % of the last term we need to compute.
  fun = @(m) m.*(log(1/1048576) + 1) - log(m).*(m + 1/2) + ...
    (sum(NDig)+1)*log(10) - 1/2*(log(2*pi) - 1);
  
  mterms = termsbisector(fun);
  
  % run the Taylor series in reverse now. This is why I wanted to know how
  % many terms would be necessary. By running the loop backwards, I
  % avoid divisions, and a divide is more expensive than a multiply for
  % an hpf number.
  E1 = hpf('1',NDig);
  Fact = E1;
  H = waitbar(0,['Computing exp(1) to ',num2str(maxdigits),' digits']);
  for m = mterms:-1:1
    % update a waitbar
    if mod(m,100) == 0
      waitbar((mterms - m)/mterms,H)
    end
    
    Fact = Fact*m;
    E1 = z*E1 + Fact;
  end
  delete(H)
  
  % because we ran the loop backwards, in the end we need to
  % divide by factorial(mterms)
  E1 = E1./Fact;
  
  % raise E1 to the 1048576'th power, to get exp(1). Squaring E1 20 times
  % will suffice. Along the way, compute a good approximation to
  % exp(log(10)) in binary fractional powers of e
  E1 = E1.*E1;  % 2^-19'th power
  E1 = E1.*E1;  approx10 = E1; % 18
  E1 = E1.*E1;  % 17
  E1 = E1.*E1;  % 16
  E1 = E1.*E1;  approx10 = approx10.*E1; % 15
  E1 = E1.*E1;  approx10 = approx10.*E1; % 14
  E1 = E1.*E1;  % 13
  E1 = E1.*E1;  approx10 = approx10.*E1; % 12
  E1 = E1.*E1;  approx10 = approx10.*E1; % 11
  E1 = E1.*E1;  approx10 = approx10.*E1; % 10
  E1 = E1.*E1;  % 9
  E1 = E1.*E1;  approx10 = approx10.*E1; % 8
  E1 = E1.*E1;  % 7
  E1 = E1.*E1;  approx10 = approx10.*E1; % 6
  E1 = E1.*E1;  approx10 = approx10.*E1; % 5
  E1 = E1.*E1;  % 4
  E1 = E1.*E1;  % 3
  E1 = E1.*E1;  approx10 = approx10.*E1; % 2
  E1 = E1.*E1;  % 1
  E1 = E1.*E1;  approx10 = approx10.*E1.*E1;% 0
  
  % extract the digits as a uint8 vector to be saved back into
  % the special numbers mat file.
  E1dig = mantissa(E1);
  edigits = uint8(E1dig(1:maxdigits));
  
end

% update pi
if ismember('pi',constantname) && (maxdigits > numel(pidigits))
  NDig = [maxdigits,50];
  H = waitbar(0,['Computing pi to ',num2str(maxdigits),' digits (quadratic convergence)']);
  
  a0 = hpf('1',NDig);
  b0 = sqrt(hpf('0.5',NDig));
  t0 = hpf('0.25',NDig);
  p0 = hpf('1',NDig);
  
  notdone = true;
  while notdone
    a1 = (a0 + b0)/2;
    b1 = sqrt(a0.*b0);
    t1 = t0 - p0.*(a0 - a1).^2;
    p1 = 2 .*p0;
    
    % has our estimate converged?
    k = find(t0.Migits ~= t1.Migits,1,'first');
    if isempty(k)
      % t has converged to all of its digits, so we can generate
      % the final estimate for pi, and terminate the loop
      piest = (a1 + b1).^2./(4 .*t1);
      
      notdone = false;
    else
      % update the waitbar...
      waitbar(k/maxdigits,H)
      
      a0 = a1;
      b0 = b1;
      t0 = t1;
      p0 = p1;
    end
  end
  % zap the waitbar
  delete(H)
  piestdig = mantissa(piest);
  pidigits = uint8(piestdig(1:maxdigits));
  
end

% update log(10)
if ismember('ln10',constantname) && (maxdigits > numel(ln10digits))
  % use the current value for log(10) as an offset to make the
  % natural log series converge like a bat out of a very warm place.
  ln10 = hpf('ln10',[numel(ln10digits),0]);
  NDig = [maxdigits,50];
  ln10 = augmentdigits(ln10,NDig);
  
  % first compute exp(ln10) for the old approximation, to the
  % desired number of digits, so it will not be exactly 10.
  % I can't use exp to do so directly, since it uses the old
  % value of ln10 inside and exp uses ln10.
  z = ln10 - 2 - 1/4 - 1/32 - 1/64 - 1/256 - 1/1024 - 1/2048 - 1/4096 - 1/16384 - 1/32768 - 1/262144;
  
  fun = @(m) m.*(log(double(z)) + 1) - log(m).*(m + 1/2) + ...
    (sum(NDig)+1)*log(10) - 1/2*(log(2*pi) - 1);
  
  mterms = termsbisector(fun);
  
  % again, run the exponential Taylor series in reverse.
  Ez = hpf('1',NDig);
  Fact = Ez;
  H = waitbar(0,['Computing log(10) to ',num2str(maxdigits),' digits']);
  for m = mterms:-1:1
    % update a waitbar
    if mod(m,100) == 0
      waitbar((mterms - m)/mterms,H)
    end
    
    Fact = Fact*m;
    Ez = z*Ez + Fact;
  end
  delete(H)
  
  % because we ran the loop backwards, in the end we need to
  % divide by factorial(mterms)
  Ez = Ez./Fact;
  
  % multiply by our approximation to exp(log(10)) from before
  Ez = Ez.*approx10;
  
  % Ez is now exp(ln10). We can use this to bootstrap our computation
  % log(10) to the desired number of digits. Thus divide 10 by Ez. That
  % will be a number fairly close to 1, but slightly above or below 1.
  z = reciprocal(Ez).*10;
  
  % transform z. this transformation creates y VERY near zero, since z
  % was incredibly close to 1.
  y = (z - 1)./(z + 1);
  
  % we can find the log(z) with a series in y, that will take vanishingly
  % few terms, since we already had a very good approximation for log(x).
  logz = 2.*y;
  ysq = y.*y;
  yterm = logz;
  notdone = true;
  n = 1;
  while notdone
    n = n + 2;
    yterm = yterm.*ysq;
    logz = logz + yterm./n;
    if yterm.Exponent < -sum(NDig)
      notdone = false;
    end
  end
  
  % we have found an approximation for ln10 that is accurate to the
  % requested number of digits
  ln10 = ln10 + logz;
  
  % extract those digits as uint8
  ln10dig = mantissa(ln10);
  ln10digits = uint8(ln10dig(1:maxdigits));
  
end

% save out the new file to be used for future computations
save(snpath,'edigits','pidigits','ln10digits')

% we need to cause the persistent variables in hpf to be reinitialized
clear hpf

% =============================================================
%      end mainline hpf, begin subfunctions
% =============================================================

function mterms = termsbisector(fun)
% mterms is the last term in the Taylor series that we will
% need to get a good approximation from that series. Always take
% at least 5 terms in the series.

if fun(4) <= 0
  % always take at least 5 terms in the series
  mterms = 5;
else
  % bracket the root
  ab = [4 8];
  fab = fun(ab);
  while fab(2) > 0
    ab = [ab(2),2*ab(2)];
    fab = [fab(2),fun(ab(2))];
  end
  % simple bisection step
  while diff(ab) > 1
    c = mean(ab);
    fc = fun(c);
    if fc <= 0
      ab = [ab(1),c];
      fab = [fab(1),fc];
    else
      ab = [c,ab(2)];
      fab = [fc,fab(2)];
    end
  end
  
  % pick the upper end point.
  mterms = ab(2);
end % function termsbisector

Contact us