[C Programming] 擴充Fibonacci數-延伸

請參見擴充Fibonnaci數。
延伸問題：
1. 請問再給了n之後這個程式一共用了多少個乘法，多少個加法。
2. 請用陣列的功能改寫這個程式，請問效率會有所改善嗎？
    請證實您的想法；如果沒有，用陣列豈不沒多大意義？
3. 這個題目不過是Fibonacci數的擴充，而用矩陣可以得到很高效率的程式，
    那麼有沒有辦法用矩陣的方式來解這個題目呢？
    提示：請考慮下式並化簡。
    An+1 = X*An-1 + Y*(An-Fn) + Fn + Fn+1
              = X*An-1 + Y*An + (1-Y)Fn + Fn+1
    這個解法來自Jan Tijmen Udding。
4. 請仔細研究一下我們的程式，或使用的式子，找出可以省略下來的運算，
    並且改良此地的程式(這是一個極為簡單的練習，您應該要試一試)。

個人解答：
1.
從i = 2開始到i = n，一共是n-1個迴圈。
計算fn+1是2乘法1加法，計算An+1是2乘法3加法。
如果不算那個減法，也不把i++當成加法的話，
一共使用了(n-1)(4m+4p)=4(n-1)(m+p)的運算，m表乘法p表加法。
2.
基本上如果我們從小到大計算的話，計算次數仍然會相同，
但會省略掉遞移值的動作，理論上會速度稍微快上一點點，
但記憶體使用空間也算是效率的一環，所以方法不變的話，
只是拿兩者來交換而已。
但還有一個好處是我們可以記錄下前面的Fn和An值，
這會使得重複查詢可以反覆利用。
3.
我們把式子後面兩項再化簡一下：
   (1-Y)Fn + Fn+1
= (1-Y)Fn + Y*Fn + X*Fn-1
= Fn + X*Fn-1
我們可以把式子寫成矩陣的關係式，
左邊的四項為An+1, Fn+1, An   , Fn，
右邊的四項為An    , Fn    , An-1, Fn-1，
則中間的四階矩陣為：
y 1 x x
0 y 0 x
1 0 0 0
0 1 0 0
利用這個矩陣的次方來做像之前那樣子的運算。
值得注意的是，這當中有4個被標藍色的位置無論如何都是0。
有興趣的話可以自行做一次乘法試試看。
最後求出來以後，將第一行的值相乘加總就可以得到An了！
4.
這裡最簡單的省略就是將迴圈結束條件設為i <= n-1。
然後在其後我們算出tmpfn1和tmpAn1後return tmpAn1即可，
這樣我們就不用做最後一次值的遞移。

[C Programming] 擴充Fibonacci數

我們定義一組叫做擴充的Fibonacci 數如下：
已知X與Y兩個數，擴充Fibonacci數Fi為
          1                             i =  0
Fi  =   1                              i  =  i
          X*Fi-2 + Y * Fi-1      i  > 1

因此，當X=Y=1時，這一組擴充的Fibonacci 數就變成一般的Fibonacci 數。
現在請寫一個函數，接受一個n值，不用陣列與遞迴的手法算出下面的結果：

F(0)*F(n)+F(1)*F(n-1)+...+F(i)*F(n-1)+...+F(n-1)*F(1)+F(n)*F(0)

本題出自Jan L.A.van de Snepscheut.

答：
注意這邊原題定義上有些問題，一般來說應該會把X係數給F(n-1)，
但此處剛好相反。原Code裡面按其定義，
則其有一行程式碼應改為temp1 = x * f0 + y * f1;。
我們下面的解答全部都以X係數是在F(n) = X * F(n-2) + Y * F(n-1)為準。

Code：

unsigned long ext_fibonacci(int n, int x, int y)
{
 // Fn+1, Fn-1, Fn , An+1, An-1, An
 unsigned long tmpfn1, fn_1, fn, tmpAn1, an_1, an;
 int i;
 if( n == 0 )
  return 1;
 else if( n == 1 )
  return 2;
 else
  for(an_1=fn=fn_1=1, an=2, i=2; i<=n; i++)
  {
   // Fn+1 = x * Fn-1 + y * Fn, 
   // An+1 = x * An-1 + y * (An - Fn) + Fn + Fn+1
   tmpfn1 = x * fn_1 + y * fn;
   tmpAn1 = x * an_1 + y * (an - fn) + fn + tmpfn1;
   
   // shift the value to next
   fn_1 = fn; fn = tmpfn1;
   an_1 = an; an = tmpan1;
  }
 return an;
}

[C Programming] 快速Fibonacci數算法-延伸(下)

延伸問題：
3. 在上一題與本題的解答中分別有兩個計算Fibonacci數的程式，
    一個很短，很簡單，但另一個卻很長、很複雜，
    請用比較大的n， 並且透過反覆計算1000次的方式，
    看看哪一個程式比較快。

答：
我們利用time.h的函數來取結尾和開頭的時間相減，
分別在兩個主程式中加入以下的程式碼：

Code：

int i;
clock_t start_time, end_time;
float total_time = 0;
start_time = clock();
for( i = 1; i <= 1000; i++) // 做1000次
 fibonacci(8000000);
end_time = clock();
total_time = (float) (end_time - start_time)/CLOCKS_PER_SEC;
printf("Time : %f sec \n", total_time);

結果如下：
Traditional: 32.514999 sec
Matrix: 0.001000 sec

另外測試的狀況，使用矩陣方式數量再多2個層級，
我們才能取到有實際代表性的執行時間。
但再多兩個層級的話一般的做法就會花很久很久的時間了XD

結論：
在n很大的時候，二階的Fibonacci相當節省時間，
這跟前面理論上推導的狀況是一樣的。

[C Programming] 快速Fibonacci數算法-延伸(中)

延伸問題：
2. (續上題) 請評估一下上一個習題程式的時間值，
    當k愈來愈大時，計算時間會以什麼方式增加；
    這個方法與上一個問題中習題3的做法相比又如何呢？

答：
我們每次算k階矩陣的相乘時，
每個元素需要用到k個乘法，k-1個加法。
一共有k^2個元素，如果乘法花費m單位時間，加法為p單位時間，
總時間最多則為(2 log n)*(k^2)((m+p)k-p)。 (底數為2)
在k越大的狀況下這個時間會以k^3的樣子增長。
比較起上一個問題的k階算法，
概念是每次往後推一直推到f(n,k)，
每次要算k項的加總為k-1次的加法，
全部要做的次數為k+1~n，為n-k次，
總共為(n-k)(k-1)p次。

結論：
當計算時是n dominates的時候，現在這個快速算法，
會比前面使用加法來得吃香。
但當k dominates的時候，這個算法的時間就會花的比之前的方法兇。

[C Programming] 快速Fibonacci數算法-延伸(上)

請參見快速Fibonacci數算法。
延伸問題：
1. 請擴充這個方法，讓它可以算出 k 階 Fibonacci 數。
2. (續上題) 請評估一下上一個習題程式的時間值，
    當k愈來愈大時，計算時間會以什麼方式增加；
    這個方法與上一個問題中習題3的做法相比又如何呢？
3. 在上一題與本題的解答中分別有兩個計算Fibonacci數的程式，
    一個很短，很簡單，但另一個卻很長、很複雜，請用比較大的n，
    並且透過反覆計算1000次的方式，看看哪一個程式比較快。

我們今天先解答第一題，剩下兩小題留待明天。

個人解答：
1.
我們在k階的Fibonacci中，
如同原本的方法一樣，將fn....f(n-k+1)和f(n-1)...f(n-k)的關係寫出來。
我們寫出來的矩陣應該會像：
1 1 1 1 ... 1 1
1 0 0 0 ... 0 0
0 1 0 0 ... 0 0
0 0 1 0 ... 0 0
.....................
0 0 0 0 ... 1 0

所以最後化簡到fk ... f1的時候，
我們需要將這個矩陣的(n-k)次方，
而且加總是把第一行的所有值共k個加起來。
這裡我們的m矩陣就是最開始的矩陣，
xm矩陣則是用來承載中間過度的矩陣，
跟原本2階的xa~xd是一樣的概念；
mm則是乘法完成時的矩陣。
由於全部使用陣列，相當於傳遞指標，
所以就不用擔心reference的問題了XD!
原則上下面在做的就是將原本的四行乘法，
擴充到k階的矩陣乘法而已。

別忘了要初始化m矩陣(unik的二維陣列)，
當然如果一次執行只想求一個f(n,k)值，
可以直接令為{0}後，將1的部分填入就好，
但為了檢驗是否中間都有算對，
我們還是得乖一點每次都確實初始化才行~

Code：

#include <stdio.h>
#include <stdlib.h>


#define MAX_K 20
unsigned long unik[MAX_K+1][MAX_K+1];

/**                         (n-2)
 +--------+    +--------+
 |        |    |        |
 |   mm   | =  |   m    |   (m == unik)
 |        |    |        |
 +--------+    +--------+

*/


