diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,339 @@ + GNU GENERAL PUBLIC LICENSE + Version 2, June 1991 + + Copyright (C) 1989, 1991 Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The licenses for most software are designed to take away your +freedom to share and change it. By contrast, the GNU General Public +License is intended to guarantee your freedom to share and change free +software--to make sure the software is free for all its users. This +General Public License applies to most of the Free Software +Foundation's software and to any other program whose authors commit to +using it. (Some other Free Software Foundation software is covered by +the GNU Lesser General Public License instead.) You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +this service if you wish), that you receive source code or can get it +if you want it, that you can change the software or use pieces of it +in new free programs; and that you know you can do these things. + + To protect your rights, we need to make restrictions that forbid +anyone to deny you these rights or to ask you to surrender the rights. +These restrictions translate to certain responsibilities for you if you +distribute copies of the software, or if you modify it. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must give the recipients all the rights that +you have. You must make sure that they, too, receive or can get the +source code. And you must show them these terms so they know their +rights. + + We protect your rights with two steps: (1) copyright the software, and +(2) offer you this license which gives you legal permission to copy, +distribute and/or modify the software. + + Also, for each author's protection and ours, we want to make certain +that everyone understands that there is no warranty for this free +software. If the software is modified by someone else and passed on, we +want its recipients to know that what they have is not the original, so +that any problems introduced by others will not reflect on the original +authors' reputations. + + Finally, any free program is threatened constantly by software +patents. We wish to avoid the danger that redistributors of a free +program will individually obtain patent licenses, in effect making the +program proprietary. To prevent this, we have made it clear that any +patent must be licensed for everyone's free use or not licensed at all. + + The precise terms and conditions for copying, distribution and +modification follow. + + GNU GENERAL PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. This License applies to any program or other work which contains +a notice placed by the copyright holder saying it may be distributed +under the terms of this General Public License. The "Program", below, +refers to any such program or work, and a "work based on the Program" +means either the Program or any derivative work under copyright law: +that is to say, a work containing the Program or a portion of it, +either verbatim or with modifications and/or translated into another +language. (Hereinafter, translation is included without limitation in +the term "modification".) Each licensee is addressed as "you". + +Activities other than copying, distribution and modification are not +covered by this License; they are outside its scope. The act of +running the Program is not restricted, and the output from the Program +is covered only if its contents constitute a work based on the +Program (independent of having been made by running the Program). +Whether that is true depends on what the Program does. + + 1. You may copy and distribute verbatim copies of the Program's +source code as you receive it, in any medium, provided that you +conspicuously and appropriately publish on each copy an appropriate +copyright notice and disclaimer of warranty; keep intact all the +notices that refer to this License and to the absence of any warranty; +and give any other recipients of the Program a copy of this License +along with the Program. + +You may charge a fee for the physical act of transferring a copy, and +you may at your option offer warranty protection in exchange for a fee. + + 2. You may modify your copy or copies of the Program or any portion +of it, thus forming a work based on the Program, and copy and +distribute such modifications or work under the terms of Section 1 +above, provided that you also meet all of these conditions: + + a) You must cause the modified files to carry prominent notices + stating that you changed the files and the date of any change. + + b) You must cause any work that you distribute or publish, that in + whole or in part contains or is derived from the Program or any + part thereof, to be licensed as a whole at no charge to all third + parties under the terms of this License. + + c) If the modified program normally reads commands interactively + when run, you must cause it, when started running for such + interactive use in the most ordinary way, to print or display an + announcement including an appropriate copyright notice and a + notice that there is no warranty (or else, saying that you provide + a warranty) and that users may redistribute the program under + these conditions, and telling the user how to view a copy of this + License. (Exception: if the Program itself is interactive but + does not normally print such an announcement, your work based on + the Program is not required to print an announcement.) + +These requirements apply to the modified work as a whole. If +identifiable sections of that work are not derived from the Program, +and can be reasonably considered independent and separate works in +themselves, then this License, and its terms, do not apply to those +sections when you distribute them as separate works. But when you +distribute the same sections as part of a whole which is a work based +on the Program, the distribution of the whole must be on the terms of +this License, whose permissions for other licensees extend to the +entire whole, and thus to each and every part regardless of who wrote it. + +Thus, it is not the intent of this section to claim rights or contest +your rights to work written entirely by you; rather, the intent is to +exercise the right to control the distribution of derivative or +collective works based on the Program. + +In addition, mere aggregation of another work not based on the Program +with the Program (or with a work based on the Program) on a volume of +a storage or distribution medium does not bring the other work under +the scope of this License. + + 3. You may copy and distribute the Program (or a work based on it, +under Section 2) in object code or executable form under the terms of +Sections 1 and 2 above provided that you also do one of the following: + + a) Accompany it with the complete corresponding machine-readable + source code, which must be distributed under the terms of Sections + 1 and 2 above on a medium customarily used for software interchange; or, + + b) Accompany it with a written offer, valid for at least three + years, to give any third party, for a charge no more than your + cost of physically performing source distribution, a complete + machine-readable copy of the corresponding source code, to be + distributed under the terms of Sections 1 and 2 above on a medium + customarily used for software interchange; or, + + c) Accompany it with the information you received as to the offer + to distribute corresponding source code. (This alternative is + allowed only for noncommercial distribution and only if you + received the program in object code or executable form with such + an offer, in accord with Subsection b above.) + +The source code for a work means the preferred form of the work for +making modifications to it. For an executable work, complete source +code means all the source code for all modules it contains, plus any +associated interface definition files, plus the scripts used to +control compilation and installation of the executable. However, as a +special exception, the source code distributed need not include +anything that is normally distributed (in either source or binary +form) with the major components (compiler, kernel, and so on) of the +operating system on which the executable runs, unless that component +itself accompanies the executable. + +If distribution of executable or object code is made by offering +access to copy from a designated place, then offering equivalent +access to copy the source code from the same place counts as +distribution of the source code, even though third parties are not +compelled to copy the source along with the object code. + + 4. You may not copy, modify, sublicense, or distribute the Program +except as expressly provided under this License. Any attempt +otherwise to copy, modify, sublicense or distribute the Program is +void, and will automatically terminate your rights under this License. +However, parties who have received copies, or rights, from you under +this License will not have their licenses terminated so long as such +parties remain in full compliance. + + 5. You are not required to accept this License, since you have not +signed it. However, nothing else grants you permission to modify or +distribute the Program or its derivative works. These actions are +prohibited by law if you do not accept this License. Therefore, by +modifying or distributing the Program (or any work based on the +Program), you indicate your acceptance of this License to do so, and +all its terms and conditions for copying, distributing or modifying +the Program or works based on it. + + 6. Each time you redistribute the Program (or any work based on the +Program), the recipient automatically receives a license from the +original licensor to copy, distribute or modify the Program subject to +these terms and conditions. You may not impose any further +restrictions on the recipients' exercise of the rights granted herein. +You are not responsible for enforcing compliance by third parties to +this License. + + 7. If, as a consequence of a court judgment or allegation of patent +infringement or for any other reason (not limited to patent issues), +conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot +distribute so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you +may not distribute the Program at all. For example, if a patent +license would not permit royalty-free redistribution of the Program by +all those who receive copies directly or indirectly through you, then +the only way you could satisfy both it and this License would be to +refrain entirely from distribution of the Program. + +If any portion of this section is held invalid or unenforceable under +any particular circumstance, the balance of the section is intended to +apply and the section as a whole is intended to apply in other +circumstances. + +It is not the purpose of this section to induce you to infringe any +patents or other property right claims or to contest validity of any +such claims; this section has the sole purpose of protecting the +integrity of the free software distribution system, which is +implemented by public license practices. Many people have made +generous contributions to the wide range of software distributed +through that system in reliance on consistent application of that +system; it is up to the author/donor to decide if he or she is willing +to distribute software through any other system and a licensee cannot +impose that choice. + +This section is intended to make thoroughly clear what is believed to +be a consequence of the rest of this License. + + 8. If the distribution and/or use of the Program is restricted in +certain countries either by patents or by copyrighted interfaces, the +original copyright holder who places the Program under this License +may add an explicit geographical distribution limitation excluding +those countries, so that distribution is permitted only in or among +countries not thus excluded. In such case, this License incorporates +the limitation as if written in the body of this License. + + 9. The Free Software Foundation may publish revised and/or new versions +of the General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + +Each version is given a distinguishing version number. If the Program +specifies a version number of this License which applies to it and "any +later version", you have the option of following the terms and conditions +either of that version or of any later version published by the Free +Software Foundation. If the Program does not specify a version number of +this License, you may choose any version ever published by the Free Software +Foundation. + + 10. If you wish to incorporate parts of the Program into other free +programs whose distribution conditions are different, write to the author +to ask for permission. For software which is copyrighted by the Free +Software Foundation, write to the Free Software Foundation; we sometimes +make exceptions for this. Our decision will be guided by the two goals +of preserving the free status of all derivatives of our free software and +of promoting the sharing and reuse of software generally. + + NO WARRANTY + + 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY +FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN +OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES +PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED +OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS +TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE +PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, +REPAIR OR CORRECTION. + + 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR +REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, +INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING +OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED +TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY +YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER +PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE +POSSIBILITY OF SUCH DAMAGES. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +convey the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +Also add information on how to contact you by electronic and paper mail. + +If the program is interactive, make it output a short notice like this +when it starts in an interactive mode: + + Gnomovision version 69, Copyright (C) year name of author + Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, the commands you use may +be called something other than `show w' and `show c'; they could even be +mouse-clicks or menu items--whatever suits your program. + +You should also get your employer (if you work as a programmer) or your +school, if any, to sign a "copyright disclaimer" for the program, if +necessary. Here is a sample; alter the names: + + Yoyodyne, Inc., hereby disclaims all copyright interest in the program + `Gnomovision' (which makes passes at compilers) written by James Hacker. + + , 1 April 1989 + Ty Coon, President of Vice + +This General Public License does not permit incorporating your program into +proprietary programs. If your program is a subroutine library, you may +consider it more useful to permit linking proprietary applications with the +library. If this is what you want to do, use the GNU Lesser General +Public License instead of this License. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/gmndl.cabal b/gmndl.cabal new file mode 100644 --- /dev/null +++ b/gmndl.cabal @@ -0,0 +1,32 @@ +Name: gmndl +Version: 0.1 +Synopsis: Mandelbrot Set explorer using GTK + +Description: + A Mandelbrot Set explorer, with (currently) only one + control (left-click to center and zoom in). Multple + render threads use higher precision maths at higher + zoom levels. Suggested usage: + . + @gmndl +RTS -N -qa -RTS --width=640 --height=480@ + +Cabal-version: >=1.4 +License: GPL-2 +License-file: LICENSE +Author: Claude Heiland-Allen +Maintainer: claudiusmaximus@goto10.org +Category: Graphics +Build-type: Simple + +Executable gmndl + Main-is: gmndl.hs + Build-depends: base >= 4 && < 5, + array >= 0.3 && < 0.4, + gtk >= 0.11 && < 0.12, + gtkglext >= 0.11 && < 0.12, + mtl, + OpenGL >= 2.4 && < 2.5, + priority-queue >= 0.2.1 && < 0.3, + qd >= 0.2 && < 0.3 + GHC-options: -O2 -Wall -threaded -fno-excess-precision -funbox-strict-fields + GHC-prof-options: -O2 -Wall -threaded -fno-excess-precision -funbox-strict-fields -prof -auto-all -caf-all diff --git a/gmndl.hs b/gmndl.hs new file mode 100644 --- /dev/null +++ b/gmndl.hs @@ -0,0 +1,346 @@ +{-# LANGUAGE BangPatterns #-} +module Main (main) where + +import Control.Concurrent (killThread) +import Control.Monad (when) +import Control.Monad.Trans (liftIO) +import Data.Array.IO (IOUArray, newArray, readArray, writeArray, inRange) +import Data.IORef (newIORef, readIORef, writeIORef) +import Data.List (isPrefixOf) +import Data.PriorityQueue (PriorityQueue, newPriorityQueue, enqueue, enqueueBatch, dequeue) +import Foreign (mallocBytes, nullPtr, plusPtr, pokeArray, pokeByteOff, Word8) +import Foreign.C (CDouble) +import GHC.Conc (forkOnIO, numCapabilities) +import Graphics.UI.Gtk +import Graphics.UI.Gtk.OpenGL +import qualified Graphics.Rendering.OpenGL as GL +import Graphics.Rendering.OpenGL (($=), GLfloat) +import Numeric.QD.DoubleDouble (DoubleDouble(DoubleDouble)) +import Numeric.QD.QuadDouble (QuadDouble(QuadDouble)) +import qualified Numeric.QD.DoubleDouble as DD +import qualified Numeric.QD.QuadDouble as QD +import Numeric.QD.FPU.Raw (fpu_fix_start) +import Unsafe.Coerce (unsafeCoerce) + +type B = Word8 +type N = Int +type R = Double + +convert :: (Real a, Fractional b) => a -> b +convert = realToFrac +convertDouble2CDouble :: Double -> CDouble +convertDouble2CDouble !x = unsafeCoerce x +convertCDouble2Double :: CDouble -> Double +convertCDouble2Double !x = unsafeCoerce x +convertDouble2DoubleDouble :: Double -> DoubleDouble +convertDouble2DoubleDouble !x = convertCDouble2DoubleDouble . convertDouble2CDouble $ x +convertCDouble2DoubleDouble :: CDouble -> DoubleDouble +convertCDouble2DoubleDouble !x = DoubleDouble x 0 +convertDoubleDouble2Double :: DoubleDouble -> Double +convertDoubleDouble2Double !(DoubleDouble x _) = convertCDouble2Double x +convertDoubleDouble2CDouble :: DoubleDouble -> CDouble +convertDoubleDouble2CDouble !(DoubleDouble x _) = x +{-# RULES "convert/Double2CDouble" convert = convertDouble2CDouble #-} +{-# RULES "convert/CDouble2Double" convert = convertCDouble2Double #-} +{-# RULES "convert/Double2DoubleDouble" convert = convertDouble2DoubleDouble #-} +{-# RULES "convert/CDouble2DoubleDouble" convert = convertCDouble2DoubleDouble #-} +{-# RULES "convert/DoubleDouble2Double" convert = convertDoubleDouble2Double #-} +{-# RULES "convert/DoubleDouble2CDouble" convert = convertDoubleDouble2CDouble #-} + +data Complex c = {-# UNPACK #-} !c :+ {-# UNPACK #-} !c deriving (Read, Show, Eq) + +instance Num c => Num (Complex c) where + {-# SPECIALIZE instance Num (Complex Float) #-} + {-# SPECIALIZE instance Num (Complex Double) #-} + {-# SPECIALIZE instance Num (Complex DoubleDouble) #-} + {-# SPECIALIZE instance Num (Complex QuadDouble) #-} + (!(a :+ b)) + (!(c :+ d)) = {-# SCC "C+" #-} ((a + c) :+ (b + d)) + (!(a :+ b)) - (!(c :+ d)) = {-# SCC "C-" #-} ((a - c) :+ (b - d)) + (!(a :+ b)) * (!(c :+ d)) = {-# SCC "C*" #-} ((a * c - b * d) :+ (a * d + b * c)) + negate !(a :+ b) = (-a) :+ (-b) + abs x = error $ "Cx.abs: " ++ show x + signum x = error $ "Cx.signum: " ++ show x + fromInteger !x = fromInteger x :+ 0 + +class Num c => Turbo c where + sqr :: c -> c + sqr !x = x * x + twice :: c -> c + twice !x = x + x + +instance Turbo Float where +instance Turbo Double where +instance Turbo CDouble where + +instance Turbo c => Turbo (Complex c) where + {-# SPECIALIZE instance Turbo (Complex Float) #-} + {-# SPECIALIZE instance Turbo (Complex Double) #-} + {-# SPECIALIZE instance Turbo (Complex DoubleDouble) #-} + {-# SPECIALIZE instance Turbo (Complex QuadDouble) #-} + sqr !(r :+ i) = (sqr r - sqr i) :+ (twice (r * i)) + twice !(r :+ i) = (twice r) :+ (twice i) + +instance Turbo DoubleDouble where + sqr !x = DD.sqr x + twice !(DoubleDouble a b) = DoubleDouble (twice a) (twice b) + +instance Turbo QuadDouble where + sqr !x = QD.sqr x + twice !(QuadDouble a b c d) = QuadDouble (twice a) (twice b) (twice c) (twice d) + +hsv2rgb :: R -> R -> R -> (R, R, R) +hsv2rgb !h !s !v + | s == 0 = (v, v, v) + | h == 1 = hsv2rgb 0 s v + | otherwise = + let !i = floor (h * 6) `mod` 6 :: N + !f = (h * 6) - fromIntegral i + !p = v * (1 - s) + !q = v * (1 - s * f) + !t = v * (1 - s * (1 - f)) + in case i of + 0 -> (v, t, p) + 1 -> (q, v, p) + 2 -> (p, v, t) + 3 -> (p, q, v) + 4 -> (t, p, v) + 5 -> (v, p, q) + _ -> (0, 0, 0) + +colour :: Complex Double -> Complex Double -> N -> (B, B, B) +colour !(zr:+zi) !(dzr:+dzi) !n = + let !il2 = 1 / log 2 + !zd2 = sqr zr + sqr zi + !dzd2 = sqr dzr + sqr dzi + !d = (fromIntegral n :: R) - log (log zd2 / log escapeR2) * il2 + !dwell = fromIntegral (floor d :: N) + !finala = atan2 zi zr + !de = (log zd2 * il2) * sqrt zd2 / sqrt dzd2 + !dscale = log de * il2 + 32 + !hue = log d * il2 / 3 + !saturation = 0 `max` (log d * il2 / 8) `min` 1 + !value = 0 `max` (1 - dscale / 48) `min` 1 + !h = hue - fromIntegral (floor hue :: N) + !k = dwell / 2 + !satf = if k - fromIntegral (floor k :: N) >= (0.5 :: R) then 0.9 else 1 + !valf = if finala < 0 then 0.9 else 1 + (!r, !g, !b) = hsv2rgb h (satf * saturation) (valf * value) + !rr = floor $ 0 `max` 255 * r `min` 255 + !gg = floor $ 0 `max` 255 * g `min` 255 + !bb = floor $ 0 `max` 255 * b `min` 255 + in (rr, gg, bb) + +data Job c = Job !N !N !(Complex c) !(Complex c) !(Complex c) !N + +priority :: Job c -> N +priority !(Job _ _ _ _ _ n) = n + +addJob :: RealFloat c => N -> N -> Complex c -> N -> PriorityQueue IO (Job c) -> IOUArray (N,N) Bool -> N -> N -> IO () +{-# SPECIALIZE addJob :: N -> N -> Complex Float -> N -> PriorityQueue IO (Job Float) -> IOUArray (N,N) Bool -> N -> N -> IO () #-} +{-# SPECIALIZE addJob :: N -> N -> Complex Double -> N -> PriorityQueue IO (Job Double) -> IOUArray (N,N) Bool -> N -> N -> IO () #-} +{-# SPECIALIZE addJob :: N -> N -> Complex DoubleDouble -> N -> PriorityQueue IO (Job DoubleDouble) -> IOUArray (N,N) Bool -> N -> N -> IO () #-} +{-# SPECIALIZE addJob :: N -> N -> Complex QuadDouble -> N -> PriorityQueue IO (Job QuadDouble) -> IOUArray (N,N) Bool -> N -> N -> IO () #-} +addJob !w !h !c !zoom todo sync !i !j = do + already <- readArray sync (j, i) + when (not already) $ do + writeArray sync (j, i) True + enqueue todo $! Job i j (coords w h c zoom i j) 0 0 0 + +renderer :: (Turbo c, RealFloat c) => ((N,N),(N,N)) -> (N -> N -> B -> B -> B -> IO ()) -> Complex c -> N -> IO (IO ()) +{-# SPECIALIZE renderer :: ((N,N),(N,N)) -> (N -> N -> B -> B -> B -> IO ()) -> Complex Float -> N -> IO (IO ()) #-} +{-# SPECIALIZE renderer :: ((N,N),(N,N)) -> (N -> N -> B -> B -> B -> IO ()) -> Complex Double -> N -> IO (IO ()) #-} +{-# SPECIALIZE renderer :: ((N,N),(N,N)) -> (N -> N -> B -> B -> B -> IO ()) -> Complex DoubleDouble -> N -> IO (IO ()) #-} +{-# SPECIALIZE renderer :: ((N,N),(N,N)) -> (N -> N -> B -> B -> B -> IO ()) -> Complex QuadDouble -> N -> IO (IO ()) #-} +renderer rng output !c !zoom = do + workerts <- mapM (\w -> forkOnIO w $ worker rng c zoom output w) [ 0 .. workers - 1 ] + return $ do + mapM_ killThread workerts + +coords :: RealFloat c => N -> N -> Complex c -> N -> N -> N -> Complex c +{-# SPECIALIZE coords :: N -> N -> Complex Float -> N -> N -> N -> Complex Float #-} +{-# SPECIALIZE coords :: N -> N -> Complex Double -> N -> N -> N -> Complex Double #-} +{-# SPECIALIZE coords :: N -> N -> Complex DoubleDouble -> N -> N -> N -> Complex DoubleDouble #-} +{-# SPECIALIZE coords :: N -> N -> Complex QuadDouble -> N -> N -> N -> Complex QuadDouble #-} +coords !w !h !c !zoom !i !j = c + ( (fromIntegral (i - w`div`2) * k) + :+(fromIntegral (h`div`2 - j) * k)) + where !k = convert (1/2^^zoom :: Double) + +border :: N -> N -> [(N, N)] +border !w !h = concat $ + [ [ (j, i) | i <- [ 0 .. w - 1 ], j <- [ 0 .. workers - 1 ] ] + , [ (j, i) | j <- [ 0 .. h - 1 ], i <- [ 0 .. workers - 1 ] ] + , [ (j, i) | j <- [ 0 .. h - 1 ], i <- [ w - workers .. w - 1 ] ] + , [ (j, i) | i <- [ 0 .. w - 1 ], j <- [ h - workers .. h - 1 ] ] + ] + +worker :: (Turbo c, RealFloat c) => ((N,N),(N,N)) -> Complex c -> N -> (N -> N -> B -> B -> B -> IO ()) -> N -> IO () +{-# SPECIALIZE worker :: ((N,N),(N,N)) -> Complex Float -> N -> (N -> N -> B -> B -> B -> IO ()) -> N -> IO () #-} +{-# SPECIALIZE worker :: ((N,N),(N,N)) -> Complex Double -> N -> (N -> N -> B -> B -> B -> IO ()) -> N -> IO () #-} +{-# SPECIALIZE worker :: ((N,N),(N,N)) -> Complex DoubleDouble -> N -> (N -> N -> B -> B -> B -> IO ()) -> N -> IO () #-} +{-# SPECIALIZE worker :: ((N,N),(N,N)) -> Complex QuadDouble -> N -> (N -> N -> B -> B -> B -> IO ()) -> N -> IO () #-} +worker rng@((y0,x0),(y1,x1)) !c !zoom output !me = do + sync <- newArray rng False + queue <- newPriorityQueue priority + let addJ = addJob w h c zoom queue sync + js = filter mine (border w h) + w = x1 - x0 + 1 + h = y1 - y0 + 1 + mapM_ (flip (writeArray sync) True) js + enqueueBatch queue (map (\(j,i) -> Job i j (coords w h c zoom i j) 0 0 0) js) + compute rng addJ output queue + where + mine (j, _) = j `mod` workers == me + +compute :: (Turbo c, RealFloat c) => ((N,N),(N,N)) -> (N -> N -> IO ()) -> (N -> N -> B -> B -> B -> IO ()) -> PriorityQueue IO (Job c) -> IO () +{-# SPECIALIZE compute :: ((N,N),(N,N)) -> (N -> N -> IO ()) -> (N -> N -> B -> B -> B -> IO ()) -> PriorityQueue IO (Job Float) -> IO () #-} +{-# SPECIALIZE compute :: ((N,N),(N,N)) -> (N -> N -> IO ()) -> (N -> N -> B -> B -> B -> IO ()) -> PriorityQueue IO (Job Double) -> IO () #-} +{-# SPECIALIZE compute :: ((N,N),(N,N)) -> (N -> N -> IO ()) -> (N -> N -> B -> B -> B -> IO ()) -> PriorityQueue IO (Job DoubleDouble) -> IO () #-} +{-# SPECIALIZE compute :: ((N,N),(N,N)) -> (N -> N -> IO ()) -> (N -> N -> B -> B -> B -> IO ()) -> PriorityQueue IO (Job QuadDouble) -> IO () #-} +compute rng addJ output queue = do + mjob <- dequeue queue + case mjob of + Just (Job i j c z dz n) -> do + let done' !(zr:+zi) !(dzr:+dzi) !n' = do + let (r, g, b) = colour (convert zr :+ convert zi) (convert dzr :+ convert dzi) n' + output i j r g b + sequence_ + [ addJ x y + | u <- spreadX + , v <- spreadY + , let x = i + u + , let y = j + v + , inRange rng (y, x) + ] + todo' z' dz' n' = enqueue queue $! Job i j c z' dz' n' + calculate c limit z dz n done' todo' + compute rng addJ output queue + Nothing -> return () + +calculate :: (Turbo c, RealFloat c) => Complex c -> N -> Complex c -> Complex c -> N -> (Complex c -> Complex c -> N -> IO ()) -> (Complex c -> Complex c -> N -> IO ()) -> IO () +{-# SPECIALIZE calculate :: Complex Float -> N -> Complex Float -> Complex Float -> N -> (Complex Float -> Complex Float -> N -> IO ()) -> (Complex Float -> Complex Float -> N -> IO ()) -> IO () #-} +{-# SPECIALIZE calculate :: Complex Double -> N -> Complex Double -> Complex Double -> N -> (Complex Double -> Complex Double -> N -> IO ()) -> (Complex Double -> Complex Double -> N -> IO ()) -> IO () #-} +{-# SPECIALIZE calculate :: Complex DoubleDouble -> N -> Complex DoubleDouble -> Complex DoubleDouble -> N -> (Complex DoubleDouble -> Complex DoubleDouble -> N -> IO ()) -> (Complex DoubleDouble -> Complex DoubleDouble -> N -> IO ()) -> IO () #-} +{-# SPECIALIZE calculate :: Complex QuadDouble -> N -> Complex QuadDouble -> Complex QuadDouble -> N -> (Complex QuadDouble -> Complex QuadDouble -> N -> IO ()) -> (Complex QuadDouble -> Complex QuadDouble -> N -> IO ()) -> IO () #-} +calculate !c !m0 !z0 !dz0 !n0 done todo = go m0 z0 dz0 n0 + where + go !m !z@(zr:+zi) !dz !n + | not (sqr zr + sqr zi < er2) = done z dz n + | m <= 0 = todo z dz n + | otherwise = go (m - 1) (sqr z + c) (let !zdz = z * dz in twice zdz + 1) (n + 1) + !er2 = convert escapeR2 + +renderer' :: Real c => ((N,N),(N,N)) -> (N -> N -> B -> B -> B -> IO ()) -> Complex c -> N -> IO (IO ()) +renderer' rng output !c !zoom + | zoom < 20 = renderer rng output (f c :: Complex Float ) zoom + | zoom < 50 = renderer rng output (f c :: Complex Double ) zoom + | zoom < 100 = renderer rng output (f c :: Complex DoubleDouble) zoom + | otherwise = renderer rng output (f c :: Complex QuadDouble ) zoom + where f !(cr :+ ci) = convert cr :+ convert ci + +data Args = Args{ aWidth :: N, aHeight :: N } + +defaultArgs :: Args +defaultArgs = Args{ aWidth = 788, aHeight = 576 } + +combineArgs :: Args -> String -> Args +combineArgs a0 s + | "--width=" `isPrefixOf` s = a0{ aWidth = read $ "--width=" `dropPrefix` s } + | "--height=" `isPrefixOf` s = a0{ aHeight = read $ "--height=" `dropPrefix` s } + | "-w=" `isPrefixOf` s = a0{ aWidth = read $ "-w=" `dropPrefix` s } + | "-h=" `isPrefixOf` s = a0{ aHeight = read $ "-h=" `dropPrefix` s } + | otherwise = a0 + +dropPrefix :: String -> String -> String +dropPrefix p s = drop (length p) s + +roundUp2 :: N -> N +roundUp2 n = head . dropWhile (< n) . iterate (2*) $ 1 + +main :: IO () +main = do + args <- foldl combineArgs defaultArgs `fmap` unsafeInitGUIForThreadedRTS + let width = aWidth args + height = aHeight args + size = roundUp2 (width `max` height) + rng = ((0, 0), (height - 1, width - 1)) + _ <- initGL + glconfig <- glConfigNew [ GLModeRGBA, GLModeDouble ] + canvas <- glDrawingAreaNew glconfig + widgetSetSizeRequest canvas width height + imgdata <- mallocBytes $ width * height * 3 + pokeArray imgdata (replicate (height * width * 3) (255 :: B)) + let output x y r g b = do + let p = imgdata `plusPtr` ((y * width + x) * 3) + pokeByteOff p 0 r + pokeByteOff p 1 g + pokeByteOff p 2 b + window <- windowNew + eventb <- eventBoxNew + set window [ containerBorderWidth := 0, containerChild := eventb,windowResizable := False ] + set eventb [ containerBorderWidth := 0, containerChild := canvas ] + mapM_ (flip forkOnIO $ fpu_fix_start nullPtr) [ 0 .. numCapabilities - 1 ] + stop0 <- renderer' rng output c0 zoom0 + sR <- newIORef (c0, zoom0, stop0) + _ <- eventb `on` buttonPressEvent $ {-# SCC "cb/event" #-} tryEvent $ do + LeftButton <- eventButton + (x, y) <- eventCoordinates + liftIO $ do + (c, zoom, stop) <- readIORef sR + stop + pokeArray imgdata (replicate (height * width * 3) (255 :: B)) + let c' = c + ((convert x :+ convert (-y)) - (fromIntegral width :+ fromIntegral (-height)) * (0.5 :+ 0)) * ((1/2^^zoom) :+ 0) + zoom' = zoom + 1 + stop' <- renderer' rng output c' zoom' + writeIORef sR (c', zoom', stop') + print (c', zoom') -- FIXME replace with GUI widgets + _ <- onRealize canvas $ {-# SCC "cb/realize" #-}withGLDrawingArea canvas $ \_ -> do + GL.clearColor $= (GL.Color4 0.0 0.0 0.0 0.0) + GL.matrixMode $= GL.Projection + GL.loadIdentity + GL.ortho 0.0 1.0 0.0 1.0 (-1.0) 1.0 + GL.drawBuffer $= GL.BackBuffers + [tex] <- GL.genObjectNames 1 + GL.texture GL.Texture2D $= GL.Enabled + GL.textureBinding GL.Texture2D $= Just tex + GL.texImage2D Nothing GL.NoProxy 0 GL.RGB' (GL.TextureSize2D (fromIntegral size) (fromIntegral size)) 0 (GL.PixelData GL.RGB GL.UnsignedByte nullPtr) + GL.textureFilter GL.Texture2D $= ((GL.Nearest, Nothing), GL.Nearest) + GL.textureWrapMode GL.Texture2D GL.S $= (GL.Repeated, GL.ClampToEdge) + GL.textureWrapMode GL.Texture2D GL.T $= (GL.Repeated, GL.ClampToEdge) + _ <- onExpose canvas $ {-# SCC "cb/expose" #-} \_ -> do + withGLDrawingArea canvas $ \glwindow -> do + let v :: GLfloat -> GLfloat -> GLfloat -> GLfloat -> IO () + v tx ty vx vy = GL.texCoord (GL.TexCoord2 tx ty) >> GL.vertex (GL.Vertex2 vx vy) + w = fromIntegral width + h = fromIntegral height + sx = fromIntegral width / fromIntegral size + sy = fromIntegral height / fromIntegral size + GL.clear [ GL.ColorBuffer ] + GL.texSubImage2D Nothing 0 (GL.TexturePosition2D 0 0) (GL.TextureSize2D w h) (GL.PixelData GL.RGB GL.UnsignedByte imgdata) + GL.renderPrimitive GL.Quads $ do + v 0 sy 0 0 >> v 0 0 0 1 >> v sx 0 1 1 >> v sx sy 1 0 + glDrawableSwapBuffers glwindow + return True + _ <- onDestroy window mainQuit + _ <- timeoutAdd (widgetQueueDraw canvas >> return True) 200 + widgetShowAll window + mainGUI + +spreadX, spreadY :: [ N ] +spreadX = [ -workers, 0, workers ] +spreadY = [ -1, 0, 1 ] + +workers :: N +workers = numCapabilities + +limit :: N +limit = (2^(11::N)-1) + +c0 :: Complex QuadDouble +c0 = 0 + +zoom0 :: N +zoom0 = 6 + +escapeR, escapeR2 :: R +escapeR = 65536 +escapeR2 = escapeR * escapeR