Şifreli matematik formülü kodda

Python'da normal dağılımı hesaplamak için bir işleve sahibim:

def norm_cdf(z):
  """ Use the norm distribution functions as of Gale-Church (1993) srcfile. """
  # Equation 26.2.17 from Abramowitz and Stegun (1964:p.932)

  t = 1.0/(1+0.2316419*z) # t = 1/(1+pz) , p=0.2316419
  probdist = 1 - 0.3989423*math.exp(-z*z/2) * ((0.319381530 * t)+ \
                                         (-0.356563782* math.pow(t,2))+ \
                                         (1.781477937 * math.pow(t,3)) + \
                                         (-1.821255978* math.pow(t,4)) + \
                                         (1.330274429 * math.pow(t,5)))
  return probdist

Fakat PEP8 ve 80 kar marjına bağlı kalmak zorundayım, bu yüzden çirkin \ s. Kodumu başka nasıl doğrulamalıyım?

Matematiksel olarak

$$\begin{align*} \textrm{norm_cdf}(z) = 1 - 0.3989423 e^\frac{-z^2}{2} (&1.330274429 t^5 - 1.821255978 t^4 \\ &+ 1.781477937 t^3 - 0.356563782 t^2 + 0.319381530 t ) \end{align*}$$

nerede

$$ t = \ kırık {1} {1 + 0,2316419 z} $$

42
Python'daki Jeez yorumları çirkindir ... başka bir notta, matematiksel sabitleriniz için C tarzı typedef gibi bir şey yapabilir misiniz? Kodu çok daha okunaklı kılacaktı, imo. typdef 0.319381530 NORMALITY gibi bir şey. Ama biliyorsun, daha iyi isimler seç. Python için, NORMALITY = 0.319381530 gibi bir yöntemle tüm sabitlerinizi yöntemin üstünde (okunabilirlik için) kullanabilirsiniz. Zaten matematiksel fonksiyonları organize etmek için tercih ettiğim yöntem bu.
katma yazar tai, kaynak
@Davidmh Ah, tamam. Bu benim yanlış anlaşmamdı :)
katma yazar tai, kaynak
@Davidmh Python'da typedef olmadığını biliyorum, bu yüzden NORMALITY = 0.319381530 gibi bir "sabit" tanımlamayı önerdim. Şahsen, @ Davidmh’in cevap bölümünde önerildiği gibi liste yaklaşımını beğenmedim: coeff = [1 , 0.319381530, -0.356563782, 1.781477937, -1.821255978, 1.330274429] Bir çok sabit varsa, katsayı [0] , katsayı [ 1] vb. Uzun bir fonksiyonda, özellikle sabitleri defalarca kullanan fonksiyonlarda okunamaz hale gelmeye başlar. Yine, sadece benim görüşüm.
katma yazar tai, kaynak
@ChrisCirefice Er, typedef 'in C’de çalışması böyle değil.
katma yazar Nic Hartley, kaynak
@ 200_success formül görüntüsü için seviyor !!
katma yazar hippytea, kaynak
@ChrisCirefice Bir liste bir döngüde kullanışlıdır, katsayılar temalar tarafından özel bir anlama sahip değildir ve yalnızca polinomdan türetilen ifadelerde (türevler, integraller, vb.) Tekrar kullanılabilirler.
katma yazar Davidmh, kaynak
@ChrisCirefice Python'da typedefs gibi bir şey yoktur. Katsayıları bir listede kaydetmek daha okunaklı bir yaklaşımdır.
katma yazar Davidmh, kaynak

8 cevap

Harika kitabı C ++ Nümerik Tarifler (ama aynı zamanda uygulanabilir) alıntı yapalım:

Bir polinomu bu şekilde değerlendirmek için asla asla bilmediğinizi biliyoruz:

<�Ön> p = C [0] + c [1] * x + c [2] * x * x + c [3] * x * x * x + c [4] * x * x * x * x;      

veya (daha da kötüsü!),

<�Ön> p = C [0] + c [1] * x + c [2] * pow (x, 2.0) + c [3] * pow (x, 3.0) + c [4] * pow (x, 4.0);      

Gelip (bilgisayar) devrimi, bu tür suç davranışlarından suçlu bulunan tüm kişiler özetle uygulanacak ve programları olmayacak!

