1.0k
views6
comments

All-

Here is a Plouffe-Bellard Pi calculator I wrote for an upcoming article. I copied the workings from a C# application and as such, it's slower than the C# app. Anyone have any suggestions on speeding this up? You can download the full test solution from [link:64.32.243.109]

Here's the source:

module pBopen Systemlet rec is_prime(n:float) =
	let mutable ret = true in
	if (n%2.0) = 0.0 then ret <- false
	else begin
		let r = Int32.of_float(sqrt(n)) in
		let mutable i = 3 in            
		while i <= r && ret do
			if (n % Float.of_int(i)) = 0.0 then ret <- false;
			i <- i + 2
		done;
	end;
	ret
        
let rec next_prime (a:float) =
	let mutable n = a in
	n <- n + 1.0;
	while not(is_prime(n))  do
		n <- n + 1.0;
	done;
	n
        
let rec inv_mod( x, y )  =
	let mutable q = 0 in
	let mutable u = 0 in
	let mutable v = 0 in
	let mutable a = 0 in 
	let mutable c = 1 in
	let mutable t = 0 in
	
	u <- Int32.of_float(x);       
	v <- Int32.of_float(y);
	
	q <- v/u;

	t <- c;
	c <- a-q*c;
	a <- t;        t <- u;
	u <- v-q*u;
	v <- t;                       
	while not(u=0) do
		q <- v/u;

		t <- c;
		c <- a-q*c;
		a <- t;                        t <- u;
		u <- v-q*u;
		v <- t;
	done;                
	a <- a%Int32.of_float(y);
	if (a<0) then a<-Int32.of_float(y)+a;
	Float.of_int32(a)       
    
let rec mul_mod( a, b, c ) = (a*b) % clet rec pow_mod( (a:float), (inb:int), (m:float)) =   
	let mutable r = 1.0 in
	let mutable aa = 0.0 in 
	let mutable run = true in
	let mutable b = 0 in        
	b <- inb;
	aa <- a;        
	while run do
		if not((b &&& 1) = 0) then r <- mul_mod(r, aa, m);                     b <- b >>> 1;            if (b = 0) then run <- false
		else aa <- mul_mod(aa, aa, m)
	done;        
	r    
                    
let rec CalculatePiDigits (n:int) = 
		let mutable av = 0.0 in
		let mutable vmax = 0.0 in
		let mutable sum = 0.0 in
						
		let N = floor((Float.of_int(n)+20.0) * (log(10.0) / log(2.0))) in
		let mutable a = 3.0 in
		let mutable i = 0.0 in
		let mutable s = 0.0 in
		let mutable t = 0.0 in
		let mutable num = 0.0 in
		let mutable den = 0.0 in        
		//System.Diagnostics.Trace.WriteLine("N=" + N.ToString());
		while a <= (2.0 * N) do
			//System.Diagnostics.Trace.WriteLine("a=" + a.ToString());
			vmax <- floor(log(2.0 * N) / log(a));
			//System.Diagnostics.Trace.WriteLine("vmax=" + vmax.ToString());
			av <- 1.0;
			i <- 0.;
			while i < vmax do
							av <- av * a;
							i <- i + 1.0
			done;
			//System.Diagnostics.Trace.WriteLine("av=" + av.ToString());
			s <- 0.0;
			num <- 1.0;
			den <- 1.0;
			let mutable v = 0.0 in
			let mutable k = 1.0 in
			let mutable kq = 1.0 in
			let mutable kq2 = 1.0 in
			while k <= N do
					t <- k;
					if (kq >= a) then begin
									t <- t / a;
									v <- v - 1.0;
									while ((t % a) = 0.0) do
													t <- t / a;
													v <- v - 1.0;
									done;                
					kq <- 0.0
					end;
					kq <- kq + 1.;
					//System.Diagnostics.Trace.WriteLine("num in=" + num.ToString());
					//System.Diagnostics.Trace.WriteLine("t=" + t.ToString());
					//System.Diagnostics.Trace.WriteLine("av=" + av.ToString());
					num <- mul_mod(num, t, av);
					//System.Diagnostics.Trace.WriteLine("num out=" + num.ToString());
					
					t <- 2. * k - 1.;
					
					if (kq2 >= a) then begin
						if (kq2 = a) then begin
								t <- t / a;
								v <- v + 1.0;
								while ((t % a) = 0.0) do
												t <- t / a;
												v <- v + 1.0;
								done;       
						end;
						kq2 <- kq2 - a;
					end;
					den <- mul_mod(den, t, av);
					kq2 <- kq2 + 2.;
					
					if (v > 0.) then begin
							//System.Diagnostics.Trace.WriteLine("den=" + den.ToString() + " av=" + av.ToString());
							t <- inv_mod(den, av);
							//System.Diagnostics.Trace.WriteLine("t=" + t.ToString());
							t <- mul_mod(t, num, av);
							t <- mul_mod(t, k, av);                

							i <- v;
							while i < vmax do
											t <- mul_mod(t, a, av);
											i <- i + 1.;
							done;                                s <- s + t;                                if (s >= av) then s <- (s - av);                                                
					end;    
					k <- k + 1.;            
		done;
		
		//System.Diagnostics.Trace.WriteLine("av=" + av.ToString());
		t <- pow_mod(10., n - 1, av);
		//System.Diagnostics.Trace.WriteLine("t=" + t.ToString());
		//System.Diagnostics.Trace.WriteLine("s in=" + s.ToString());
		s <- mul_mod(s, t, av);
		//System.Diagnostics.Trace.WriteLine("s out=" + s.ToString());
		sum <- (sum + s / av) % 1.0;
		//System.Diagnostics.Trace.WriteLine("sum=" + sum.ToString());
		
		// next 
		a <- next_prime(a)
	done;
	let Result = Int32.of_float(sum * 1e9) in
	string.Format("{0:D9}", Result)

