#!/usr/bin/env python
# coding: utf-8

# # [“กฎของ 72”](https://en.wikipedia.org/wiki/Rule_of_72)
# 
# เมื่อเราเอาเงินไปลงทุนที่อัตราตอบแทนทบต้น $r$ เปอร์เซ็นต์ต่อปี “กฎของ 72” บอกว่าจะใช้เวลา $\frac{72}{r}$ ปีที่จะทำให้เงินงอกเงยเป็นสองเท่า เช่นถ้าผลตอบแทนเท่ากับ 3% ต่อปีจะใช้เวลา $\frac{72}{3} = 24$ ปี ถ้าผลตอบแทนเท่ากับ 8% ต่อปีจะใช้เวลา$\frac{72}{8} = 9$ ปี
# 
# (กฎนี้ไม่ใช่กฎที่ไอน์สไตน์คิดนะครับ น่าจะมีมาก่อนไอน์สไตน์หลายร้อยปีแล้ว แต่ไอน์สไตน์เคยพูดเรื่องดอกเบี้ยทบต้นไว้ว่า “Compound interest is the eighth wonder of the world. He who understands it earns it … he who doesn’t … pays it.” เลยมีคนเข้าใจผิดว่าไอน์สไตน์เป็นคนคิด)
# 
# กฎนี้เป็นการประมาณเท่านั้น แต่ประมาณได้ใกล้เคียงกับความเป็นจริงมากเมื่อผลตอบแทนไม่สูงเกินไป (เช่นต่ำกว่า 30% ต่อปีลงมา) บางครั้งเราจะใช้เลข 70 แทน 72 ก็ได้ ขึ้นอยู่กับว่าหาร 70 หรือ 72 ง่ายกว่า
# 
# สำหรับเด็กๆที่รู้จัก logarithm เราสามารถคำนวณเวลาที่ทำให้เงินงอกเป็นสองเท่าได้เป๊ะๆเท่ากับ $\frac{log(2)}{log(1+\frac{r}{100})}$ 
# 
# ถ้าเด็กรู้จัก logarithm และ calculus มีคำอธิบายอ่านได้ที่  [กฎของ "70" (หรือ "72")](http://mpec.sc.mahidol.ac.th/forums/index.php/topic,1807.0.html) ครับ เป็นผลลัพธ์จากการกระจาย Taylor's series ของ $\frac{log(2)}{log(1+\frac{r}{100})}$ 
# 
# แต่ถ้าไม่รู้จักคณิตศาสตร์เหล่านั้น เราสามารถใช้วิธี bisection ที่เราเคยใช้ในอดีตมาคำนวณเวลาที่ทำให้เงินเพิ่มเป็นสองเท่าก็ได้
# 
# ก่อนอื่นทดลองเปรียบเทียบเวลาที่เงินกลายเป็นสองเท่าที่คำนวณจากสูตร $\frac{log(2)}{log(1+\frac{r}{100})}$ กับเวลาที่คำนวณจากกฎของ 72 และ 70 ดู จะพบว่าถ้าผลตอบแทนน้อยกว่า 30% ลงมา คำตอบประมาณจากกฎของ 72 และ 70 จะผิดเพี้ยนน้อย:
# 
# 
# 

# In[1]:


from math import log  #ใช้ฟังก์ชั่น math.log ในการคำนวณ logarithm

print("%\tคำตอบ\tกฎ 72\tกฎ 70")
for r in (1, 2, 5, 8, 10, 15, 20, 25, 30, 50, 100):
    ans = log(2)/(log(1+r/100)) # คำตอบเป๊ะๆ
    ans_72 = 72/r               # คำตอบประมาณจาก 72/r
    ans_70 = 70/r               # คำตอบประมาณจาก 70/r
    print(f"{r}\t{ans:.1f}\t{ans_72:.1f}\t{ans_70:.1f}") #ใช้ f-string แสดงทศนิยม 1 ตำแหน่ง


# สมมุติว่าเราไม่รู้จัก logarithm และ calculus เราจะใช้วิธี bisection มาหาคำตอบว่าใช้เวลากี่ปีเงินถึงจะเพิ่มเป็นสองเท่าได้ ก่อนอื่นเราเอาฟังก์ชั่น bisection4 มาจากการเรียนในอดีตที่ [วิทย์ม.ต้น: เขียนโปรแกรมไพธอนหาว่าพาราโบลาตัดแกน x ที่ไหน, วิธีหาคำตอบโดย Bisection Method](https://witpoko.com/?p=6595) 

# In[2]:


