Забавна математика с Python

Спомените ми от първите години в Техническия университет (освен разбиващите купони и хроничното студентско безпаричие) са свързани с писането на огромно количество протоколи от измервания и сума ти свързани с това изчисления. Понякога се изкушавах да си напиша кратка програма за да си сметна и начертая графиките от измерените резултати вместо налагания от някои преподаватели методи с калкулатор и милиметрова хартия. Проблемът беше, че сметките често не бяха елементарни, а Паскал (и по-късно Java) не са от най-лесните за ползване езици за бързи сметки от мързеливи студенти като мен. По-късно се налагаше да правим преобрзувания на уравнения и да ги опростяваме, а единствения софтуер който донякъде помагаше беше MathLab. Въоръжени с него и с теоремите на Нагласаяев и Натъмъняев успявахме да открием нови клонове в инженерната математика. Наистина се изискваше много усилия от предварително дадения ни за сравнение резултат да приложим реверсивен инженериг и да скалъпим обратно оригиналното задание на проекта, но историята познава и по-големи героични постъпки (като например свързването на амперметър като волтметър или паралелното свързванео на електролитни кондензатори към променливотокови трансформаторни вериги).

От тези времена ми остана един респект към многостъпковите изчисления (и към електролитните кондензатори) и често след това съм се чудил дали няма наистина удобен инструмент точно за такива уморителни, но необходими математически гимнастики. Изискванията ми за подобен софтуер са скромни – да е софтуер с отворен код, преносим (Linux/Windows), да е простичък за инсталация и употреба и по възможност да се разширява лесно. По едно време се бях отказал да търся (а и не ми трябваше, честно казано) докато преди месец попаднах на SymPy и mpmath. Първата е библиотека занимаваща се със символна алгебра (от типа колко е (ax2 + by3)2 * (x + y2)2 в разгъната форма) изчисляваща също интеграли и диференциали от символни уравнения (такива с неизвестни като x, y и z например). Втората библиотека и за смятане на реални и комплексни числа с произволна точност. Точно така – с произволна. И двете библиотеки са писани на Python които освен че е страшно лесен за учене и експериментиране е също и много мощен скриптов език. И двете библиотеки са с отворен код и понеже са на Python са достъпни на всички софтуерни платформи поддържани от езика. Ала един пример говори повече от сто реклами, така че ето един пример:

>>>from sympy import *
>>>x = Symbol('x')
>>>y = Symbol('y')
>>>a = ((x**2 + y**3)**2 * (x+y**2)**2)
>>>b.expand()
x**6 + y**10 + x**2*y**6 + x**4*y**4 + 2*x*y**8 + 2*x**2*y**7 + 2*x**4*y**3 + 2*x**5*y**2 + 4*x**3*y**5
>>>

Трите символа “>” са от промпт-а на Python интепретатора. Последния ред е всъщност разгънатата форма на по-горното уравнение. Звздичката е знак за умножение а двойната звезда – повдигане на степен. За домашно – сметнете уравнението на степен 3 🙂
Това далеч не е всичко. Ето друг пример (взет от ръководствотото на SymPy):

>>> from sympy import *
>>> x=Symbol("x")
>>> limit(sin(x)/x, x, 0)
1
>>> limit(x, x, oo)
oo # това е знака за безкрайност :)
>>> limit((5**x+3**x)**(1/x), x, oo)
5

Някой спомня ли си как се смятаха границите (лимеси) от училище? Е, с няколко реда код няма да ви се налага да ги смятате.

Впечатлих ли ви? А ето какво може да прави mpmath:

>>> from mpmath import *
>>> pi
mpf('3.1415926535897932384626433832793')

Хм, дотук – нищо впечатляващо. Ала нека да променим прецизността след десетичната точка (която по подразбиране е 30 знака):

>>> from mpmath import *
>>> mpf.dps = 100 # задаваме броя на знаците след десетичната точка
>>> pi
mpf('3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798')

Това е числото Пи със сто знака след десетичната точка (или запетая – кое е по-правилното?). А ето го със 1000 знака:

3.1415926535897932384626433832795028841971693993751058209749445923078164062862
089986280348253421170679821480865132823066470938446095505822317253594081284811
174502841027019385211055596446229489549303819644288109756659334461284756482337
867831652712019091456485669234603486104543266482133936072602491412737245870066
063155881748815209209628292540917153643678925903600113305305488204665213841469
519415116094330572703657595919530921861173819326117931051185480744623799627495
673518857527248912279381830119491298336733624406566430860213949463952247371907
021798609437027705392171762931767523846748184676694051320005681271452635608277
857713427577896091736371787214684409012249534301465495853710507922796892589235
420199561121290219608640344181598136297747713099605187072113499999983729780499
510597317328160963185950244594553469083026425223082533446850352619311881710100
031378387528865875332083814206171776691473035982534904287554687311595628638823
53787593751957781857780532171226806613001927876611195909216420199