An interesting sanity check for F# code generation. Certainly there is no reason why this will run slower in principle - there may be a glitch or two in IL code gen or inlining but that will be easy to fix - we would expect to generate essentially the same IL as C# in this case. Can you send the perf results you've observed? Also, have you inspected the IL produced by F# (with -O3) and with C#? And have you used any profiling tools on the code?

BTW let's move this to a public facing folder (optionScalper - can you do that for us? Would there be a way for me to do that is fsAdmin role?)

Thanks!

Don

By dsyme on 4/1/2006 10:04 AM (permalink)

Don-

Thanks for the feedback! I'm really looking for better ways of accomplishing the PB algorithm in F#. Are there native F# functions I should use instead of the pow_mod(), inv_mod(), and mul_mod() functions? I'm making frequent use of mutable due to the fact I'm copying C# code, how should some of this be done in F#? What I've looked around for (and been unable to find) is a reference that discusses algorithm implementation in F# vs. C#. Can you make any suggestions? I guess I'm looking for F# best practices (which may lead to improved performance) rather than just improving performance.

I should also mention that Fabrice Bellard's algorithm is discussed here.

Thanks for any additional feedback!

-Chad

By calbrecht on 4/2/2006 7:29 AM (permalink)

Another note. I converted a majority of the F# module to do integer math. This is more in line with the C# code. This increases the performance a bit (291 ms. vs. 484 ms.). The C# class can still run 100 digits in 62 ms. whereas the F# code take 291 ms.

By calbrecht on 4/2/2006 4:28 PM (permalink)

Hi Chad,

Here's my quick translation of the first two functions - you will see that the use of recursive functions makes for a very different style.

open System
// note is_prime is buggy on 2!
let is_prime n = 
  (n%2 > 0) && 
  let r = Float.to_int(Math.Sqrt(float(n))) in
  let rec check i = (i > r) or (n % i > 0) && check (i+2) in 
  check 3

let rec next_prime a =
  let n = a+1 in 
  if is_prime n then n else next_prime(n);;

The key things here are (a) recursive functions (b) and, or and if-then-else expressions to simplify control-flow and (c) local functions - check is a function that utilizes the value 'r' .

Re performance - we will be taking a good look at the generated IL tomorrow. We've actually been doing a bit of work on the optimizer recently in conjunction with optimiing the F# BigInteger implementation that will make a lot of difference to the code you first showed. I would also expect the code above to run about as fast as the C# equivalent (if not, it would be considered a bug in F#). Both languages have strongly typed input, so full optimizations down to IL are possible. Make sure you're using -O3, of course.

By dsyme on 4/5/2006 7:12 PM (permalink)

As another point of reference, here is the C# class I copied.

 

