So a little while ago, I was trying to write a nice, clear version of the language shootout n-body problem, using the ST monad, and STUArrays. Though I ran into a snag; I was getting segmentation faults, for very small values of n (< 23), and not knowing enough about such things, I have no idea where to begin in diagnosing the problem. So, I shall post my code here, and hope that someone might have some ideas. I’ve talked to Don Stewart about it, but he hasn’t got back to me yet. (All the MUArray instance stuff is from Don’s work, which I nicked) So, here we have it:
> {-# OPTIONS -O2 -funbox-strict-fields #-} > {-# LANGUAGE BangPatterns #-} > {-# LANGUAGE UnboxedTuples #-} > {-# LANGUAGE MagicHash #-} > {-# LANGUAGE MultiParamTypeClasses #-} > > module Main where > > import System.Environment > import Data.Array.Base > import GHC.Exts > import GHC.ST > import GHC.Arr > import Control.Monad > import Control.Monad.ST > import Data.List > import Foreign.Storable > import Text.Printf > import Data.STRef > > > data Vec3 = V !Double !Double !Double deriving (Show) > data Body = B { > pos :: {-# UNPACK #-} !Vec3, > vel :: {-# UNPACK #-} !Vec3, > mass :: !Double > } deriving (Show) > > > main = do > n <- getArgs >>= \xs -> if null xs then return 1000 else (readIO.head) xs > let (e,e') = runST $ do > ps <- planets > offset ps > e1 <- energy ps > advn n ps > e2 <- energy ps > return (e1,e2) > prtnum e > prtnum e' > > > offset ps = do > v <- newSTRef (V 0 0 0) > (B p ve m) <- readArray ps 0 > forM_ [0..numPs] $ \i -> do > p <- readArray ps i > modifySTRef v (.+. momentum p) > modifySTRef v (./ solar_mass) > v' <- readSTRef v > writeArray ps 0 (B p (ve.-. v') m) > > > > advn 0 ps = do > return ps > advn n ps = do > advance deltat ps > -- e <- energy ps > -- unsafeIOToST (prtnum e) > advn (n-1) ps > > advance dt ps = do > forM_ [0..numPs] $ \i -> do > (B p v m) <- readArray ps i > vvar <- newSTRef v > forM_ [i+1..numPs] $ \j -> do > (B p' v' m') <- readArray ps j > let dp = p .-. p' > dst = mag dp > magnit = dt / (dst dst dst) > change = magnit . dp > sca1 = m' . change > sca2 = m . change > -- writeArray ps I (B p (v .-. sca1) m) > modifySTRef vvar (.-. sca1) > writeArray ps j (B p' (v' .+. sca2) m') > vvar' <- readSTRef vvar > writeArray ps I (move dt (B p vvar' m)) > > > energy ps = do > e <- newSTRef (0 :: Double) > > forM_ [0..numPs] $ \i -> do > b1 <- readArray ps i > modifySTRef e ((+) (kinetic b)1) > > forM_ [i+1..numPs] $ \j -> do > b2 <- readArray ps j > let x = potential b1 b2 > modifySTRef e (subtract x) > e' <- readSTRef e > return e' > > > > move dt (B p v m) = B (p .+. dt . v) v m > > > > > ----------------------- > solar_mass, deltat, days_per_year :: Double > days_per_year = 365.24 > dp = days_per_year > solar_mass = 4 pi pi > deltat = 0.01 > numPs = length planets' - 1 > > planets'' = newArray_ (0,numPs) :: ST s (STUArray s Int Body) > -- (B (V 0 0 0) (V 0 0 0) 1) > > planets = do > arr <- planets'' > forM_ (zip [0..] planets') $ \(i,p) -> do > writeArray arr I p > return arr > > planets' = > [B (V 0.0 > 0.0 > 0.0) > (V 0.0 > 0.0 > 0.0) > solar_mass, > > B (V 4.84143144246472090e+00 > (-1.16032004402742839e+00) > (-1.03622044471123109e-01)) > (V (1.66007664274403694e-03dp) > (7.69901118419740425e-03dp) > ((-6.90460016972063023e-05)dp)) > (9.54791938424326609e-04 solar_mass), > > B (V 8.34336671824457987e+00 > 4.12479856412430479e+00 > (-4.03523417114321381e-01)) > (V ((-2.76742510726862411e-03)dp) > (4.99852801234917238e-03dp) > (2.30417297573763929e-05dp)) > (2.85885980666130812e-04 solar_mass), > > B (V 1.28943695621391310e+01 > (-1.51111514016986312e+01) > (-2.23307578892655734e-01)) > (V (2.96460137564761618e-03dp) > (2.37847173959480950e-03dp) > ((-2.96589568540237556e-05)dp)) > (4.36624404335156298e-05 solar_mass), > > B (V 1.53796971148509165e+01 > (-2.59193146099879641e+01) > 1.79258772950371181e-01) > (V (2.68067772490389322e-03dp) > (1.62824170038242295e-03dp) > ((-9.51592254519715870e-05)dp)) > (5.15138902046611451e-05 solar_mass) > ] > > > ----------------------- > > infix 4 .-. > infix 4 .+. > -- infix 5 .. > infixr 5 . > infixl 5 ./ > > {-# INLINE (.+.) #-} > {-# INLINE (.-.) #-} > -- {-# INLINE (..) #-} > V x y z .+. V x' y' z' = V (x+x') (y+y') (z+z') > V x y z .-. V x' y' z' = V (x-x') (y-y') (z-z') > -- V x y z .. V x' y' z' = V (xx') (yy') (zz') > > {-# INLINE (.) #-} > n . (V x y z) = V (nx) (ny) (nz) > {-# INLINE (./) #-} > V x y z ./ n = V (x/n) (y/n) (z/n) > > {-# INLINE momentum #-} > momentum :: Body -> Vec3 > momentum (B _ v m) = m . v > > {-# INLINE kinetic #-} > kinetic x = 0.5 mass x (mag2.vel) x > > > potential a b = > let dst = dist' a b > in if dst <= 0 > then 0 > else (mass a mass b) / dist' a b > > {-# INLINE mag #-} > mag :: Vec3 -> Double > mag = sqrt . mag2 > > {-# INLINE mag2 #-} > mag2 (V x y z) = xx + yy + zz > > {-# INLINE dist #-} > dist v1 v2 = sqrt $ dist2 v1 v2 > > {-# INLINE dist2 #-} > dist2 v1 v2 = mag2 $ v1 .-. v2 > > {-# INLINE dist' #-} > dist' b1 b2 = dist (pos b1) (pos b2) > > prtnum x = putStrLn $ printf "%.9f" x > > -------------------- > bodyScale :: Int# -> Int# > bodyScale n = (scale # 7#) *# n > where I# scale = sizeOf (undefined :: Double) > > instance MArray (STUArray s) Body (ST s) where > > {-# INLINE getBounds #-} > getBounds (STUArray l h _ _) = return (l,h) > > {-# INLINE getNumElements #-} > getNumElements (STUArray _ _ n _) = return n > > {-# INLINE unsafeNewArray_ #-} > unsafeNewArray_ (l,u) = unsafeNewArraySTUArray_ (l,u) bodyScale -- bytes > > {-# INLINE newArray_ #-} > newArray_ bounds = newArray bounds (B (V 0 0 0) (V 0 0 0) 1) > > > {-# INLINE unsafeRead #-} > unsafeRead (STUArray _ _ _ arr) (I# i) = ST $ \s1 -> > let j = bodyScale i in > case readDoubleArray# arr j s1 of { (# s2, xd #) -> > case readDoubleArray# arr (j +# 1#) s2 of { (# s3, yd #) -> > case readDoubleArray# arr (j +# 2#) s3 of { (# s4, zd #) -> > case readDoubleArray# arr (j +# 3#) s4 of { (# s5, vxd #) -> > case readDoubleArray# arr (j +# 4#) s5 of { (# s6, vyd #) -> > case readDoubleArray# arr (j +# 5#) s6 of { (# s7, vzd #) -> > case readDoubleArray# arr (j +# 6#) s7 of { (# s8, massd #) -> > (# s8, B { pos = V (D# xd) (D# yd) (D# zd) > , vel = V (D# vxd) (D# vyd) (D# vzd) > , mass = D# massd } #) }}}}}}} > > {-# INLINE unsafeWrite #-} > unsafeWrite (STUArray _ _ _ arr) (I# i) > (B (V (D# xd) (D# yd) (D# zd)) (V (D# vxd) (D# vyd) (D# vzd)) (D# massd)) > = ST $ \s1 -> > let j = bodyScale i in > case writeDoubleArray# arr j xd s1 of {s2 -> > case writeDoubleArray# arr (j +# 1#) yd s2 of {s3 -> > case writeDoubleArray# arr (j +# 2#) zd s3 of {s4 -> > case writeDoubleArray# arr (j +# 3#) vxd s4 of {s5 -> > case writeDoubleArray# arr (j +# 4#) vyd s5 of {s6 -> > case writeDoubleArray# arr (j +# 5#) vzd s6 of {s7 -> > case writeDoubleArray# arr (j +# 6#) massd s7 of {s8 -> > (# s8, () #) }}}}}}} >
I’ve just read Marcus S. Zarra’s post about software licencing, piracy and stupidity (to paraphrase), and I think he’s hit the nail right on the head.
Treating potential customers like potential thieves makes them feel more inclined to follow these feeling. I’ll admit that I’ve pirated software, but I’ve also bought a lot of it too. The stuff I’ve pirated is uaually something that I feel is far overpriced for what they offer. Like Marcus says, software is something you make once, and sell thousands of times; it’s making money from nothing. The things I’ve bought are usually pretty cheap, and that makes me want to buy them, I want to help the developer and reward them for making something which is useful to me, and reasonably priced, or even a bargain.
I do not like being forced to buy your software. I do not like being forced to buy your software to use it for something that is useful, but not exactly (or at all) worth paying $30 for. Nor do I like being told that I have not bought your software. Give me the full thing to try, without restrictions, at all, give me a month and see if I feel I need it. If after that month, (and that’s a month of use, not a month from the date of download), I feel I need your software, then I will buy it. If I never hit that month, then no, I will not be buying it, it has not been useful enough for me to consider it. Doesn’t mean it’s not great, or that it’s buggy and annoying (I will tell you if it is), it just means I have no ude for it.
I’ll take an example, VMware Fusion. Now, this is made by a very large company, a big, corporate business, that makes expensive software for some rather important tasks. But what happened when Fusion v2.0 came out? they gave it to us for free! I, and I’m sure many others would have paid for this update, they did a lot of work for it, added some very nice new features. But instead of them saying “hey! they’re getting new stuff, they HAVE to pay!”, they seemed to take the stance that these new features were missing features, and that v2.0 was heading towards making the program more complete, rather than expanding it.
This is an excellent example of how I feel software licencing should work. Be nice to your customers, give them treats, and they will do more than stay quiet about your products, they will tell their friends, or anyone they know who might be interested. They feel an affinity with you and your software, and will go the point of arguing and fighting over which is better, your product, or someone else’s. Loyalty will always bring you more than fear, and being nice will bight you in the arse far less often than being a dick.
I don’t think I have summed this up as well as Marcus managed to (I’ve just written this from start to finish, and can’t be bothered to re-read it [pan's labyrinth awaits]), but I did feel like showing my support for his views, because he got how I feel so right.
Anyway, more later
~ Axman
So, we’ve just finished setting up Hendo’s new bacon blog. Here’s what she has to say about it:
Me: So do you like you new blog? Hendo: Yup, it’s pretty! Me: Yes it is… -_-
Yep, she sure is cute. Anyway, hopefully she’ll keep up with it (and so will I). So check it out!
Sp for my birthday last November, my girlfriend bought me a very lovely new pen. It is a lovely Acme pen, the Blueprint one on that page (my studying Engineering/IT ‘n all). For christmas, she bought me the fountain pen nib replacement, and far out, I love it! I’ve never used a fountain pen before, but they are just wonderful. I can’t say what it is about it that I like, but it feels so much more natural than a ballpoint.
This, added with the shear weight of the thing makes it the perfect writing tool, it’s quite a leap from my normal Lamp Spirit pencil, which I still love also.
now i’m not one for writing my hand much, my hand writing is pathetically unreadable most of the time, but I do enjoy fine stationary, and this seems like it could be something that turns into somewhat of an addiction…
Welcome to my new blog powered by Movable Type. This is the first post on my blog and was created for me automatically when I finished the installation process. But that is ok, because I will soon be creating posts of my own!
So last night, my girlfriend’s dad gave me an old XP machine they didn’t have a need for. Turns out it has fairly decent specs, so I got to work installing FreeBSD on it. In the next few days, i shall be turning it into my new web server, and possibly putting my blog back on it (although it’ll be slower to serve it there, i can at least install plugins, which I’m still pissed off about with this blog). In the mean time though, i have been installing pfSense on the machine that was my web server. the specs of the old one are as follows:
and the new one:
so basically it’s a much faster machine (though now I come to think of it, maybe swapping the disks wouldn’t be a bad idea). So what you say? Well, what I wanted to tell you is that pfSense is awesomely easy to use, it’s wonderful. the webGUI is well laid out and works very well, installing new packages like nmap, spamd, bandwidthd, OpenBGPD etc. is as easy as visiting the web interface and clicking a button. You get some nice graphs about most aspects of system performance (though they haven’t been working for me, not sure why, getting no data).
The other thing I was going to mention is how nice Seagate’s warranty returns site is. Our 250GB drive died recently (not sure when, it was hooked up to the web server and i never used it), so I checked to see if it was still under warranty. Just sick in the model number and serial number, and they’ll tell you exactly when the warranty ends. If you need to return it, you can go through two or three pages, choose the nearest point of contact for their repairs, and then print out a label with the address and your return address to stick on your parcel. They even give you wrapping instructions which is nice (though they seem to insist that you use a “SeaShell” plastic case to put the drive in, which of course I don’t own).
So, in the next few days, hopefully I’ll have a nice new web server, running my blog and being awesome!
So i was talking to my girlfriend today about my blog, after I showed it to her last night. She was not pleased. I had neglected to mention her at all, and she did not appreciate it. So now I am making an effort to make it up to my lovely, beautiful, amazing girlfriend of a year and a month, Sarah. I love you Sarah! In other news, Happy Obama Day. About time the US did something right.
I assert though that the GPL does not give us freedom, merely different restrictions. I even assert that it does not give us all of the freedoms listed above. In this post, I will argue that the GPL is inherently flawed when it comes to freeing software.
My good IRC friend Beelsebob decided independently to start a new blog, on the exact same say I did, but hopefully I Beelsebob’d him this time!* Anyway, he’s made a better start than I have by writing a nice article on why the GPL is in fact not free, but just restricting you in a different way. So go ahead and read the article
Ok, so my plans to get things like markdown working have been shot. I had no idea Wordpress.com was so lame, I thought they might at least allow you to add some plugins, but nooooo. oh well, I’ll make do.
So, I’m wanting to put more examples on the Monad/ST page on the Haskell wiki, but I’m not sure what else I can do with it that’s simple. I’ve been using Haskell for too long to remember what sort of problems I’d write statefully! So, if you have any ideas, either let me know, or just post them on the wiki! My first example is:
import Control.Monad.ST import Data.STRef import Control.Monad sumST :: Num a => [a] -> a sumST xs = runST $ do -- runST takes out stateful code and makes it pure again. n <- newSTRef 0 -- Create an STRef (place in memory to store values) forM_ xs $ x -> do -- For each element of xs .. -- (forM_ is just flip mapM_) modifySTRef n (+x) -- add it to what we have in n. readSTRef n -- read the value of n, and return it.
which is a version of the sum function, that uses the ST monad to update a variable in place, but also run in a purely functional way (notice the lack of ST in the type signature) And foldl and Fibonacci using ST:
foldlST :: (a -> b -> a) -> a -> [b] -> a foldlST f acc xs = runST $ do acc' <- newSTRef acc -- Create a variable for the accumulator forM_ xs $ x -> do -- For each x in xs... a <- readSTRef acc' -- read the accumulator writeSTRef acc' (f a x) -- apply f to the accumulator and x readSTRef acc' -- and finally read the result fibST :: Integer -> Integer fibST n = if n < 2 then n else runST $ do x <- newSTRef 0 y <- newSTRef 1 fibST' n x y where fibST' 0 x _ = readSTRef x fibST' n x y = do x' <- readSTRef x y' <- readSTRef y writeSTRef x y' writeSTRef y (x'+y') fibST' (n-1) x y
*Output created by HsColor