void kmatrix_power(unsigned long m[MAX_K+1][MAX_K+1], int n, int k,
       unsigned long mm[MAX_K+1][MAX_K+1]        )
{
 unsigned long xm[MAX_K+1][MAX_K+1];
 int i, j, l;
 if(n==1)
 {
  for(i=1;i<=k;i++)
   for(j=1;j<=k;j++)
    mm[i][j] = m[i][j];
 }
 else if(n & 0x01 == 1) // n = 2k+1, break into m*m^2k
 {
  // compute m^2k, then multiply m
  kmatrix_power(m, n-1, k, xm);
  for(i=1;i<=k;i++)
   for(j=1;j<=k;j++)
   {
    mm[i][j] = m[i][1] * xm[1][j];
    for(l=2;l<=k;l++)
     mm[i][j] += m[i][l] * xm[l][j];
   }
 }
 else // n = 2k, break into (m^k) * (m^k)
 {
  kmatrix_power(m, n>>1, k, xm);
  for(i=1;i<=k;i++)
   for(j=1;j<=k;j++)
   {
    mm[i][j] = xm[i][1] * xm[1][j];
    for(l=2;l<=k;l++)
     mm[i][j] += xm[i][l] * xm[l][j];
   }
 }
}

unsigned long fibonacci(int n, int k)
{
 int i;
 unsigned long mm[MAX_K+1][MAX_K+1];
 if( n <= k )
  return 1;
 else
 {
  kmatrix_power(unik, n-k, k, mm);
  unsigned long tmp;
  for(i = 2, tmp = mm[1][1]; i <= k; i++)
   tmp += mm[1][i];
  return tmp;
 }
}

int main(int argc, char *argv[])
{ 
 int n, i, j, k = 5;
 for(n = 1; n <= 30; n++)
 {
  for(i = 1; i <= k; i++)   // first row is 1 1 1 1 ... 1 
   unik[1][i] = 1UL;
  for(i = 2; i <= k; i++)      // then 1 0 0 0 0 ... 0,
  {        // 0 1 0 0 0 ... 0, ...
   for(j = 1; j <= k; j++){
    unik[i][j] = 0UL;
   }
   unik[i][i-1] = 1UL;
  } 
  printf("f(%d) = %ld\n", n, fibonacci(n, k));
 }
 system("PAUSE");
 return 0;
}

[C Programming] 快速Fibonacci數算法

續Fibonacci數列問題。在非遞迴的Fibonacci程式中，
我們在算fn時最多不超過n-2個加法，
請寫作一個速度更快的程式，或許您可以用乘法。
如果每一個乘法用m單位時間，每一個加法用p單位時間，
於是非遞迴的寫法因為最多有n-2個加法，因此最多用(n-2)p時間。
請問，您的改良程式要用多少時間？
假設只考慮加法與乘法而已。

答：

令上面的矩陣為M，在算M^(n-2)時，
我們可以依照n的奇偶數來像之前算m^n一樣來化簡。
           單位矩陣   r=0
M^r =  (M^k)^2     r=2k
           M*(M^2k)  r=2k+1

因此使用遞迴，我們用2 log n次矩陣相乘就可以得到M^n了。(log底數為2)
2x2的矩陣相乘一次需要8個乘法，4個加法，
所以最多總時間不超過(2 log n)*(8m+4p)。
當n夠大的時候，n > log n，其他的m和p和n無關所以可以視為常數。
所以當n非常大了，此方法所用的時間就小於(n-2)p，因此會比較快。

Code：

#include <stdio.h>
#include <stdlib.h>

/**                         (n-2)
 +--------+    +--------+
 | aa  bb |    | a    b |
 |        | =  |        |
 | cc  dd |    | c    d |
 +--------+    +--------+

*/
void matrix_power(unsigned long a, unsigned long b,
      unsigned long c, unsigned long d, int n,
      unsigned long *aa, unsigned long *bb,
      unsigned long *cc, unsigned long *dd    )
{
 unsigned long xa, xb, xc, xd;
 if(n==1)
  *aa = a, *bb = b, *cc = c, *dd = d;
 else if(n & 0x01 == 1) // n = 2k+1, break into m*m^2k
 {
  // compute m^2k, then multiply m
  matrix_power(a, b, c, d, n-1, &xa, &xb, &xc, &xd);
  *aa = a*xa + b*xc;
  *bb = a*xb + b*xd;
  *cc = c*xa + d*xc;
  *dd = c*xb + d*xd;
 }
 else // n = 2k, break into (m^k) * (m^k)
 {
  matrix_power(a, b, c, d, n>>1, &xa, &xb, &xc, &xd);
  *aa = xa*xa + xb*xc;
  *bb = xa*xb + xb*xd;
  *cc = xc*xa + xd*xc;
  *dd = xc*xb + xd*xd;
 }

}

unsigned long fibonacci(int n)
{
 unsigned long a, b, c, d;
 if( n <= 2 )
  return 1;
 else
 {
  matrix_power(1UL, 1UL, 1UL, 0UL, n-2, &a, &b, &c, &d);
  return a + b;
 }

}

int main(int argc, char *argv[])
{ 
 int n;
 for(n = 1; n <= 20; n++)
  printf("f(%d) = %ld\n", n, fibonacci(n));
 system("PAUSE");
 return 0;
}

[C Programming] Fibonacci數非遞迴解-延伸

請參見Fibonacci數非遞迴解。
延伸問題：
1. 使用遞迴的解法會使用return f(n-1)+f(n-2);型態的方式來進行計算，
    這種方式會產生許多重複。
    試以手算的方式來計算f(13)會有多少次的重複計算。
2. (續上題) 請問，算fn時一共用了多少個加法，請導出一條公式。
3. k階的Fibonacci數f1k, f2k, ... , fik, ...是這樣定義的(k是上標)，
    f1k=f2k=...=fkk=1
    fik = f(i-1)k + f(i-2)k + ... + f(n-k+1)k  ( i > k )。
   換言之，fjk是前k個fik的和，因此當k=2時就是一般的Fibonacci數列，
   請不用遞迴(但或許可以用陣列)，寫一個函數fibK(n, k)，
   接收n與k算出fnk。
4. 為什麼函數值要用unsigned long 而不用int呢？請說出你的看法。
個人解答：
1. f(13)  = f(12) + f(11) = 2f(11) + f(10) = 3 f(10)+2 f(9)
             =5 f(9) + 3 f(8) = 8 f(8) + 5 f(7) = 13 f(7) + 8 f(6)
             = 21 f(6) + 13 f(5) = 34 f(5) + 21 f(4) = 55 f(4) + 34 f(3)
             = 89 f(3) + 55 f(2) = 144 f(2) + 89 f(1) = 233
   這當中由於最後都會分解到f(2)和f(1)，
   所以實際分解成233個項次，扣掉f(2)和f(1)，一共重複了231次。
2. Tfn =  fn - 1   , n > 0
3. 利用原題的方法，每次將前面的數加總後遞移往下，
    只是我們需要限定一個陣列最大值K，沒有用到的就不要用。

Code：

#include <stdio.h>
#include <stdlib.h>

#define MAX_K 20
unsigned long fnk[MAX_K+1];

unsigned long fibK(int n, int k)
{
 unsigned long tmp;
 int i, j;
 if( n <= k )
  return 1;
 else
 {
  for(i=1; i<=k; i++)
   fnk[i] = 1;
  for(i=k+1; i<=n; i++)
  {
   tmp  = fnk[k];     // tmp = fnk[1]+...+fnk[k]
   for(j=1; j<k; j++)
    tmp += fnk[j];

   // fnk[1] = fnk[2], fnk[2] = fnk[3], ... , fnk[k] = tmp
   for(j=1; j<k; j++)    
    fnk[j] = fnk[j+1];
   fnk[k] = tmp;
  }
 }
 return tmp;
}

int main(int argc, char *argv[])
{ 
 int n;
 for(n = 1; n <= 20; n++)
  printf("f(%d) = %ld\n", n, fibK(n,4));
 system("PAUSE");
 return 0;
}

4. unsigned long比int大阿廢話，
    簡單說可以放更大的可能值，同時unsigned可以多出一個bit放值，
    因為我們在這個狀況下不需要用到負值。

[C Programming] 數值自乘非遞迴解-延伸1

請參見數值自乘非遞迴解。
*1. 這個程式的乘法數目還可以再降低，
     其實我們保持完整的m, m^2, m^4, m^8, ...是沒多大需要的，
     請以此為出發點來改良程式。
感謝Jennya Chang的建議。
以下來試寫一個尚未驗證能否確實改進效率的方法。
這個方法的核心觀念在於：
如果像15=1111很多1的會很難作，
那麼我們是不是可以算16=10000，回過頭來再扣回1(除以m)就可以了？
整體而言，我們的演算法是：
如果~n的bitsum比較小，就算m^(n+~n+1)/m^(~n+1)，
這裡的~n是以從n的最高位以下的所有bit作反轉。

等式例如： 101011 = 1000000 - ( 010100 + 1)

