title: Is it dark outside? Draw a sun map!

text:
# Question
Title: Is it dark outside? Draw a sun map!
Our closest star, the sun, is quite fidgety. The times it rises and sets depend on where you are, and whether it is winter or not.

We would like to be able to deduce if the sun is shining outside without having to leave the comforts of our basements - which is why we need an up-to-date sun map (a.k.a. daylight map). You are the one to write a program that generates just that!

**Rules:** Your program should output an image (in a known format) or an ASCII art representation of our planet, showing (an approximation of) which parts are currently lit by the sun. Your program must be *original* and *self-contained*: you are not allowed to copy, use, include or call any code except your programming language's standard libraries.

If you still don't have a clue what I'm talking about, here's an example from Wikipedia:

This is a *popularity contest*. You should note in your answer which of the following you're trying to achieve (multiple choices possible):

* Correctness. Note that the rules say 'an approximation of' - the better your approximation, the more points in this category. You can check your implementation against Wolfram Alpha's, Time and Date's or die.net's.
* Functionality. For example, what about interactivity? Marking specific locations? Mapping other planets?
* Aesthetics. Drawing continents? Bonus points. Textured continents? Bonus points. On a 3D earth? With clouds? Stars? *Correct* stars? Massive bonus points. And so on.
* Using uncommon, old, or just plain wrong technology. Sure, you could whip this up in Mathematica, but have you considered using `m4`? SQL? Forth? x86 assembly?
* Fun. Want to use a Dymaxion-projection map? Go ahead!
* Short code. This is the Code Golf SE, after all.

Have fun!

# Answer
> 89 votes
# Haskell - low quality code

I was extremely tired when I wrote this.

I might have gone too far with projections idea, anyway, here's the projection the program uses. Basically like projecting earth onto a cube and then unfolding it. Besides, in this projection, the shadow is made of straight lines.
The program uses current date/time, and outputs a PPM file on stdout.

```
import Data.Time.Clock
import Data.Time.Calendar
import Control.Applicative
import Data.Fixed
import Data.Maybe

earth :: [[Int]]
earth = [[256],[256],[256],[256],[64,1,1,2,1,5,14,16,152],[56,19,3,27,1,6,50,1,2,1,90],[53,6,1,11,2,36,26,1,2,1,16,2,1,1,2,1,24,4,66],[47,2,5,14,4,35,22,7,54,2,1,3,60],[38,1,2,2,3,1,6,1,2,1,2,7,6,1,1,33,24,3,3,1,56,2,60],[34,2,1,4,2,1,3,1,1,3,3,2,15,3,3,29,57,5,19,1,2,11,17,1,1,1,34],[40,3,10,2,1,8,16,27,54,3,18,19,18,1,36],[33,6,5,3,2,3,1,3,2,2,1,5,16,21,1,2,53,2,10,1,6,19,1,7,4,3,9,2,33],[32,4,1,7,1,2,3,2,1,1,3,11,14,23,53,2,10,3,1,4,2,33,7,7,29],[8,5,25,10,5,3,2,14,10,2,1,18,1,2,31,6,18,1,7,4,1,60,22],[5,18,2,12,3,5,1,3,2,2,1,3,4,2,3,8,11,18,30,13,9,2,7,3,2,72,1,6,8],[4,36,2,1,1,4,3,7,1,4,3,9,8,15,34,18,2,2,2,17,1,78,4],[4,1,1,27,3,1,1,24,6,3,1,1,1,3,6,13,13,1,20,15,1,4,1,104,1],[3,31,1,24,1,2,4,8,10,9,12,6,18,7,3,7,1,1,2,99,3,2,2],[7,50,2,2,2,1,2,1,3,2,1,2,10,7,15,1,20,7,2,111,7,1],[4,35,1,15,9,1,1,3,4,1,12,5,34,8,3,110,10],[4,9,1,2,1,37,12,6,16,3,34,8,3,96,5,6,13],[6,6,1,1,8,32,12,6,3,1,49,9,4,2,1,86,1,3,4,2,19],[9,2,1,1,11,31,11,11,40,1,8,1,2,4,5,83,12,3,20],[8,1,16,33,9,11,39,2,8,1,2,3,3,83,13,5,19],[28,33,5,12,40,2,7,3,6,62,1,19,13,5,20],[27,36,2,15,34,3,2,2,6,71,1,22,11,2,22],[30,21,1,11,2,16,33,3,1,4,2,72,1,24,1,1,9,1,23],[31,21,1,26,39,4,1,98,1,1,33],[31,42,7,1,40,100,1,1,33],[33,25,2,15,4,4,35,102,36],[33,23,2,1,2,14,8,1,36,27,1,9,1,61,3,1,33],[33,26,5,14,42,10,1,11,2,2,2,7,3,5,1,9,1,44,38],[33,26,1,2,1,9,2,1,45,7,1,2,2,9,8,6,2,6,1,53,4,2,33],[33,26,1,4,1,6,44,8,6,2,3,7,9,5,3,56,1,1,4,3,33],[33,37,45,8,7,2,3,6,2,4,3,6,4,53,43],[33,36,46,6,6,1,4,1,2,2,3,16,3,47,1,5,8,2,34],[34,34,46,7,11,1,3,2,2,16,3,45,6,2,8,1,35],[34,33,48,5,11,1,4,1,4,16,2,49,3,2,6,2,35],[35,32,54,8,17,60,5,2,4,4,35],[36,30,50,12,18,60,8,2,1,1,38],[38,27,50,15,16,61,6,2,41],[38,25,51,18,3,4,6,62,6,1,42],[39,1,1,17,2,3,51,93,49],[40,1,1,11,9,2,49,31,1,10,2,50,49],[40,1,2,9,10,2,48,33,1,10,2,49,49],[41,1,2,8,11,1,47,34,2,10,5,44,50],[42,1,2,7,58,36,1,11,2,1,8,36,51],[46,6,58,36,2,15,7,34,2,1,49],[46,6,12,2,43,38,2,14,7,2,1,12,1,15,55],[46,6,5,2,7,2,41,38,2,14,10,10,4,10,59],[47,6,3,3,10,3,38,37,3,12,11,8,6,9,2,1,57],[49,10,51,38,3,9,13,7,8,9,9,2,48],[51,7,51,40,2,7,15,6,9,1,1,8,8,2,48],[55,7,47,41,1,6,17,4,12,8,8,1,49],[57,5,47,42,1,2,20,4,13,8,9,1,47],[59,3,8,1,38,43,22,4,13,1,2,4,10,2,46],[60,2,6,5,38,41,1,4,18,3,17,3,10,2,46],[61,2,1,1,2,3,1,7,34,45,18,2,18,1,60],[63,1,2,13,33,44,22,1,12,1,16,3,45],[66,14,33,43,22,1,13,1,14,1,1,1,46],[66,18,30,4,1,1,5,30,34,1,2,2,9,3,50],[66,19,43,27,34,2,2,1,7,3,52],[65,20,43,26,36,2,1,2,5,5,51],[65,21,42,24,39,3,4,7,2,1,1,1,1,1,44],[56,1,7,23,41,16,1,6,41,2,4,6,7,1,44],[64,25,39,16,1,5,42,3,4,5,2,1,8,1,2,1,37],[64,29,35,22,43,3,1,1,2,3,2,1,1,1,2,1,1,2,1,7,6,1,27],[63,31,35,20,45,2,11,1,9,7,4,2,26],[64,32,34,19,67,1,2,6,1,2,28],[65,31,34,12,1,6,48,4,18,6,31],[65,31,34,19,54,2,1,2,2,1,10,2,2,1,30],[66,29,36,14,1,3,57,1,19,2,28],[66,29,36,14,1,4,63,1,42],[67,27,36,15,1,4,63,5,3,2,33],[67,26,37,20,5,2,53,2,1,4,4,2,33],[68,25,37,20,4,3,52,9,3,3,32],[70,23,36,20,3,4,53,11,1,4,31],[71,22,37,17,5,4,51,18,31],[71,22,37,16,7,3,50,20,30],[71,21,39,15,6,3,5,1,42,24,29],[71,20,40,15,6,3,47,26,28],[71,17,43,15,6,3,46,28,27],[71,16,45,13,8,1,48,27,27],[71,16,45,12,58,28,26],[71,16,45,12,58,28,26],[70,16,47,10,59,28,26],[70,15,49,9,60,27,26],[70,14,50,7,62,7,6,13,27],[70,13,51,6,63,6,8,1,1,9,28],[70,10,138,10,28],[69,12,139,7,29],[69,11,141,5,19,3,8],[69,8,167,3,9],[69,8,166,1,1,1,10],[70,5,149,2,16,2,12],[69,6,166,3,12],[68,6,166,2,14],[68,5,166,3,14],[68,6,182],[67,6,183],[68,4,184],[68,4,6,2,176],[69,4,183],[70,5,20,1,160],[256],[256],[256],[256],[256],[256],[78,1,1,1,109,1,65],[75,2,115,1,23,1,39],[72,3,80,1,1,5,20,42,32],[74,1,70,1,4,21,5,52,2,1,25],[67,1,2,2,1,4,64,28,4,62,21],[69,9,34,1,1,1,1,1,1,1,2,48,3,69,15],[50,1,5,1,16,5,34,130,14],[32,1,1,2,4,1,3,1,4,29,32,128,18],[20,1,1,54,32,128,20],[17,49,34,137,19],[9,1,2,54,20,4,6,143,17],[16,51,18,5,10,135,21],[11,1,4,54,25,140,21],[12,66,4,155,19],[12,231,13],[0,6,9,5,2,234],[0,256],[0,256]]
main = do
    header
    mapM_ line [0..299]
    where
        header = do
            putStrLn "P3"
            putStrLn "# Some PPM readers expect a comment here"
            putStrLn "400 300"
            putStrLn "2"
        line y = mapM_ (\x -> pixel x y >>= draw) [0..399]
            where
                draw (r, g, b) = putStrLn $ (show r) ++ " " ++ (show g) ++ " " ++ (show b)
                pixel x y = fromMaybe (return (1, 1, 1)) $
                    mapRegion (\x y -> (50, -x, y)) (x - 50) (y - 50)
                    <|> mapRegion (\x y -> (-x, -50, y)) (x - 150) (y - 50)
                    <|> mapRegion (\x y -> (-x, y, 50)) (x - 150) (y - 150)
                    <|> mapRegion (\x y -> (-50, y, -x)) (x - 250) (y - 150)
                    <|> mapRegion (\x y -> (y, 50, -x)) (x - 250) (y - 250)
                    <|> mapRegion (\x y -> (y, -x, -50)) (x - 350) (y - 250)
                    where
                        mapRegion f x y = if x >= -50 && y >= -50 && x < 50 && y < 50 then
                            Just $ fmap (worldMap . shade) getCurrentTime
                            else Nothing
                                where
                                    t (x, y, z) = (atan2 y z) / pi
                                    p (x, y, z) = asin (x / (sqrt $ x*x+y*y+z*z)) / pi * 2
                                    rotate o (x, y, z) = (x, y * cos o + z * sin o, z * cos o - y * sin o)
                                    tilt o (x, y, z) = (x * cos o - y * sin o, x * sin o + y * cos o, z)
                                    shade c = ((t $ rotate yearAngle $ tilt 0.366 $ rotate (dayAngle - yearAngle) $ f x y)) `mod'` 2 > 1
                                        where
                                            dayAngle = fromIntegral (fromEnum $ utctDayTime c) / 43200000000000000 * pi + pi / 2
                                            yearAngle = (fromIntegral $ toModifiedJulianDay $ utctDay c) / 182.624 * pi + 2.5311
                                    worldMap c = case (c, index (t $ f x y) (p $ f x y)) of
                                            (False, False) -> (0, 0, 0)
                                            (False, True) -> (0, 0, 1)
                                            (True, False) -> (2, 1, 0)
                                            (True, True) -> (0, 1, 2)
                                            where
                                                index x y = index' (earth !! (floor $ (y + 1) * 63)) (floor $ (x + 1) * 127) True
                                                    where
                                                        index' [] _ p = False
                                                        index' (x:d) n p
                                                            | n < x = p
                                                            | otherwise = index' d (n - x) (not p)