(Basımınızda sayfayı analitik indeksinde "özellikle kötü şeyler" girişi altında bulabilirsiniz. Bu kitabı çok seviyorum.)

Bunu yapmamak için iki neden var: doğruluk ve performans. Polinomu değerlendirmenin doğru yolu şudur:

-t * (0.319381530  +  t * (-0.356563782 + t * (1.781477937 + t * (-1.821255978 + 1.330274429 * t))))

And you can, of course, split at your convenience, as the newlines inside parenthesis are ignored. Remember the PEP: " The preferred place to break around a binary operator is after the operator, not before it."

-t * (0.319381530  +  t * (-0.356563782 +
    t * (1.781477937 + t * (-1.821255978 + 1.330274429 * t))))

Diğer bir seçenek de katsayıları bir listeye kaydetmektir:

coeff = [0, 0.319381530, -0.356563782, 1.781477937, -1.821255978, 1.330274429]
poly = coeff[-1]
for c in coeff[-2::-1]:
    poly *= x
    poly += c

Bellek ayırmaktan ve dağıtmaktan kaçınmak için işlemler yapıyorum, ancak bu yalnızca x bir NumPy dizisi ise geçerlidir. Tek bir sayı üzerinde değerlendirme yapıyorsanız, bunun yerine daha güzel olan ifadeyi kullanabilirsiniz:

poly = poly * x + coeff[i]

Ama ilkine sadık kalacağım çünkü daha genel.

Tabii ki, sonuç prefactor ile çarpılmalıdır:

return 1 - 0.3989423*math.exp(-z*z/2) * poly

Ya da yerinde yapmak istiyorsanız:

z2 = z * z # Be careful not to modify your input!
z2 *= 0.5  # Multiplication is faster than division.
np.exp(z2, out=z2)

probd = z2 * poly
probd *= -0.3989423
probd += 1
return probd

Bonus izleme:

Bu işlevi büyük dizilere (binden fazla sayıya) uyguluyorsanız, ilk tekniği numexpr'de kullanmaktan yararlanabilirsiniz:

expr += '1 - 0.3989423* exp(-z * z/2) * '
expr += '(-t * (0.319381530  +  t * (-0.356563782 +  t * '
expr += '(1.781477937 + t * (-1.821255978 + 1.330274429 * t)))))'
ne.evaluate(expr)

Bu, ifadeyi sizin için derleyecektir ve sahip olduğunuz kadar çekirdeği şeffaf bir şekilde kullanacaktır.

40
katma
@EmilioMBumachar Polinomun durumuna göre değişir. Ancak pratikte kullanılan birçok polinom için naif yaklaşım, mutlak hatayı arttıran benzer büyüklükteki değerlerin çıkarılmasıyla sonuçlanacaktır;
katma yazar Ben Crowell, kaynak
Uygulamanın ağır performans cezalarını cesaretini kırılmış yollardan birinde görebiliyorum, ancak doğruluk cezalarını göremiyorum. Neden doğruluk için kötü? SO'dan sordum: stackoverflow.com/q/25203873/174365
katma yazar Czyrek, kaynak
Bu polinomları değerlendirme yöntemi "Horner'ın şeması" (veya nerede yaşadığınıza bağlı yöntem) olarak bilinir.
katma yazar mikeazo, kaynak
Belki bir bağlamı özlüyorum, ancak hız uğruna daha uzun, daha çirkin bir kod çıkarıyor gibi görünüyorsunuz, bu da OP için bir endişe gibi görünmüyor.
katma yazar jas663, kaynak
Bu ifade için "sahip olduğunuz sayıda çekirdek" kullanmanın bir faydası olduğunu sanmıyorum. Etkili bir şekilde paralel hale gelmek çok önemsiz.
katma yazar Yahya Uddin, kaynak
@Ruslan "büyük diziler için". Aynı işlevi uyguladığınızda 1e4 sayıları paralel yürütmeden faydalanır.
katma yazar Davidmh, kaynak
@FranciscoCouzo mükemmel bir gelişmedir! Düzenledim.
katma yazar Davidmh, kaynak
@ raptortech97 Her iki seçenek de veriyorum. Neyse, önemli olan, polinom değerlendirmesinin doğru olması için doğru yapılması gerektiğidir.
katma yazar Davidmh, kaynak
xrange (len (katsayı) - 2, -1, -1) kullanmak ve ardından katsayıi ile indekslemek yerine, katsayısı [-2 :: - 1] ile yinelenir
katma yazar Francisco Couzo, kaynak