我們先算n在二進位中有多少個digit，假設為d，
接著將1左移d位就是n+~n+1了!我們把它叫作x。
那麼~n就等於(x-1)和n作XOR。
同時我們就可以計算n的二進位含1的個數bit。
如果bit < d-bit 就直接算m^n，(也就是原來的方法)
否則就算m^x / m^(~n+1)。

這個方法關鍵點在於bitsum計算的速度，
以及最後會多出一個除法。
只要這兩件事情比起少做的乘法數還要快，
程式就會比較節省時間。
等有時間在來計算跑多次的執行時間列出來XD~

Code：

int bitsum(unsigned long n, int *d)
{
 int bitsum = 0;
 *d = 0;
 while(n) // while n > 0
 {
  *d = *d + 1;
  if ( n & 0x01UL == 1 )
   bitsum++;
  n >>= 1;
 }
 return bitsum;
}
unsigned long iterative_power(unsigned long m, unsigned long n)
{
 unsigned long m_x, n_c1, temp = 1, x = 1;
 int bit, *d;
 d = (int*)malloc(sizeof(int));
 bit = bitsum(n, d);
 if( bit < *d - bit)
  while(n > 0)
  {
   if ( n & 0x01UL == 1) // last bit is 1
   temp *= m;
   m *= m;      // m -> m^2 -> m^4 -> m^8 ......
   n >>= 1;     // throw the last bit
  }
 else
 {
  x <<= *d;        // 100000....0
  n_c1 = ((x-1)^n) + 1;   // n_c1 equals ((x-1) XOR n) + 1
  m_x = m;
  while(*d)
  {
    *d  = *d - 1;
   m_x *= m_x;
  }
  while(n_c1)     // while n_c1 > 0
  {
   if ( n_c1 & 0x01UL == 1) // last bit is 1
   temp *= m;
   m *= m;        // m -> m^2 -> m^4 -> m^8 ......
   n_c1 >>= 1;    // throw the last bit
  }
  temp = m_x / temp; // m^x / m^(n_c + 1)
 }
 return temp;
}

[C Programming] Fibonacci數非遞迴解

Fibonacci數列f1,f2,...,fn的定義是這樣的：
        1                        n = 1 or 2
fn =  f(n-1) + f(n-2)    n > 2

請用不用遞迴的方法，也不用陣列，寫一個函數，它接收n的值，傳回fn。

解：
我們只要記錄下fn-1和fn-2，就可以算出fn。
相對來說只要從1往上記錄，每次記錄下前2個，
最後就可以達到fn了!
不用遞迴的好處是節省重複運算，
當然也可以用陣列來記錄遞迴，
這個方法在之前[Algorithm] 以空間換取時間(以Fibonacci為例)裡有提過。
這裡不使用陣列的好處是節省空間，但相對而言，
最後就不會記錄下n-3以前的數據(會全部被蓋掉)。

Code：

unsigned long fibonacci(int n)
{
 unsigned long tmp, fn_2, fn_1;
 int i;
 if( n <= 2 )
  return 1;
 else
  for(fn_2=fn_1=1, i=3; i<=n; i++)
  {
   tmp  = fn_1 + fn_2;
   fn_2 = fn_1;
   fn_1 = tmp;
  }
 return tmp;
}

[C Programming] 數值自乘非遞迴解-延伸

請參見數值自乘非遞迴解。
*1. 這個程式的乘法數目還可以再降低，
     其實我們保持完整的m, m^2, m^4, m^8, ...是沒多大需要的，
     請以此為出發點來改良程式。
 2. 請用手算方式列出計算2^13, 2^15, 2^16的過程，
     是不是這個程式也會得到算2^15時要多幾個乘法呢？
　 請問一般而言，當n是多少時會比較慢？
 3. 如果您用的硬體系統會偵測出整數溢位(Overflow)，
     那麼這個程式可能就會有問題，因為求出m^45時會同時多算了m^64，
     那麼m^64就可能溢位。請修改程式克服這一點問題。

個人解答：
1. 這題不太確定題目的意思。
    因為不論如何還是會經過m的偶數次方。
    如果是指想要快速經過中間二進位的0值的話，
    大概再將3.的答案改成如下：

Code：

while(1)
{
 while( n & 0x01UL == 0 )
  if((n >>= 1) > 0)
   m *= m;
 if ( n & 0x01UL == 1 ) // last bit is 1
  temp *= m;
 if((n >>= 1) > 0)
  m *= m;
 else
  break;
}

這並不是好作法，而且會降低程式可讀性。
如果有人能理解這題到底是想改進什麼的話，
請不吝指教XD~

2. 換成二進位來看時，13 = 1101 , 15 = 1111, 16 = 10000
    我們以只考慮temp*=m和m*=m這兩個乘法作的總次數為前提來看的話，
    1101會做7次，1111會做8次，10000會做6次，
    簡而言之，乘法的總次數相當於(二進位值的總位數+二進位值裡的1的總數)。
    所以在二進位位數同等狀況下，裡面含有越多1，程式就做越慢。
    這點跟之前遞迴解時，有些n會導致乘法次數變多的原因是一樣的。
3. 由於n為正整數，我們可以將判斷迴圈繼續的條件搬到迴圈裡面，
    先判斷n丟掉最右邊的bit後是否還大於0，
    還大於0的話才表示有繼續將m往上計算的價值，
    否則的話我們就可以直接break。

Code：

while(1)
{
 if ( n & 0x01UL == 1) // last bit is 1
  temp *= m;
 if((n >>= 1) > 0)
  m *= m;
 else
  break;
}

[C Programming] 數值自乘非遞迴解

續求m^n問題(m與n是正整數)。前面的提示會得到一個遞迴的程式，
請發展一個非遞迴，但效率相同的程式。

答：
我們先以求m^45為例，45的二進位為101101，
所以m^45 = m^(1*2^5 + 0*2^4 +1*2^3 + 1*2^2 + 0*2^1 + 1*2^0)
我們可以用一個tmp來記錄從最低位算到最高位的值，
先令tmp=1, 當碰到這個位數是1時，就乘上m，
每次位數檢查完以後都讓m自乘，
顯然剛好m -> m^2 -> m^4 -> m^8 ->...，
會符合每一位所代表其為m的正確的次方數。
每次檢查完以後，我們就將n的尾數丟棄(n>>=1)，
記錄到最後n==0時，就離開迴圈輸出值。

Code：

#include <stdio.h>
#include <stdlib.h>

unsigned long iterative_power(unsigned long m, unsigned long n)
{
 unsigned long temp = 1;
 while(n > 0)
 {
  if ( n & 0x01UL == 1) // last bit is 1
   temp *= m;
  m *= m;      // m -> m^2 -> m^4 -> m^8 ......
  n >>= 1;     // throw the last bit
 }
 return temp;
}

int main(int argc, char *argv[])
{ 
 unsigned long m, n;
 m = 4;
 n = 13;
 printf("\n%ld^%ld = %ld\n",m , n, iterative_power(m,n));
 system("PAUSE");
 return 0;
}

[C Programming] 數值自乘遞迴解-延伸

請參見數值自乘遞迴解。
延伸問題：
1. 請用手算模擬求2^13，2^15，2^17。
2. 請算一算程式計算2^7, 2^8, 2^9,...., 2^14所花的時間。
    在時間上有沒有什麼異常的情況？能解釋問題點在哪嗎？
3. (續上題)請以求2^11, 2^13, 2^15(手算)的步驟來解答上一題的異常情形。
*4. 將程式改成不用遞迴的寫法；且不要用堆疊來模擬。
5. 每次把遞迴用的n值印出來或許有助於第四題。
    但是這些值與n的二進位表示法有什麼關聯？
6. 請擴充這個程式，使得m與n可以是負整數或0。
    如何處理那些不正常的情況？

個人解答：
1. 13->1+6*2->1+3*2*2->1+(1+1+1)*2*2
    15->1+7*2->1+(1+3*2)*2->1+(1+(1+1+1)*2)*2
    17->1+8*2->1+4*2*2->1+2*2*2*2->1+(1+1)*2*2*2
    拆分的時候，如果是奇數就要多拆出1，剩下拆成等份，
    偶數則直接拆成等份。
    如此，每一個+號在程式裡都是一個乘法，
    每一個乘號也是。計算總符號個數就知道做了多少次乘法運算。
2. 以1.的方法，
        7   -> 1+(1+1+1)*2         4次
        8   -> (1+1)*2*2             3次
        9   -> 1+(1+1)*2*2         4次

      10   -> (1+(1+1)*2)*2      4次

      11   -> 1+(1+(1+1)*2)*2  5次

      12   -> (1+1+1)*2*2         4次

      13   -> 1+(1+1+1)*2*2     5次

      14   -> (1+(1+1+1)*2)*2  5次
    11、13、14會是花最多時間的，最少的則為8，
    它會比算7來的快。異常情形在於時間上不是單純數字較大就較慢。
    解釋的話，請見3.。

3. 15   -> 1+(1+(1+1+1)*2)*2  6次   

簡單來說當每次拆分無法整除時，只能多花一次乘法去乘以m。

因此這種狀況越多，使用的乘法就會越多。

