ACM數論之旅11---淺談指數與對數（長篇）（今天休息，不學太難的數論＞ 3＜）

時間 2019-12-10

標籤 acm 數論之旅淺談指數對數長篇今天休息太難欄目應用數學简体版

原文原文鏈接

c/c++語言中，關於指數，對數的函數我也就知道那麼多php

exp(),pow(),sqrt(),log(),log10(),linux

exp(x)就是計算e的x次方，sqrt(x)就是對x開根號ios

pow()函數但是十分強大的(￣ε￣)c++

pow(a, b)能夠算a的b次方，可是b不限於整數，小數也能夠算法

因此pow(x, 0.5)至關於sqrt(x)ide

　　pow(M_E, x)至關於exp(x) (M_E就是e)函數

這是我在math.h發現的能夠直接用post

 1 #ifdef __STRICT_ANSI__
 2 #undef __STRICT_ANSI__
 3 #endif
 4 #include<cstdio>
 5 #include<cmath>
 6 int main(){
 7     printf("%.20f\n", M_E);
 8     printf("%.20f\n", M_PI);
 9 }
10 /*
11 在math頭文件的前面，最好加上最上面那三行
12 由於編譯器不一樣，系統不一樣，都有可能致使用不了M_E
13 好比codeforces，不加前三行，沒法識別M_E
14 */

可是函數總有他存在的意義，要否則你們都用pow(x, 0.5)，沒人用sqrt了測試

因此我認爲，後者比前者要快，或者可能更精確(oﾟωﾟo)優化

好比sqrt，我來給你們講一個鬼故事(ΘˍΘ=)，有一個比sqrt還要快計算出根號的函數

一個關於被稱做「魔數」0x5F3759DF的故事

如下摘自http://www.guokr.com/post/90718/

（不看能夠跳過）

http://www.douban.com/note/93460299/

Quake-III Arena (雷神之錘3)是90年代的經典遊戲之一。該系列的遊戲不但畫面和內容不錯，並且即便計算機配置低，也能極其流暢地運行。這要歸功於它3D引擎的開發者約翰-卡馬克（John Carmack）。事實上早在90年代初DOS時代，只要能在PC上搞個小動畫都能讓人驚歎一番的時候，John Carmack就推出了石破天驚的Castle Wolfstein, 而後再接再勵，doom, doomII, Quake...每次都把3-D技術推到極致。他的3D引擎代碼資極度高效，幾乎是在壓榨PC機的每條運算指令。當初MS的Direct3D也得聽取他的意見，修改了很多API。

最近，QUAKE的開發商ID SOFTWARE 遵照GPL協議，公開了QUAKE-III的原代碼，讓世人有幸目擊Carmack傳奇的3D引擎的原碼。
這是QUAKE-III原代碼的下載地址：
http://www.fileshack.com/file.x?fid=7547