```

That's right - triangular `where`-code, nested `case`s, invalid IO usage.

# Answer
> 62 votes
# Haskell, in the 'because it's there' category

I was curious so I wrote one. The formulas are reasonably accurate\[1\], but then I go and use some ascii art instead of a proper Plate Carrée map, because it looked nicer (the way I convert pixels to lat/long only works correctly for Plate Carrée)

```
import Data.Time
d=pi/180
tau=2*pi
m0=UTCTime(fromGregorian 2000 1 1)(secondsToDiffTime(12*60*60))
dark lat long now =
  let
    time=(realToFrac$diffUTCTime now m0)/(60*60*24)
    hour=(realToFrac$utctDayTime now)/(60*60)
    mnlong=280.460+0.9856474*time
    mnanom=(357.528+0.9856003*time)*d
    eclong=(mnlong+1.915*sin(mnanom)+0.020*sin(2*mnanom))*d
    oblqec=(23.439-0.0000004*time)*d
    ra=let num=cos(oblqec)*sin(eclong)
           den=cos(eclong) in
       if den<0 then atan(num/den)+pi else atan(num/den)
    dec=asin(sin(oblqec)*sin(eclong))
    gmst =6.697375+0.0657098242*time+hour
    lmst=(gmst*15*d)+long
    ha=(lmst-ra)
    el=asin(sin(dec)*sin(lat)+cos(dec)*cos(lat)*cos(ha))
  in
  el<=0

td x = fromIntegral x :: Double
keep="NSEW"++['0'..'9']
pixel p dk=if dk && p`notElem`keep then if p==' ' then '#' else '%' else p
showMap t= do
  let w=length(worldmap!!0)
      h=length worldmap
  putStrLn (worldmap!!0)
  putStrLn (worldmap!!1)
  mapM_(\y->do
           mapM_(\x->let
                    lat=(0.5-td y/td h)*pi
                    long=(0.5-td x/td w)*tau
                    in
                     putStr [pixel ((worldmap!!(y+2))!!x) (dark lat long t)]) [0..(w-1)]
           putStrLn "") [0..(h-4)]
  putStrLn (last worldmap)

main = do {t<-getCurrentTime; showMap t}

worldmap=[
 "180 150W  120W  90W   60W   30W  000   30E   60E   90E   120E  150E 180",
 "|    |     |     |     |     |    |     |     |     |     |     |     |",
 "+90N-+-----+-----+-----+-----+----+-----+-----+-----+-----+-----+-----+",
 "|          . _..::__:  ,-\"-\"._       |7       ,     _,.__             |",
 "|  _.___ _ _<_>`!(._`.`-.    /        _._     `_ ,_/  '  '-._.---.-.__|",
 "|.{     \" \" `-==,',._\\{  \\  / {)     / _ \">_,-' `                mt-2_|",
 "+ \\_.:--.       `._ )`^-. \"'      , [_/(                       __,/-' +",
 "|'\"'     \\         \"    _L       oD_,--'                )     /. (|   |",
 "|         |           ,'         _)_.\\\\._<> 6              _,' /  '   |",
 "|         `.         /          [_/_'` `\"(                <'}  )      |",
 "+30N       \\\\    .-. )          /   `-'\"..' `:._          _)  '       +",
 "|   `        \\  (  `(          /         `:\\  > \\  ,-^.  /' '         |",
 "|             `._,   \"\"        |           \\`'   \\|   ?_)  {\\         |",
 "|                `=.---.       `._._       ,'     \"`  |' ,- '.        |",
 "+000               |    `-._        |     /          `:`<_|h--._      +",
 "|                  (        >       .     | ,          `=.__.`-'\\     |",
 "|                   `.     /        |     |{|              ,-.,\\     .|",
 "|                    |   ,'          \\   / `'            ,\"     \\     |",
 "+30S                 |  /             |_'                |  __  /     +",
 "|                    | |                                 '-'  `-'   \\.|",
 "|                    |/                                        \"    / |",
 "|                    \\.                                            '  |",
 "+60S                                                                  +",
 "|                     ,/           ______._.--._ _..---.---------._   |",
 "|    ,-----\"-..?----_/ )      _,-'\"             \"                  (  |",
 "|.._(                  `-----'                                      `-|",
 "+90S-+-----+-----+-----+-----+----+-----+-----+-----+-----+-----+-----+",
 "Map 1998 Matthew Thomas. Freely usable as long as this line is included"]

```

Example output, from a more interesting time of year (we're near the equinox, so Wander Nauta's rectangular blobs are fairly accurate :) ) - this is for Jan 16 13:55:51 UTC 2014:

```
180 150W  120W  90W   60W   30W  000   30E   60E   90E   120E  150E 180
|    |     |     |     |     |    |     |     |     |     |     |     |
%90N%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%##########%#%%%%%%%%##%%%%%%%#######%7#######%#####%%%%%#############%
%##%%%%%#%#%%%%%%%%%%%%%%####%########%%%#####%%#%%%##%##%%%%%%%%%%%%%%
%%%#####%#%#%%%%%%%%%%%##%##%#%%#####%#%#%%%%%%#%################%%%2%%
%#%%%%%%%#######%%%#%%%%%#%%######, [_/(         ##############%%%%%%#%
%%%%#####%#########%####%%#####  oD_,--'            ####%#####%%#%%###%
%#########%###########%%#####    _)_.\\._<> 6        ######%%%#%##%###%
%#########%%#########%######    [_/_'` `"(             ###%%%##%######%
%30N#######%%####%%%#%#####     /   `-'"..' `:._       ###%%##%#######%
%###%########%##%##%%#####     /         `:\  > \  ,-^. #%%#%#########%
%#############%%%%###%%###     |           \`'   \|   ?_)##%%#########%
%################%%%%%%%#      `._._       ,'     "`  |' %%#%%########%
%000###############%####`-._        |     /          `:`<_%%%%%%######%
%##################%####    >       .     | ,          `=.%%%%%%%#####%
%###################%%#    /        |     |{|              %%%%%#####%%
%####################%#  ,'          \   / `'            ,"#####%#####%
%30S#################%  /             |_'                |  %%##%#####%
%####################% |                                 '-'##%%%###%%%
%####################|/                                      ##%####%#%
%####################\.                                       #####%##%
%60S################                                          ########%
%##################   ,/           ______._.--._ _..---.-------%%%%###%
%####%%%%%%%%%%%%%--_/ )      _,-'"             "                ##%##%
%%%%%###########       `-----'                                    ##%%%
%90S%%%%%%%%%----+-----+-----+----+-----+-----+-----+-----+-----+----%%
Map 1998 Matthew Thomas. Freely usable as long as this line is included

```

\[1\] they're the same as you'll find elsewhere, except without the extra work to keep degrees between 0 and 360, hours between 0 and 24, and radians between 0 and 2pi. I think those are holdovers from the days we used slide rules; trig functions work just fine outside those ranges...

# Answer
> 39 votes
.

# Bash, 882\* characters

This is my second entry, this time in the **Aesthetics**, **Weird tech**, **Fun** and **Short code** categories. It's inspired by Ram Narasimhan's entry and Peter Taylor's comment.

The script first generates a low-res texture of the world, bundled as base64-encoded data. It then generates 24 PovRay-scenes containing a sphere with that texture, each one rotated to 'face the sun'. Finally, the frames are combined into a GIF animation using ImageMagick. *This means you'll have to have both PovRay and ImageMagick installed for the script to work - feel free to ignore this entry if you think that should disqualify it.*

Like Ram's entry, and my first entry, this does not account for seasonal change, which means it's not very precise. It is, however, shorter, prettier and more precise than my first entry - and unlike Ram's entry, the map data and the code for generating the GIF animation are included.

```
                               echo '
                    iVBO  Rw0KGgoAAAA       NS
              UhE  U g      AAAEgAAAA                     kAQMAAAAQFe4lAAAABlB
    MVEUAFFwAbxKgAD63 AAAA   AWJLR0                  QAiAUdSAAAAAlwSFlzAAALEwAACx
 MB AJqcGAAAAAd0SU1FB9  4DE  hUWI   op      Fp5MAAADDSURBVBhXrcYhTsNQGADgr3ShE4Qi
    h4BeYQFBgqAJN8Lh    +r                jBb rArIJHPobgAgkzgeSQkVHT7MWThAHzq44
           /j/jezy6jSH  M6fB           gd  9T Nbxdl99R4Q+XpdNRISj4dlFRCz
            oI11FxIpup4uIRDe5           fokp0Y2W25jQFDfrGNGsDNsoqBaGj34D2
             bA7TcAwnmRoDZM             5tLkePUJb6uIT2rEq7hKaUhUHCXWpv7Q
             PqEv1rsuoc7X              RbV Bn2d   kGTYKMQ3C7H8z2+wc/eMd S
              QW39v8kAAA               AA      SUVOR K5CYII='|base64 \
               -di>t;for                X in     {0..23};do R=$((90-(\
                $X*15)                )); echo "camera{location <0,
                 0,                   -5> angle 38 }    light_source{
                  <0,0,               -1000> rgb < 2,2,   2>} sphere
                    {<0              ,0,0> 1 pigment      {
                      /**/            image_map{\"t\"        map_type
                        1}}                rotate           <0,$R,0>
                        }">s               ;povray             +Is +H300\
                        +Of$X.png          +W400
                        mogrify            -fill                     white    \
                        -annotate           +0+10                    "$X:00" \
                         -gravity           south                    f$X.png
                         done;              convert                -delay     \
                         100                -loop                 0 $(ls f*  \
                         |sort               -V)                  ani.gif
                        exit;

