Undersökning av primtal

Imorse steg jag upp som vanligt, duschade, gjorde kaffe, och slog mig ner i min favoritfåtölj för lite morgonläsning. En artikel, om primtal, i facebook-flödet fångade min uppmärksamhet. Primtal har lite geek-status, primtal är tacksamma, alla människor begriper vad ett primtal är. Alla tänker på primtal. Och för de som gillar mer abstrakt tänkande erbjuder primtalsöknarna en lagom utmaning att gymnastisera lilla hjärnan med. Jag tror de icke-matematiskt lagda vädrar morgonluft, en reva i matematikens rationalitet. Det finns ännu ingen primtalslag. Det tror jag är kärnan i primtalens lockelse.

Nåja, artikeln handlade om en bristande slumpegenskap hos primtalen. Om man har ett primtal som slutar på 1, t ex, hur ser då fördelningen av slutsiffror ut i det följande primtalet? Om primtalet slutar på 1, hur ofta slutar nästa primtal på 1, 3, 7 eller 9? Man skulle gissa att det borde vara ungefär lika ofta en 1:a, 3:a, 7:a eller 9:a. Men, det är det inte, säger artikeln.

Eftersom det där är en otroligt enkel sak att räkna på själv, gjorde jag det. Såklart. Det går inte att låta bli. Nedan min lilla matlab-snutt för att räkna på fördelningen efter en 1:a. Ursäkta den klumpiga programmeringen, jag orkade inte tänka ut en mer slimmad kod. Dessutom är jag inte programmerare heller. I alla fall. Jag fick, ungefär, samma fördelning som artikeln sa. 18 % ettor, 30 % treor, 31 % sjuor och 21 % nior. På 100 miljoner tal. Man kan lägga till motsvarande rader för att kolla fördelningen efter en trea, sjua och nia. Jag gjorde det, den blev inte lika som för en etta, men den blev inte heller jämn.

Ja, det var väl kul?

—-
function [Nd,j1] = findtheprimes(ANTAL)
V = primes(ANTAL);
Vr = rem(V,10);
N = hist(Vr,9);
Nd = N([1 3 7 9]);
L = length(Vr);
j1 = [0 0 0 0];
for i = 1:L-1
   if Vr(i) == 1
      if Vr(i+1) == 1
         j1(1) = j1(1)+1;
      elseif Vr(i+1) == 3
         j1(2) = j1(2)+1;
      elseif Vr(i+1) == 7
         j1(3) = j1(3)+1;
      elseif Vr(i+1) == 9
         j1(4) = j1(4)+1;
      end
   end
end

—-

Ett litet tillägg: Doktor T tyckte att jag på kvinnligt vis nedvärderade och bad om ursäkt för mina programmerings-skills. Det tycker inte jag. Och till mitt försvar vill jag dessutom anföra att koden faktiskt är okommenterad. Som jag brukar göra det. Wild and Crazy.

En orgie i statistik

Antag att man i ett anfall av ostyrbar nyfikenhet får för sig att man ska väga alla godisbitar i sin fredagsgodispåse, för att man t ex har en ny våg, och för att man bara måste få veta. Man väger bitarna och skriver upp resultatet, bit för bit, i en tabell. Och så sitter man där med en massa siffror.

Då tänker man förstås att man skulle vilja veta lite mer om egenskapen hos mängden godis. T ex vill man kanske veta vad medelgodisen väger, och vilken spridning det är på godisarna, dvs vilken standardavvikelse man har i godispåsen. Och den absolut brännande frågan är förstås är godismassor normalfördelade, eller kommer de från någon annan fördelning? Exponential? Weibull? Rektangelfördelning?

Det funderade jag på idag. Att hitta fördelningen för en datamängd. Frågan kom nämligen upp på jobbet. Man kan ju göra en hypotestest och kolla om man har en normalfördelning t ex. Men sådär lite ingenjörsmässigt göra en anpassning, fit, till data. Jag har ingen funktion för det. Hade inte.

