Haskell, plötslig prestandapuckel vid större beräkning.

2017-05-08, 22:06 #1

Moderator

Reg: Jul 2010

Inlägg: 12 961

Inspirerad av (FB) Hur lär man sig stora kvadrater och deras rötter? testade jag att koda https://en.wikipedia.org/wiki/Method..._roots#Example vilket resulterade i följande kod:

Citat:

sqrt' :: (Floating a, Ord a) => a -> (a,Int)
sqrt' a = f a (1,0)
where
f a (s,n)
| (((s+a/s)/2)**2-a) < 0.000000000001 = (s,n+1)
| otherwise = f a (((s+a/s)/2),n+1)

Den lilla tupeln, (a,Int) finns med för att Int-värdet visar hur många varv beräkningen behöver snurra.
I jämförelse med den inbyggda funktionen sqrt räknar min funktion tillräckligt bra,
sqrt' 12345678900987654321 = (3.5136418513220444e9,36)
sqrt 12345678900987654321 = 3.513641828785008e9
de första 7 decimalerna är lika och det tar inte någon av mig uppfattbar tid.

Men...

När jag kör
sqrt' 1000000000000 så kommer svaret, (1000000.0000000054,25), direkt, men när jag kör
sqrt' 10000000000000 så kommer svaret aldrig (fast jag har som mest väntat någon minut).

Vad kan det bero på?

GHC 7.10.3. Samma sak oavsett om jag kompilerar eller kör i GHCI.

Citera

2017-05-08, 22:11 #2

Medlem

Reg: Nov 2015

Inlägg: 1 237

Citat:

Ursprungligen postat av Neksnor

Inspirerad av (FB) Hur lär man sig stora kvadrater och deras rötter? testade jag att koda https://en.wikipedia.org/wiki/Method..._roots#Example vilket resulterade i följande kod:

Den lilla tupeln, (a,Int) finns med för att Int-värdet visar hur många varv beräkningen behöver snurra.
I jämförelse med den inbyggda funktionen sqrt räknar min funktion tillräckligt bra,
sqrt' 12345678900987654321 = (3.5136418513220444e9,36)
sqrt 12345678900987654321 = 3.513641828785008e9
de första 7 decimalerna är lika och det tar inte någon av mig uppfattbar tid.

Men...

När jag kör
sqrt' 1000000000000 så kommer svaret, (1000000.0000000054,25), direkt, men när jag kör
sqrt' 10000000000000 så kommer svaret aldrig (fast jag har som mest väntat någon minut).

Vad kan det bero på?

GHC 7.10.3. Samma sak oavsett om jag kompilerar eller kör i GHCI.

Villkoret < 0.00...1 är kanske för strängt? Har du testat att printa var 1000 loop eller något? För att se vad som händer?
Testa mindre noggrannhet.

Citera

2017-05-08, 22:16 #3

Medlem

Reg: Feb 2011

Inlägg: 3 618

Prova med double istället för float.

Citera

2017-05-08, 22:19 #4

Moderator

Reg: Jul 2010

Inlägg: 12 961

Citat:

Ursprungligen postat av bevarass

Villkoret < 0.00...1 är kanske för strängt? Har du testat att printa var 1000 loop eller något? För att se vad som händer?
Testa mindre noggrannhet.

Sträng får man tydligen inte vara, inget händer om jag plockar bort en nolla eller två...
Eller tre.
Men ändrar jag värdet till 0,01 så funkar det. 0,001 funkar inte.

Citera

2017-05-08, 22:26 #5

Moderator

Reg: Jul 2010

Inlägg: 12 961

Citat:

Ursprungligen postat av Motherofgod

Prova med double istället för float.

Jag använder inte Float, jag använder Floating (typeclass). Både Float och Double är Floating. Men även om jag slänger in
sqrt' (10000000000000 :: Double)
så hänger den sig.

Citera

2017-05-08, 22:27 #6

Medlem

Reg: Feb 2011

Inlägg: 3 618

Citat:

Ursprungligen postat av Neksnor

Jag använder inte Float, jag använder Floating (typeclass). Både Float och Double är Floating. Men även om jag slänger in
sqrt' (10000000000000 :: Double)
så hänger den sig.

Ah, sry tänkte fel.

Citera

2017-05-08, 22:43 #7

Medlem