```

As a **bonus**, here's a GIF that uses NASA's Blue Marble image instead of the space-saving 1-bit texture, i.e. what the result would have looked like without any size restriction: https://i.stack.imgur.com/euOZf.gif

\*: 882 characters not counting decorative whitespace, 1872 characters total.

# Answer
> 25 votes
I decided to kick off the contest with an entry of my own, in the **short code** category. It's 923 characters long, not counting newlines.

# C: 923 characters

Here's the code:

```
i;j;w=160;f=40;t;b;p;s;e;k;d=86400;q=599;
char* m="M('+z EDz :!#\"!*!8S$[\"!$!#\"\")\"!3R)V$'!!()1M./!F)\"!!!!)'/GE5@\"\"!&%.3&,Y$D\"!!%$)5i\"\"\"F\"%&&6%!e'A#!#!!#&$5&!f&A'$*\"5&!c-#'3''8\"$!!#\"U'\"=5$'8#$$\"S(#=7!*5\"!\"#['!A@6#!^H=!#6bH;!!!\"6_!!I;<&!&\"!!$\"F\"!I8;&\"#\"$&#\"C#\"I7<%#!\"/\"BP5=$*,\"=#\"$!L4A%&\"\"G\"\"\"#M1@)*F\"%P/@,!N#!S(E;!@W'E=!!!<Y&D7!&!\"$7\\$D8!)$4_$C8!('&#&!!a&@9!&(%$&g$>9!$*#(%h\">:!!-\"(%&!b!$&5:!\"+\"(!!#$!!!c+5<-!'!'!#!e)5:.!(!&!\"\"e,:25!!!\"!\"\"h-;07#\"$h.9/:\"\"$!!#\"a17-;'!\"$!!\"$!X46,<\"%\"&$\\45,>#&!$$#!W45,C!!!'!\"!$!V26,H\"#!$!\"!\"!S17-#!A!!#\"!_07,\"#A&!\"`.7+#\"A*.!Q.7*$\">/^-9)$\"=0^*<)$!>1]*<(D1])>&E2\\)>&F&!)\\)@#G$%(\\'w%]'x#,\"P%z .\"P%z .!R$z -\"S$z b#z c#z d#z 3";
main(){
t=(time(0)%d*160)/d;
printf("P2\n%d 62\n5\n",w);
for(;i<q;i++){
for(j=m[i]-' ';j>0;j--){
p=k%w,s=(t-f),e=(t+f);
printf("%c ","1324"[b*2+((p>s&&p<e)||(p>s+w&&p<e+w)||(p>s-w&&p<e-w))]);
k++;
}
b=!b;
}
}

```

Here's how it works:

> A crude bitmap of the world* is run-length encoded as a string. Every character in the string represents a run of either land or sea pixels. Long runs of sea are split into a run of sea, then 0 land pixels, then another run of sea, to avoid including unprintable characters in the string. The Python script I wrote to convert PBM files into this format is here.
>
> I then use time() to find out how many seconds have passed in Greenwich since midnight, 1 January 1970. I modulo that to find out how many seconds have passed there today, using that information to position the light portion of the map more-or-less accordingly (I hope).
>
> Correctness is a joke. There's no math at all. The code assumes the earth is a cylinder (block-shaped day/night), that the sun is directly above the equator (no summer/winter), and that you like the color gray (no color).
>
> On the plus side, I do draw continents.
>
> The output is in Portable Graymap (PGM) format, which can then be converted to PNG by something like ImageMagick or the GIMP.

Here's an example output, converted to PNG (larger version):

\*: The entire world except Antarctica, but who lives there anyway...

# Answer
> 24 votes
# Haskell - it's Hammer time.

I did another one. Adapted from my previous version, this one uses an oblique Hammer projection to show both poles at the same time (in fact you're seeing the whole earth in each frame). Just for added weirdness, instead of using a bitmap directly, I sampled the earth along a spiral to give approximately equal area coverage; this is what allows me to distort the earth and rotate it easily. Hammer's projection is equal area too; my idea was that pairing these two things would lead to less distortion when I fill in the gaps. I display a graticule on the projection too, using Bresenham's algorithm to draw the lines. The globe and terminator line both move over the course of the day.

Edited to change the image to a higher resolution but (deliberately) coarser underlying bitmap, so you can see the effect of the spiral. This uses 5000 points (sampled from ~260000), equivalent to a 50x100 bitmap, but giving more resolution to the equator than to the poles.

To use the code, compile with ghc, run with an optional numeric parameter which is the hour offset; files are generated like 'earth0.pgm', 'earth1.pgm'.

```
import System.Environment
import Data.List (intercalate,unfoldr)
import qualified Data.Set as Set
import Data.List.Split
import Data.List
import Data.Maybe (catMaybes)
import qualified Data.Map as Map
import Data.Time
import Debug.Trace
d=pi/180
tau=2*pi
m0=UTCTime(fromGregorian 2000 1 1)(secondsToDiffTime(12*60*60))
dark::Double->Double->UTCTime->Bool
dark lat long now =
  let
    time=(realToFrac$diffUTCTime now m0)/(60*60*24)
    hour=(realToFrac$utctDayTime now)/(60*60)
    mnlong=280.460+0.9856474*time
    mnanom=(357.528+0.9856003*time)*d
    eclong=(mnlong+1.915*sin(mnanom)+0.020*sin(2*mnanom))*d
    oblqec=(23.439-0.0000004*time)*d
    ra=let num=cos(oblqec)*sin(eclong)
           den=cos(eclong) in
       if den<0 then atan(num/den)+pi else atan(num/den)
    dec=asin(sin(oblqec)*sin(eclong))
    gmst =6.697375+0.0657098242*time+hour
    lmst=(gmst*15*d)+long
    ha=(lmst-ra)
    el=asin(sin(dec)*sin(lat)+cos(dec)*cos(lat)*cos(ha))
  in
  el<=0
infill(open, known)=
  if null open then known else infill gen
  where
    neighbours (x,y)=catMaybes $ map ((flip Map.lookup) known) [(x+1,y),(x-1,y),(x,y+1),(x,y-1),(x+1,y+1),(x-1,y+1),(x-1,y-1),(x-1,y-1)]
    vote a= if null a then Nothing
             else Just ((sum a)`div`(length a))
    fill x (open',  known')=
      case vote (neighbours x) of
        Nothing->(x:open',known')
        Just c->(open',(x,c):known')
    gen=(\(o,k)->(o,Map.fromList k))$foldr fill ([], Map.toList known) open