(下面是官方的下載網址，搜索「quake3-1.32b-source.zip」能夠找到一大堆中文網頁的
ftp://ftp.idsoftware.com/idstuff/source/quake3-1.32b-source.zip)

咱們知道，越底層的函數，調用越頻繁。3D引擎歸根到底仍是數學運算。那麼找到最底層的數學運算函數（在game/code/q_math.c），必然是精心編寫的。裏面有不少有趣的函數，不少都使人驚奇，估計咱們幾年時間都學不完。

在game/code/q_math.c裏發現了這樣一段代碼。它的做用是將一個數開平方並取倒，經測試這段代碼比(float)(1.0/sqrt(x))快4倍：
float Q_rsqrt( float number )
{
long i;
float x2, y;
const float threehalfs = 1.5F;

x2 = number * 0.5F;
y = number;
i = * ( long * ) &y; // evil floating point bit level hacking
i = 0x5f3759df - ( i >> 1 ); // what the fuck?
y = * ( float * ) &i;
y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
// y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed

#ifndef Q3_VM
#ifdef __linux__
assert( !isnan(y) ); // bk010122 - FPE?
#endif
#endif
return y;
}

函數返回1/sqrt(x)，這個函數在圖像處理中比sqrt(x)更有用。
注意到這個函數只用了一次疊代！（其實就是根本沒用疊代，直接運算）。編譯，實驗，這個函數不只工做的很好，並且比標準的sqrt()函數快4倍！要知道，編譯器自帶的函數，但是通過嚴格仔細的彙編優化的啊！

這個簡潔的函數，最核心，也是最讓人費解的，就是標註了「what the fuck?」的一句
i = 0x5f3759df - ( i >> 1 );

再加上y = y * ( threehalfs - ( x2 * y * y ) );
兩句話就完成了開方運算！並且注意到，核心那句是定點移位運算，速度極快！特別在不少沒有乘法指令的RISC結構CPU上，這樣作是極其高效的。

算法的原理其實不復雜,就是牛頓迭代法,用x-f(x)/f'(x)來不斷的逼近f(x)=a的根。

簡單來講好比求平方根,f(x)=x^2=a ,f'(x)= 2*x,f(x)/f'(x)=x/2,把f(x)代入

x-f(x)/f'(x)後有(x+a/x)/2，如今咱們選a=5,選一個猜想值好比2，
那麼咱們能夠這麼算
5/2 = 2.5; (2.5+2)/2 = 2.25; 5/2.25 = xxx; (2.25+xxx)/2 = xxxx ...
這樣反覆迭代下去，結果一定收斂於sqrt(5)，沒錯，通常的求平方根都是這麼算的
可是卡馬克(quake3做者)真正牛B的地方是他選擇了一個神祕的常數0x5f3759df 來計算那個猜想值
就是咱們加註釋的那一行,那一行算出的值很是接近1/sqrt(n),這樣咱們只須要2次牛頓迭代就能夠達到咱們所須要的精度.
好吧若是這個還不算NB,接着看:

普渡大學的數學家Chris Lomont看了之後以爲有趣，決定要研究一下卡馬克弄出來的
這個猜想值有什麼奧祕。Lomont也是個牛人，在精心研究以後從理論上也推導出一個
最佳猜想值，和卡馬克的數字很是接近, 0x5f37642f。卡馬克真牛，他是外星人嗎？

傳奇並無在這裏結束。Lomont計算出結果之後很是滿意，因而拿本身計算出的起始
值和卡馬克的神祕數字作比賽，看看誰的數字可以更快更精確的求得平方根。結果是
卡馬克贏了... 誰也不知道卡馬克是怎麼找到這個數字的。

最後Lomont怒了，採用暴力方法一個數字一個數字試過來，終於找到一個比卡馬克數
字要好上那麼一丁點的數字，雖然實際上這兩個數字所產生的結果很是近似，這個暴
力得出的數字是0x5f375a86。

Lomont爲此寫下一篇論文，"Fast Inverse Square Root"。

論文下載地址：
http://www.math.purdue.edu/~clomont/Math/Papers/2003/InvSqrt.pdf
http://www.matrix67.com/data/InvSqrt.pdf

參考：<IEEE Standard 754 for Binary Floating-Point Arithmetic><FAST INVERSE SQUARE ROOT>

最後，給出最精簡的1/sqrt()函數：
float InvSqrt(float x)
{
float xhalf = 0.5f*x;
int i = *(int*)&x; // get bits for floating VALUE
i = 0x5f375a86- (i>>1); // gives initial guess y0
x = *(float*)&i; // convert bits BACK to float
x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
return x;
}
你們能夠嘗試在PC機、5一、AVR、430、ARM、上面編譯並實驗，驚訝一下它的工做效率。

百度百科給的Carmack的sqrt()函數

static float CarmackSqrt ( float x)

{

float xhalf = 0.5f * x;

int i = *( int *)&x; // get bits for floating VALUE

i = 0x5f3759df - (i>>1); // gives initial guess y0

x = *( float *)&i; // convert bits BACK to float

x = x*(1.5f - xhalf*x*x); // Newton step, repeating increases accuracy

return (1 / x);

}

看完故事是否是以爲技巧很重要

故事告一段落，要不要來練練手(￣ｏ￣)

poj 2109

http://poj.org/problem?id=2109

題目大意： K ^ N = P, 給N 和 P，求K。數據規模：1<=n<= 200, 1<=p<10¹⁰¹ 並且保證存在 k, 1<=k<=10⁹ 。

正常不就是二分+高精度算法嗎？

AC代碼：

1 #include<cstdio>
2 #include<cmath>
3 double n, p;
4 int main(){
5     while(scanf("%lf%lf", &n, &p) != EOF){
6         printf("%.0f\n", pow(p, 1/n));
7     }
8 }

View Code

哇哈哈，看到沒有，看到沒有，這就是技巧(*ﾟ▽ﾟ*)

double雖然精度只有16位左右，可是咱們只要知道前16位就夠了，後面任憑他用科學計數法去表示吧，反正咱們不須要。

由於當n錯一位，K的值前16位都會變化很大，因此這樣計算出來的K確定是惟一的。

下面來講說對數：

C語言中，有兩個log函數，分別爲log和log10函數

log()是以e爲底數的，數學上咱們是寫做ln(x)的

log10()是以10爲底數的

那若是我想以2做爲底數怎麼辦

這麼寫 log(x) / log(2) 數學公式，還記得嗎<（￣︶￣）>

定義類型：double log(double x);

　　　　　double log10(double x);

固然咱們通常用double的，它不僅能接受double

double log (double x);
      float log (float x);
long double log (long double x);
     double log (T x);           // additional overloads for integral types

最後一句模板T類型只有c++11支持，基本你不會本身去重載因此用不上
而後，從c++98開始，就支持  and <complex><valarray>兩個類型了



待會我會講講<complex>頭文件，這是複數類

在比較a^b和c^d次方，若是b和d很是大怎麼辦

好比這題：hdu 5170

http://acm.hdu.edu.cn/showproblem.php?pid=5170

告訴你a,b,c,d，要你比較a^b和c^d，輸出"<"，">"或"="

1≤a,b,c,d≤1000

因此直接用log的性質

log(a^b) = b * log(a)

若是兩邊同時log一下再比較，那就方便多了（注意log有精度偏差）

完整性質：

AC代碼：

 1 #include<cstdio>
 2 #include<cmath>
 3 int main(){
 4     int a, b, c, d;
 5     double l, r;
 6     while(~scanf("%d%d%d%d", &a, &b, &c, &d)){
 7         l = b * log(a);
 8         r = d * log(c);
 9         if(fabs(l - r) < 1e-6){//精度偏差，通常小於0.000001能夠認爲相等 
10             puts("=");
11         }else if(l < r){
12             puts("<");
13         }else{
14             puts(">");
15         }
16     }
17 }
18 //關於1e-6,有人寫1e-7,1e-8,連1e-10都有,看喜愛咯

View Code

或許有比較數學化的比較方法，可是精度用的好的人真是無敵了(☆ﾟ∀ﾟ)

有沒有想過遇到x^y^z怎麼辦

cf 621D

http://codeforces.com/problemset/problem/621/D

給你三個數x,y,z，比較這12個式子，問你哪一個式子最大

0.1 ≤ x, y, z ≤ 200.0

x^(y^z)

這個式子log一下

變成

原式 = y^z*log(x)

再log一下變成

= log(y^z*log(x))

= log(y^z) + log(log(x))

= z * log(y) + log(log(x))

原本這樣就能夠比較了

但是題目的範圍是0.1

log()小數會產生負數

log負數就沒意義了

因此對於log(log(x))這麼寫不行

那怎麼辦

哼哼，技巧

double範圍 -1.7*10(-308)～1.7*10(308)
long double範圍 128 18-19 -1.2*10(-4932)～1.2*10(4932)

雖然他們兩精度都是16位，可是200的200次方long double居然存的下

因此只要一次log就行了

而後愉快的寫代碼吧

AC代碼：

 1 #include<cstdio>
 2 #include<cmath>
 3 #include<iostream>
 4 using namespace std;
 5 char str[12][10] = {
 6     "x^y^z",
 7     "x^z^y",
 8     "(x^y)^z",
 9     "(x^z)^y",
10     "y^x^z",
11     "y^z^x",
12     "(y^x)^z",
13     "(y^z)^x",
14     "z^x^y",
15     "z^y^x",
16     "(z^x)^y",
17     "(z^y)^x",
18 };
19 long double x, y, z;
20 long double mx, t;
21 int pos;
22 void judge(int x){
23     //printf("t = %llf\n", t);
24     if(fabs(mx - t) <= 1e-6) return ;
25     else if(mx < t){
26         pos = x;
27         mx = t;
28     }
29 }
30 int main(){
31     cin >> x >> y >> z;
32     pos = 0;
33     mx = pow(y, z)*log(x);
34     t = pow(z, y)*log(x);
35     judge(1);
36     t = z*log(pow(x, y));
37     judge(2);
38     t = y*log(pow(x, z));
39     judge(3);
40     
41     t = pow(x, z)*log(y);
42     judge(4);
43     t = pow(z, x)*log(y);
44     judge(5);
45     t = z*log(pow(y, x));
46     judge(6);
47     t = x*log(pow(y, z));
48     judge(7);
49     
50     t = pow(x, y)*log(z);
51     judge(8);
52     t = pow(y, x)*log(z);
53     judge(9);
54     t = y*log(pow(z, x));
55     judge(10);
56     t = x*log(pow(z, y));
57     judge(11);
58     
59     printf("%s\n", str[pos]);
60 }

View Code

其實log()一個負數是能夠解的

還記得當年大明湖畔的歐拉公式嗎

e^iπ = -1

由於e的i∏次方等於-1

因此log(-1) = i∏

因此負數迎刃而解

log(-2) = log(-1 * 2) = log(-1) + log(2)

那log(i)呢

根號-1等於i

因此log(i) = log( -1^(1/2) ) = 1/2 * log(-1) = 1/2 * i∏

那log(a + bi)

歐拉原公式寫做

e^ix= cosx + isinx

那麼

因此說嘛，年輕人就應該拿一本複變函數去看去(,,• ₃ •,,)

附上剛剛那題用複數計算的AC代碼

 1 #include <iostream>
 2 #include <complex>
 3 #include <string>
 4 using namespace std;
 5 bool bigger (complex<long double> a, complex<long double> b) {
 6     if (imag(a) == 0 && imag(b) == 0) {//沒有虛部 
 7         return real(a) > real(b);//比較實部 
 8     } else if (imag(a) == 0 && imag(b) != 0) { //有虛部的確定小 
 9         return true;
10     } else if (imag(a) != 0 && imag(b) == 0) {
11         return false;
12     } else if (imag(a) != 0 && imag(b) != 0) {//都有虛部，按實部反過來比 
13         return real(a) < real(b);
14     }
15 }
16 
17 int main () {
18     long double ax, ay, az;
19     cin >> ax >> ay >> az;
20     
21     complex<long double> x (ax, 0.0L);
22     complex<long double> y (ay, 0.0L);
23     complex<long double> z (az, 0.0L);
24     
25     complex<long double> cmaz (3, 3);
26     string ans = "xd";
27     
28     if (bigger(z * log(y) + log(log(x)), cmaz)) {
29         cmaz = z * log(y) + log(log(x));
30         ans = "x^y^z";
31     }
32     if (bigger(y * log(z) + log(log(x)), cmaz)) {
33         cmaz = y * log(z) + log(log(x));
34         ans = "x^z^y";
35     }
36     if (bigger(log(y * z) + log(log(x)), cmaz)) {
37         cmaz = log(y * z) + log(log(x));
38         ans = "(x^y)^z";
39     }
40     
41     if (bigger(z * log(x) + log(log(y)), cmaz)) {
42         cmaz = z * log(x) + log(log(y));
43         ans = "y^x^z";
44     }
45     if (bigger(x * log(z) + log(log(y)), cmaz)) {
46         cmaz = x * log(z) + log(log(y));
47         ans = "y^z^x";
48     }
49     if (bigger(log(x * z) + log(log(y)), cmaz)) {
50         cmaz = log(x * z) + log(log(y));
51         ans = "(y^x)^z";
52     }
53     
54     if (bigger(y * log(x) + log(log(z)), cmaz)) {
55         cmaz = y * log(x) + log(log(z));
56         ans = "z^x^y";
57     }
58     if (bigger(x * log(y) + log(log(z)), cmaz)) {
59         cmaz = x * log(y) + log(log(z));
60         ans = "z^y^x";
61     }
62     if (bigger(log(x * y) + log(log(z)), cmaz)) {
63         cmaz = log(x * y) + log(log(z));
64         ans = "(z^x)^y";
65     }
66     
67     cout << ans << endl;
68 }