Reg: Nov 2015

Inlägg: 1 237

Citat:

Ursprungligen postat av Neksnor

Sträng får man tydligen inte vara, inget händer om jag plockar bort en nolla eller två...
Eller tre.
Men ändrar jag värdet till 0,01 så funkar det. 0,001 funkar inte.

Och Haskells egen funktion? Den har inte samma problem?
Finns det inget smart tricks för rötter man kan använda eller det är iteration som gäller? DU använder ju iofs en beräkning i beräkningen.

Citera

2017-05-09, 01:22 #8

Moderator

Reg: Jul 2010

Inlägg: 12 961

Citat:

Ursprungligen postat av bevarass

Och Haskells egen funktion? Den har inte samma problem?
Finns det inget smart tricks för rötter man kan använda eller det är iteration som gäller? DU använder ju iofs en beräkning i beräkningen.

Finns olika algoritmer, jag länkade lite i trådstarten. Tyckte att den algoritmen så lämplig ut för implementering i Haskell. Har ingen aning om vilken algoritm som används i Haskells egna funktion, sqrt, men den implementeringen verkar inte använda operatorn < eftersom funktionen, till skillnad från min, inte kräver att a är en Ord. En annan lustig detalj är att min funktion funkar med det avsevärt större talet 12345678900987654321. Kan det vara någon flyttalsrepresentation som spökar?

Citera

2017-05-09, 06:39 #9

Moderator

Reg: Jul 2010

Inlägg: 12 961

Roade mig med att köra beräkningen manuellt, alltså ett steg i taget, kopiera, klistra in...
Efter 27 eller 28 varv kom min kod fram till att kvadratroten ur 10000000000000 = 3162277.6601683795 vilket råkar vara samma svar som den inbyggda funktionen sqrt ger.
Men min kod kvarderar svaret i en kontroll

Kod:

      | (((s+a/s)/2)**2-a) < 0.000000000001 = (s,n+1)

och 3162277.6601683795^2 - 10000000000000 = 1.953125e-3 ≈0.002
vilket är mer än 0.000000000001.
Då går samma tal in i algoritmen igen, och igen...
La in en ny guard som "nödstopp", bör ge tillräckligt bra resultat för duktigt stora tal:

Kod:

sqrt' :: (Floating a, Ord a) => a -> (a,Int)
sqrt' a = f a (1,0)
  where
    f a (s,n)
      | (((s+a/s)/2)**2-a) < 0.000000000001 = (s,n+1)
      | n > 100                                              = (s,n+1)
      | otherwise                                            = f a (((s+a/s)/2),n+1)

sqrt' 10000000000000 = (3162277.6601683795,102)

Citera

2017-05-10, 03:37 #10

Medlem

Reg: Feb 2011

Inlägg: 14 549

Citat:

Ursprungligen postat av Neksnor

Roade mig med att köra beräkningen manuellt, alltså ett steg i taget, kopiera, klistra in...
Efter 27 eller 28 varv kom min kod fram till att kvadratroten ur 10000000000000 = 3162277.6601683795 vilket råkar vara samma svar som den inbyggda funktionen sqrt ger.
Men min kod kvarderar svaret i en kontroll

Kod:

      | (((s+a/s)/2)**2-a) < 0.000000000001 = (s,n+1)

och 3162277.6601683795^2 - 10000000000000 = 1.953125e-3 ≈0.002
vilket är mer än 0.000000000001.
Då går samma tal in i algoritmen igen, och igen...
La in en ny guard som "nödstopp", bör ge tillräckligt bra resultat för duktigt stora tal:

Kod:

sqrt' :: (Floating a, Ord a) => a -> (a,Int)
sqrt' a = f a (1,0)
  where
    f a (s,n)
      | (((s+a/s)/2)**2-a) < 0.000000000001 = (s,n+1)
      | n > 100                                              = (s,n+1)
      | otherwise                                            = f a (((s+a/s)/2),n+1)

sqrt' 10000000000000 = (3162277.6601683795,102)

Är inte felet att du försöker avbryta när det absoluta felet är tillräckligt litet?

Om du t.ex. beräknar roten ur 1,0E+101 så är 3,16227766016838E+50 ett bra resultat.
Om felet är mindre än 1,0E+37 så är det helt OK, men att få ner felet till t.ex. 1 i det fallet går inte med vanlig "datoraritmetik" motsvarande float eller double.