mpoint (a,b)=case a of Nothing->Nothing;Just c->Just(c,b)
grid w h n g lut= map (\y->map (\x->if Set.member (x,y) g then 3 else case Map.lookup (x,y) lut of Nothing->7;Just c->c) [1..w]) [1..h]
unknowns w h lut=concatMap (\y->concatMap (\x->let z=1-(2*x//w-1)^2-(2*y//h-1)^2 in case Map.lookup (x,y) lut of Nothing->if z<0 then [] else [(x,y)];_->[]) [1..w]) [1..h]
main=do
  args <- getArgs
  let off = if null args then 0 else read(args!!0)
  actual <- getCurrentTime
  let now=((fromIntegral off)*60*60) `addUTCTime` actual
  let tod=realToFrac(utctDayTime now)/86400+0.4
  let s=5000
  let w=800
  let h=400
  let n=6
  -- pbm <- readFile "earth.pbm"
  -- let bits=ungrid s$parsepbm pbm
  let bits=[0,23,4,9,1,3,1,2,6,10,1,10,4,1,3,7,10,7,4,2,2,1,2,6,12,1,1,2,1,5,4,1,8,1,3,
            1,21,7,2,2,35,1,4,3,2,2,2,2,16,1,25,1,2,8,1,4,1,2,13,3,2,1,26,1,1,10,3,3,8,
            2,3,6,1,3,25,2,1,10,15,5,1,6,2,3,30,10,15,19,32,11,16,20,35,11,1,2,14,22,27,
            1,8,14,16,22,2,1,22,1,1,2,1,1,2,1,2,1,3,16,14,25,1,2,21,1,6,1,2,1,1,2,3,17,
            14,26,1,2,1,1,26,1,1,3,3,1,1,19,13,28,4,1,26,6,6,21,11,35,40,21,11,37,41,20,
            2,4,4,1,1,39,19,1,6,1,16,19,2,4,5,40,18,2,7,1,17,19,1,1,1,1,1,2,3,46,7,1,5,
            4,25,16,3,1,1,3,5,44,1,4,5,4,3,6,4,1,19,22,5,46,2,3,4,6,2,9,22,22,2,50,1,5,
            2,1,1,6,1,8,24,15,5,1,2,51,2,5,1,1,1,5,1,10,23,14,9,55,1,4,2,17,16,1,4,14,9,
            57,4,1,3,17,13,20,11,54,2,1,3,1,2,20,12,18,13,47,4,3,8,21,10,17,15,44,5,1,1,
            4,1,3,2,22,10,15,16,46,4,3,1,2,2,25,9,17,15,47,1,1,3,30,9,18,13,46,2,1,4,25,
            2,1,11,16,13,46,8,24,2,2,9,16,11,45,12,22,1,3,7,17,10,45,12,21,1,3,7,19,8,
            43,12,25,6,19,8,41,12,25,5,20,7,40,11,25,4,20,6,40,5,3,2,48,6,38,3,54,4,30,
            1,6,2,55,2,29,1,5,1,53,3,28,1,55,3,49,1,30,2,76,1,284,3,4,1,15,1,17,10,1,9,
            7,1,13,21,4,4,1,2,6,17,2,8,3,63]
  let t(phi,lambda)=unitsphere$rx (-pi/4)$rz (-tod*4*pi)$sphereunit(phi, lambda)
  let hmr=(fmap (\(x,y)->(floor((fl w)*(x+4)/8),floor((fl h)*(y+2)/4)))).hammer.t
  let g=graticule hmr n
  let lut = Map.fromList$ catMaybes $map mpoint$map (\((lat,long),bit)->(hmr(lat,long),bit*4+2-if dark lat long now then 2 else 0))  $zip (spiral s) (rld bits)
  -- let lut = Map.fromList$ catMaybes $map mpoint$map (\((lat,long),bit)->(hmr(lat,long),bit))$zip (spiral s) (rld bits)
  let lut' = infill ((unknowns w h lut), lut)
  let pgm = "P2\n"++((show w)++" "++(show h)++" 7\n")++(intercalate "\n" $ map (intercalate " ")$chunksOf 35 $ map show(concat$grid w h n g lut'))++"\n"
  writeFile ("earth"++(show off)++".pgm") pgm

fl=fromIntegral
spiral::Int->[(Double,Double)]
spiral n=map (\k-> let phi=acos(((2*(fl k))-1)/(fl n)-1) in rerange(pi/2-phi,sqrt((fl n)*pi)*phi)) [1..n]
rld::[Int]->[Int]
rld bits=concat$rld' (head bits) (tail bits)
  where
   rld' bit []=[]
   rld' bit (run:xs) = (replicate run bit):(rld' (case bit of 1->0;_->1) xs)
rle::[Int]->[Int]
rle bits=(head bits):(map length$group bits)
sample::Int->Int->Int->[(Int,Int)]
sample n w h = map (\(phi, theta)->((floor((fl h)*((phi-(pi/2))/pi)))`mod`h, (floor((fl w)*(theta-pi)/(tau)))`mod`w )) $ spiral n
ungrid::Int->[[Int]]->[Int]
ungrid n g = rle $ map (\(y, x)->(g!!y)!!x) (sample n w h)
  where w = length$head g
        h = length g
parsepbm::[Char]->[[Int]]
parsepbm pbm=
    let header = lines pbm
        format = head header
        [width, height] = map read$words (head$drop 1 header)
        rest = drop 2 header
        d = ((map read).concat.(map words)) rest
    in chunksOf width d
rerange(phi,lambda)
 | abs(phi)>pi = rerange(phi - signum(phi)*tau, lambda)
 | abs(phi)>pi/2 = rerange(phi-signum(phi)*pi, lambda+pi)
 | abs(lambda)>pi = rerange(phi, lambda - signum(lambda)*tau)
 | otherwise = (phi, lambda)
laea(phi,lambda)=if isInfinite(z) then Nothing else Just (z*cos(phi)*sin(lambda),z*sin(phi)) where z=4/sqrt(1+cos(phi)*cos(lambda))
hammer(phi,lambda)=case laea(phi, lambda/2) of Nothing->Nothing; Just(x,y)->Just (x, y/2)
bresenham :: (Int, Int)->(Int, Int)->[(Int, Int)]
bresenham p0@(x0,y0) p1@(x1,y1)
  | abs(dx)>50||abs(dy)>50=[]
  | x0>x1 = map h$ bresenham (h p0) (h p1)
  | y0>y1 = map v$ bresenham (v p0) (v p1)
  | (x1-x0) < (y1-y0) = map f$ bresenham (f p0) (f p1)
  | otherwise = unfoldr (\(x,y,d)->if x>x1 then Nothing else Just((x,y),(if 2*(d+dy)<dx then(x+1,y,d+dy)else(x+1,y+1,d+dy-dx)))) (x0,y0,0)
      where
        h(x,y)=(-x,y)
        v(x,y)=(x,-y)
        f(x,y)=(y,x)
        dx=x1-x0
        dy=y1-y0
globe n k=
  (concatMap (\m->map (meridian m) [k*(1-n)..k*(n-1)]) [k*(1-2*n),k*(2-2*n)..k*2*n])
  ++(concatMap (\p->map (parallel p) [k*(-2*n)..k*2*n]) [k*(1-n),k*(2-n)..k*(n-1)])
  where
  meridian m p=(radians(p,m),radians(p+1,m))
  parallel p m=(radians(p,m),radians(p,m+1))
  radians(p,m)=rerange((p//(k*n))*pi/2,(m//(k*n))*pi/2)
graticule f n=Set.fromList $ concatMap (\(a,b)->case (f a,f b) of (Nothing,_)->[];(_,Nothing)->[];(Just c,Just d)->bresenham c d) (globe n 4)
rx theta (x,y,z) = (x, y*(cos theta)-z*(sin theta), y*(sin theta)+z*(cos theta))
ry theta (x,y,z) = (z*(sin theta)+x*(cos theta), y, z*(cos theta)-x*(sin theta))
rz theta (x,y,z) = (x*(cos theta)-y*(sin theta), x*(sin theta)+y*(cos theta), z)
sphereunit (phi, theta) = (rz theta (ry (-phi) (1,0,0)))
unitsphere (x,y,z) = (asin z, atan2 y x)
x//y=(fromIntegral x)/(fromIntegral y)

```

# Answer
> 21 votes
## C, using pnm images

Late answer, focusing on **correctness** and **aesthetics**. The output is a blend of two input images (day.pnm and night.pnm), including a stripe of twilight. I'm using images based on NASAs blue marble here.

The code uses my own img.h for clarity (just imagine it being included verbatim in the .c for strict rule compliance...). Everything in there is implemented via C macros. The animations are built with imagemagicks convert from multiple frames - the program itself will only output static images. Code is below.

**Now:** (Aug 13, ~13:00 CEST)

**One day:** (Jan 1st)

**One year:** (12:00 UTC)

**sun.c**

```
  #include <math.h>
  #include <time.h>

  #include "img.h"

  #ifndef M_PI
  #define M_PI 3.14159265359
  #endif

  double deg2rad(double x) {return x / 180.0 * M_PI;}
  double rad2deg(double x) {return x * 180.0 / M_PI;}

  double  sind(double x) {return  sin(deg2rad(x));}
  double  cosd(double x) {return  cos(deg2rad(x));}
  double asind(double x) {return rad2deg(asin(x));}

  double elevation(double latitude, double longitude, int yday, int hour, int min, int sec)
  {
     double fd = (hour + (min + sec / 60.0) / 60.0) / 24.0;
     double fyd = 360.0 * (yday + fd) / 366.0;

     double m = fyd - 3.943;
     double ta = -1.914 * sind(m) + 2.468 * sind(2 * m + 205.6);
     double hourangle = (fd - 0.5) * 360.0 + longitude + ta;
     double decl = 0.396 - 22.913 * cosd(fyd) + 4.025 * sind(fyd) - 0.387 * cosd(2 * fyd) + 0.052 * sind(2 * fyd) - 0.155 * cosd(3 * fyd) + 0.085 * sind(3 * fyd);

     return asind(cosd(hourangle) * cosd(decl) * cosd(latitude) + sind(decl) * sind(latitude));
  }

  int main(int argc, char* argv[])
  {
     Image day, night, out;
     int x, y;
     time_t t = time(0);
     struct tm* utc = gmtime(&t);
     int yday = utc->tm_yday, hour = utc->tm_hour, min = utc->tm_min, sec = utc->tm_sec;

     imgLoad(day, "day.pnm");
     imgLoad(night, "night.pnm");
     imgLoad(out, "day.pnm");
     for(y = 0; y < day.height; ++y)
     {
        double latitude = 90.0 - 180.0 * (y + 0.5) / day.height;
        for(x = 0; x < day.width; ++x)
        {
           double longitude = -180.0 + 360.0 * (x + 0.5) / day.width;
           double elev = elevation(latitude, longitude, yday, hour, min, sec);
           double nf = elev > -0.8 ? 0.0 : elev > -6.0 ? 0.5 : 1.0;
           double df = 1.0 - nf;
           Color dc = imgGetColor(day, x, y);
           Color nc = imgGetColor(night, x, y);
           imgDotC3(out, x, y, df * dc.r + nf * nc.r, df * dc.g + nf * nc.g, df * dc.b + nf * nc.b);
        }
     }
     imgSave(out, "out.pnm");
  }

```

**img.h**

```
  #include <stdlib.h>
  #include <stdio.h>
  #include <string.h>

  typedef struct
  {
     unsigned char r;
     unsigned char g;
     unsigned char b;
  } Color;

  typedef struct
  {
     Color* data;
     int width;
     int height;
     Color c;
  } Image;

  #define imgCreate(img, w, h)           {\
                                            int length;\
                                            (img).width = (w);\
                                            (img).height = (h);\
                                            length = (img).width * (img).height * sizeof(Color);\
                                            (img).data = malloc(length);\
                                            memset((img).data, 0, length);\
                                            (img).c.r = (img).c.g = (img).c.b = 0;\
                                         }

  #define imgDestroy(img)                {\
                                            free((img).data);\
                                            (img).width = 0;\
                                            (img).height = 0;\
                                            (img).c.r = (img).c.g = (img).c.b = 0;\
                                         }

  #define imgSetColor(img, ur, ug, ub)   {\
                                            (img).c.r = (ur);\
                                            (img).c.g = (ug);\
                                            (img).c.b = (ub);\
                                         }

  #define imgDot(img, x, y)              {\
                                            (img).data[(int)(x) + (int)(y) * (img).width] = (img).c;\
                                         }

  #define imgDotC3(img, x, y, ur, ug, ub) {\
                                            (img).data[(int)(x) + (int)(y) * (img).width].r = (ur);\
                                            (img).data[(int)(x) + (int)(y) * (img).width].g = (ug);\
                                            (img).data[(int)(x) + (int)(y) * (img).width].b = (ub);\
                                         }

  #define imgDotC(img, x, y, c)          {\
                                            (img).data[(int)(x) + (int)(y) * (img).width] = (c);\
                                         }

  #define imgGetColor(img, x, y)         ((img).data[(int)(x) + (int)(y) * (img).width])

  #define imgLine(img, x, y, xx, yy)     {\
                                            int x0 = (x), y0 = (y), x1 = (xx), y1 = (yy);\
                                            int dx =  abs(x1 - x0), sx = x0 < x1 ? 1 : -1;\
                                            int dy = -abs(y1 - y0), sy = y0 < y1 ? 1 : -1;\
                                            int err = dx + dy, e2;\
                                            \
                                            for(;;)\
                                            {\
                                               imgDot((img), x0, y0);\
                                               if (x0 == x1 && y0 == y1) break;\
                                               e2 = 2 * err;\
                                               if (e2 >= dy) {err += dy; x0 += sx;}\
                                               if (e2 <= dx) {err += dx; y0 += sy;}\
                                            }\
                                         }

  #define imgSave(img, fname)            {\
                                            FILE* f = fopen((fname), "wb");\
                                            fprintf(f, "P6 %d %d 255\n", (img).width, (img).height);\
                                            fwrite((img).data, sizeof(Color), (img).width * (img).height, f);\
                                            fclose(f);\
                                         }

  #define imgLoad(img, fname)            {\
                                            FILE* f = fopen((fname), "rb");\
                                            char buffer[16];\
                                            int index = 0;\
                                            int field = 0;\
                                            int isP5 = 0;\
                                            unsigned char c = ' ';\
                                            while(field < 4)\
                                            {\
                                               do\
                                               {\
                                                  if(c == '#') while(c = fgetc(f), c != '\n');\
                                               } while(c = fgetc(f), isspace(c) || c == '#');\
                                               index = 0;\
                                               do\
                                               {\
                                                  buffer[index++] = c;\
                                               } while(c = fgetc(f), !isspace(c) && c != '#' && index < 16);\
                                               buffer[index] = 0;\
                                               switch(field)\
                                               {\
                                                  case 0:\
                                                     if (strcmp(buffer, "P5") == 0) isP5 = 1;\
                                                     else if (strcmp(buffer, "P6") == 0) isP5 = 0;\
                                                     else fprintf(stderr, "image format \"%s\" unsupported (not P5 or P6)\n", buffer), exit(1);\
                                                     break;\
                                                  case 1:\
                                                     (img).width = atoi(buffer);\
                                                     break;\
                                                  case 2:\
                                                     (img).height = atoi(buffer);\
                                                     break;\
                                                  case 3:\
                                                     index = atoi(buffer);\
                                                     if (index != 255) fprintf(stderr, "image format unsupported (not 255 values per channel)\n"), exit(1);\
                                                     break;\
                                               }\
                                               field++;\
                                            }\
                                            imgCreate((img), (img).width, (img).height);\
                                            if (isP5)\
                                            {\
                                               int length = (img).width * (img).height;\
                                               for(index = 0; index < length; ++index)\
                                               {\
                                                  (img).data[index].r = (img).data[index].g = (img).data[index].b = fgetc(f);\
                                               }\
                                            }\
                                            else\
                                            {\
                                               fread((img).data, sizeof(Color), (img).width * (img).height, f);\
                                            }\
                                            fclose(f);\
                                         }