Nu använde jag bara godisexemplet för att anknyta till en vardagstillämpning, istället för abstrakta datapunkter från god knows where. En vardagstillämpningar som mina barn väl känner till, eftersom deras fredagsgodis får väga 400 gram, med en avvikelse uppåt på 10 %. Som i ett svagt ögonblick förvandlades till 10% + 6g. Deras fredagsgodis väger alltså 446 gram. Exakt. De vet vad alla godissorter väger. De plockar och anpassar. Vår fredagsgodisshopping är en orgie i matematik, eller snarare i statistik.

Hursomhelst hittade jag funktionen allfitdist på matlab central. Jag testade den på enkelt vis.

for i = 1:1000
E(i) = random(‘exp’,2);
N(i) = random(‘norm’,2,1);
end

[DE,PDE]=allfitdist(E,’PDF’);
[DN,PDN]=allfitdist(N,’PDF’);

normal

För exponentialfördelningen, E, fick jag tillbaka förslaget Inverse Gaussian. Inte helt nöjd med det.

För normalfördelningen fick jag förslaget Normalfördelning. Med rätt medelvärde och standardavvikelse. Bra.

I badrumsskåpet hittade jag de sista resterna av ett badskum från Lush. Frågan är nu om jag ska bada skumbad, eller om jag ska undersöka allfitdist lite noggrannare?

lush

Och det får väl illustrera hur matematiken kommer in i vardagen på ett sätt som man knappt lägger märke till.

The magic of numbers. Eller siffermagi. Men det låter inte lika coolt.

Idag hade jag en deadline. Inte superskarp kanske, men jag hade lovat att ha en text klar till klockan två. När mina amerikanska kollegor, i alla fall de på östkusten, startade dagen på kontoret. Och det man lovar vill man hålla. Jag hade bestämt att från klockan elva skulle jag vara tvungen att fokusera på texten för att få den klar. Det var sista iterationen, jag har hållit på några dagar från och till. Det är en text som beskriver en teknisk lösning som vi jobbar på.

Klockan elva släppte jag fokus på de analyser jag höll på med, och började skriva det sista.

Tjugo över elva fick jag för mig att jag skulle kolla om 1713 är ett primtal. Jag fick en rent oändlig lust. Så jag googlade. Hittade snabbt en tabell. 1713 är inget primtal. Då undrar man förstås vilka tal det egentligen är delbart med. Eller hur? Man måste få veta.

Matlab har en funktion, isprime(X), för att kolla om ett tal är ett primtal. Men man får bara svar True eller False (1 eller 0). Men jag ville veta faktorerna i talet. Jag tänkte att jag skriver en liten snutt i matlab som kollar det där. Så det gjorde jag. Den här snutten:

function J = isprime(P);
%
% find out if P is a prime number
% J = isprime(P)
%
K = zeros(1,P);
for i = 1:P
if rem(P,i) ~0;
K(i) = 1;
end
end
I = find(K==0);
if length(I)>2
J = I(2:length(I)-1);
else
disp([num2str(P) ‘ is a prime’])
J = 0;
return
end

1713 är delbart med 3 och 571.

Jag hann klart med min text i tid.

Tvångsprogrammering

Det är väl klart att jag läst programmering. Och programmerat. Grundligt och mycket. Fortran, Pascal och ve och fasa, Assembler (kursboken, den om M 68000 var utmärkt att lägga under spjälsängens huvudända när barnen var förkylda). Hur kul är det? Inte kul alls. Numeriska metoder, lösa matriser och indexering till dödagar. Ordo av första och andra ordningen. Matlab använder jag varje dag, har gjort i 25 år sådär. Shell-script och awk måste jag, motvilligt, ge mig i kast med ibland för att fixa smidig hantering av stora mängder beräkningsdata. Men jag tycker egentligen inte att själva programmeringen är kul. Bara praktisk.

Men när jag var mammaledig andra gången, med tvillingarna, var jag så svältfödd på intellektuell stimulans, då när jag vaknade ur amningsdimman, att till och med programmering kändes som en rolig aktivitet. Så jag lärde mig lite databashantering i MySQL och en del php. Inte världens mesta, men jag skrev några tillämpningar i alla fall. Installerade Apache-server på datorn hemma och stod i.

I vilket fall som helst, har all denna erfarenhet nu lett fram till kvällens bragd: Jag har lagt till en liten rad i mitt WP-tema, så att taggarna syns under posterna. Jag är mycket nöjd.