Du borde IMO ange ett fel som relaterar till resultatet i stället.

Kanske först beräkna ett närmevärde och sedan fortsätta med en miljardel (t.ex.) av närmevärdet som felgräns?

Enklare kanske att bara ta argumentet dela exponenten med två och sedan resultatet gånger vad du vill ha som "delta" ("relativt fel"), kanske 1,0E-09, t.ex.?

–

__________________
Senast redigerad av lasternassumma 2017-05-10 kl. 03:43.

Citera

2017-05-10, 14:24 #11

Moderator

Reg: Jul 2010

Inlägg: 12 961

Citat:

Ursprungligen postat av lasternassumma

Är inte felet att du försöker avbryta när det absoluta felet är tillräckligt litet?

Om du t.ex. beräknar roten ur 1,0E+101 så är 3,16227766016838E+50 ett bra resultat.
Om felet är mindre än 1,0E+37 så är det helt OK, men att få ner felet till t.ex. 1 i det fallet går inte med vanlig "datoraritmetik" motsvarande float eller double.

Du borde IMO ange ett fel som relaterar till resultatet i stället.

Kanske först beräkna ett närmevärde och sedan fortsätta med en miljardel (t.ex.) av närmevärdet som felgräns?

Enklare kanske att bara ta argumentet dela exponenten med två och sedan resultatet gånger vad du vill ha som "delta" ("relativt fel"), kanske 1,0E-09, t.ex.?

–

Det är väl en del i problemet, men den stora grejen var att med 1.0e13 så hamnade man i ett läge där algoritmen plötsligt returnerar samma tal som den matas med, trots att svaret är ganska långt från korrekt. Talet 1.0e101 stannar automatiskt efter 173 varv.
Man skulle väl kunna ändra den översta guarden till

Kod:

| (((s+a/s)/2)**2-a) < a*0.000000000001 = (s,n+1)

då blir felet beroende på talets storlek.

Gjorde en lite kodsnutt om någon vill se vad resultatet blir efter n iterationer.

Kod:

f :: Floating a => a -> Int -> a
f a n = f' a 1 n
  where
    f' a s n
      | n == 1    = ((s+a/s)/2)
      | otherwise = f' a ((s+a/s)/2) (n-1)

f är ett jävligt bra funktionsnamn.
Här kan man se att 1.0e101 inte ändras mer efter 173 iterationer.

Citera

2017-05-10, 15:02 #12

Medlem

Reg: Feb 2011

Inlägg: 14 549

Kanske blir det lite off-topic att fundera för mycket över algoritmen, men...
Talet som vi ska beräkna roten ur kan vi skriva som mantissa x bas ^ exponent
Antag t.ex. att basen internt är 2. Då är mantissan mellan 0,5 och 1, utom då argumentet är noll.
Justera mantissan, vid behov, så att exponenten blir jämn (öka ev. med ett). Nu har vi en mantissa som är mellan 0,25 och 1 att beräkna roten ur. Exponenten är bara att dela med två.
Roten ur [0,25-1[ finns i intervallet [0,5-1[
Men ett startvärde mitt i det intevallet konvergerar det nog snabbt.

I Haskell? Varför inte?

Citat:

Ursprungligen postat av Neksnor

Det är väl en del i problemet, men den stora grejen var att med 1.0e13 så hamnade man i ett läge där algoritmen plötsligt returnerar samma tal som den matas med, trots att svaret är ganska långt från korrekt. Talet 1.0e101 stannar automatiskt efter 173 varv.
Man skulle väl kunna ändra den översta guarden till

Kod:

| (((s+a/s)/2)**2-a) < a*0.000000000001 = (s,n+1)

då blir felet beroende på talets storlek.

Gjorde en lite kodsnutt om någon vill se vad resultatet blir efter n iterationer.

Kod:

f :: Floating a => a -> Int -> a
f a n = f' a 1 n
  where
    f' a s n
      | n == 1    = ((s+a/s)/2)
      | otherwise = f' a ((s+a/s)/2) (n-1)

f är ett jävligt bra funktionsnamn.
Här kan man se att 1.0e101 inte ändras mer efter 173 iterationer.

Citera

Haskell, plötslig prestandapuckel vid större beräkning.

Stöd Flashback