# bisection เวอร์ชั่น 4 ฟังก์ชั่นไม่เรียกตัวเองแล้ว (non-recursive)
# ไม่พิมพ์ดูว่า xmin, xmax, f(xmin)*f(xmax) เป็นเท่าไรแล้ว
# คำตอบที่ได้จะเทียบบัญญัติไตรยางค์ระหว่างจุด (xmin, f(xmin)) และ จุด (xmax, f(xmax))
# ซื่งน่าจะได้คำตอบที่ใกล้ความจริงมากกว่าแบบตอบตรงกลางระหว่าง xmin และ xmax

def bisection4(f,xmin, xmax):
    "พยายามหาค่า x ที่ทำให้ f(x)==0, xmin และ xmax คือช่วงที่เดาว่าคำตอบอยู่ในนั้น"
    
    tolerance = 1e-6  #ตั้งค่า tolerance ไว้ให้หยุดทำงาน
                      #จะหยุดทำงานเมื่อ xmin ห่างจาก xmax น้อยกว่า tolerance
    
    if xmin > xmax:   # จัดการให้ xmin น้อยกว่า xmax เสมอ
        xmin, xmax = xmax, xmin
    
    if f(xmin)*f(xmax) > 0:
        print("Bad xmin, xmax") #ถ้าฟังก์ชั่น f ไม่เปลี่ยนเครื่องหมายจากบวกไปลบหรือลบไปบวก
                              #ระหว่าง xmin และ xmax ก็ควรไปเดาใหม่ว่า xmin, xmax คืออะไร
        return None
    
    
    while xmax-xmin > tolerance: # ทำตรงนี้วนๆไปตราบใดท่ี xmin และ xmax ยังห่างกันกว่า tolerance
        
        #print(f"xmin, xmax = {xmin}, {xmax}\t f(xmin)*f(xmax) = {f(xmin)*f(xmax):.5f}")
        
        if f(xmin) == 0:
            return xmin
        
        if f(xmax) == 0:
            return xmax

        xmid = (xmin + xmax)/2    #หาจุดกลางระหว่าง xmin และ xmax
        if f(xmin)*f(xmid) < 0:   #ถ้าฟังก์ชั่นเปลี่ยนเครื่องหมายระหว่าง xmin กับจุดกลาง ก็ไปหาต่อในช่วงนี้
            xmin, xmax = xmin, xmid
        else:                     #ไม่งั้นก็หาต่อในช่วงจุดกลางถึง xmax
            xmin, xmax = xmid, xmax
            
    slope = (f(xmax)-f(xmin))/(xmax-xmin) 
    x = xmin + (0-f(xmin))/slope  # คำตอบที่ได้จะเทียบบัญญัติไตรยางค์ระหว่างจุด (xmin, f(xmin)) และ จุด (xmax, f(xmax))
    
    return x


# In[3]:


#เราสามารถตั้งชื่อฟังก์ชั่น bisection4 ใหม่ให้เป็นชื่อที่ง่ายขึ้นก็ได้ดังนี้
solve = bisection4


# In[4]:


# ทดลองหาเวลาที่ทำให้การลงทุนทบต้นที่ผลตอบแทน 10% ต่อปีเติบโตเป็นสองเท่าด้วยวิธี bisection

# กำหนดฟังก์ชั่นที่รับตัวแปร y คิดเป็นปีเข้าไป 
# แล้วจะให้ค่าเท่ากับ 0 เมื่อ (1+10/100) ** y เท่ากับ 2
def f10(y):
    return (1+10/100)**y - 2 

ymin = 0  #เดาว่าคำตอบอยู่ระหว่าง 0-10 ปี
ymax = 10

#แก้สมการด้วย bisection (ซึ่งเราตั้งชื่อใหม่ว่า solve แล้วจากข้างบน)
ans_bisect = solve(f10, ymin, ymax)

#คำตอบเป๊ะจากสูตร log(2)/log(1+10/100)
ans = log(2)/log(1+10/100)

#เปรียบเทียบกัน
print(f"คำตอบจาก bisection:\t{ans_bisect}")
print(f"คำตอบจากสูตร log():\t{ans}")


# ลองเปรียบเทียบคำตอบจากกฎ 72, 70, และคำตอบที่หาจากวิธี bisection ดู
# 
# จะพบว่าถ้าเรารู้จักใช้ bisection เราก็สามารถหาคำตอบได้ถูกต้องเหมือนๆสูตรจาก log 

# In[5]:


print("%\tคำตอบ\tกฎ 72\tกฎ 70\tbisect")
for r in (1, 2, 5, 8, 10, 15, 20, 25, 30, 50, 100):
    ans = log(2)/(log(1+r/100)) # คำตอบเป๊ะๆ
    ans_72 = 72/r               # คำตอบประมาณจาก 72/r
    ans_70 = 70/r               # คำตอบประมาณจาก 70/r
    def f(y):                   # ฟังก์ชั่นนี้รับตัวแปร y เป็นปีเข้าไป
        return (1+r/100)**y - 2 # จะให้ค่าเป็น 0 เมื่อทบต้น y ปีแล้วเงินเพิ่มเป็น 2 เท่า
    ans_bisect = solve(f,0,100) # คำตอบจากวิธี bisection
    print(f"{r}\t{ans:.1f}\t{ans_72:.1f}\t{ans_70:.1f}\t{ans_bisect:.1f}") #ใช้ f-string แสดงทศนิยม 1 ตำแหน่ง