4. 在下一題的時候會提到。

5. 分解時，相當於在將n轉為二進位，

    只是這個方法並沒有留存n的記錄。

6. 0^0->無意義

    m^0, m!=0 -> 1

    0^n, n!=0 -> 0

    m為負整數的做法和正整數的做法沒有區別，

    剩下n為負整數的部分可先求m^(-n)，

    求完以後再行倒數，注意使用可支援小數點的資料型態如double。

[C Programming] 數值自乘遞迴解

如果n與m是正整數，那麼m^n就是把m連乘n次，
這是一個很沒有效率的方法。請寫一個更有效率的程式，
並且分析您的程式中一共用了多少個乘法。
你應該以少於n-1個乘法做為設計準則。

解：
對於遞迴關係而言，最重要的就是分析n的值的狀況，
然後做出相對應的動作，且要有終止的條件。
就此例而言，我們可以將n每次分成2份，
如果n是偶數的話就可以變成(m^(n/2))*(m^(n/2))，
奇數的話則會多出一個m，
最後分到沒得分(0)的時候就傳回m^0=1。
如此一來遞迴做到最後我們只需要負擔最多2*log n次的乘法(底數為2)，
最少則是log n次。(在前面每次拆分都是偶數的情況下。)

Code：

#include <stdio.h>
#include <stdlib.h>

unsigned long recursive_power(unsigned long m, unsigned long n)
{
 unsigned long temp;
 if(n == 0)       // m ^ 0 = 1
  return 1;
 else if ( n & 0x01UL == 0) // even
 { 
  temp = recursive_power(m, n >> 1); // m, n/2
  return temp * temp;
 }
 else              // odd
  return m * recursive_power(m, n-1);
}

int main(int argc, char *argv[])
{ 
 unsigned long m, n;
 m = 5;
 n = 11;
 printf("\n%ld^%ld = %ld\n",m , n, recursive_power(m,n));
 system("PAUSE");
 return 0;
}

[C Programming] 因數分解-延伸

請參照因數分解。
延伸問題：
1. 程式中都有work>1的測試，究竟這有沒有必要？
    請說明你的理由。
2. 在第二個for迴圈中，我們分別用k=3,5,7,9,11,13,...去除，
    但事實上用3去除時已經去掉所有3的倍數了，
    因此再用9、15等數去除自然不可能再整除，
    所以我們等於做了多餘的事情。請修改這個程式，
    使得只會用質數去除，因此9,15,...等就不會出現了。
3. 請問，這個程式所輸出的因數為什麼都是質數？請證明它。
4. 在程式中，其實factor[]與exps[]，
    這兩個用來存放因數與對應的指數的陣列是不必要的。
    為什麼？請修改程式，不用這兩個陣列，但效果完全相同。
個人解答：
1. 沒有必要。work不斷去除以k的時候，
    最後一定會讓work變成1，從而不符合work%k == 0的迴圈繼續條件。
2. 我們可以利用前面的篩法的概念，先將合成數全部篩掉。
    首先開出陣列且初始化為0，
    簡單的使用倍數的方式，由3開始，將所有合成數都標記為1。
    如此一來，在陣列中所有的奇數且為合成數的都會被標記為1。
    在求因數的時候，先檢查標記是否為0，不是的情況下再行做除法。
    下面Code範例是求出從2到50的因數分解。
Code：

#include <stdio.h>
#include <stdlib.h>

#define MAXSIZE 50
void factor(unsigned long n, unsigned long *fact)
{
 unsigned long i, j, work;
 printf("\n%ld = ", n);
 // Deal with 2 First
 for(i=0, work=n; (work & 0x01UL)==0 ; work>>=1, i++) ;
  
 if(i) // NOT zero
  printf("%ld(%ld)",2,i);
 for(j=3; j<=work; j+=2){
  if(fact[j]==0)
   for(i=0; work%j==0 ; work/=j, i++)
   ;
  if(i)
   printf("%ld(%ld)",j,i);
 }
}

int main(int argc, char *argv[])
{ 
 unsigned long i, j, n = MAXSIZE, fact[MAXSIZE+1]={0};
 for(i=3;i<=MAXSIZE;i+=2)
  if(fact[i]==0)
   for(j=i<<1;j<=MAXSIZE;j+=i)
    if(fact[j]!=1)
     fact[j] = 1; 
 for(i=2;i<=n;i++)
  factor(i, fact);
 printf("\n");
 system("PAUSE");
 return 0;
}

3. 如果輸出的因數具有合成數q=p*r,其中p為質數：
    那麼q>p，所以p在迴圈中一定會先被拿來做為work/=k中的k來使用。
    既然迴圈只要在work%k==0的狀況下就會持續，
    那麼所有存在p為因數的數，在經過k=p的迴圈時，p的因子皆會被除盡，
    因此q再拿來當作k來嘗試時，不可能除盡，因此輸出中不會出現，矛盾。
    故不存在合成數能作為輸出的因數。
4. 在基本題的用法就是不用陣列的做法。
    明顯的優點是節省空間，缺點就是只能輸出而沒有記錄。

[C Programming] 因數分解

對於一個大於2的正整數而言，
我們一定可以將其分解成一或多個質因數的乘積和。
請寫一個程式，將輸入的數分解，
輸出格式以括號表示該因數的次方數，
例如：輸入:240  輸出:2(4)3(1)5(1)     (即240=(2^4)*3*5)

解：
最基本的想法，是我們從2開始，嘗試將每一個數都去除它，
此時會有幾種狀況：
1.餘數等於0 => 整除，繼續使用同一個數嘗試去除。
2.餘數不等於0 => 不整除，將嘗試的數遞增。
持續上面的動作一直到值被除到等於1的時候，
就表示分解完畢可以輸出結果。

另一方面，由於2的倍數除了2以外都是質數，
我們可以讓嘗試的數以2、3、5、7、9...的方式來嘗試。

以這個方法而言，我們可以保證得出來的因數分解，
絕對是質因數分解。
何以見得？
(下面所提的所有數全部都是正整數)
首先我們不用考慮2的倍數會不會出現，因為迴圈中本來就不使用。
再者，如果合成數p在迴圈中被使用，
且有作為因數分解到，
設p = r*q ，r為質數，q為1以外的正整數。
既然如此，那麼因為r < p，所以迴圈中一定會先使用到r來除。
如此一來，假設所求數n = (r^i) * k， ( i >= 1, k>=1)，
那麼r^i的部分一定會被r除盡。
對於p來說，它並沒有辦法整除剩下的k；
因此產生矛盾，所以我們得證在因數分解中，
這個方法最後輸出的所有因數均為質因數。

在程式中，我們可以選擇每次找到一項質因數及其平方，
就進行輸出，如此一來，我們就不需要使用陣列來記錄。
優點當然是會省去記錄質因數及次方數的陣列，
缺點則是輸出過以後沒有任何留存。
我們直接進行輸出，這是延伸問題的第4題所要求的，
如果需要記錄陣列的話，可以比照類似昨天線性篩法附帶的因式分解來處理。
請見線性篩法-延伸(下)。
不過，在延伸問題2的時候有可能會用到需要記錄的架構，
到那時候再用那個方式處理吧XD

Code：

#include <stdio.h>
#include <stdlib.h>

void factor(unsigned long n)
{
 unsigned long i, j, work;
 
 printf("\n%ld = ", n);
 // Deal with 2 First
 // work和0x01UL(也就是1)進行&運算也就是work的最後一位為0時這項才為真(偶數)。
for(i=0, work=n; (work & 0x01UL)==0 && work > 1; work>>=1, i++) ;
 if(i) // NOT zero
  printf("%ld(%ld)",2,i);
 for(j=3; j<=work; j+=2){
  for(i=0; work%j==0 && work > 1; work/=j, i++)
  ; // 這個分號是為了避免if被當作是內層for所要做的事情。
  if(i) // 有分解到的狀況就輸出這項的分解
   printf("%ld(%ld)",j,i);
 }
}

int main(int argc, char *argv[])
{ 
 int i, n = 500;
 for(i=2;i<=n;i++)
  factor(i);
 system("PAUSE");
 return 0;
}

[C Programming] 線性篩法-延伸(下)

3. 在這裡我們來看看如何修改本問題，
    而同時求出2到n之間每一個數的因數分解，
    這是Dexter Kozen提出來的方法。
    因為每個合成數r都可以寫成p^i * q的形式，
    p為r最小的質因數。
    當要刪除r時，我們在程式中就會知道p、i、q的值，
    我們可以準備三個陣列來存P[r]、EXP[r]、Q[r]。
    只要先將P[]初始化成0，在程式執行完後，
    P[r]=0的r肯定是質數，其他的狀況則為合成數。
    那我們知道r^i次方的部分，接著再往下找P[Q[r]]，
    一路找下去就可以將r給分解完成。

個人解答：
我們將輸出的格式定為如下：
60 = 2(2)3(1)5   (括號內的數為次方項)
在這種狀況下，我們為了要求每一個數的因數分解，
那麼就沒辦法無視2的倍數了。
所以前面只用next方法依舊是可以用，
但初始化時要改為每次+1，且LEFT(x) 改為x-1。
這裡就直接使用最原本的線性篩法來解答。
在三重迴圈中，除了相應的陣列外，
我們多加了一個pcnt，每次從新的q開始時就要重置，
而在第三層迴圈是每次乘上pcnt，自然要進行遞增。
當我們要進行刪去時，
有兩種情況：
1. prime == q
這種情況下p^pcnt*q相當於p^(pcnt+1)，
因此給值P[mult]=prime,EXP[mult]=pcnt+1,Q[mult]=1
2. prime  <  q
給值P[mult]=prime,EXP[mult]=pcnt,Q[mult]=q

然後在讀取每一個數的時候，分下列三種狀況：
1. P[i] == 0
主程式中由於我們會先用迴圈對陣列初始化為0，
因此當碰到0的狀況時，表示它已經是質數了，
因此就輸出這個數本身即可。
2. Q[i] == 1
這個狀況代表說這個數為p的EXP[mult]，
那麼只要輸出P[mult](EXP[mult])即可。
3. else
除此以外的狀況，表示在將p^i的部分表示出來以後，
q還可以被分解，那麼就用遞迴的方式，
將Q[mult]做為參數傳遞進去，最後就可以完全分解。



Code：

void L_Sieve(unsigned long n,unsigned long *P,unsigned long *Q,unsigned long *EXP)
{
 unsigned long prev[MAXSIZE+1], next[MAXSIZE+1];
 unsigned long prime, q, i, mult, cnt = 0;
 unsigned long pcnt;
 INITIAL(n);
  
 for(prime=2; prime*prime <= n; NEXT(prime))
  for(q=prime, pcnt=1; prime*q <= n; NEXT(q), pcnt=1)
   for(mult=prime*q; mult <= n; mult*=prime, pcnt++)
   { 
    if(prime==q){
     P[mult] = prime; EXP[mult] = pcnt+1; Q[mult] = 1;
    }
    else{
     P[mult] = prime; EXP[mult] = pcnt; Q[mult] = q;
    }
    REMOVE(mult);
   }
 
 for(i = 2; i!= NULL; NEXT(i))
 {
  if(cnt % 8 == 0) printf("\n");
  printf("%6ld", i);
  cnt++;
 }
 printf("\n\nTotal Prime Count: %ld \n", cnt);
}
void factor(unsigned long i, unsigned long *P, unsigned long *Q, unsigned long *EXP){
  if(P[i]==0) 
   printf("%ld", i);
  else if(Q[i] == 1){
   printf("%ld(%ld)", P[i], EXP[i]);
  }
  else{
   printf("%ld(%ld)", P[i], EXP[i]);
   factor(Q[i],P,Q,EXP);
  } 
}

int main(int argc, char *argv[])
{
 // Initializing.
 unsigned long P[MAXSIZE+1], EXP[MAXSIZE+1], Q[MAXSIZE+1];
 unsigned long i;
 for(i=2;i<=MAXSIZE+1;i++){
  P[i]=EXP[i]=Q[i]=0;
 }
 
 L_Sieve(MAXSIZE,P,Q,EXP);
 for(i=2;i<=MAXSIZE;i++){
  printf("%ld = ",i); 
  factor(i,P,Q,EXP);
  printf("\n");
 }
 system("PAUSE");
 return 0;
}

[C Programming] 線性篩法-延伸(中)

(延伸問題)
2. 此地的記憶體用量明顯比上一題的傳統篩法大許多。
    可否將它作節省呢？
　首先我們不要prev[]只留下next[]，那對於next[i]來說，
    如果i-1沒有被刪除，那i的上一個是i-1，但若i-1已經被刪除了，
    我們用next[i-1]來存原來的prev[i]，換句話說，
    i的上一個就是i-1與next[i-1]中值較少的那一個。
    請問，還有什麼邊界條件要考慮，並且將此概念寫成程式。
    這個技巧在第一個問題之下能否適用呢？

個人解答：
以下的解法非常非常非常(沒錯，就是三個非常)容易搞混，
所以請務必仔細的讀每一句話，並且留意：
"值"和"索引值"的差別！
首先，在第一個延伸問題的前提下，這個技巧仍然可以適用，
只不過，index的基礎移動變為2了。
原本的的核心概念是：
當自己被刪掉時，告訴自己的上一個數:"你的下一個變成XX"，
再告訴自己的下一個數:"你的上一個變成OO"。
這樣子就結束了。
但在現在，我們只剩下next可以放"下一個"，
所以在上一個數的存放的時候，我們要利用已被"刪除"的數的空間。
整體而言，當一個數被刪掉的時候，我們應該要做三件事：

1. 把next[x]值，存給前一個index的next[]。
2. 在index為next[x]的位置的前一個(-2)，存入前一個位置的index。
3. 將自己的位置，也存入前一個位置的index。

在這個機制下，要如何判斷一個數x的前面一個是哪位的方式，
就是先找往前2單位的index，
看它裡面的值是否比index大，較大代表未刪去，較小代表已刪去。
前者的狀況表示x-2即是x的前一個，後者的狀況，
x-2裡面存放的，就是x的前一個數的index。

很難理解是嗎？
我們舉幾個例子：
考慮3、5、7、9、11，一開始是
index  3  5  7  9    11
value  5  7  9  11  13
9被刪去的時候，我們做上述的三件事：
1. 告訴前一個數，你的next[]要改存11了。
    9-2是7，next[7]=9，表示它目前是沒被刪除的，
    所以我們把11存到next[7]裡面。
2. 在9的下一個位置next[9]的前一個(-2)位置，存入前一個位置(7)的index。
    也就是在11-2=9的地方存入7。
3. 將自己的位置，也存入前一個位置的index。(也就是在9的地方存入7)
這當中其實3可以不用做結果也一樣，
會這麼做只是因為在定義上來說，我們希望只要是已刪去的數，
儲存的值都要比自己小。
做完以後的樣子會變成：

index  3  5  7    9    11
value  5  7  11  7    13

再舉個例子，考慮31、33、35、37、39，
篩去的順序會是33->39->35，依序做完的狀況如下：

index    31   33   35   37   39
value    33   35   37   39   41

after 1  35   31   37   39   41  
after 2  35   31   37   41   37
after 3  37   31   31   41   37
過程不再贅述，就是仍然按照上面三個步驟進行就是了。
或許有人會說，看起來2和3是做相同的事情呀！
你每次繞來繞去，最後不就是把自己改掉就好了，
那我們不做2直接做3就好啦！

2和3兩個情況，在大多數是相同的，但還是有例外：
以下的刪去順序是先27再25，因為3比較早篩。

index    23   25   27   29
value    25   27   29   31

after 1  25   29   25   31
after 2  29   23   23   31
發現了嗎？因為29要知道自己的前面是23的話，
我們必須告訴27，要27裡面放23進去，
這時候29-2是27，而不是25，如果我們只改掉25裡面的值，
就會產生問題。

除此以外還有一個很重要的邊界條件：
在原問題中，我們使用NULL做為邊界。
在C語言中，NULL其實就是指0的意思，
但記得我們的迴圈是怎麼判斷繼續的嗎？
沒錯，就是乘積<n。
使用NULL的狀況下，
在我們這樣子把值指定給"刪去"後的index時，
會發生我們會讀到NULL的狀況，
這時候mult會一直等於0，也就一直跳不出去，
而且偵測到NULL就break的方式，
也會受限於0的值會傳遞到奇怪的地方而有時會出錯。
怎麼辦呢？
既然用0不行，那我們將其設成MAXSIZE+1如何？
如此一來，當跑到最後面，迴圈肯定能夠跳出來。
當然缺點是能跑的值會更加受限於資料型態的最大範圍，
但這個方法本來就會受限於乘法隨數字增大，
而令計算時間變長的問題，其實也不能跑多大的數字XD

判斷前一個值的時候，我們使用
if(LEFT(x) < next[LEFT(x)])   當中LEFT(x)可以定義成x-2或x-1，
端看今天我們有沒有先行將2的倍數處理掉。
如果前一位(i-2)的index比所存的value小，代表沒被刪掉，
所以i-2就是i的前一個。
否則，就是i-2所存的value值，是i的前一個的index。

Code的部分，我們將prev相關的通通除去，
而在最後的走訪print出所有質數的時候，
將條件改成：for(i = 2; i!= MAXSIZE+1; NEXT(i))即可。
為了方便一點閱讀，多設了一個變數tmp，
讓if和else執行的狀況看起來更有條理。

Code：

/* Coded by Desolve(雨蕭) 2013.2.9
   REMOVE先判斷前一個的位置後，做三件事情：
   1. 把next[x]值，存給前一個index的next[]。
   2. 在index為next[x]的位置的前一個(-2)，存入前一個位置的index。
   3. 將自己的位置，也存入前一個位置的index。
   第3點可以不用做，這只是為了概念上較容易懂；
   但注意2和3不完全等義，所以不能只做3而不做2。
*/
#define REMOVE(x) {          \
 if(LEFT(x) < next[LEFT(x)]){     \
  next[LEFT(x)] = next[x]; \
  next[LEFT(next[x])] = LEFT(x); \
  next[x] = LEFT(x);  \
 }else{    \
  tmp = next[LEFT(x)];   \
  next[tmp] = next[x]; \
  next[LEFT(next[x])] = tmp; \
  next[x] = tmp;    \
 }\
}
#define LEFT(x) x-2

