uniform sampling av en cirkel

Jag har en cirkel med centrum i origon och radien r, och jag vill ha slumpmässiga punkter som är jämnt fördelade över denna cirkel. Först tänkte jag att jag bara kan slumpa fram ett tal φ = [0, 2pi] och ett annat θ = [0, r] och sedan ta punkten θ(cos(φ), sin(φ)). Efter lite funderande kom jag på att jag kommer väl få mycket fler punkter nära cirkelns mitt än vad jag får i dess utkant? Tex kommer hälften ligga innanför cirkeln med radien r/2, även om denna bara utgör en fjärdedel av den totala arean. Ett sätt att lösa det på är att slumpa två värden α och β i intervallen [-r, r], och sedan kasta bort resultatet och göra det en gång till om α^2 + β^2 > r^2. Sannolikheten att det händer borde vara pi/4 (arean av cirkeln delat med arean av kvadraten), vilket gör det inte helt osannolikt att ganska ofta måste köras 2-4 gånger innan den hittar ett bra värde. Värt att nämna är att den här punkten ska slumpas fram 1024*1024*300 = 314 572 800 gånger, så det ska gärna gå snabbt.

Fler småfrågor: hur hur tung är sin och cos-beräkningar av en dator? Kan man anta att det görs någon polynomutveckling på 5-10 termer och alltså innebär 5-10 multiplikationer och additioner? Hur förhåller det sig till antar klockcykler för att generera ett slumpmässigt tal? Är eventuellt den andra metoden snabbare trots att den ibland måste köra flera gånger?
Om sannolikheten att den körs n gånger är 1 - (pi/4)^n, för n >= 1, hur många gånger kör den i genomsnitt (jag borde egentligen kunna räkna ut det här själv, men jag har lite tidspress och mycket annat att göra). tacksam för svar.

Om din slumptalsgenerering inte är exceptionellt långsam bör det inte vara något som helst problem att dra (x, y) uniformt och loopa tills de ligger inom cirkeln.

Exempel: med taus2-generatorn från GSL spottar min dator ut 1e9 slumptal på drygt 5 sekunder. Att generera 300*1024**2 koordinater i en cirkel så som du vill tar 7.1 sekunder.

Tack för svaret. Om det tex skulle krävas i genomsnitt fyra slumptal i stället för två skulle det alltså utgöra cirka 7 sekunder av renderingstiden, vilket är väldigt lite med tanke på att den totala tiden är på över en timme.
Tar gärna emot svar på de andra frågorna också, även om loop-metoden förmodligen ligger hyfsat nära det optimala.

Du skulle kunna: slumpa t likformigt i [0, r²], låta ρ = sqrt(t), slumpa φ likformigt i [0, 2π] och ta punkten (ρ cos φ, ρ sin φ). Detta ger dig en likformig fördelning på cirkelskivan. Jag tror inte detta går avsevärt snabbare än att slumpa punkter i kvadraten och slänga svaren som hamnar utanför, men du kanske skulle kunna testa empiriskt vad som är snabbast.

Sagt och gjort:
Output:
Square sampling... 1 billion samples completed in 16172 ms
Circle sampling... 1 billion samples completed in 30297 ms
sine/cosine sampling... 1 billion samples completed in 12094 ms

Den första metoden är alltså dubbelt så snabbt. Den tredje är ännu lite snabbare men tyvärr felaktig.
EDIT: glömde multiplikation samt sin och cos i andra metoden. Fick nu tiderna 16266/89922/12063 vilket väldigt mycket långsammare.

Det där ger dig alltid exakt 0. Om du nödvändigtvis måsta använda rand() (som inte garanteras vara lämplig till något, allra minst simuleringar, även om vissa moderna implementationer är sunda), så måste du göra din cast till double innan du dividerar för annars blir det ju bara en heltalsdivision. Dessutom får du inte glömma att multiplicera med 2 och subtrahera 1 om du vill täcka hela cirkeln och inte bara första kvadranten.

Att multiplicera med radien är förresten rätt onödigt - gör i stället det *efter* att du hittat koordinater som ligger innanför enhetscirkeln, så sparar du ett (egentligen ganska försumbart) antal multiplikationer.

Vad gäller implementation av trigonometriska funktioner: Sådant varierar säkerligen med processorarkitektur och kompilator/mattebibliotek, men du kan ju se om du hittar källkoden till den implementation du använder. På x86 finns instruktionerna fsin och fcos och hur många klockcykler de kan tänkas ta man du nog leta efter i dokumentationen från Intel eller AMD osv.

När de gäller slumptalsgeneratorer finns det rätt många att välja på, med rätt olika egenskaper. Några böcker som tar upp ämnet nämns på http://en.wikipedia.org/wiki/Pseudo-...mber_generator . Det lilla jag kan har jag nog främst lärt mig av Numerical Recipes (vars kod är horribel) och av GSL-dokumentationen.

Tack för tipsen, jag såg felet i att jag alltid fick noll, samt att jag samplade bara första kvadranten precis efter att jag postat inlägget, men jag struntade i att rätta det eftersom det började kännas som jag diskuterade ämnet med mig själv.
Angående sin: googlade och hittade en c-implementation med bland annat den här kommentaren: "sin(x) is approximated by a polynomial of degree 13 on [0,pi/4]" så det stämmer väl hysat med hur jag trodde den beräknades.

Tack, det var den typen av optimeringar jag var ute efter, även om den är väldigt liten.

Jag tycker det verkar väldigt vanligt att folk klagar på random-funktionerna i standardbiblioteken, är den här kritiken verkligen befogad? Även om fördelningsfunktionen inte är exakt 1 i [0,1] (vilken i princip borde vara en omöjlighet, även om man tar hänsyn till att man bara kan representera ett ändligt antal tal med 23/64 bitar) är den inte tillräcklig för ändamålet?

Förresten, en sak till angående avrundningen till noll. Jag märkte att jag det blev noll med
så jag ändrade till
vilket gjorde att exekveringen tog mycket längre tid. Jag testade sen:
och det fungerar (x och y fördelas mellan -2 och 2)

I värsta fall är rand() en linjär kongruensgenerator med kortare period än antalet slumptal man behöver, och problemet man försöker lösa dessutom känsligt för korrelerade (eller triviellt icke-oberoende) slumptal. Men visst, om man inte gör simuleringar där man samlar statistik över lång tid så spelar det knappast någon roll att slumpgeneratorn inte är lämpad för just det.

Något man dock bör tänka på (särskilt om det finns en risk att generatorn man använder är just en linjär kongruensgenerator) är att inte generera små heltal genom att slänga de höga bitarna av slumptalen. Det är särskilt olämpligt att göra något i stil med
if(rand() & 1) ...
eftersom vissa slumptalsgeneratorer gör att villkoret blir sant precis varannan gång.