As it turns out, a similar question was asked recently on Math.SE. Rather than , take advantage of built-in functionality in Python.

norm_cdf (z) 'nuz sadece sayısal bir yaklaşımdır.

$$ P (z) = \ kırık {1} {\ sqrt {2 \ pi}} \ int _ {- \ infty} ^ {z} e ^ {- t ^ 2/2} \ dt = \ int _ {- \ infty} ^ {z} Z (t) \ dt = \ frac {1} {2} \ left (1 + \ mathrm {erf} \ left (\ frac {z} {\ sqrt {2}} \ sağ) \ sağ) = \ frac {1} {2} \ mathrm {erfc} \ left (- \, \ frac {z} {\ sqrt {2}} \ sağ) $$

Bu nedenle, yalnızca math.erfc() (Python 2.7'den beri kullanılabilir) ve daha doğru bir sonuç elde edin (özellikle \ $ z \ $ negatif değerleri için).

import math

def norm_cdf(z):
    return 0.5 * math.erfc(-x/math.sqrt(2))

Daha da iyisi, sadece scipy.stats.norm. cdf() !

22
katma
@GraniteRobert Aksine, sorudaki norm_cdf (z) , cruder yaklaşımıdır. norm_cdf (0) , 0.5 yerine 0.499999975962'dir. Herhangi bir z ≤ -2.0 negatif sonuç verir. 0.5 * math.erfc (-x/math.sqrt (2)) bu sorunlardan hiçbirine sahip değildir.
katma yazar wjv, kaynak
Veya bu erf() genel uygulamasını çalmanız yeterlidir.
katma yazar wjv, kaynak
Python ≤ 2.6 ile çalışırken neden SciPy ürününü gerektirmiyor veya yeniden dağıtmıyorsunuz?
katma yazar wjv, kaynak
@alvas Tüm verileriniz ortalamanın üstünde mi?
katma yazar wjv, kaynak
aslında norm_cdf() , her zaman başka bir maliyet fonksiyonunun mutlak bir değerini alırsa, tamamdır.
katma yazar hippytea, kaynak
bunun nedeni, herhangi bir cümle uzunluğu için alt sınırın 0 olması nedeniyle girişlerin hiçbir zaman negatif olamamasıdır.
katma yazar hippytea, kaynak
scipy.stats.norm.cdf() içeri alınırken bir try-hariç tam koda sahibim, sadece bu durumda, kullanıcılar scipy yüklememişlerdi ama math.erfc iyi olurdu.
katma yazar hippytea, kaynak
ama sanırım <= py 2.6 kullanıcısı ...
katma yazar hippytea, kaynak
"Aynı" hesaplamayı yapmak için farklı bir algoritmaya geçmeden önce dikkatli olun. Normal dağılım fonksiyonu için Abramowitz ve Stegun (A&S) 'deki 8 farklı yaklaşımdan biri bu özel yaklaşımı seçti. Seçim için çok iyi bir sebep olabilirdi. Hesaplamanın çirkin detaylarını anlamak için zaman ayırmazsanız, kullanıma hazır bir kütüphane fonksiyonunun eldeki problem için yeterli doğruluk sağlayacağını varsaymayın. A & S'de dolaşan insanlar sıradan kullanıcılar olma eğilimindedir; doğruluk ve hata sınırları çok önemli olabilir.
katma yazar atos, kaynak
@ 200_success: A & S 26.2.17, geçerli girdilerin negatif olmadığını açıkça söylüyor. Kod negatif değerlerle çağrılmamalı ve geçersiz aramaların muhtemel olduğu bir ortamda, belki bir kontrol içermelidir. Hata, 7.5 * 10 ** -8 olarak listelenmiştir. Doğruluk bir şekilde math.erfc için belirtilmiş mi? Herhangi bir algoritmanın daha iyi ya da daha kötü olduğunu öne sürmüyorum, sadece değişimlerin çözülmekte olan sorun bağlamında anlaşılması gerektiğini söylüyorum.
katma yazar atos, kaynak

Sadece birkaç gözden geçiricinin bahsettiği “sihirli sayılar” ı ele alacağım.

Bazen, saf matematikte çalışırken, ilk bakışta göründüğü gibi görünen şey "sihirli sayı" aslında değil . Rakamların kendileri de sorun ifadesinin bir parçası olabilir. Sanırım soru şuna bağlı: Numaradan daha açıklayıcı bir isim bulabilir misiniz? İyi bir isim varsa, muhtemelen kullanmalısın.

İlk bakışta, numaralarınızın sorunun doğal bir parçası olduğunu düşündüm. Fakat Abramowitz ve Stegun'a baktığımda, başvurulan formülün zaten çirkin görünümlü sabitlerinizi adlandırdığını gördüm. İsimler p (bir yorumda bahsettiğiniz) ve b1 ile b5 . Bu adları kodda kullanmalısınız, çünkü orijinal formül tanımına çok açık bir bağlantı oluştururlar.

p = 0.2316419 yorumunu eklemenin iyi bir fikir olduğuna karar verdiğinde, basılması gereken sayının adlandırılması çok güçlü bir kanıt oldu. (Ve kod , p = 0.2316419 dediğinde, yorum kaldırılmalıdır.)

Bu arada, yoruma tam olarak Abramowitz ve Stegun referansını dahil etmeniz çok iyi oldu.

18
katma

math.pow yerine, yerleşik ** operatörünü kullanın. EOL'de \ s öğesine ihtiyacınız yoktur, çünkü ifadeyi çevreleyen parantez içinde satır satırlarının birden fazla satır örtülmesine izin verilir. Bu iki değişikliği de yaptıktan sonra:

def norm_cdf(z):
    """ Use the norm distribution functions as of Gale-Church (1993) srcfile. """
    # Equation 26.2.17 from Abramowitz and Stegun (1964:p.932)

    t = 1.0/(1 + 0.2316419*z) # t = 1/(1+pz) , p=0.2316419
    probdist = 1.0 - (   (0.319381530  * t)
                       + (-0.356563782 * t**2)
                       + (1.781477937  * t**3)
                       + (-1.821255978 * t**4)
                       + (1.330274429  * t**5)) * 0.3989423*math.exp(-z*z/2)
    return probdist

Ben de önceliği biraz daha açık ve okunaklı kılmak için çarpmalardan birini yeniden düzenlemiştim.

Bu değişikliklerden sonra, hala görebildiğim tek sorun tüm sihirli sayılar. Bu sabitlere nasıl ulaşıldığını bilmiyorum, ancak sabitlere anlamlı isimler verilebilirse okunabilirliğe yardımcı olabilir. Bazen formüllerle, gerçekte verilecek anlamlı bir isim yoktur.

16
katma

Bunun işe yarayıp yaramadığından emin değilsiniz ancak bir polinomun değerini belirli bir pozisyonda değerlendirmek için kolayca bir fonksiyon tanımlayabilirsiniz.

def evaluate_polynom(pol, x):
    return sum(a * math.pow(x, i) for i, a in enumerate(pol))

Sonra

(0.319381530 * t) + (-0.356563782* math.pow(t,2)) + (1.781477937 * math.pow(t,3)) + (-1.821255978* math.pow(t,4)) + (1.330274429 * math.pow(t,5))

olur:

evaluate_polynom([0, 0.319381530, -0.356563782, 1.781477937, -1.821255978, 1.330274429], t)
15
katma
Açıkçası bunu anlamlı bir ismi olan bir fonksiyona koymak için yeterli bilgiye sahip değilim, ancak bunun gerekli olduğundan bile emin değilim. Katılan matematiksel nesneler gibi görünmesini sağlamak için gerekli soyutlama düzeyini getirdiğini düşünüyorum ve bunu farklı bir yere taşımak matematiksel formülü parçalara ayırabilir.
katma yazar Mike, kaynak
Sahip olduğum tek nitpick, daha açık hale getirmek için pol ile katsayı veya katsayılarını yeniden adlandırmak veya muhtemelen yerleşimi aşağıda açıklandığı şekilde kullanmak olabilir.
katma yazar jmagnusson, kaynak
Eğer bu doğruysa, tanımladığınız kod pasajını çağıran başka bir işlevi tanımlayabilirsiniz, bir tür yardımcı yöntem.
katma yazar Pimgd, kaynak

Denetleyicileri (ayrıca) kızdırma riski altında, Nümerik Tarifler kitaplarının 1 herhangi birinden kaynaklanan nümeriklerle ilgili herhangi bir tavsiyede bulunmaya karşı tavsiyede bulunacağım.

Horner'ın yöntemi bazı durumlarda yardımcı olabilirken, her derde deva uzun bir yoldur. Değerleri toplarken, yuvarlama sorununu düzeltmek yerine, bu sorunu önlemeye çalışır. Maalesef, bunu yaparken sadece kısmen başarılı. Polinomların çoğu için, sonuçların göreceli olarak düşük olduğu, hatta en iyi ihtimalle, girdiler olacaktır.

Polinom terimlerinin ürettiği değerler toplanırken sayısal kararsızlığa yol açabilirse, her birini ayrı ayrı oluşturmayı düşünebilirsiniz, sonra Bu terimleri toplamak için Kahan toplamı . Eğer umursuyorsanız, bu size toplamın kendisiyle birlikte bir hata payı da verebilir.

Muhtemelen daha iyi hala Langlois, et al telafi edilmiş Horner'ın planı . En azından en son dikkatlice baktığımda, bu polinomların değerlendirilmesinde en son teknolojinin durumu gibi görünüyordu. Horner'ın şemasını, hassasiyetin iki katı olan kayan nokta sayılarını kullanarak elde ettiğinizde elde edeceğiniz sonuç kabaca aynı kesinliği korur (örneğin, 64-bit bir çift kullanarak, Horner'ın şemasıyla aynı şekilde 128-bit dörtlü) hassas kayan nokta, ancak normalde taşıyacak neredeyse hız cezası olmadan). Kahan toplamı gibi, bu da sonuçta bağlı bir hatanın hesaplanmasını destekler.


1. Başvuru için aşağıdaki eleştirilere bakın:
http : //www.stat.uchicago.edu/~lekheng/courses/302/wnnr/nr.html
http://www.lysator.liu.se/c/ num-tarifleri-in-c.html
Bence "harika" ifadesinden daha kesin bir özet: " Nümerik Tarifler 'in bir insanın başını derde sokmasına yetecek kadar bilgi sağladığını buldum, çünkü NR'yi okuduktan sonra bir kişinin anlayacağını düşünüyor neler oluyor."

6
katma

buna karşılık:

a*x**3 + b*x**2 + c*x = ((a*x + b)*x + c)*x

h/t @ "Horner" referansı için Emily L.:

https://en.wikipedia.org/wiki/Horner%27s_method

ve h/t @ Davidmh, bu yöntemin hesaplama hızında/hassasiyetinde gelişmeler olduğunu belirtti.

gale-church, 1990'da şöyle belirtti:

import math

def pnorm(z):

    t = 1/(1 + 0.2316419 * z)
    pd = (1 - 0.3989423 *
      math.exp(-z * z/2) *
        ((((1.330274429 * t - 1.821255978) * t
           + 1.781477937) * t - 0.356563782) * t + 0.319381530) * t)

    return pd

Bu yöntem rahatça t ^ n sorununu önler.

atıf:

enter image description here

enter image description here

kaynak:

http://www.aclweb.org/anthology/J93-1004

sayfa 28/28 pdf

Derginin 95. sayfası Hesaplamalı Dilbilim Cilt 19, Sayı 1

Ben "güzelleştirmek" olabilir:

def pnorm(z):

t = 1/(1 + 0.2316419 * z)
pd = (1 - 0.3989423 * math.exp(-z * z/2) *
      ((((1.330274429 * t - 
          1.821255978) * t + 
          1.781477937) * t - 
          0.356563782) * t + 
          0.319381530) * t )

return pd

eğer kontrol edersen

Abromowitz ve Stegun, Matematiksel Fonksiyonlar El Kitabı

sayfa 932 denklemi 26.2.17

atıf:

http://people.math.sfu.ca/~cbm/aands/ page_932.htm

bunu bulacaksınız:

enter image description here

bize bir tablo yaratabiliriz:

def pnorm(z):

    p  =  0.2316419
    b1 =  0.319381530
    b2 = -0.356563782
    b3 =  1.781477937
    b4 = -1.821255978
    b5 =  1.330274429
    t = 1/(1 + p * z)
    pd = (1 - 0.3989423 * math.exp(-z * z/2) *
          ((((b5 * t + b4) * t + b3) * t + b2) * t + b1) * t )

    return pd

Sonra önceki sayfadan; 931 bulacaksınız:

enter image description here

Zx = (1/√(2* π))*e(-z*z/2)

python içinde:

Zx = (1/math.sqrt(2* math.pi))*math.exp(-z*z/2)

ve şunu bulduk (1/√ (2 * π)) = 0.3989423

ayrıca, sanırım bu hoşuma gitti:

t * (b1 + t * (b2 + t * (b3 + t * (b4 + t * b5))))

daha iyi:

(((b5 * t + b4) * t + b3) * t + b2) * t + b1) * t 

öyleyse, nihayet:

import math

def pnorm(z):

    p  =  0.2316419
    b1 =  0.319381530
    b2 = -0.356563782
    b3 =  1.781477937
    b4 = -1.821255978
    b5 =  1.330274429
    t  = 1/(1 + p * z)
    Zx = (1/math.sqrt(2 * math.pi)) * math.exp(-z * z/2)
    pd = Zx *  t * (b1 + t * (b2 - t * (b3 + t * (b4 - t * b5))))

    return (1 - pd) 

çalışmalarımın operasyona karşı kontrol edilmesi

import matplotlib.pyplot as plt
import numpy as np
import math




def norm_cdf(z):
  """ Use the norm distribution functions as of Gale-Church (1993) srcfile. """
  # Equation 26.2.17 from Abramowitz and Stegun (1964:p.932)

  t = 1.0/(1+0.2316419*z) # t = 1/(1+pz) , p=0.2316419
  probdist = 1 - 0.3989423*math.exp(-z*z/2) * ((0.319381530 * t)+ \
                                         (-0.356563782* math.pow(t,2))+ \
                                         (1.781477937 * math.pow(t,3)) + \
                                         (-1.821255978* math.pow(t,4)) + \
                                         (1.330274429 * math.pow(t,5)))

  return probdist

for z in np.arange (-3,3,0.01):
    zf = pnorm(z)
    plt.plot(z,zf, c='red', marker = '.', ms=1)

for z in np.arange (-3,3,0.01):
    zf = norm_cdf(z)+0.1 #offset 0.1
    plt.plot(z,zf, c='blue', marker = '.', ms=1)

plt.show()
plt.pause(0.1)

enter image description here

Horner yönteminin daha hızlı olmasını bekliyordum, bu yüzden yerine geçerek bir zaman testi yaptım:

#Zx = (1.0/math.sqrt(2.0 * math.pi)) * math.exp(-z * z/2.0)
Zx = 0.3989423* math.exp(-z * z/2.0)

adil kılmak ve np.arrange kararını 0.0001’e yükseltmek için:

t0 = time.time()
for z in np.arange (-3,3,0.0001):
    zf = pnorm(z)
t1 = time.time()
for z in np.arange (-3,3,0.0001):
    zf = norm_cdf(z)
t2 = time.time()

print ('pnorm time    : %s' % (t1-t0))
print ('norm_cdf time : %s' % (t2-t1))

ve dört çekirdekli AMD 7950 FM2 + w/16G’yi döndürerek elde ettiğim sonuçlar oldukça zor (bazı diğer uygulamaları da olsa da)

>>>
pnorm time    : 81.4725670815
norm_cdf time : 80.7865998745

Horner yöntemi daha hızlı değildi

2
katma
Ayrıca MathJax ile formülleri geliştirmeyi de tavsiye ederim. Bunu nasıl yapacağını bilmiyorsan, daha sonra çalışabilirim.
katma yazar Jamal, kaynak
StackExchange Kod İncelemesine Hoş Geldiniz! Harika bir tarih, paylaşım için teşekkürler.
katma yazar Stephen Rauch, kaynak
hmm Mathjax ile daha önce hiç çalışmadım, sonucu görmek isterdim Jamal; hoşgeldin Stephen Rauch için teşekkürler
katma yazar litepresence, kaynak

Burada dürüst olacağım: Tamamen pitona yatarım.

Bununla birlikte, bu sayısal değişmezlerin bir kısmını sabit olarak bildirmek mümkün mü? Bu, kodunuzu temizler ve kodun kendisini biraz daha açıklığa kavuşturur.

1
katma