#define INITIAL(n) { unsigned long i;    \
    for(i=3; i<=n; i+=2)      \
     next[i]=i+2;    \
    next[2]=3, next[n]=MAXSIZE+1;  \
  } 

[C Programming] 線性篩法-延伸(上)

請參照線性篩法。
由於2的部分我還需要花點時間確認其正確性，
因此將其分開來解答。
敬請期待XD~

延伸問題：
1. 這個程式花了一半的時間將2到n之間的偶數刪去，
    但這其實是多餘的，因為除了2以外，
    所有的偶數都不是質數。
    請修改此地的程式，讓它不處理2以外的偶數，
    因此省下一半的時間。(2與n之間一半元素是偶數)
    可參照前面使用一般篩法的想法。
2. 此地的記憶體用量明顯比上一題的傳統篩法大許多。
    可否將它作節省呢？
　首先我們不要prev[]只留下next[]，那對於next[i]來說，
    如果i-1沒有被刪除，那i的上一個是i-1，但若i-1已經被刪除了，
    我們用next[i-1]來存原來的prev[i]，換句話說，
    i的上一個就是i-1與next[i-1]中值較少的那一個。
    請問，還有什麼邊界條件要考慮，並且將此概念寫成程式。
    這個技巧在第一個問題之下能否適用呢？
3. 在這裡我們來看看如何修改本問題，
    而同時求出2到n之間每一個數的因數分解，
    這是Dexter Kozen提出來的方法。
    因為每個合成數r都可以寫成p^i * q的形式，
    p為r最小的質因數。
    當要刪除r時，我們在程式中就會知道p、i、q的值，
    我們可以準備三個陣列來存P[r]、EXP[r]、Q[r]。
    只要先將P[]初始化成0，在程式執行完後，
    P[r]=0的r肯定是質數，其他的狀況則為合成數。
    那我們知道r^i次方的部分，接著再往下找P[Q[r]]，
    一路找下去就可以將r給分解完成。

個人解答：
1. 因為目的是省時間，我們只要在一開始的時候，
    將每一個next和prev跳著接就好了。
    也就是說，每次的next[i]=i+2, prev[i]=i-2;
    除了要注意next[2]=3, prev[3]=2,
    而且最後的n要視其是否為偶數來定位prev。
    這樣相當於我們完全直接跳過了偶數的部分(直接刪掉。)
Code：
首先我們在INITIAL(n)之前先加入此式使其成為奇數：

n -= (n-1)%2;  // make sure odd

接下來，把prime=2的初始條件改為3，並且更改INITIAL：

#define INITIAL(n) { unsigned long i;    \
    for(i=3; i<=n; i+=2)      \
     prev[i]=i-2, next[i]=i+2; \
    prev[2]=next[n]=NULL;   \
    next[2]=3, prev[3]=2;   \
  }

[C Programming] 線性篩法

在做篩法求質數的問題時，我們知道在刪除非質數的資料時，
有很多是重複刪除的；比如說如果有一個數是3*7*17*23，
那麼在刪除3的倍數時會刪到它，刪除7、17、23的倍數時也會刪它，
當質數範圍很大的時候，重複的狀況就會越來越多，
因而程式效率就會大打折扣。
基於篩法的觀念，請寫一個程式，在刪除非質數時「絕對」不做重複的工作。

解答：
我們使用所謂的線性篩法，
基本想法是先找第一個質數(2)，刪除2^2、2^3、2^4...
接著刪除2*3、2^2 *3、2^3 * 3、.......
接著是2*5、2^2 * 5、2^3 * 5、......以此類推
每次我們都先找一個質數p，接著找下一個比這個質數大，
但還沒有被刪去的數q(注意q不一定是質數!)，
刪去p^i *q (i=1,2,...)。由於每次刪去的數都為當時的q*(質數p的i次方)，
所以當找到下一個沒被刪去的q'時，不存在q* p^i = q' * p^i'的可能性。
因為假如等號成立，q' = q* p^(i-i')，i != i'，所以q' = q*(p^j)，
如此一來，在以q作刪去的時候就會刪掉q'了，矛盾。
依此推論，這個篩法不會去篩重複的數。

對於所有合成數來說應該都可以寫成如下的形式：
r = p^i * q，p為r的最小的質因數。
由於我們的q是使用所有大於p且沒有被刪去的數，
也就是說對於所有含有質因數p在內的數，應該都會被篩掉。
因為從頭開始篩到x，因為不論如何這之前的數都會進篩選中，
也就是所有x之前的質因數都被選為p給篩過去了，
所以小於x的質數全部都進行篩選過，
那麼現在x只能是質數，否則它就會被篩掉。
換句話說，這樣一個一個掃過去，選到的p必定是質數。

這是一個很不嚴謹的說明而非證明，
只是嘗試將線性篩法的概念表達出來。
想要了解更完備的證明，請參照這篇paper：
A Linear Sieve Algorithm for Finding Prime Numbers
在實作上的部分來說，
要找到2到n的所有質數，那麼p^2最大不能超過n，
同理，p^i*q也不能超過n，
篩選部分的函式大概可以寫成如下的三層迴圈：
for(p=2; p*p<=n ; NEXT(p))
     for(q=p; p*q<=n; NEXT(q))
         for(m=p*q; m<=n; m=p*m)
             REMOVE(m);

我們可以用Double Linked List來做數列的連結和刪去。
這當中往前找上一個值的prev[i]的值為i-1，
往後找下一個值的next[i]的值則為i+1，
prev[2]和next[n]則到頭尾了，設為NULL。
刪除某個數m呢?
我們就把prev[m]的next指向m的下一個，
即next[prev[m]] = next[m]；
next[m]的prev指向m的上一個，
即prev[next[m]] = prev[m]。
最後，要抓下一個待篩的元素，
就令 x = next[x]，就可以得到x的下一個元素了。

Code：

#include <stdio.h>
#include <stdlib.h>

#define MAXSIZE 1000
#define NEXT(x) x = next[x]
#define REMOVE(x) { next[prev[x]]=next[x]; prev[next[x]]=prev[x];}
#define INITIAL(n) { unsigned long i;    \
    for(i=2; i<=n; i++)      \
     prev[i]=i-1, next[i]=i+1; \
    prev[2]=next[n]=NULL;   \
  } 

void L_Sieve(unsigned long n)
{
 unsigned long prev[MAXSIZE+1], next[MAXSIZE+1];
 unsigned long prime, q, i, mult, cnt = 0;
 INITIAL(n);
 for(prime=2; prime*prime <= n; NEXT(prime))
  for(q=prime; prime*q <= n; NEXT(q))
   for(mult=prime*q; mult <= n; mult *= prime)
    REMOVE(mult);
 for(i = 2; i!= NULL; NEXT(i))
 {
  if(cnt % 8 == 0) printf("\n");
  printf("%6ld", i);
  cnt++;
 }
 printf("\n\nTotal Prime Count: %ld \n", cnt);
}

int main(int argc, char *argv[])
{
 L_Sieve(1000);
 system("PAUSE");
 return 0;
}

[C Programming] 篩法-延伸

請參照篩法問題。
延伸：
1. 為什麼在程式工作過程中，如果sieve[i]為KEPT時，
　2i+3就是一個質數？
2. 請分析一下程式一共查過了多少個元素。
3. 連3的倍數都不處理的話，每6個數剩下的就只有6n+1, 6n+5而已，
    換句話說只要處理原來的三分之一的元素就行了。
    請把這個觀念寫成程式。

額外(個人討論)：
4.將程式改寫成只篩到N/2就結束，同樣要輸出質數及總數目。

個人解答：
1. 在查到第i個的時候，表示已經用所有前面的數篩過了，
    既然到這裡沒有被篩掉，那麼就表示其沒有質因數，
    所以2i+3就是一個質數。
2. 首先程式要每一個數經過(共n次)，
    再來每次濾掉2i+3的倍數，
    就是n/3 + n/5 + n/7 + ...也就是n*(n以內質數的倒數和-1/2)。
    兩者相加就剛好是n*(n以內的質數的倒數和)。
    n越大的時候其實它會發散就是了= =......
3.
對於前一個程式而言，x[i]=2i+1，
我們每一次移動1單位的index跳的值是2。
如果我們要直接去用6i+1和6i+5來記錄，
單位的移動會變得很難掌握。
我們可以改成，
一開始將全部都設為KEPT，
接著把所有3的倍數再設為DELETED。
除3以外剩下為KEPT的是5、7、11、13、17、...
也就是從x[1]開始每次值是2、4、2、4的移動，
index則是1、2、1、2的移動。
換句話說，
在x[]中的初始樣態應該長這樣：
index:     0  1  2  3  4   5   6   7   8   9  10  11
代表值:  3  5  7  9  11 13 15 17 19 21 23 25
K/D:      K  K  K  D  K   K  D   K  K   D   K   K