```

# Answer
> 20 votes
# R: Using ggplot2 and Map projection

Inspired by @mniip's post, I decided to try using R's mapproj package, wherein we can orient the globe by specifying where the North Pole should be when computing the projection.

Based on current GMT time, I compute the longitude where is it currently noon and make that point the map's center. We are looking at Earth from the "Sun's point of view" so whatever is visible is in daylight.

Much of the code is just aesthetics. The only part I had to figure out was to compute the "noon Longitude," that is the longitude value where it was noon given any GMT time.

```
library(ggplot2);library(maps);library(ggmap)
world <- map_data("world")# a lat-long dataframe from the maps package
worldmap <- ggplot(world, aes(x=long, y=lat, group=group)) +
  geom_path(color="orange") +
  theme(panel.background= element_rect("black"),
        axis.text.y=element_blank(),
        axis.ticks=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        panel.grid.major = element_line(colour="blue", size=0.75),
        panel.grid.minor = element_line(colour="blue")
  )

#Create a function that takes in the current GMT time
print_3d_coordmap <- function (current_gmt_time) {
  curr_gmt_mins <- as.POSIXlt(current_gmt_time)$hour*60 + as.POSIXlt(current_gmt_time)$min
  noon_longitude <- 180 - (curr_gmt_mins * 360/1440)
  #centered at wherever longitude where it is Noon now on (lat:equator)
  worldmap + coord_map("ortho", orientation=c(0, noon_longitude, 0))
}

#test it out
print_3d_coordmap(Sys.time() + 7*60*60) # my location is 7 hours behind UTC

```

I then used the R animation package to generate 24 images and stitched them into one GIF.

# Answer
> 9 votes
## JavaScript – by Martin Kleppe (http://aem1k.com/)

I want to stress that this is not my work, but the work of Martin Kleppe. I just think it fits so perfectly that it should not be missing here:

Online Demo (or just paste it into the console)

```
eval(z='p="<"+"pre>"/*        ######## */;for(y in n="zw24l6k\
4e3t4jnt4qj24xh2 x/*    *############### */42kty24wrt413n243n\
9h243pdxt41csb yz/*  #################### */43iyb6k43pk7243nm\
r24".split(4)){/*     *#################*   */for(a in t=pars\
eInt(n[y],36)+/*          ###############*   */(e=x=r=[]))for\
(r=!r,i=0;t[a/*               ############*   */]>i;i+=.05)wi\
th(Math)x-= /*                #############    */.05,0<cos(o=\
new Date/1e3/*                #########*       */-x/PI)&&(e[~\
~(32*sin(o)*/*                     ####*       */sin(.5+y/7))\
+60] =-~ r);/*                         *###    */for(x=0;122>\
x;)p+="   *#"/*                        #####  */[e[x++]+e[x++\
]]||(S=("eval"/*                      *##### */+"(z=\'"+z.spl\
it(B = "\\\\")./*      ###*           ####  */join(B+B).split\
(Q="\'").join(B+Q/*                  ###* */)+Q+")//m1k")[x/2\
+61*y-1]).fontcolor/*               ##   */(/\\w/.test(S)&&"#\
03B");document.body.innerHTML=p+=B+"\\n"}setTimeout(z)')//m1k\