# # ตัวอย่างการประยุกต์หาคำตอบด้วย bisection เพิ่มเติม
# 
# สมมุติว่าเราต้องการทราบว่าเราต้องออมเงินปีละเท่าไรเป็นเวลา 40 ปี ด้วยผลตอบแทนทบต้นปีละ 8% ถึงจะมีเงินเมื่อสิ้นปีที่ 40 เท่ากับ 40,000,000 บาท ถ้าเราไม่รู้สูตรอะไรเลยเราก็สามารถใช้ bisection หาคำตอบให้ได้เหมือนกัน
# 
# ก่อนอื่นเราเอาฟังก์ชั่นคำนวณเงินอนาคตจาก [วิทย์ม.ต้น: การเติบโตของเงินออมเป็นประจำ, หัดใช้ MU-EDITOR](https://witpoko.com/?p=6673) มาใช้เพื่อความสะดวกก่อน:

# In[6]:


def future_value(present_value, interest, periods):
    """คำนวณมูลค่าเงินในอนาคต (future value) 
    ถ้าเรารู้ค่าเงินปัจจุบัน (present value), 
    อัตราผลตอบแทนต่อปีหรือต่อเดือน (interest), 
    และจำนวนปีหรือเดือน (periods) ที่ออมเงินไว้ให้งอกเงย
    """
    return present_value * (1+interest) ** periods

def total_future_value(deposit, interest, n_periods):
    """คำนวณค่าเงินในอนาคตทั้งหมด ถ้าออมเงินเท่าๆกันทุกปี = deposit
    ด้วยผลตอบแทนต่อปี = interest (เช่น 5% = 0.05)
    เป็นจำนวนปีทั้งหมด = n_periods
    """
    sum = 0.0
    for i in range(1,n_periods+1):
        sum = sum + future_value(deposit, interest, i)
    return sum


# แล้วเราก็ทำการคำนวณคำตอบแบบนี้ก็ได้:

# In[7]:


# สร้างฟังก์ชั่นที่รับตัวแปร deposit (คือเงินฝากต่อปี เป็นตัวแปรที่เราต้องการหาคำตอบ) 
# ฟังก์ชั่นนี้ให้ค่าเท่ากับ 0 ถ้าเงินฝากปีละ deposit 
# ด้วยผลตอบแทนทบต้น 8% ต่อปีเป็นเวลา 40ปี แล้วมีเงินทั้งหมด 40 ล้าน
def fd(deposit):  
    return total_future_value(deposit, 0.08, 40) - 40_000_000

# แก้ด้วยวิธี bisection เดาว่าคำตอบอยู่ระหว่าง 100,000-1,000,000 
ans = solve(fd, 100_000, 1_000_000)
print(f"ต้องออมเงินปีละ {ans:.0f} บาท ที่ผลตอบแทนทบต้น 8% เป็นเวลา 40 ปีถึงจะมีเงิน 40 ล้านบาท")


# สมมุติว่าเราต้องการออมเงินทุกปีเป็นเวลา 20 ปี ปีละ 100,000 บาท แล้วเราต้องการให้เงินตอนจบเป็น 10 ล้านบาท เราต้องหาผลตอบแทนทบต้นปีละกี่เปอร์เซ็นต์ เราก็สามารถคำนวณดังนี้ได้:

# In[8]:


# สร้างฟังก์ชั่นที่รับตัวแปร r (อัตราผลตอบแทนต่อปี เป็นตัวแปรที่เราต้องการหาคำตอบ)
# ฟังก์ชั่นนี้จะให้ค่าเท่ากับ 0 ถ้าฝากปีละ 100,000 บาท
# ทบต้นที่อัตราผลตอบแทน r ต่อปีเป็นเวลา 20 ปี แล้วมีเงินทั้งหมด 10 ล้าน
def fr(r):
    return total_future_value(100_000, r, 20) - 10_000_000

# แก้ด้วยวิธี bisection เดาว่าผลตอบแทนอยู่ระหว่าง 0% และ 20%
ans = solve(fr,0,0.2)
print(f"ต้องหาผลตอบแทนทบต้น {ans*100:.1f}% ต่อปี ฝากปีละ 100,000 เป็นเวลา 20 ปี ถึงจะมีเงิน 10 ล้านบาท")


# In[ ]:





# In[ ]:





# In[ ]: