Sieve of Eratosthene

Eratosthenes เป็นนักคณิตศาสตร์ชาวกรีก ที่คิดค้นสิ่งต่าง ๆ ไว้มากพอสมควร แต่ วันนี้เราจะมาพูดกันถึง การหาค่า prime หรือจำนวนเฉพาะ ที่ใช้วิธีตัดออกทีละชุด (screen หรือ sieve นั่นเอง)
วิธีการของเค้าก็ง่าย ๆ ครับ คือ ในตัวเลขจำกัดจำนวนนึงเราสามารถหาค่า จำนวนเฉพาะได้โดยใช้ ตัวเลขเริ่มต้น (2) เพื่อที่จะนำไป หาค่าที่เป็นตัวคูณเพื่อ screen ตัวที่เป็น ค่าที่เกิดจากการคูณออกจากความน่าจะเป็นที่จะเป็น จำนวนเฉพาะ โดยตัวที่จะนำมาสกรีน ก็นำมาจถึงแค่ ตัวเลขที่ น้อยกว่าหรือเท่ากับ sqrt ของ ค่าสูงสุดที่ต้องการหา (ตรงนี้พิสูจน์ไม่ยาก)
เช่น หาจำนวนเฉพาะ จนถึง 20
2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
ใช้ 2 screen (2^2 = 4 < 20)
2, 3, 5, 7, 9, 11, 13, 15, 17, 19
ใช้ 3 screen (3^2 = 9 < 20)
2, 3, 5, 7, 11, 13, 17, 19
5 ไม่ต้องนำมาสกรีนอีก เพราะ 5^2 > 20
เพราะฉะนั้นค่าที่ไม่ถูกตัดทิ้งทั้งหมดก็คือ เลขจำนวนเฉพาะ
เขียนด้วย C++ ได้ดังนี้

vector<bool> GetPrimeBitSets(int from, int to) 
{
    vector<bool> bs (to + 1); //Create a bool vector to hold a primality
    for(int i = 0; i < to + 1; i++) //Initialize the primality vector
    {
        bs[i] = bool(i & 1);        //even number will be set as false
    }
    int j = 0;
    from = (from <= 3)? 3 : from; 
    bs[0] = false;
    bs[1] = false;
    bs[2] = true;                   //2 is prime
    bs[3] = true;                   //3 is prime
    for(int i = 3; i * i <= to; i += 2) //All odds numbers are considering
    {
        if(!bs[i]) continue;        //if the number is not a prime then no need to do
        for(j = i + i; j <= to; j += i) // populate all the multiply of i
        {       
            bs[j] = false;
        }       
    }
    return bs;                      //when done, the left number are prime
}

Binary GCD algorithm

จากที่เห็นกันมาแล้วว่า Euclidean algorithm ง่ายต่อการเข้าใจแล้วก็ง่ายต่อการเขียน อีกทั้งยังเร็วกว่าวิธีที่เราคิดกันขึ้นมาเอง แต่ว่า ความเร็วของ Euclidean algorithm ก็ยังไม่สามารถเทียบเท่ากับ binary GCD algorithm ข้อดีอีกอย่างของ algorithm นี้คือสามารถนำมาเขียนบน computer programming ได้ง่ายเช่นเดียวกับ Euclidean algorithm
Binary GCD algorithm มีวิธีการดังนี้

1. ถ้า a = 0 จะทำให้ gcd(a, b) = b เพราะอะไรก็ตามนำมาหาร 0 ก็ได้ 0 แน่นอนอยู่แล้ว และ b ก็เป็นตัวนึง
2. ถ้าทั้ง a และ b เป็นเลขคู่ gcd(a, b) = 2*gcd(a/2, b/2)
3. ถ้า a เป็นเลขคู่ และ b เป็นเลขคี่ gcd(a, b) = gcd(a/2, b) เพราะว่ายังไงซะ 2 ก็ไม่ใช่ ตัวประกอบนึงแน่ ๆ
4. ถ้า a เป็นเลขคี่ และ b เป็นเลขคู่ gcd(a, b) = gcd(a, b/2) เหตุผลเหมือนข้อ 3
5. ถ้าทั้ง a และ b เป็นเลขคี่ gcd(a, b) = gcd(abs((a – b)/2), b) = gcd(a, abs((a – b)/2)) เหตุผลมาจาก Euclidean algorithm รวมกะ ข้อ 3 และ 4

เขียนเป็น C/C++ ได้ดังนี้

unsigned int gcd(unsigned int a, unsigned int b)
{
  unsigned int c = 0;
  while((a & 1) == 0 && (b & 1) == 0) // both are even number
  {
      a >>= 1;    // divided by 2
      b >>= 1;    // divided by 2
      c++;        // power of 2 to multiply back
  }

  do
  {
      if((a & 1) == 0)
          a >>= 1;
      else if((b & 1) == 0)
          b >>= 1;
      else if(a >= b)
          a = (a - b) >> 1;
      else  
          b = (b - a) >> 1;
  } while (a > 0);
  return b << c;
}

เขียนแบบ Recursive ได้เป็น

unsigned int gcd(unsigned int a, unsigned int b)
{
  if(a == 0)
      return b;
  else if((a & 1) == 0 && (b & 1) == 0)
      return 2 * gcd(a >> 1, b >> 1);
  else if((a & 1) == 0)
      return gcd(a >> 1, b);
  else if((b & 1) == 0)
      return gcd(a, b >> 1);
  else if(a >= b)
      return gcd((a - b) >> 1, b);
  else
      return gcd(a, (b - a) >> 1);
}

ที่มา และข้อมูลเพิ่มเติม

Binary GCD algorithm at answers.com

Eucledian algorithm at my blog or at my blogger.com

Greatest Common Divisor

วันนี้กลับไปอ่านหนังสือ discrete math สมัยเรียนมหาลัย ไปเจอบทนึงพูดถึงเรื่อง GCD หรือ Great Common Divisor หรือในไทยก็ ห.ร.ม. หรือ หารร่วมมาก

วิธีการหาหารร่วมมากโดยปรกติ ก็ คือ แยก factor ของตัวเลข 2 ตัวที่ต้องการหา GCD แล้วก็นำมาหารกันจนเหลือตัวที่ไม่สามารถหารกันลงแล้ว

เช่น

GCD(14, 56) = 2x7 and 2x2x2x7 = 1 and 2x7 = 14

วิธีเขียน Program ก็อาจจะทำได้โดย

int GCD(int a, int b)
{
    int temp = a;
    int temp2 = b;

    int gcd = 1;
    while((temp % 2 == 0) && (temp2 % 2 == 0))
    {
        temp /= 2;
        temp2 /= 2;
        gcd *= 2;
    }
    for(int i = 3; i <= a && i <= b; i+=2)
    {
        while((temp  % i == 0) && (temp2 % i == 0))
        {       
            temp /= i;
            temp2 /= i;
            gcd *= i;
        }       
    }
    return gcd;
}

แต่มีอีกหลายวิธีที่ทำได้รวดเร็วกว่าวิธีที่เสนอมานี้เยอะมาก วิธีนึงก็คือ วิธี Euclidean วิธีให้ความคิดไว้ว่า GCD ของ a กับ b เท่ากับ GCD ของ b กับ abs(a-b) หรือสามารถเขียนเป็น program ได้ดังนี้

int GCD(int a, int b)
{
   if(a % b == 0)
   {
       return b;
   }
   else   
   {
       return GCD(b, abs(a - b));
   }
}

จะเห็นได้ว่า วิธี Euclidean ทำให้ program ง่ายขึ้นและรวดเร็วขึ้นด้วย