public class Plouffe_Bellard
    {
        public Plouffe_Bellard() {}

        private static int mul_mod(int a, int b, int m) 
        {
            return (int) (((long) a * (long) b) % m);
        }

        /* return the inverse of x mod y */
        private static int inv_mod(int x, int y) 
        {
            int q,u,v,a,c,t;

            u=x;
            v=y;
            c=1;
            a=0;

            do 
            {
                q=v/u;
    
                t=c;
                c=a-q*c;
                a=t;
    
                t=u;
                u=v-q*u;
                v=t;
            } while (u!=0);

            a=a%y;

            if (a<0) 
            {
                a=y+a;
            }

            return a;
        }

        /* return (a^b) mod m */
        private static int pow_mod(int a, int b, int m)
        {
            int r, aa;
   
            r=1;
            aa=a;

            while (true) 
            {
                if ((b & 1) != 0)
                {
                    r = mul_mod(r, aa, m);
                }

                b = b >> 1;

                if (b == 0)
                {
                    break;
                }

                aa = mul_mod(aa, aa, m);
            }

            return r;
        }

        /* return true if n is prime */
        private static bool is_prime(int n)
        {
            if ((n % 2) == 0) 
            {
                return false;
            }

            int r = (int) Math.Sqrt(n);

            for (int i = 3; i <= r; i += 2)
            {
                if ((n % i) == 0) 
                {
                    return false;
                }
            }

            return true;
        }

        /* return the prime number immediatly after n */
        private static int next_prime(int n)
        {
            do 
            {
                n++;
            } while (!is_prime(n));

            return n;
        }

        public String CalculatePiDigits(int n)
        {
            int av, vmax, num, den, s, t;
  
            int N = (int) ((n + 20) * Math.Log(10) / Math.Log(2));
            //System.Diagnostics.Trace.WriteLine("N=" + N.ToString());
            double sum = 0;

            for (int a = 3; a <= (2 * N); a = next_prime(a)) 
            {
                //System.Diagnostics.Trace.WriteLine("a=" + a.ToString());
                vmax = (int) (Math.Log(2 * N) / Math.Log(a));
                //System.Diagnostics.Trace.WriteLine("vmax=" + vmax.ToString());

                av = 1;

                for (int i = 0; i < vmax; i++) 
                {
                    av = av * a;
                }
                //System.Diagnostics.Trace.WriteLine("av=" + av.ToString());

                s = 0;
                num = 1;
                den = 1;
                int v = 0;
                int kq = 1;
                int kq2 = 1;

                for (int k = 1; k <= N; k++) 
                {

                    t = k;

                    if (kq >= a) 
                    {
                        do 
                        {
                            t = t / a;
                            v--;
                        } while ((t % a) == 0);

                        kq = 0;
                    }
                    kq++;
                    //System.Diagnostics.Trace.WriteLine("num in=" + num.ToString());
                    //System.Diagnostics.Trace.WriteLine("t=" + t.ToString());
                    //System.Diagnostics.Trace.WriteLine("av=" + av.ToString());
                    num = mul_mod(num, t, av);
                    //System.Diagnostics.Trace.WriteLine("num out=" + num.ToString());

                    t = 2 * k - 1;

                    if (kq2 >= a) 
                    {
                        if (kq2 == a) 
                        {
                            do 
                            {
                                t = t / a;
                                v++;
                            } while ((t % a) == 0);
                        }
                        kq2 -= a;
                    }
                    den = mul_mod(den, t, av);
                    kq2 += 2;
      
                    if (v > 0) 
                    {
                        //System.Diagnostics.Trace.WriteLine("den=" + den.ToString() + " av=" + av.ToString());
                        t = inv_mod(den, av);
                        //System.Diagnostics.Trace.WriteLine("t=" + t.ToString());
                        t = mul_mod(t, num, av);
                        t = mul_mod(t, k, av);

                        for (int i = v; i < vmax; i++)
                        {
                            t = mul_mod(t, a, av);
                        }

                        s += t;

                        if (s >= av)
                        {
                            s -= av;
                        }
                    }
      
                }

                //System.Diagnostics.Trace.WriteLine("av=" + av.ToString());
                t = pow_mod(10, n - 1, av);
                //System.Diagnostics.Trace.WriteLine("t=" + t.ToString());
                //System.Diagnostics.Trace.WriteLine("s in=" + s.ToString());
                s = mul_mod(s, t, av);
                //System.Diagnostics.Trace.WriteLine("s out=" + s.ToString());
                sum = (sum + (double) s / (double) av) % 1.0;
                //System.Diagnostics.Trace.WriteLine("sum=" + sum.ToString());
            }
   
            int Result = (int) (sum * 1e9);

            String StringResult = String.Format("{0:D9}", Result);

            return StringResult;
        }

        public int DigitsReturned()
        {
            return 9;
        }
    }

By calbrecht on 4/2/2006 4:22 PM (permalink)

Don,

I advised Chad to put this in THCF as he's publishing an article on this topic. The main idea here was to decide what was useful in the code and then to "pretty it up" into an article. At the time it didn't make sense to post these details as they should lead into a polished article (and I'm doing the same on my random article).

You have moderator rights. So on any thread that you want to moderate, just hit the "Moderate" button and select "Move". You will be presented with a list of forums that can be the target of the post. Of course, I can move this as well for you. Let me know what you prefer, i.e. leave it here, try moving it yourself, or I can move it. :-)

---O

By optionsScalper on 4/1/2006 7:08 PM (permalink)

Topic tags

Built with WebSharper

Home

Answers

Events

Courses

Groups and Conferences

Blogs

Jobs

Developers

Topic tags