From: "royaltm (Rafał Michalski)" <royaltm75@...>
Date: 2012-08-11T23:46:20+09:00
Subject: [ruby-core:47130] [ruby-trunk - Bug #6857][Open] bigdecimal/math BigMath.E/BigMath.exp R. P. Feynman inspired optimization


Issue #6857 has been reported by royaltm (Rafa�� Michalski).

----------------------------------------
Bug #6857: bigdecimal/math BigMath.E/BigMath.exp R. P. Feynman inspired optimization
https://bugs.ruby-lang.org/issues/6857

Author: royaltm (Rafa�� Michalski)
Status: Open
Priority: Normal
Assignee: 
Category: 
Target version: 
ruby -v: ruby 1.9.3p194 (2012-04-20) [x86_64-linux]


The algorythms to calculate E and exp programmed in BigMath module are the very straightforward interpretation of the series 1 + x + x^2/2! + 
x^3/3! + ....
Therefore they are slow.

Try it yourself:

      require 'bigdecimal/math'

      def timer; s=Time.now; yield; puts Time.now-s; end

      timer { BigMath.E(1000) }   #->  0.038848
      timer { BigMath.E(10000) }  #-> 16.526972
      timer { BigMath.E(100000) } #-> lost patience

That's because every iteration divides 1 by n! and the dividend grows extremely fast.

In "Surely You're Joking, Mr. Feynman!" (great book, you should read it if you didn't already) R. P. Feynman said: 

"One day at Princeton I was sitting in the lounge and overheard some mathematicians talking about the series for e^x, which is 1 + x + x^2/2! + 
x^3/3! Each term you get by multiplying the preceding term by x and dividing by the next number. For example, to get the next term after x^4/4! you 
multiply that term by x and divide by 5. It's very simple."

Yes it's very simple indeed. Why it's not been applied in such a great, modern and popular language? Is it because people just forget about simple solutions today?

Here is a Feynman's optimized version of BigMath.E:

      def E(prec)
        raise ArgumentError, "Zero or negative precision for E" if prec <= 0
        n = prec + BigDecimal.double_fig
        y = d = i = one = BigDecimal('1')
        while d.nonzero? && (m = n - (y.exponent - d.exponent).abs) > 0
          m = BigDecimal.double_fig if m < BigDecimal.double_fig
          d = d.div(i, m)
          i += one
          y += d
        end
        y
      end

Now, let's put it to the test:

      (1..1000).all? {|n| BigMath.E(n).round(n) == E(n).round(n) }
      => true
      BigMath.E(10000).round(10000) == E(10000).round(10000)
      => true

What about the speed then?

      timer { E(1_000) }     #-> 0.003832 ~ 10 times faster
      timer { E(10_000) }    #-> 0.139862 ~ 100 times faster
      timer { E(100_000) }   #-> 8.787411 ~ dunno?
      timer { E(1_000_000) } #-> ~11 minutes

The same simple rule might be applied to BigDecimal.exp() which originally uses the same straightforward interpretation of power series.
Feynman's pure ruby version of BigMath.exp (the ext version seems now pointless anyway):

      def exp(x, prec)
        raise ArgumentError, "Zero or negative precision for exp" if prec <= 0
        x = case x
        when Float
          BigDecimal(x, prec && prec <= Float::DIG ? prec : Float::DIG + 1)
        else
          BigDecimal(x, prec)
        end
        one = BigDecimal('1', prec)
        case x.sign
        when BigDecimal::SIGN_NaN
          return BigDecimal::NaN
        when BigDecimal::SIGN_POSITIVE_ZERO, BigDecimal::SIGN_NEGATIVE_ZERO
          return one
        when BigDecimal::SIGN_NEGATIVE_FINITE
          x = -x
          inv = true
        when BigDecimal::SIGN_POSITIVE_INFINITE
          return BigDecimal::INFINITY
        when BigDecimal::SIGN_NEGATIVE_INFINITE
          return BigDecimal.new('0')
        end
        n = prec + BigDecimal.double_fig
        if x.round(prec) == one
          y = E(prec)
        else
          y = d = i = one
          while d.nonzero? && (m = n - (y.exponent - d.exponent).abs) > 0
            m = BigDecimal.double_fig if m < BigDecimal.double_fig
            d = d.mult(x, m).div(i, m)
            i += one
            y += d
          end
        end
        y = one.div(y, n) if inv
        y.round(prec - y.exponent)
      end

      (1..1000).all? {|n| exp(E(n),n) == BigMath.exp(BigMath.E(n),n) }
      # => true
      (1..1000).all? {|n| exp(-E(n),n) == BigMath.exp(-BigMath.E(n),n) }
      # => true
      (-10000..10000).all? {|n| exp(BigDecimal(n)/1000,100) == BigMath.exp(BigDecimal(n)/1000,100) }
      # => true
      (1..1000).all? {|n| exp(BigMath.PI(n),n) == BigMath.exp(BigMath.PI(n),n) }
      # => true

      timer { BigMath.exp(BigDecimal('1').div(3, 10), 100) }    #-> 0.000496
      timer { exp(BigDecimal('1').div(3, 10), 100) }            #-> 0.000406 faster but not that really

      timer { BigMath.exp(BigDecimal('1').div(3, 10), 1_000) }  #-> 0.029231
      timer { exp(BigDecimal('1').div(3, 10), 1_000) }          #-> 0.004554 here we go... 

      timer { BigMath.exp(BigDecimal('1').div(3, 10), 10_000) } #-> 12.554197
      timer { exp(BigDecimal('1').div(3, 10), 10_000) }         #->  0.189462 oops :)

      timer { exp(BigDecimal('1').div(3, 10), 100_000) }        #-> 11.914613 who has the patience to compare?

Arguments with large mantissa should slow down the results of course:

      timer { BigMath.exp(BigDecimal('1').div(3, 1_000), 1_000) }   #->  0.119048
      timer { exp(BigDecimal('1').div(3, 1_000), 1_000) }           #->  0.066177

      timer { BigMath.exp(BigDecimal('1').div(3, 10_000), 10_000) } #-> 68.083222
      timer { exp(BigDecimal('1').div(3, 10_000), 10_000) }         #-> 29.439336

Though still two times faster than the ext version.

It seems Dick Feynman was not such a joker after all. I think he was a master in treating lightly "serious" things and treating very seriously things that didn't matter to anybody else.

I'd write a patch for ext version if you are with me. Just let me know.



-- 
http://bugs.ruby-lang.org/