接下來，在迴圈中的移動，
我們就直接跳著來做篩選就好，
比照之前的dist，這次初始的i是1，
dist是1，做完後要以dist = 3-dist來改變位移量。
同時參照4.來控制範圍，節省的更多。
Code部分的變化:

 int range = n/2, dist = 1, cnt = 1, KEPT = 1, DELETED = 0;
 for(i = 0; i <= n; i++)
  x[i] = KEPT;
 for(i = 3; i <= n; i+=3)
  x[i] = DELETED;
 i = 1;
 while( i <= range )
 {
  if(x[i])  // If KEPT -> start to sieve
  {
   tmp = i+i+3; // Each time index moves 2i+3
   for(j = i+tmp; j <= n ; j += tmp) 
    x[j] = DELETED;
  }
  i += dist;
  dist = 3 - dist;
 }
   

4. 宣告一個變數range = n/2，並且將for迴圈條件改成for(i = 0; i <= range; i++)。

    篩完以後，再使用迴圈加總count。
    相較於原程式，多了一次走訪全陣列，少去了篩大於n/2的倍數，
    整體而言，當n越大，這個方法就越節省時間。
額外增加的Code：

 for(i = 0; i <= n; i++)
  if(x[i]) 
   cnt++;

[C Programming] 篩法

關於求質數，早在2000年前，
人們就知道了一個方法可以不必用除法，
就可找出從2到N之間的所有質數。
假設我們有一個很神奇的篩子，我們可以給它一個數i，
它可以把i的所有倍數去掉。
請用這個方法求出2到N之間的所有質數。

注意：
要求2到N的質數的話，最多篩到N/2就可以停下來了。
因為大過N/2的數都不可能整除N。
(因為存在這樣的數的話，表示有一個質數小於2且整除N)
程式不可使用乘法或除法，只能用加或減，以求加快速度。

這個方法叫做Eratosthenes(人名)的篩法(Sieve Method)。

解：
首先我們建立陣列，x[0] = 3、x[1] = 5、x[2] = 7......
(因為2的倍數可以直接不用考慮，或者說一開始就篩掉了)
所以x[i] = 2i+3 = i + i + 3。
我們可以用標記來表示x[i]是KEPT(保留，定為1)或DELETED(篩去，定為0)。
當一個一個查過去，沒有被篩掉的必然是質數，
因為其之前不存在任何數可以篩掉這個數。


接下來我們就要考慮碰到x[i]的時候，要篩掉2i+3的倍數：
1(2i+3)、3(2i+3)、5(2i+3)、7(2i+3)、9(2i+3)....(2n+1)(2i+3)  n=0,1,2,...
偶數的部分本來就不在裡面，
又因為自己不可以篩掉。
所以每次要增加2(2i+3)的值，
即陣列的索引值每次要增加2i+3。
我們依此來標記所有對應的索引值為篩去。
同時，如果只是要把範圍內的合成數篩掉而不計數，
理論上來說不需要篩到最後一個。這個部分留待延伸討論。
Code中其實可以用把KEPT和DELETED放在#define裡，
不過我不太喜歡這種用法，所以就算了XD

Code：

int Sieve(int* x, int n)
{
 // cnt = 1 because we don't have 2 in the array.
 int i, j, tmp, cnt = 1, KEPT = 1, DELETED = 0;
 for(i = 0; i <= n; i++)
  x[i] = KEPT;
 for(i = 0; i <= n; i++)
 {
  if(x[i])  // If KEPT -> start to sieve
  {
   tmp = i+i+3; // Each time index moves 2i+3
   cnt++;   // prime count
   for(j = i+tmp; j <= n ; j += tmp) 
    x[j] = DELETED;
  }
 }
 return cnt;
}

主程式部分：

int main(int argc, char *argv[])
{
 int i, k, count, n = 500;
 int *x = (int*)malloc(sizeof(int)*(n+1));
 count = Sieve(x, n);
 printf("Total: %d prime(s)", count);
 printf("\n%6d", 2);
 for(i = 0, k = 2; i <= n; i++)
 {
  if(x[i]) // KEPT
  {
   if(k > 10)
   {
    printf("\n");
    k = 1;
   }
   printf("%6d", 2*i+3);
   k++;
  } 
 }
 printf("\n");
 system("PAUSE");
 free(x); x = 0;
 return 0;
}

[C Programming] 求質數-延伸

請參見求質數問題。
延伸問題：
1. 跳過2、3、5的倍數也是可能的，試構想作法並說明其複雜程度，
    以及是否值得寫到程式中？
2. 程式中把5放到prime[]中，是否為多此一舉？為什麼？
    不這樣做的話，有沒有什麼好的替代方案？
3. 有人說，is_prime這個變數是多餘的，其實也是如此，
    你有沒有辦法修改這個程式讓它不需要用is_prime呢？
    程式會比較容易讀嗎？
4. 修改程式讓它可以求出1~N(使用者輸入)中的所有質數。
    請問如果程式再改成讀入M與N(M<N)，
    求出M到N之間的所有質數時，程式的工作量可以少一點嗎？
    為什麼？

個人解答：
1. 考慮到5的倍數的話，
    為: 30n+1、30n+7、30n+11、30n+13、30n+17、30n+19、30n+23、30n+29
    間隔差距為：6 4 2 4 2 4 6 2 的循環，我們可以定一個dist的陣列：
    dist[8] = {6,4,2,4,2,4,6,2};
    前面的數為2、3、5、7，distCnt = 1，
    接下來每次遞增distCnt，用除以8的餘數來達到間距切換的目的。
    也就是每30個數裡面要判斷8個，
    這和原先的每6個數裡面要判斷2個比較起來，
    每經過30個數，我們就可以少判斷2個。

Code：

int* nAdvPrime(int n)
{
 int *prime = (int*)malloc(sizeof(int)*n);
 int i, isPrime, distCnt = 1, totalCnt = 4, primeCand = 7;
 prime[0] = 2; prime[1] = 3; prime[2] = 5; prime[3] = 7;
 int dist[8] = {6,4,2,4,2,4,6,2};
 // isPrime: 1 -> TRUE, 0 -> FALSE.
 while(totalCnt < n)
 {
  primeCand += dist[(distCnt++)%8];
  isPrime = 1;
  for( i = 0; prime[i]*prime[i]<=primeCand; i++ )
   if(primeCand % prime[i] == 0)
   {
    isPrime = 0; break;
   }
  if(isPrime) 
   prime[totalCnt++] = primeCand;

 }
 return prime;
}

2. 5先放到prime[]裡可以節省一點時間，認真說起來，
    我們是可以從5開始檢查，但比較浪費而已。
3. is_prime可以視作是多餘的，我們可以在迴圈結束的時候，
    再次以if(!primeCand % prime[i] == 0)來檢查，成立的話就是質數，
    這樣就不需要用到is_prime。每輪會少掉兩個指定is_prime值的操作，
    但會多出一次求餘數的檢查。(而且也比較不好理解)
    或者，我們可以把for迴圈的部分寫成一個小函式，
    在確認為合成數的時候直接return 0，或者迴圈跑到最後發現是質數的話，
    就return 1。這樣我們可以直接寫成if(funcIsPrime(...))的形式，
    程式碼會更好讀，唯一不確定的就是這樣再度呼叫函式的速度會不會受影響。
4. 求出1~N(使用者輸入)中的所有質數時，只要把while迴圈的條件改成：
    primeCand <= N，然後陣列的大小開成(N*8)%30 + 2即可，
    因為以30的倍數來看，除第一個循環外，每30個數最多有8個數有可能是質數，
    而第一輪有10個質數(2,3,5,7,11,13,17,19,23,29)，所以多加上2。

    求出M到N之間的所有質數時，程式的工作量沒辦法變少，
    因為還是要把這之前的質數找出來。
   

[C Programming] 求質數

請寫出一個程式，找出前N個(比如說200)個質數。
請使用除法的方法用最快的辦法來寫作程式。

解：
1. 質數的基本定義為除了1以及自己以外沒有其他的因數的數，
    就可以稱為質數。
2. 沒有其他因數的條件其實等同於沒有其他質因數，
    因為合成數可以拆成質數的乘積。
3. 另外，當某數x是合成數的話，那麼就會有x=i*j這樣的式子，
     i、j為x的因數；
     i、j要嘛就是一個大於根號x，一個小於根號x，
    不然就是兩個都等於根號x，因為兩者相乘必須要等於x，
    所以不能同時小於或同時大於根號x。
    所以我們只要找小於等於根號x的正整數中是否有x的因數，
    就可以確信x是否為質數。
4. 既然所有合成數都可以拆成質因數的乘積，
    當搜尋一個數是否為質數時，只要看它前面存在的質數，
    是否是它的因數，就可以判定了。