За домашна – пробвайте със 10000 знака 🙂 Повече примери и идеи за ползване на mpmath вижте нейната документация.

Забележка: За момента двете библиотеки дефинират числото Пи по свой начин, който не е съвместим, затова рестартирайте комания интерпретатор на Python между тестовете за да получите смислени резултати.

Стана ли ви интересно? Аз със сигурност съм заинтригуван! Сега се надявам по-малко ученици и студенти да четат това, че иначе много домашни ще станат безумно лесни за решаване с няколко редова програмка 😀

13 thoughts on “Забавна математика с Python

  1. На мен ми беше препоръчана Максима като компютърна алгебра система. Написана е на Lisp, разпространява се под GPL и не е много трудна за ползване.

    Ето и отговора на домашната:
    (%i6) block([fpprec:10000],bfloat(%pi));

    Програмата ми беше препоръчана от един от професорите ми и мога да кажа, че е доста добра.

    Ето за теб една домашна:
    коe е 10000 фибоначи число? 🙂

  2. Признавам, отне ми малко време да си припомня за Фибоначи. Иначе сметнах резултата. Числото е със 2090 цифри като първите и последните десет са: 2079360823…1238230626. Така става като смятам на ръка 🙂

    Иначе ето и програмката на Python за смятане на произволно число от редицата на Фибоначи (тъй като Python стандартно си работи с неограничено големи ЦЕЛИ числа, даже нямаше нужда да зареждам от специализирана библиотека):


    def fibonacci(number):
        a,  b = 0, 1
        counter = 1 # -- WRONG: must be 0. (Thanks to George) See bellow blog comments. 
        while counter < number:
            a,  b = b,  a + b
            counter = counter + 1
        return a
    
    fiboNumber = 10000
    f = fibonacci(fiboNumber)
    
    # -- remove the comment on the next line if you want to print off the whole number
    #print ("Fibonacci number ",  fiboNumber,  ": ",  f)
    
    # -- The following code just print the first and last 10 digits
    fst = str (f)
    st = fst[:10] + "..." + fst[-10:]
    print ("Formated fibonacci number ", fiboNumber, ": ", st)
    print ("Digits in the Fibonacci number: ",  len(fst))

    Formated fibonacci number 10000 : 2079360823…1238230626
    Digits in the Fibonacci number: 2090


    Ако някой има желание да отпечата цялото число, трябва само да премахне коментара от първия print().

    А сега още една домашна за такива прилежни студенти като теб. Кажи ми кое е 10000-ното просто число? 😎

  3. Относно други програми за алгебрични и символни изчисления – на страницата на SymPy има линкове към поне няколко такива. Проблемът е обаче че имам непоносимост към всякакви компютърни (и не само) езици, чиито програми изглеждат като че ли котка се е разхождала по клавиатурата. Затова може би Python ми допадна – при него празното пространство ИМА значение и програмата ТРЯБВА да е структурирана и визуално за да работи правилно. За сравнение – LISP има твърде много скоби за вкуса ми, а Perl използва твърде много специални знаци (не една, ами цяло племе котки са играли ритуални танци по клавиатурата). И за финал – пробвай следната команда под Python:

    >>> import this

    Резултатът е много… Дзен.

  4. Това изглежда добре, но е малко много код за изчисляване на фибоначи. Под максима има вградена функция fib(), но ако искаш сам да си я дефинираш можеш да ползваш кода:
    fib[0]:0; fib[1]:1; fib[n]:=fib[n-1]+fib[n-2];
    и после да извикаш примерно
    fib[150];

    Колкото до самият отговор, трябва до уточним, че това е числото при n=999 за функцията Фибоначи, но е в действителност 10000-ното число ( 0 е първото).

    За твоят въпрос трябваше да си напиша моя функция на максима:

    s:1; for i:1 while i<=10000 do s:next_prime(s); done; s;

    За да ми даде резултат:
    104729

    Предполагам е мой ред да поискам 10000-ното просто Фибоначи, но това ще ти отнеме поне няколко минути 🙂

    Поздрави

  5. Жоро, благодаря за забележката. Наистина първото число е нула (с индекс нула) а не едно както съм написал. Ето и коригирания алгоритъм:


    def fibonacci(number):
        a,  b = 0, 1
        counter = 0 # corrected!!!
        while counter < number:
            (a,  b) = (b,  a + b)
            counter = counter + 1 
        return a
    

    Formated fibonacci number 10000 : 3364476487...9947366875
    Digits in the Fibonacci number: 2090


    Съответно резултатът също е променен. Мисля по втората част на задачката...

  6. Първо някои пояснения и дефиниции:

    1. Просто число е това това положително цяло число, което се дели само на едно и на самото себе си, т.е не се разлага на множители.

    2. Алгоритъмът за намиране на прости числа е съвсем… прост – първо се написват всички числа от едно до безкрайност 😉 След това се започва от 3 (1 и 2 са прости) и се пробва да се раздели на всички преди него. Ако остатъка е цяло число и е нула, то значи числото не е просто (т.е. е съставно и може да се разложи на множители). После преминаваме към следващото число – 4 и така нататък. Проблемът при тази дефиниция е, че ако искаме да разберем дали числото N е просто, то трябва да определим дали всички числа до N-1 са също прости. С малко еквилибристики (като например премахването на всички четни числа, които се делят на 2 с остатък нула и следователно са съставни) се определя че е достаъчно да знаем дали числата до квадратен корен от N са прости. Все пак за много големи числа се появява проблем не само с пресмятането им, а също и със съхраняването. За достатъчно големи прости числа може да се окаже че компютърната памет е недостатъчна да съхрани всички прости числа в интервала 1..sqrt(N). Повече по въпроса мове да се намери в Уикипедия.

    3. За да се избегне горепосочения проблем се използват математическо “сито”, което определя с голяма вероятност (над 99.9999 %) дали N е просто или не. За практически цели това е достатъчно, макар че може да даде и лъжливи резултати (все пак грешка от 0.0001 процент при 1000000 числа дава 100% вероятност за грешка, т.е. едно сигурно лъжливо попадение). Използваните практически алгоритми в използваната по-долу програма (по-точно използваната библиотеката) реализира точно такова сито. Предимството на ситото е неговото скорост и липсата на изискването да се изчисляват и пазят огромен брой прости числа. На този блог постинг може да се намери много интересно описание на алгоритмите за намиране на просто число.

    След това уточнение ето и програмката. Този път съм използвал библиотеката GMPY, която е написана на C и съответно скороста на изчисления е многократно по-голяма отколкото на “чист” Python.


    import gmpy
    import time
    
    def genPrimeFib(amount):
        primeFibList =[]
        counter = 0
        while len(primeFibList) < amount:
            f = gmpy.fib(counter) # generate Fibonacci number
            if gmpy.is_prime(f):
                primeFibList.append(f)
                # Show the progress
                print "Found prime Fibonacci number ", len(primeFibList)
            counter +=1
        return primeFibList
    

    А ето и какво получаваме за проверка на първите 20 прости Фибоначи числа:


    fibList = genPrimeFib(20)
    
    for f in fibList:
        print "Fibonacci prime number : ", f
    
    

    Found prime Fibonacci number 1
    Found prime Fibonacci number 2
    Found prime Fibonacci number 3
    Found prime Fibonacci number 4
    Found prime Fibonacci number 5
    Found prime Fibonacci number 6
    Found prime Fibonacci number 7
    Found prime Fibonacci number 8
    Found prime Fibonacci number 9
    Found prime Fibonacci number 10
    Found prime Fibonacci number 11
    Found prime Fibonacci number 12
    Found prime Fibonacci number 13
    Found prime Fibonacci number 14
    Found prime Fibonacci number 15
    Found prime Fibonacci number 16
    Found prime Fibonacci number 17
    Found prime Fibonacci number 18
    Found prime Fibonacci number 19
    Found prime Fibonacci number 20
    Fibonacci prime number : 2
    Fibonacci prime number : 3
    Fibonacci prime number : 5
    Fibonacci prime number : 13
    Fibonacci prime number : 89
    Fibonacci prime number : 233
    Fibonacci prime number : 1597
    Fibonacci prime number : 28657
    Fibonacci prime number : 514229
    Fibonacci prime number : 433494437
    Fibonacci prime number : 2971215073
    Fibonacci prime number : 99194853094755497
    Fibonacci prime number : 1066340417491710595814572169
    Fibonacci prime number : 19134702400093278081449423917
    Fibonacci prime number : 475420437734698220747368027166749382927701417016557193662268716376935476241
    Fibonacci prime number : 529892711006095621792039556787784670197112759029534506620905162834769955134424689676262369
    Fibonacci prime number : 1387277127804783827114186103186246392258450358171783690079918032136025225954602593712568353
    Fibonacci prime number : 3061719992484545030554313848083717208111285432353738497131674799321571238149015933442805665949
    Fibonacci prime number : 10597999265301490732599643671505003412515860435409421932560009680142974347195483140293254396195769876129909
    Fibonacci prime number : 36684474316080978061473613646275630451100586901195229815270242868417768061193560857904335017879540515228143777781065869
    

    Причината да не пиша за 10000-ното просто число е… че не разполагам със суперкомпютър 🙁 . По горния алгоритъм на машина с Pentium 4 на 3.0Ghz отне 1117 секунди (близо 19 минути) , за да сметне първите 25 прости числа от редицата на Фибоначи. Приемам идеи за подобряване на алгоритъма.

  7. Жоро – имам още една задачка за тебе 🙂

    Ето модифицирана версия на предходния скрипт, която измерва и времето за изчисление на всяко просто Фибоначи число:


    import gmpy
    import time
    
    def genPrimeFib(amount):
        primeFibList =[] # will contain tuples (time_taken, prime_Fibonacci_number)
        counter = 0
        start = time.time()
        while len(primeFibList) < amount:
            f = gmpy.fib(counter) # generate Fibonacci number
            if gmpy.is_prime(f):
                timeTaken = time.time() - start # capture the time it took for the computation
                primeFibList.append((timeTaken,f))
                print "Found prime Fibonacci number ", len(primeFibList), " for ", timeTaken, " seconds"
                start = time.time() # reset start time
            counter +=1
        return primeFibList
    
    totalStartTime = time.time() 
    
    # Compute the first 23 prime Fibonacci numbers
    a = genPrimeFib(23)
    
    totalTime = time.time() - totalStartTime
    for t,n in a:
        print "Time: ", t, " seconds"
        print "Number : ", n
       
    print "Total time in seconds: ", totalTime
    

    А ето и резултата от намирането на първите 23 прости числа на Фибоначи и времето което отне намирането им на машина Pentium 4 на 3.0GHz:


    Found prime Fibonacci number 1 for 0.0 seconds
    Found prime Fibonacci number 2 for 0.0 seconds
    Found prime Fibonacci number 3 for 0.0 seconds
    Found prime Fibonacci number 4 for 0.0 seconds
    Found prime Fibonacci number 5 for 0.0 seconds
    Found prime Fibonacci number 6 for 0.0 seconds
    Found prime Fibonacci number 7 for 0.0 seconds
    Found prime Fibonacci number 8 for 0.0 seconds
    Found prime Fibonacci number 9 for 0.0 seconds
    Found prime Fibonacci number 10 for 0.0 seconds
    Found prime Fibonacci number 11 for 0.0 seconds
    Found prime Fibonacci number 12 for 0.0 seconds
    Found prime Fibonacci number 13 for 0.0 seconds
    Found prime Fibonacci number 14 for 0.0 seconds
    Found prime Fibonacci number 15 for 0.0160000324249 seconds
    Found prime Fibonacci number 16 for 0.0150001049042 seconds
    Found prime Fibonacci number 17 for 0.0160000324249 seconds
    Found prime Fibonacci number 18 for 0.0469999313354 seconds
    Found prime Fibonacci number 19 for 0.0469999313354 seconds
    Found prime Fibonacci number 20 for 0.0780000686646 seconds
    Found prime Fibonacci number 21 for 0.0160000324249 seconds
    Found prime Fibonacci number 22 for 18.388999939 seconds
    Found prime Fibonacci number 23 for 79.1040000916 seconds
    Time: 0.0 seconds
    Number : 2
    Time: 0.0 seconds
    Number : 3
    Time: 0.0 seconds
    Number : 5
    Time: 0.0 seconds
    Number : 13
    Time: 0.0 seconds
    Number : 89
    Time: 0.0 seconds
    Number : 233
    Time: 0.0 seconds
    Number : 1597
    Time: 0.0 seconds
    Number : 28657
    Time: 0.0 seconds
    Number : 514229
    Time: 0.0 seconds
    Number : 433494437
    Time: 0.0 seconds
    Number : 2971215073
    Time: 0.0 seconds
    Number : 99194853094755497
    Time: 0.0 seconds
    Number : 1066340417491710595814572169
    Time: 0.0 seconds
    Number : 19134702400093278081449423917
    Time: 0.0160000324249 seconds
    Number : 475420437734698220747368027166749382927701417016557193662268716376935476241
    Time: 0.0150001049042 seconds
    Number : 529892711006095621792039556787784670197112759029534506620905162834769955134424689676262369
    Time: 0.0160000324249 seconds
    Number : 138727712780478382711418610318624639225845035817178[...cut here...]5225954602593712568353
    Time: 0.0469999313354 seconds
    Number : 306171999248454503055431384808371720811128543235373[...cut here...]8149015933442805665949
    Time: 0.0469999313354 seconds
    Number : 105979992653014907325996436715050034125158604354094[...cut here...]3254396195769876129909
    Time: 0.0780000686646 seconds
    Number : 366844743160809780614736136462756304511005869011952[...cut here...]0515228143777781065869
    Time: 0.0160000324249 seconds
    Number : 960412006189225538239428833609248650261049174118770[...cut here...]4331637646183008074629
    Time: 18.388999939 seconds
    Number : 357103560641909860720907774139063454445569926582843[...cut here...]7148642001438470316229
    Time: 79.1040000916 seconds
    Number : 500195636126957292905024512596972806695803345136243[...cut here...]4302854387700053591957
    Total time in seconds: 97.728000164
    Script terminated.
    

    Задачката за тебе, Жоро е следната – като използваш посочените по-горе времена, определи формула по-която (приблизително) може да се изчисли колко време ще отнеме на същата тази машина да сметне 10000-ното число. С други думи, базирано на горните времена есктраполирай времето за намиране на 10000-ното просто Фибоначи число.

    И нека силата бъде с тебе 😛

  8. Само да добавя, че и Maxima използва същото математическо сито за генериране и проверка на прости числа както и GMPY. Така че този път силите са изравнени. За справка – документацията (PDF) на Maxima и документацията (PDF) на GMP (на чиято основа GMPY е базирана).

  9. На пръв поглед, гледайки времената, не мисля, че ще има функция, която да изчисли времето за n, но обещавам, като се върна от работа ( сега трябва да тръгвам), да помисля повече 🙂

    Ето и времената на моят компютър ( Athlon 4000+ ) под Линукс:
    Found prime Fibonacci number 16 for 0.0152499675751 seconds
    Found prime Fibonacci number 17 for 0.00510787963867 seconds
    Found prime Fibonacci number 18 for 0.00592708587646 seconds
    Found prime Fibonacci number 19 for 0.0155727863312 seconds
    Found prime Fibonacci number 20 for 0.015280008316 seconds
    Found prime Fibonacci number 21 for 0.00786304473877 seconds
    Found prime Fibonacci number 22 for 6.78547692299 seconds
    Found prime Fibonacci number 23 for 23.3749740124 seconds

  10. Не съм съгласен че домашните ще се решават дори с многоредова програмка. Никоя от домашните, които съм получавал, не е била изчислителна, а символна задача – дори неща от рода на опростяване на дроби (примерно (5 – й^4)(й^2-4)/(й^2 – 1)), доказване на разни хитри и на пръв поглед очевидни леми, които лесничко се смятат за числа от рода на 10,000,000 но пък са доста трудни да се докажат при произволно n.

    Според мен числените методи и математиката нямат нищо общо 😉 Макар че матлаб може да извърши опростяването без да извърши деленето и да излезе резултат от рода на 3.575757575~

    Относно протоколите и прочее – писал съм софт за смятането на почти всички задачи от ИГМ и разни други дисциплини, включително програмка за сплайниране на Паскал, като кода беше адаптиран от Фортран. Още си пазя книжката – “Числени методи за програмиране”, сички примери на Фортран. Много е забавно наистина. И почти толкова безполезно когато говорим за математиката която се изучава във “висшите” ни учебни заведения. Хъм, заведения е правилната дума! 🙂

  11. Някой от вас има ли интересни логически задачи, които може да се претворят визуално във Flash, ако някой има добра идея да пише

    Коментарът е редактиран – премахнати са рекламни линкове.

  12. аз имам задачка за вас (на която не знам решението, защото не я разбирам):
    да се разложи произволно цяло число на всички възможни суми на числа на Фибоначи, в които няма повторение на еднакви числа.

  13. @g.ivanov: Задачата в твоя въпрос може да се опрости по следния начин:
    1. Нека избраното произволно цяло число е N. Да се намерят всички числа на Фибоначи, които са <= N. Нека да означим тази поредица от числа F. Това вече го направихме в предишните коментари. 2. Да се намерят всички комбинации от числа в предицата F чиято сума е равна на N. Решението на тази задачка ще оставя за теб и останалите читатели 🙂

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.