```

---
Tags: popularity-contest, graphical-output
---


title: Find the craftiest prime

text:
# Question
Title: Find the craftiest prime
# Intro

Consider the process of taking some positive integer *n* in some base *b* and replacing each digit with its representation in the base of the digit to the right.

* If the digit to the right is a 0, use base *b*.
* If the digit to the right is a 1, use unary with 0's as tally marks.
* If there is no digit to the right (i.e. you are in the ones place), loop around to the most significant digit.

As an example let *n* = 160 and *b* = 10. Running the process looks like this:

```
The first digit is 1, the digit to the right is 6, 1 in base 6 is 1.
The next digit is 6, the digit to the right is 0, 0 is not a base so use b, 6 in base b is 6.
The last digit is 0, the digit to the right (looping around) is 1, 0 in base 1 is the empty string (but that's ok).

Concatenating '1', '6', and '' together gives 16, which is read in the original base b = 10.

```

The exact same procedure *but moving left instead of right* can also be done:

```
The first digit is 1, the digit to the left (looping around) is 0, 0 is not a base so use b, 1 in base b is 1.
The next digit is 6, the digit to the left is 1, 6 in base 1 is 000000.
The last digit is 0, the digit to the left is 6, 0 in base 6 is 0.

Concatenating '1', '000000', and '0' together gives 10000000, which is read in the original base b = 10.

```

Thus, we've made two numbers related to 160 (for *b* = 10): 16 and 10000000.

*We will define n to be a **crafty number** if it evenly divides at least one of the two numbers generated in this process into 2 or more parts*

In the example *n* is crafty because 160 divides 10000000 exactly 62500 times.

203 is NOT crafty because the resulting numbers are 2011 and 203 itself, which 203 cannot fit evenly into 2 or more times.

# Challenge

(For the rest of the problem we will only consider *b* = 10.)

The challenge is to write a program that finds the highest crafty number that is also prime.

The first 7 crafty primes (and all that I have found so far) are:

```
2
5
3449
6287
7589
9397
93557 <-- highest so far (I've searched to 100,000,000+)

```

I am not officially certain whether more exist, but I expect they do. If you can prove that there are (or aren't) finitely many I'll give you +200 bounty rep.

The winner will be the person who can provide the highest crafty prime, provided it is apparent that they have been active in the search and are not intentionally taking glory from others.

# Rules

* You may use any prime finding tools you want.
* You may use probabilistic prime testers.
* You may reuse other peoples code **with attribution**. This is a communal effort. Cutthroat tactics will not be tolerated.
* Your program must actively search for the prime. You may start your search at the highest known crafty prime.
* Your program should be able to compute **all** of the known crafty primes within 4 hours of Amazon EC2 t2.medium instances (either four at once or one for four hours or something in between). I will not actually be testing it on them and you certainly don't need to. This is just a benchmark.

Here is my Python 3 code I used for generating the table above: (runs in a second or two)

```
import pyprimes

def toBase(base, digit):
    a = [
            ['0', '1', '2', '3', '4', '5', '6', '7', '8', '9'],
            ['', '0', '00', '000', '0000', '00000', '000000', '0000000', '00000000', '000000000' ],
            ['0', '1', '10', '11', '100', '101', '110', '111', '1000', '1001'],
            ['0', '1', '2', '10', '11', '12', '20', '21', '22', '100'],
            ['0', '1', '2', '3', '10', '11', '12', '13', '20', '21'],
            ['0', '1', '2', '3', '4', '10', '11', '12', '13', '14'],
            ['0', '1', '2', '3', '4', '5', '10', '11', '12', '13'],
            ['0', '1', '2', '3', '4', '5', '6', '10', '11', '12'],
            ['0', '1', '2', '3', '4', '5', '6', '7', '10', '11'],
            ['0', '1', '2', '3', '4', '5', '6', '7', '8', '10']
        ]
    return a[base][digit]

def getCrafty(start=1, stop=100000):
    for p in pyprimes.primes_above(start):
        s = str(p)
        left = right = ''
        for i in range(len(s)):
            digit = int(s[i])
            left += toBase(int(s[i - 1]), digit)
            right += toBase(int(s[0 if i + 1 == len(s) else i + 1]), digit)
        left = int(left)
        right = int(right)
        if (left % p == 0 and left // p >= 2) or (right % p == 0 and right // p >= 2):
            print(p, left, right)
        if p >= stop:
            break
    print('DONE')

getCrafty()

```

# Answer
> 7 votes
## Mathematica, finds 93,557 in 0.3s (no further crafty primes below 2\*1010)

This is just a naive exhaustive search through all primes. To begin with it checks about 1,000,000 primes every 55 seconds (which is bound to get slower as the primes get larger).

I'm using a bunch of helper functions:

```
lookup = {
  {0, 1, 2, 3, 4, 5, 6, 7, 8, 9},
  {{}, 0, {0, 0}, {0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0},
   {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0}},
  {0, 1, {1, 0}, {1, 1}, {1, 0, 0}, {1, 0, 1}, {1, 1, 0}, {1, 1, 1}, {1, 0, 0, 0},
   {1, 0, 0, 1}},
  {0, 1, 2, {1, 0}, {1, 1}, {1, 2}, {2, 0}, {2, 1}, {2, 2}, {1, 0, 0}},
  {0, 1, 2, 3, {1, 0}, {1, 1}, {1, 2}, {1, 3}, {2, 0}, {2, 1}},
  {0, 1, 2, 3, 4, {1, 0}, {1, 1}, {1, 2}, {1, 3}, {1, 4}},
  {0, 1, 2, 3, 4, 5, {1, 0}, {1, 1}, {1, 2}, {1, 3}},
  {0, 1, 2, 3, 4, 5, 6, {1, 0}, {1, 1}, {1, 2}},
  {0, 1, 2, 3, 4, 5, 6, 7, {1, 0}, {1, 1}},
  {0, 1, 2, 3, 4, 5, 6, 7, 8, {1, 0}}
};
convertBase[d_, b_] := lookup[[b + 1, d + 1]];
related[n_] := (
   d = IntegerDigits[n];
   {FromDigits[Flatten[convertBase @@@ Transpose[{d, RotateRight@d}]]],
    FromDigits[Flatten[convertBase @@@ Transpose[{d, RotateLeft@d}]]]}
);
crafty[n_] := (
   {ql, qr} = related[n]/n;
   IntegerQ[ql] && ql > 1 || IntegerQ[qr] && qr > 1
);

```

And then this loop does the actual search:

```
p = 2;
start = TimeUsed[];
i = 1;
While[True,
 If[crafty[p], Print@{"CRAFTY PRIME:", p, TimeUsed[] - start}];
 p = NextPrime@p;
 If[Mod[++i, 1000000] == 0,
  Print[{"Last prime checked:", p, TimeUsed[] - start}]
 ]
]

```

I'll keep updating the post, if I find any primes or can think of optimisations.

It currently checks all primes up to 100,000,000 in about 5.5 minutes.

**Edit:** I decided to follow the OP's example and switched to a lookup table for base conversion. That gave roughly a 30% speedup.

## Crafty Numbers in general

I'm stopping my search for crafty primes now, since I'd need several days just to catch up with where the Perl answer already got. Instead, I started searching for all crafty numbers. Maybe their distribution helps to find a proof the the number of crafty primes is finite or infinite.

I define *right-crafty* numbers those which divide the related number obtained by interpreting the digit to the *right* as the new base, and *left-crafty* numbers accordingly. It will probably help to tackle these individually for a proof.

Here are all left-crafty numbers up to 2,210,000,000:

```
{2, 5, 16, 28, 68, 160, 222, 280, 555, 680, 777, 1600, 2605, 2800,
 6800, 7589, 7689, 9397, 9777, 16000, 16064, 16122, 22222, 24682,
 26050, 28000, 55555, 68000, 75890, 76890, 93557, 160000, 160640,
 161220, 247522, 254408, 260500, 280000, 680000, 758900, 768900,
 949395, 1600000, 1606400, 1612200, 2222222, 2544080, 2605000,
 2709661, 2710271, 2717529, 2800000, 3517736, 5555555, 6800000,
 7589000, 7689000, 9754696, 11350875, 16000000, 16064000, 16122000,
 25440800, 26050000, 27175290, 28000000, 28028028, 35177360, 52623721,
 68000000, 68654516, 75890000, 76890000, 113508750, 129129129, 160000000,
 160640000, 161220000, 222222222, 254408000, 260500000, 271752900,
 275836752, 280000000, 280280280, 333018547, 351773600, 370938016,
 555555555, 680000000, 758900000, 768900000, 777777777, 877827179,
 1135087500, 1291291290, 1600000000, 1606400000, 1612200000, 1944919449}

```

And here are all right-crafty numbers in that range:

```
{2, 5, 16, 28, 68, 125, 128, 175, 222, 284, 555, 777, 1575, 1625,
 1875, 3449, 5217, 6287, 9375, 14625, 16736, 19968, 22222, 52990,
 53145, 55555, 58750, 93750, 127625, 152628, 293750, 529900, 587500,
 593750, 683860, 937500, 1034375, 1340625, 1488736, 2158750, 2222222,
 2863740, 2937500, 5299000, 5555555, 5875000, 5937500, 6838600,
 7577355, 9375000, 12071125, 19325648, 21587500, 28637400, 29375000,
 29811250, 42107160, 44888540, 52990000, 58750000, 59375000, 68386000,
 71461386, 74709375, 75773550, 93750000, 100540625, 116382104,
 164371875, 197313776, 207144127, 215875000, 222222222, 226071269,
 227896480, 274106547, 284284284, 286374000, 287222080, 293750000,
 298112500, 421071600, 448885400, 529900000, 555555555, 587500000,
 593750000, 600481125, 683860000, 714613860, 747093750, 757735500,
 769456199, 777777777, 853796995, 937500000, 1371513715, 1512715127,
 1656354715, 1728817288, 1944919449, 2158750000}

```

Note that there's an infinite number of left-crafty and right-crafty numbers, because there are several ways to generate them from existing ones:

* One can append an arbitrary number of `0`s to any left-crafty number whose least significant digit is greater than its most significant digit to obtain another left-crafty number.
* Likewise, one can append an arbitrary number of `0`s to any right-crafty number whose least significant digit in *less* than its most significant digit. This (and the previous statement) is because the `0` will be appended to both the crafty number and its related number, effectively multiplying both of them by 10.
* Every odd number of `2`s or `5`s is crafty.
* Every odd number of `777`s is crafty.
* It appears that an odd number of `28` joined by `0`s, like `28028028` is always left-crafty.

Other things to note:

* There are at least four 10-digit numbers which consist of two repeated five-digit numbers (which are themselves not crafty, but there might be some pattern here anyway).
* There are a lot of right-crafty numbers that are a multiple of `125`. It might be worth investigating those to find another production rule.
* I haven't found a left-crafty number that starts with 4 or ends with 3.
* Right-crafty numbers can start with any digit but I haven't found a right-crafty number ending in 1 or 3.

I suppose this list would be more interesting if I omitted those whose existence is implied by a smaller crafty number, especially since these are never primes by the construction rules discovered so far. So here are all crafty primes that cannot constructed with one of the above rules:

```
Left-crafty:
{16, 68, 2605, 7589, 7689, 9397, 9777, 16064, 16122, 24682,
 93557, 247522, 254408, 949395, 2709661, 2710271, 2717529, 3517736,
 9754696, 11350875, 52623721, 68654516, 129129129, 275836752,
 333018547, 370938016, 877827179, 1944919449}

Right-crafty:
{16, 28, 68, 125, 128, 175, 284, 1575, 1625, 1875, 3449, 5217,
 6287, 9375, 14625, 16736, 19968, 52990, 53145, 58750, 127625,
 152628, 293750, 593750, 683860, 1034375, 1340625, 1488736, 2158750,
 2863740, 7577355, 12071125, 19325648, 29811250, 42107160, 44888540,
 71461386, 74709375, 100540625, 116382104, 164371875, 197313776,
 207144127, 226071269, 227896480, 274106547, 284284284, 287222080,
 600481125, 769456199, 853796995, 1371513715, 1512715127, 1656354715,
 1728817288, 1944919449}

```

Note also that there a few *doubly-crafty* numbers (those which appear in both lists and hence divide both related numbers):

```
{2, 5, 16, 28, 68, 222, 555, 777, 22222, 55555, 2222222, 5555555, 1944919449}

```

There exist infinitely many of these as well. ~~But as you can see, except for `16`, `28`, `68` these all consist only of a single (repeated) digit. It would also be interesting prove whether any larger numbers can be doubly crafty without having that property, but that might just drop out as a corollary of a proof for singly crafty numbers.~~ Found the counter-example `1944919449`.

# Answer
> 6 votes
# Perl (1e5 in 0.03s, 1e8 in 21s). Max 93557 to 1e11.

Very similar to the original. Changes include:

* transpose the base lookup. Small language-dependent savings.
* mod the incremented right shift instead of if. Language dependent micro-opt.
* use Math::GMPz because Perl 5 doesn't have auto-magic bigints like Python and Perl 6.
* Use 2s \<= left instead of int(left/s) \>= 2. Native integer shift vs. bigint divide.

Does first 1e8 primes in 21 seconds on my fast machine, 3.5 minutes for 1e9, 34 minutes for 1e10. I'm a little surprised it is at all faster than the Python code for small inputs. We could parallelize (Pari/GP's new `parforprime` would be nice for this). Since it is a search we can parallelize by hand I suppose (`forprimes` can take two arguments). `forprimes` is basically like Pari/GP's `forprime` -- it does segmented sieves internally and calls the block with each result. It's convenient, but for this problem I don't think it is a performance area.

```
#!/usr/bin/env perl
use warnings;
use strict;
use Math::Prime::Util qw/forprimes/;
use Math::GMPz;

my @rbase = (
  [   0,"",       0,   0,  0, 0, 0, 0, 0, 0],
  [qw/1 0         1    1   1  1  1  1  1  1/],
  [qw/2 00        10   2   2  2  2  2  2  2/],
  [qw/3 000       11   10  3  3  3  3  3  3/],
  [qw/4 0000      100  11  10 4  4  4  4  4/],
  [qw/5 00000     101  12  11 10 5  5  5  5/],
  [qw/6 000000    110  20  12 11 10 6  6  6/],
  [qw/7 0000000   111  21  13 12 11 10 7  7/],
  [qw/8 00000000  1000 22  20 13 12 11 10 8/],
  [qw/9 000000000 1001 100 21 14 13 12 11 10/],
);

my($s,$left,$right,$slen,$i,$barray);
forprimes {
  ($s,$slen,$left,$right) = ($_,length($_),'','');
  foreach $i (0 .. $slen-1) {
    $barray = $rbase[substr($s,$i,1)];
    $left  .= $barray->[substr($s,$i-1,1)];
    $right .= $barray->[substr($s,($i+1) % $slen,1)];
  }
  $left = Math::GMPz::Rmpz_init_set_str($left,10) if length($left) >= 20;
  $right = Math::GMPz::Rmpz_init_set_str($right,10) if length($right) >= 20;
  print "$s      $left $right\n" if (($s<<1) <= $left && $left % $s == 0)
                                 || (($s<<1) <= $right && $right % $s == 0);
} 1e9;

```

# Answer
> 5 votes
# C++11, with threads and GMP

Timing (on a MacBook Air):

* 4 threads
    + 10^8 in 2.18986s
    + 10^9 in 21.3829s
    + 10^10 in 421.392s
    + 10^11 in 2557.22s
* 1 thread
    + 10^8 in 3.95095s
    + 10^9 in 37.7009s

Requirements:

```
#include <vector>
#include <iostream>
#include <chrono>
#include <cmath>
#include <future>
#include <mutex>
#include <atomic>
#include "primesieve.hpp"
#include "gmpxx.h"

using namespace std;

using ull = unsigned long long;

mutex cout_mtx;
atomic<ull> prime_counter;

string ppnum(ull number) {
    if (number == 0) {
        return "0 * 10^0";
    }
    else {
        int l = floor(log10(number));
        return to_string(number / pow(10, l)) + " * 10^" + to_string(int(l));
    }
}

inline void bases(int& base, int& digit, mpz_class& sofar) {
    switch (base) {
        case 0:
            sofar *= 10;
            sofar += digit;
            break;
        case 1:
            sofar *= pow(10, digit);
            break;
        case 2:
            switch (digit) {
                case 0: sofar *= 10; break;
                case 1: sofar *= 10; sofar += 1; break;
                case 2: sofar *= 100; sofar += 10; break;
                case 3: sofar *= 100; sofar += 11; break;
                case 4: sofar *= 1000; sofar += 100; break;
                case 5: sofar *= 1000; sofar += 101; break;
                case 6: sofar *= 1000; sofar += 110; break;
                case 7: sofar *= 1000; sofar += 111; break;
                case 8: sofar *= 10000; sofar += 1000; break;
                case 9: sofar *= 10000; sofar += 1001; break;
            }
            break;
        case 3:
            switch (digit) {
                case 0: sofar *= 10; break;
                case 1: sofar *= 10; sofar += 1; break;
                case 2: sofar *= 10; sofar += 2; break;
                case 3: sofar *= 100; sofar += 10; break;
                case 4: sofar *= 100; sofar += 11; break;
                case 5: sofar *= 100; sofar += 12; break;
                case 6: sofar *= 100; sofar += 20; break;
                case 7: sofar *= 100; sofar += 21; break;
                case 8: sofar *= 100; sofar += 22; break;
                case 9: sofar *= 1000; sofar += 100; break;
            }
            break;
        case 4:
            switch (digit) {
                case 0: sofar *= 10; break;
                case 1: sofar *= 10; sofar += 1; break;
                case 2: sofar *= 10; sofar += 2; break;
                case 3: sofar *= 10; sofar += 3; break;
                case 4: sofar *= 100; sofar += 10; break;
                case 5: sofar *= 100; sofar += 11; break;
                case 6: sofar *= 100; sofar += 12; break;
                case 7: sofar *= 100; sofar += 13; break;
                case 8: sofar *= 100; sofar += 20; break;
                case 9: sofar *= 100; sofar += 21; break;
            }
            break;
        case 5:
            switch (digit) {
                case 0: sofar *= 10; break;
                case 1: sofar *= 10; sofar += 1; break;
                case 2: sofar *= 10; sofar += 2; break;
                case 3: sofar *= 10; sofar += 3; break;
                case 4: sofar *= 10; sofar += 4; break;
                case 5: sofar *= 100; sofar += 10; break;
                case 6: sofar *= 100; sofar += 11; break;
                case 7: sofar *= 100; sofar += 12; break;
                case 8: sofar *= 100; sofar += 13; break;
                case 9: sofar *= 100; sofar += 14; break;
            }
            break;
        case 6:
            switch (digit) {
                case 0: sofar *= 10; break;
                case 1: sofar *= 10; sofar += 1; break;
                case 2: sofar *= 10; sofar += 2; break;
                case 3: sofar *= 10; sofar += 3; break;
                case 4: sofar *= 10; sofar += 4; break;
                case 5: sofar *= 10; sofar += 5; break;
                case 6: sofar *= 100; sofar += 10; break;
                case 7: sofar *= 100; sofar += 11; break;
                case 8: sofar *= 100; sofar += 12; break;
                case 9: sofar *= 100; sofar += 13; break;
            }
            break;
        case 7:
            switch (digit) {
                case 0: sofar *= 10; break;
                case 1: sofar *= 10; sofar += 1; break;
                case 2: sofar *= 10; sofar += 2; break;
                case 3: sofar *= 10; sofar += 3; break;
                case 4: sofar *= 10; sofar += 4; break;
                case 5: sofar *= 10; sofar += 5; break;
                case 6: sofar *= 10; sofar += 6; break;
                case 7: sofar *= 100; sofar += 10; break;
                case 8: sofar *= 100; sofar += 11; break;
                case 9: sofar *= 100; sofar += 12; break;
            }
            break;
        case 8:
            switch (digit) {
                case 0: sofar *= 10; break;
                case 1: sofar *= 10; sofar += 1; break;
                case 2: sofar *= 10; sofar += 2; break;
                case 3: sofar *= 10; sofar += 3; break;
                case 4: sofar *= 10; sofar += 4; break;
                case 5: sofar *= 10; sofar += 5; break;
                case 6: sofar *= 10; sofar += 6; break;
                case 7: sofar *= 10; sofar += 7; break;
                case 8: sofar *= 100; sofar += 10; break;
                case 9: sofar *= 100; sofar += 11; break;
            }
            break;
        case 9:
            switch (digit) {
                case 0: sofar *= 10; break;
                case 1: sofar *= 10; sofar += 1; break;
                case 2: sofar *= 10; sofar += 2; break;
                case 3: sofar *= 10; sofar += 3; break;
                case 4: sofar *= 10; sofar += 4; break;
                case 5: sofar *= 10; sofar += 5; break;
                case 6: sofar *= 10; sofar += 6; break;
                case 7: sofar *= 10; sofar += 7; break;
                case 8: sofar *= 10; sofar += 8; break;
                case 9: sofar *= 100; sofar += 10; break;
            }
            break;
    };
}

vector<ull> crafty(ull start, ull stop) {
    cout_mtx.lock();
    cout << "Thread scanning from " << start << " to " << stop << endl;
    cout_mtx.unlock();
    vector<ull> res;

    auto prime_iter = primesieve::iterator(start);
    ull num;
    int prev, curr, next, fprev;
    int i, size;
    mpz_class left, right;
    unsigned long num_cpy;
    unsigned long* num_ptr;
    mpz_class num_mpz;

    while ((num = prime_iter.next_prime()) && num < stop) {
        ++prime_counter;
        left = 0;
        right = 0;
        size = floor(log10(num));
        i = pow(10, size);
        prev = num % 10;
        fprev = curr = num / i;
        if (i != 1) {
            i /= 10;
            next = (num / i) % 10;
        }
        else {
            next = prev;
        }
        for (size += 1; size; --size) {
            bases(prev, curr, left);
            bases(next, curr, right);
            prev = curr;
            curr = next;
            if (i > 1) {
                i /= 10;
                next = (num / i) % 10;
            }
            else {
                next = fprev;
            }
        }
        num_cpy = num;

        if (num != num_cpy) {
            num_ptr = (unsigned long *) #
            num_mpz = *num_ptr;
            num_mpz << sizeof(unsigned long) * 8;
            num_mpz += *(num_ptr + 1);
        }
        else {
            num_mpz = num_cpy;
        }
        if ((left % num_mpz == 0 && left / num_mpz >= 2) || (right % num_mpz == 0 && right / num_mpz >= 2)) {
            res.push_back(num);
        }
    }
    cout_mtx.lock();
    cout << "Thread scanning from " << start << " to " << stop << " is done." << endl;;
    cout << "Found " << res.size() << " crafty primes." << endl;
    cout_mtx.unlock();
    return res;
}

int main(int argc, char *argv[]) {
    ull start = 0, stop = 1000000000;
    int number_of_threads = 4;

    if (argc > 1) {
        start = atoll(argv[1]);
    }
    if (argc > 2) {
        stop = atoll(argv[2]);
    }
    if (argc > 3) {
        number_of_threads = atoi(argv[3]);
    }
    ull gap = stop - start;

    cout << "Start: " << ppnum(start) << ", stop: " << ppnum(stop) << endl;
    cout << "Scanning " << ppnum(gap) << " numbers" << endl;
    cout << "Number of threads: " << number_of_threads << endl;

    chrono::time_point<chrono::system_clock> tstart, tend;
    tstart = chrono::system_clock::now();

    cout << "Checking primes..." << endl;

    using ptask = packaged_task<decltype(crafty)>;
    using fur = future<vector<ull>>;

    vector<thread> threads;
    vector<fur> futures;
    for (int i = 0; i < number_of_threads; ++i) {
        auto p = ptask(crafty);
        futures.push_back(move(p.get_future()));
        auto tstop = (i + 1 == number_of_threads) ? (stop) : (start + gap / number_of_threads * (i + 1));
        threads.push_back(thread(move(p), start + gap / number_of_threads * i, tstop));
    }

    vector<ull> res;

    for (auto& thread : threads) {
        thread.join();
    }

    for (auto& fut : futures) {
        auto v = fut.get();
        res.insert(res.end(), v.begin(), v.end());
    }

    cout << "Finished checking primes..." << endl;

    tend = chrono::system_clock::now();
    chrono::duration<double> elapsed_seconds = tend - tstart;

    cout << "Number of tested primes: " << ppnum(prime_counter) << endl;
    cout << "Number of found crafty primes: " << res.size() << endl;
    cout << "Crafty primes are: ";
    for (auto iter = res.begin(); iter != res.end(); ++iter) {
        if (iter != res.begin())
            cout << ", ";
        cout << *iter;
    }
    cout << endl;
    cout << "Time taken: " << elapsed_seconds.count() << endl;
}

```

Output:

```
Start: 0 * 10^0, stop: 1.000000 * 10^11
Scanning 1.000000 * 10^11 numbers
Number of threads: 4
Checking primes...
Thread scanning from 25000000000 to 50000000000
Thread scanning from 0 to 25000000000
Thread scanning from 50000000000 to 75000000000
Thread scanning from 75000000000 to 100000000000
Thread scanning from 75000000000 to 100000000000 is done.
Found 0 crafty primes.
Thread scanning from 50000000000 to 75000000000 is done.
Found 0 crafty primes.
Thread scanning from 25000000000 to 50000000000 is done.
Found 0 crafty primes.
Thread scanning from 0 to 25000000000 is done.
Found 7 crafty primes.
Finished checking primes...
Number of tested primes: 4.118055 * 10^9
Number of found crafty primes: 7
Crafty primes are: 2, 5, 3449, 6287, 7589, 9397, 93557
Time taken: 2557.22

```

# Answer
> 1 votes
# C, with GMP, multithreaded (1e8 in 17s for 1 thread)

Similar in concept to the rest, with probably a bit of optimizations here and there.

Compile: `gcc -I/usr/local/include -Ofast crafty.c -pthread -L/usr/local/lib -lgmp && ./a.out`

Please donate your CPU power. I don't have a fast computer.
1e8 in 17 seconds with 1 thread on my macbook air.

```
#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <gmp.h>
#include <pthread.h>
#include <string.h>

#define THREAD_COUNT 1           // Number of threads
#define MAX_DIGITS   32768       // Maximum digits allocated for the string... some c stuff
#define MAX_NUMBER   "100000000" // Number in string format
#define START_INDEX  1           // Must be an odd number >= 1
#define GET_WRAP_INDEX(index, stringLength) index<0?stringLength+index:index>=stringLength?index-stringLength:index

static void huntCraftyPrime(int startIndex) {

    char lCS [MAX_DIGITS];
    char rCS [MAX_DIGITS];
    char tPS [MAX_DIGITS];

    mpz_t tP, lC, rC, max;
    mpz_init_set_ui(tP, startIndex);
    mpz_init(lC);
    mpz_init(rC);
    mpz_init_set_str(max, MAX_NUMBER, 10);

    int increment = THREAD_COUNT*2;

    if (START_INDEX < 9 && startIndex == START_INDEX) {
        printf("10 10 2\n\n");
        printf("10 10 5\n\n");
    }

    while (mpz_cmp(max, tP) > 0) {
        mpz_get_str(tPS, 10, tP);
        int tPSLength = strlen(tPS);
        int l = 0, r = 0, p = 0;
        while (p < tPSLength) {
            char lD = tPS[GET_WRAP_INDEX(p-1, tPSLength)];
            char d  = tPS[GET_WRAP_INDEX(p  , tPSLength)];
            char rD = tPS[GET_WRAP_INDEX(p+1, tPSLength)];
            if (d == '0') {
                if (lD != '1') lCS[l++] = '0';
                if (rD != '1') rCS[r++] = '0';
            } else if (d == '1') {
                lCS[l++] = (lD != '1') ? '1' : '0';
                rCS[r++] = (rD != '1') ? '1' : '0';
            } else if (d == '2') {
                if (lD == '1') {
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                } else if (lD == '2') {
                    lCS[l++] = '1';
                    lCS[l++] = '0';
                } else {
                    lCS[l++] = '2';
                }
                if (rD == '1') {
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                } else if (rD == '2') {
                    rCS[r++] = '1';
                    rCS[r++] = '0';
                } else {
                    rCS[r++] = '2';
                }
            } else if (d == '3') {
                if (lD == '1') {
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                } else if (lD == '2') {
                    lCS[l++] = '1';
                    lCS[l++] = '1';
                } else if (lD == '3') {
                    lCS[l++] = '1';
                    lCS[l++] = '0';
                } else {
                    lCS[l++] = '3';
                }
                if (rD == '1') {
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                } else if (rD == '2') {
                    rCS[r++] = '1';
                    rCS[r++] = '1';
                } else if (rD == '3') {
                    rCS[r++] = '1';
                    rCS[r++] = '0';
                } else {
                    rCS[r++] = '3';
                }
            } else if (d == '4') {
                if (lD == '1') {
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                } else if (lD == '2') {
                    lCS[l++] = '1';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                } else if (lD == '3') {
                    lCS[l++] = '1';
                    lCS[l++] = '1';
                } else if (lD == '4') {
                    lCS[l++] = '1';
                    lCS[l++] = '0';
                } else {
                    lCS[l++] = '4';
                }
                if (rD == '1') {
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                } else if (rD == '2') {
                    rCS[r++] = '1';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                } else if (rD == '3') {
                    rCS[r++] = '1';
                    rCS[r++] = '1';
                } else if (rD == '4') {
                    rCS[r++] = '1';
                    rCS[r++] = '0';
                } else {
                    rCS[r++] = '4';
                }
            } else if (d == '5') {
                if (lD == '1') {
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                } else if (lD == '2') {
                    lCS[l++] = '1';
                    lCS[l++] = '0';
                    lCS[l++] = '1';
                } else if (lD == '3') {
                    lCS[l++] = '1';
                    lCS[l++] = '2';
                } else if (lD == '4') {
                    lCS[l++] = '1';
                    lCS[l++] = '1';
                } else if (lD == '5') {
                    lCS[l++] = '1';
                    lCS[l++] = '0';
                } else {
                    lCS[l++] = '5';
                }
                if (rD == '1') {
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                } else if (rD == '2') {
                    rCS[r++] = '1';
                    rCS[r++] = '0';
                    rCS[r++] = '1';
                } else if (rD == '3') {
                    rCS[r++] = '1';
                    rCS[r++] = '2';
                } else if (rD == '4') {
                    rCS[r++] = '1';
                    rCS[r++] = '1';
                } else if (rD == '5') {
                    rCS[r++] = '1';
                    rCS[r++] = '0';
                } else {
                    rCS[r++] = '5';
                }
            } else if (d == '6') {
                if (lD == '1') {
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                } else if (lD == '2') {
                    lCS[l++] = '1';
                    lCS[l++] = '1';
                    lCS[l++] = '0';
                } else if (lD == '3') {
                    lCS[l++] = '2';
                    lCS[l++] = '0';
                } else if (lD == '4') {
                    lCS[l++] = '1';
                    lCS[l++] = '2';
                } else if (lD == '5') {
                    lCS[l++] = '1';
                    lCS[l++] = '1';
                } else if (lD == '6') {
                    lCS[l++] = '1';
                    lCS[l++] = '0';
                } else {
                    lCS[l++] = '6';
                }
                if (rD == '1') {
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                } else if (rD == '2') {
                    rCS[r++] = '1';
                    rCS[r++] = '1';
                    rCS[r++] = '0';
                } else if (rD == '3') {
                    rCS[r++] = '2';
                    rCS[r++] = '0';
                } else if (rD == '4') {
                    rCS[r++] = '1';
                    rCS[r++] = '2';
                } else if (rD == '5') {
                    rCS[r++] = '1';
                    rCS[r++] = '1';
                } else if (rD == '6') {
                    rCS[r++] = '1';
                    rCS[r++] = '0';
                } else {
                    rCS[r++] = '6';
                }
            } else if (d == '7') {
                if (lD == '1') {
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                } else if (lD == '2') {
                    lCS[l++] = '1';
                    lCS[l++] = '1';
                    lCS[l++] = '1';
                } else if (lD == '3') {
                    lCS[l++] = '2';
                    lCS[l++] = '1';
                } else if (lD == '4') {
                    lCS[l++] = '1';
                    lCS[l++] = '3';
                } else if (lD == '5') {
                    lCS[l++] = '1';
                    lCS[l++] = '2';
                } else if (lD == '6') {
                    lCS[l++] = '1';
                    lCS[l++] = '1';
                } else if (lD == '7') {
                    lCS[l++] = '1';
                    lCS[l++] = '0';
                } else {
                    lCS[l++] = '7';
                }
                if (rD == '1') {
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                } else if (rD == '2') {
                    rCS[r++] = '1';
                    rCS[r++] = '1';
                    rCS[r++] = '1';
                } else if (rD == '3') {
                    rCS[r++] = '2';
                    rCS[r++] = '1';
                } else if (rD == '4') {
                    rCS[r++] = '1';
                    rCS[r++] = '3';
                } else if (rD == '5') {
                    rCS[r++] = '1';
                    rCS[r++] = '2';
                } else if (rD == '6') {
                    rCS[r++] = '1';
                    rCS[r++] = '1';
                } else if (rD == '7') {
                    rCS[r++] = '1';
                    rCS[r++] = '0';
                } else {
                    rCS[r++] = '7';
                }
            } else if (d == '8') {
                if (lD == '1') {
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                } else if (lD == '2') {
                    lCS[l++] = '1';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                } else if (lD == '3') {
                    lCS[l++] = '2';
                    lCS[l++] = '2';
                } else if (lD == '4') {
                    lCS[l++] = '2';
                    lCS[l++] = '0';
                } else if (lD == '5') {
                    lCS[l++] = '1';
                    lCS[l++] = '3';
                } else if (lD == '6') {
                    lCS[l++] = '1';
                    lCS[l++] = '2';
                } else if (lD == '7') {
                    lCS[l++] = '1';
                    lCS[l++] = '1';
                } else if (lD == '8') {
                    lCS[l++] = '1';
                    lCS[l++] = '0';
                } else {
                    lCS[l++] = '8';
                }
                if (rD == '1') {
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                } else if (rD == '2') {
                    rCS[r++] = '1';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                } else if (rD == '3') {
                    rCS[r++] = '2';
                    rCS[r++] = '2';
                } else if (rD == '4') {
                    rCS[r++] = '2';
                    rCS[r++] = '0';
                } else if (rD == '5') {
                    rCS[r++] = '1';
                    rCS[r++] = '3';
                } else if (rD == '6') {
                    rCS[r++] = '1';
                    rCS[r++] = '2';
                } else if (rD == '7') {
                    rCS[r++] = '1';
                    rCS[r++] = '1';
                } else if (rD == '8') {
                    rCS[r++] = '1';
                    rCS[r++] = '0';
                } else {
                    rCS[r++] = '8';
                }
            } else if (d == '9') {
                if (lD == '1') {
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                } else if (lD == '2') {
                    lCS[l++] = '1';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                    lCS[l++] = '1';
                } else if (lD == '3') {
                    lCS[l++] = '1';
                    lCS[l++] = '0';
                    lCS[l++] = '0';
                } else if (lD == '4') {
                    lCS[l++] = '2';
                    lCS[l++] = '1';
                } else if (lD == '5') {
                    lCS[l++] = '1';
                    lCS[l++] = '4';
                } else if (lD == '6') {
                    lCS[l++] = '1';
                    lCS[l++] = '3';
                } else if (lD == '7') {
                    lCS[l++] = '1';
                    lCS[l++] = '2';
                } else if (lD == '8') {
                    lCS[l++] = '1';
                    lCS[l++] = '1';
                } else if (lD == '9') {
                    lCS[l++] = '1';
                    lCS[l++] = '0';
                } else {
                    lCS[l++] = '9';
                }
                if (rD == '1') {
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                } else if (rD == '2') {
                    rCS[r++] = '1';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                    rCS[r++] = '1';
                } else if (rD == '3') {
                    rCS[r++] = '1';
                    rCS[r++] = '0';
                    rCS[r++] = '0';
                } else if (rD == '4') {
                    rCS[r++] = '2';
                    rCS[r++] = '1';
                } else if (rD == '5') {
                    rCS[r++] = '1';
                    rCS[r++] = '4';
                } else if (rD == '6') {
                    rCS[r++] = '1';
                    rCS[r++] = '3';
                } else if (rD == '7') {
                    rCS[r++] = '1';
                    rCS[r++] = '2';
                } else if (rD == '8') {
                    rCS[r++] = '1';
                    rCS[r++] = '1';
                } else if (rD == '9') {
                    rCS[r++] = '1';
                    rCS[r++] = '0';
                } else {
                    rCS[r++] = '9';
                }
            }
            ++p;
        }
        lCS[l] = '\0';
        rCS[r] = '\0';

        mpz_set_str(lC, lCS, 10);
        mpz_set_str(rC, rCS, 10);

        if ((mpz_divisible_p(lC, tP) && mpz_cmp(lC, tP) > 0) || (mpz_divisible_p(rC, tP) && mpz_cmp(rC, tP) > 0)){
            if (mpz_millerrabin(tP, 25)) {
                gmp_printf("%Zd %Zd %Zd\n\n", lC, rC, tP);
            }
        }
        mpz_add_ui(tP, tP, increment);
    }
}

static void *huntCraftyPrimeThread(void *p) {
    int* startIndex = (int*) p;
    huntCraftyPrime(*startIndex);
    pthread_exit(NULL);
}

int main(int argc, char *argv[]) {

    struct timeval time_started, time_now, time_diff;
    gettimeofday(&time_started, NULL);

    int  startIndexes[THREAD_COUNT];
    pthread_t threads[THREAD_COUNT];

    int startIndex = START_INDEX;
    for (int i = 0; i < THREAD_COUNT; ++i) {
        for (;startIndex % 2 == 0; ++startIndex);
        startIndexes[i] = startIndex;
        int rc = pthread_create(&threads[i], NULL, huntCraftyPrimeThread, (void*)&startIndexes[i]);
        if (rc) {
            printf("ERROR; return code from pthread_create() is %d\n", rc);
            exit(-1);
        }
        ++startIndex;
    }

    for (int i = 0; i < THREAD_COUNT; ++i) {
        void * status;
        int rc = pthread_join(threads[i], &status);
        if (rc) {
            printf("ERROR: return code from pthread_join() is %d\n", rc);
            exit(-1);
        }
    }

    gettimeofday(&time_now, NULL);
    timersub(&time_now, &time_started, &time_diff);
    printf("Time taken,%ld.%.6d s\n", time_diff.tv_sec, time_diff.tv_usec);

    pthread_exit(NULL);
    return 0;
}

```

# Answer
> 0 votes
# Python, finds 93557 in 0.28s

Very similar to OP's code in that it also uses `pyprimes`. I did write this myself though xD

```
import pyprimes, time

d = time.clock()

def to_base(base, n):
    if base == 1:
        return '0'*n
    s = ""
    while n:
        s = str(n % base) + s
        n //= base
    return s

def crafty(n):
    digits = str(n)
    l, r = "", ""
    for i in range(len(digits)):
        t = int(digits[i])
        base = int(digits[i-1])
        l += to_base(base, t) if base else digits[i]
        base = int(digits[(i+1)%len(digits)])
        r += to_base(base, t) if base else digits[i]
    l, r = int(l) if l else 0, int(r) if r else 0
    if (l%n==0 and 2 <= l/n) or (r%n==0 and 2 <= r/n):
        print(n, l, r, time.clock()-d)

for i in pyprimes.primes_above(1):
    crafty(i)

```

It also prints out the time since start that it finds a number.

Output:

```
2 10 10 3.156656792490237e-05
5 10 10 0.0006756015452219958
3449 3111021 3104100 0.012881854420378145
6287 6210007 11021111 0.022036544076745254
7589 751311 125812 0.026288406792971432
9397 1231007 1003127 0.03185028207808106
93557 123121012 10031057 0.27897531840850603

```

Format is `number left right time`. As a comparison, OP's code finds 93557 in around `0.37`.

---
Tags: code-challenge, primes
---