5. 有一些判斷是可以直接省略的，
    比如2的倍數，除了2以外，絕對不會是質數；
    同理，3的倍數也是相同，5的倍數也是。
    當考慮不是2的倍數和3的倍數時，
    剩下6n+1和6n+5兩種。
    考慮到5的倍數的話，為:
    30n+1、30n+7、30n+11、30n+13、30n+17、30n+19、30n+23、30n+29
    間隔差距為：6 4 2 4 2 4 6 2 的循環，但這比較麻煩，
    我們就先做6n+1和6n+5的循環，其他留待延伸題解決。
    這種循環的間隔為2、4、2、4......。
    觀察前面的質數：2 3 5 7 11不難發現可以從5到7這裡開始迴圈。
    由於只有兩種間隔，我們設定dist=2、增加到待檢數字後，
    用dist = 6-dist來變換，這樣就可以不停切換兩種間隔。

Code：

int* nPrime(int n)
{
 int *prime = (int*)malloc(sizeof(int)*n);
 int i, isPrime, dist = 2, totalCnt = 3, primeCand = 5;
 prime[0] = 2; prime[1] = 3; prime[2] = 5;
 // isPrime: 1 -> TRUE, 0 -> FALSE.
 while(totalCnt < n)
 {
  primeCand += dist;
  dist = 6-dist;
  isPrime = 1;
  for( i = 0; prime[i]*prime[i]<=primeCand; i++ )
   if(primeCand % prime[i] == 0)
   {
    isPrime = 0; break;
   }
  if(isPrime) 
   prime[totalCnt++] = primeCand;
 }
 return prime;
}

在main中:

int main(int argc, char *argv[])
{
 int i, n, *prime;
 scanf("%d", &n);
 printf("\n");
 prime = nPrime(n);
 for(i = 0; i < n; i++)
  printf("%4d: %d\n", i+1, prime[i]);
 system("PAUSE");
 return 0;
}

[C Programming] Armstrong數-延伸

請參照Armstrong數。
延伸:
如果是要寫出輸入n的位數，求出所有n位數的數字，
等同於其每一位的n次方的和的Armstrong數呢？
可以不要用到exp(), log()等函數嗎？
資料型別要注意什麼？

解:
運用之前個人的方式，
先將0~9的n次方存起來，
接著再將每個位數一次一次拆解以後取對應的n次方值加總。
用這個方法我們也可以算從1位到n位的Armstrong數。
只要運用迴圈，就可以有效的利用0~9的n次方值。
資料型別的部分，
由於C內建最長的型態應該是unsigned long long int，
其大小為2^64 - 1 = 18446744073709551615
而最長的Armstrong數為:

No.88 11513 22190 18763 99256 50955 97973 97152 2401

有36位，所以其實以內建的長度來說只能到20位。
但......跑到10位左右其實就已經會花很長的時間了，
我相信除了閒著沒事的人應該不需要它才是Orz......

使用long long int的話，printf的參數為%lld，
unsigned long long int的話應該是%llu。(若有錯請不吝指正)

下面Code分成只算單n位的，以及從1位數算到n位的Armstrong數。

單n位的Armstrong數:

int armstrong(int n)
{
 unsigned long long int sum, num, i, ubnd=9, lbnd=1;
 int count, quoti, totalCnt = 0;
 unsigned long long int* uniE = 
   (unsigned long long int*)
     malloc(sizeof(unsigned long long int)*10);
 for(count = 0; count <= 9; count++)
  uniE[count] = count;
 for(count = 2; count <= 9; count++)
  for(i = 2;i <= n; i++)
   uniE[count] *= count;
 for(i = 2; i <= n; i++)
 {
  lbnd *= 10;   // lower_bound為1 00.....0
  ubnd += 9*lbnd;  // upper_bound為9 99.....9
 }
 for(i = lbnd; i <= ubnd; i++)
 {
  num = i; sum = 0;
  while(num != 0)
  {
   quoti = num % 10; // 每次取個位數加入指數和 
   num /= 10;
   sum += uniE[quoti];
  }
  if(sum == i)
   printf("%d:  %llu\n", ++totalCnt, i); 
 }
 free(uniE); uniE = 0;
 return totalCnt;
}

從1位算到n位的Armstrong數:

int armstrong(int n)
{
 unsigned long long int sum, num, i, ubnd=9, lbnd=1;
 int count, quoti, digitCnt = 1, totalCnt = 0;
 unsigned long long int* uniE = 
   (unsigned long long int*)
            malloc(sizeof(unsigned long long int)*10);
 uniE[0] = 0;
 for(count = 1; count <= 9; count++)
 {
  uniE[count] = count;
  printf("%2d:  %d\n", ++totalCnt, count); // 1~9直接輸出
 }
 while(++digitCnt <= n)
 {
  for(count = 2; count <= 9; count++) // 指數單位更新
   uniE[count] *= count;
  
  lbnd *= 10;   // lower_bound為1 00.....0
  ubnd += 9*lbnd;  // upper_bound為9 99.....9
  
  for(i = lbnd; i <= ubnd; i++)
  {
   num = i; sum = 0;
   while(num != 0)
   {
    quoti = num % 10; // 每次取個位數加入指數和 
    num /= 10;
    sum += uniE[quoti];
   }
   if(sum == i)
    printf("%d:  %llu\n", ++totalCnt, i); 
  }
 }
 free(uniE); uniE = 0;
 return totalCnt;
}

[Algorithm] 數值字謎-延伸

試解以下的字謎：
1.
     SEND
+  MORE
-------------
 MONEY

2. (乘法)
         RED
x       FOR
---------------
 DANGER

3.
      LPBRR
      TLSRR
+    SBRRR
----------------
      RLMPL

4.
      EIN
      EIN
      EIN
+    EIN
-------------
    VIER

5. (額外條件：SOW可構成一完全平方數。)
      HOW
          I S
+    H I S
--------------
      S  I S

解：(如果需要詳解的話請留言告訴我，
         因為有些挺長的不好打XD)
1.
     9567
+   1085
-------------
   10652

2. 我知道可以推出D=1，其他大概只能硬用程式算= =

3.
    46388
    14288
+  23888
--------------
    84564

4.
    821
    821
    821
+  821
-------------
  3284

5. SOW=961=31^2
    461
      39
+  439
-----------
    939  

下面還有更扯的題目就不做了= =........腦袋會炸開Orz

[Algorithm] 數值字謎

請寫一個程式破解如下的字謎：
     VINGT
       CINQ
 +    CINQ
------------
  TRENTE

(每一個字母代表一個0~9之中的相異的數，且首位不為0)

解:
這種問題自己找到的線索越多，程式就能越簡易，
而以這題來說其實仔細思考就可以得到答案了XD。

1. 首先I+C+C就算有進位也就是1或2，到V那位以後最多就是是進位1，
    所以T=1，而V=8或9。此時R因為進位的關係要嘛是0，要嘛是1，
    但又不能重複，所以R=0。
2. 2I+N+(從前一位的進位)後，必然有進位，而前一位進位只能是0或2。
    (如果是1的話會造成I有小數點從而不符)
    如果前一位進位是0的話，I=5(∵2I+N=10+N)，
    從第一位來看，T+2Q=E，E和T一樣應該是奇數，
    但從I+C+C+1=E來看，E為偶數，矛盾。
    ∴前一位進位只能是2。
    ∴2I+N+2 = 10+N 或 2I+N+2=20+N
       I=4 或 I=9
3. 若I=9，則V=8。(From 1.) 現在式子長這樣：
    222        (進位)
-----------
    89NG1
      C9NQ
 +   C9NQ
------------
  10EN1E


∴9+2C+2=20+E  => 2C=9+E => C=6, E=3或C=7, E=5
如果C=6, E=3 => Q=1或6 矛盾，所以C=7，E=5。
因為E=5，所以Q=2(Q=7不服) 。
因此G+2N=21=>G為奇數，但1579都被用了，所以只剩G=3的可能。
3+2N=21 => N=9，矛盾。

所以I=9的假設是錯誤的！

4. T=1，R=0，I=4，V=8 or 9
         2     (進位)
-----------
    V4NG1
      C4NQ
 +   C4NQ
------------
  10EN1E



N+4+4+2=N+10，所以進1到下一位數。
2C+5=E+10或E+20，C等於6、7、8、9其中之一。
C=7的話 => E=9、V=9 (X)
C=8的話 => E=1 (X)
C=9的話 => E=3、V=8、Q=6，但1+G+2N=21，
                     G為偶數只能是2=>N=9 (X)
所以只剩下C=6的可能性。

5.
    112     (進位)
-----------
    V4NG1
       64NQ
 +    64NQ
------------
  10EN1E


明顯可以看出V=9，E=7。 => Q=3或8
Q=8的話=>1+G+2N=21，G為偶數只能是2=>N=9 (矛盾)
∴ Q=3 => G+2N=21，G為奇數只能是5=>N=8

6.全部推完的結果: R=0、T=1、Q=3、I=4、G=5、C=6、E=7、N=8。
    11200     (進位)
-----------
    94851
      6483
 +   6483
------------
  107817