# Question
Title: Reverse-Engineer the N-Bonacci sequence[s]
**EDIT: I will be accepting an answer Monday, 2/15/2016. May the bytes be ever in your favor!**

In his "Print the N-Bonacci Sequence" challenge, @DJMcGoathem describes the N-bonacci sequences, wherein the previous *N* numbers are summed, instead of the traditional 2 of the Fibonacci sequence (said to be the "*duo*nacci sequence"). He then asked to take two inputs, X and N, then output the *X*th *N*-nacci number.

I propose the opposite.
Given a sequence, output which *N*-nacci sequence it is a subset of. I say "subset of" because:

* A) these sequences are infinite
* B) if given the start of the sequence, you could just count the number of leading 1s

**In the case that it could belong to multiple *N*-nacci sequences, chose the lowest one.**
**In the case that it does not belong to any *N-nacci* sequence, then your program may do anything other than print something that could be mistaken for output.** These behaviors include (but are not limited to): infinite loop, error, crash, delete itself (\*cough cough\* vigil \*cough cough\*), or create a black hole (as long as this black hole does not produce anything that could be mistaken for valid output).
**For the sake of this challenge, these sequences start with 1.** This means any *N*-nacci sequence starts with *N* ones. Furthermore, *N* must be a positive integer. So no **-1**-nacci, etc.

Test cases:

```
1,1,1 -> 1
49, 97 -> 7
55, 89, 144 -> 2
1 -> 1
6765 -> 2
12, 23, 45, 89 -> 12
100, 199 -> 100

```

# Answer
## Ruby, 94

I'm quite surprised how far I've been able to golf this within the same algorithm! I started with over 200!

```
->a{1.step.find{|s|x=[1]*(s+z=a.size)
x<<x[-s,s].inject(:+)while x.max<a.max&&s>1
x[-z,z]==a}}

```

Ungolfed:

```
l=->a{
    # ooh! a built-in infinite incrementer!
    1.step.find{|s|
        # if n == 1, z is important, otherwise s.
        x=[1]*(s+z=a.size)
        ## add the next nacci number until we hit or overshot max.
        ## if s == 1, we don't increment, so don't bother
        x<<x[-s,s].inject(:+)while x.max<g&&s>1
        # eval to true if there's a match, false if not
        x[-z,z]==a
    }
}

```
> 7 votes

# Answer
# APL (Dyalog Unicode), 46 bytes

```
{∧/0<⍵:⊃⊃⍸A⍷↑{⍵((+/↑),⊢)⍣{</+/¨A⍵}⍵⍴1}¨⍳⊃A←⌽⍵}

```

Try it online!

Anonymous prefix function, taking a vector of integers.

### How?

```
∧/0<⍵:   ⍝ If the given sequence contains any non-positive integers, then return nothing. Otherwise...
A←⌽⍵     ⍝ Set A to be the reverse of the given sequence
{...}¨⍳⊃ ⍝ For each n in 1 to last element of given sequence, we generate the n-bonacci sequence in reverse order:
 ⍵⍴1       ⍝ Initialize a vector of n 1s (first n terms of the n-bonacci sequence)
 ⍣         ⍝ Do:
 ⍵((+/↑),⊢)  ⍝ Prepend the sum of the most-recent (first) n elements of the vector to the rest of the vector
           ⍝ While
 {</+/¨A⍵}   ⍝ the sum of the n-bonacci sequence generated so far is less than the sum of A
               ⍝ (this doesn't terminate prematurely because:
               ⍝   For n=1, this produces k+2 1s if A consists of k 1s
               ⍝     Otherwise, worst-case is that A has 1 element
               ⍝   For n=2, F_0+F_1+F_2 ... + F_{k-2} = F_k (identity)
               ⍝     F_0 + ... F_{k-1} > F_k, and APL terminates the ⍣-loop in the next step,
               ⍝     after F_k generated
               ⍝   For n>2, F_0+F_1+F_2 ... + F_{k-2} < F_k, so an excess of terms are generated
               ⍝ )
↑        ⍝ Mix to form a matrix, where the n-th row is the n-bonacci sequence
⍸A⍷      ⍝ Find the (row, column) positions of A in the matrix
⊃⊃       ⍝ Take the row of the first one (which would have minimum n)

```
> 5 votes

# Answer
# Python 2, 176 bytes

```
def r(n,x):i,f=n,[1]*n;exec"f+=sum(f[i-n:]),;i+=1;"*(x-n);return f
l=input()
m=max(l)+len(l)
for i in range(1,m):
 if any(l==r(i,m)[b:b+len(l)]for b in range(m)):print i;break;

```

Note that this requires input in this format:

```
[1, 1, 2, 3...]

```

rather than

```
1, 1, 2, 3...

```

Fairly straightforward solution, just to get things rolling. I'll work more on golfing it down once somebody else answers. This uses a slightly modified version of the N-Bonnaci generator from @Data's answer, so props to him. Then for each N-Bonnaci in range of the input, checks if the input is a subsequence of it.
> 3 votes

# Answer
# Lua, ~~324~~ 323 Bytes

When I see other submission, I feel like there's something wrong with my code... But then, I remember that's Lua, and there's isn't all these fancy functionnalities :'(

It was a lot of fun, took me some time actually.

**Edit: Saved 1 byte with a simple trick: using a `::label::`+`goto label` instead a infinite loop done with `while''`.**

```
function f(l)c=2 y,z=table.remove,os.exit while(l[1]<2)do y(l,1)if(#l<1)then print(1)z()end end ::q:: a={}for i=1,c do a[i]=1 end b={}for i=1,#l do b[i]=l[i]end while(a[#a]<b[1])do x=0 for i=(#a-c+1>0 and #a-c+1 or 1),#a do x=x+a[i]end a[#a+1]=x if a[#a]==b[1]then y(b,1)end if #b<1 then print(c)z()end end c=c+1 goto q end

```

### Ungolfed and Explanations

Lua has no way to define a set, subset, or even check if an array/table contains a value without using its index/key. That's why i've decided to remove elements from the array I take as a parameter. that's how I keep record of which elements has already been computed, and if it matched.

```
  function f(l)
  c=2
  y,z=table.remove,os.exit           -- Create pointers on table.remove and os.exit
                                     -- saves a total of 9 bytes
  while(l[1]<2)                      -- loop used to remove leading 1
  do
    y(l,1)
    if(#l<1)                         -- we also check if it was a 1-only array
    then
      print(1)                       -- if so, we print 1 and exit
      z()
    end
  end

  ::q::                              -- label q, start of the infinite loop
    a={}for i=1,c do a[i]=1 end      -- fill an array with c 1s
    b={}for i=1,#l do b[i]=l[i]end   -- copy the sequence array
    while(a[#a]<b[1])                -- while max(a)<min(b)
    do
      x=0 for i=(#a-c+1>0            -- iterate from index a.length-c to
                    and #a-c+1       -- to a.length
                    or 1),#a
      do
        x=x+a[i]                     -- summing a's elements
      end
      a[#a+1]=x                      -- append x to a
      if a[#a]==b[1]then y(b,1)end   -- if x is equal ot a member of the sequence
                                     -- remove it
      if #b<1 then print(c)z()end    -- if b is empty, it means the subset is in a
                                     -- we print c and exit
    end                              -- else we loop again
    c=c+1                            -- with c+1
  goto q                             -- return to the start of this block
end

```

You can try Lua online, and you can copy/paste the following code sample to run some tests. As this function exit when it founds the answer (infinite loop otherwise), you will have to change the index of `test[]` used (don't forget lua is 1-indexed :)).

```
function f(l)c=2 y,z=table.remove,os.exit while(l[1]<2)do y(l,1)if(#l<1)then print(1)z()end end ::q:: a={}for i=1,c do a[i]=1 end b={}for i=1,#l do b[i]=l[i]end while(a[#a]<b[1])do x=0 for i=(#a-c+1>0 and #a-c+1 or 1),#a do x=x+a[i]end a[#a+1]=x if a[#a]==b[1]then y(b,1)end if #b<1 then print(c)z()end end c=c+1 goto q end

test={{1,1,1},
{49, 97},
{55, 89, 144},
{1},
{6765},
{12, 23, 45, 89},
{100, 199}}

print(f(test[1]))

```
> 2 votes

# Answer
# Wolfram Language (Mathematica), 71 bytes

```
(n=1;While[!SubsetQ[LinearRecurrence[s=1~Table~n,s,Max@l],l=#],n++];n)&

```

Try it online!
> 1 votes

---
Tags: code-golf, sequence, fibonacci, subsequence
---


title: Sort the pixels

text:
# Question
Title: Sort the pixels
Your task is to create a program that will, given an input image, create an output image of the same size, where all pixels are ordered by hex value.

Your program may:

* Sort the pixels from left to right and then down or first sort down in columns and then right. In any case, the top left pixel is the smallest, and the bottom right is the largest.
* Use transparency, but this is not required.
* Sort by RGB, but you may use CMY, or any other format with at least 3 values. You can choose what values to sort on. (HSV may give some nice images)
* Use any well-known image format that most computers can open.

Rules:

* Output must be written to disk or be pipeable to a file.
* Input is given as a commandline argument, in the form of a relative path to the image, or piped in from the commandline.
* This is code golf, so shortest code in bytes wins!

# Answer
# Pyth - 10 bytes

Reads image, collapses bitmap, sorts, and then splits up bitmap again, then writes.

```
.wclK'zSsK

```

Doesn't work online for obvious reasons. Takes input as relative path to image file and outputs into `o.png`.

Output from American Gothic:
> 21 votes

# Answer
# JavaScript (ES6), ~~383~~ ~~377~~ 354 bytes

```
f=s=>{d=document,i=new Image,i.src=s,i.onload=$=>{c=d.createElement`canvas`,x=c.getContext`2d`,c.width=w=i.width,c.height=h=i.height,x.drawImage(i,0,0),D=x.getImageData(0,0,w,h),t=D.data,t.set([].concat(...[...t].map((v,i,T)=>i%4?[,,,0]:T.slice(i,i+4)).sort((a,b)=>a.some((v,i)=>k=v-b[i])&&k)).slice(12*w*h)),x.putImageData(D,0,0),d.body.appendChild(c)}}

```

**Runnable demo:**

```
f=function(s) {l=[i=new Image].slice,i.crossOrigin="anonymous",i.src=s,i.onload=function($){d=document,c=d.createElement`canvas`,c.width=w=i.width,c.height=h=i.height,x=c.getContext`2d`,x.drawImage(i,0,0),D=x.getImageData(0,0,w,h),t=D.data,t.set([].concat.apply([],l.call(t).map((v,i)=>i%4?[0,0,0,0]:l.call(t,i,i+4)).sort((a,b)=>a.reduce((A,v,i)=>A||v-b[i],0))).slice(3*t.length)),x.putImageData(D,0,0),d.body.appendChild(c)}}
$("img").click(function() { f(this.src); })
```

```
canvas { padding-right: 4px; } img { cursor: pointer; }
```

```
<script src="https://ajax.googleapis.com/ajax/libs/jquery/2.1.1/jquery.min.js"></script> <img src="https://upload.wikimedia.org/wikipedia/commons/thumb/f/f4/The_Scream.jpg/114px-The_Scream.jpg"> <img src="https://upload.wikimedia.org/wikipedia/commons/thumb/6/67/Claude_Monet_La_Grenouill%C3%A9re.jpg/193px-Claude_Monet_La_Grenouill%C3%A9re.jpg"> <img src="https://upload.wikimedia.org/wikipedia/commons/thumb/c/cc/Grant_Wood_-_American_Gothic_-_Google_Art_Project.jpg/120px-Grant_Wood_-_American_Gothic_-_Google_Art_Project.jpg"><br>
```

How this code works is to use `getImageData` to get an array of the form

```
[R,G,B,A,
 R,G,B,A,
 R,G,B,A,
 ...]

```

And `map` it to an array of the form

```
[[R,G,B,A],[0,0,0,0],[0,0,0,0],[0,0,0,0],
 [R,G,B,A],[0,0,0,0],[0,0,0,0],[0,0,0,0],
 [R,G,B,A],[0,0,0,0],[0,0,0,0],[0,0,0,0],
 ...]

```

So that the R values are mapped to arrays of the RGBA set, and the B, G, and A values turn into minimum-value zero arrays. When we sort this array, all the `[0,0,0,0]` arrays sort to the bottom and the real-value arrays sort normally at the top:

```
[[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],..., [R,G,B,A],[R,G,B,A],[R,G,B,A],...]
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
               extract & flatten these sorted pixels

```

We skim the top fourth of the array (to lose the empty values we created), flatten it with `[].concat.apply`, and end up with an array of the first form again, but this time, it's sorted.

Slightly de-golfed with whitespace and comments:

```
f=s=>{
  // first, load image, then do everything else onload
  i=new Image,
  i.src = s,
  i.onload=$=>{
    // set up the canvas
    d=document,
    c=d.createElement`canvas`,
    w=c.width=i.width,
    h=c.height=i.height,
    x=c.getContext`2d`,

    // draw image to canvas and capture pixel data
    x.drawImage(i,0,0),
    D=x.getImageData(0,0,w,h),
    t=D.data,

    // set pixel data to...
    t.set(
      // the flattened array...
      [].concat(...
        // that is a mapping of the pixel data...
        [...t].map(
          // that clumps RGBA families into subarrays
          // by replacing every fourth value with [R,G,B,A]
          // and all other values to [0,0,0,0]...
          (v,i,T)=>i%4?[,,,0]:T.slice(i,i+4)
        )
        // and is then sorted...
        .sort(
          // by reducing each array to a positive, negative, or zero
          // by comparing R,G,B,& A until the first nonzero difference
          (a,b)=>a.some((v,i)=>k=v-b[i])&&k
        )
      )
      // then eliminate the low-sorted empty [0,0,0,0] values we created,
      // leaving only the top fourth, with real values
      // (note that 3*4*w*h is the same as 3*t.length)
      .slice(3*4*w*h)
    ),

    // now that `t` is set, store `D` in canvas
    x.putImageData(D,0,0),

    // show canvas
    d.body.appendChild(c)
  }
}

```

Note that most browsers may fail to run this code for large images, because it passes a huge number of arguments into `[].concat`. When the browser environment does not allow sufficient memory for all the arguments, an alternative approach is to re-map the RGBA values out of the top fourth arrays back onto the array, for a total score of **361 bytes**:

```
f=s=>{d=document,i=new Image,i.src=s,i.onload=$=>{c=d.createElement`canvas`,x=c.getContext`2d`,c.width=w=i.width,c.height=h=i.height,x.drawImage(i,0,0),D=x.getImageData(0,0,w,h),t=D.data,t.set([...t].map((v,i,T)=>i%4?[,,,0]:T.slice(i,i+4)).sort((a,b)=>a.some((v,i)=>k=v-b[i])&&k).map((v,i,A)=>A[3*w*h+(i>>2)][i%4])),x.putImageData(D,0,0),d.body.appendChild(c)}}

```

We simply replace the `[].concat(...{stuff}).slice(12*w*h)` with `{stuff}.map((v,i,A)=>A[3*w*h+(i>>2)][i%4])`.)
> 19 votes

# Answer
# Mathematica ~~86 83~~ 72 bytes

```
 f=Flatten;Image[f[Sort[f[s=ImageData@#1,1]]]~ArrayReshape~Dimensions@s]&

```

---

With 14 bytes saved thanks to @Martin Buttner.

---

**Example**

The image itself is input. Alternatively, a variable containing the image could be used.
> 13 votes

# Answer
# C (using SDL1.2), ~~333~~ ~~322~~ 315 bytes

C is most probably not the 'sharpest knife in the shelf' for this kind of work, I wanted to give it a try anyway. Tips for improving my answer are welcome. The program gets the name of the input image file as a cli argument.

```
#include <SDL.h>
#include <SDL_image.h>
#define X SDL_Surface*
#define Y const void*
C(Y a,Y b){return*(Uint32*)a-*(Uint32*)b;}main(int o,char**a){X i=IMG_Load(a[1]);X s=SDL_SetVideoMode(i->w,i->h,32,0);i=SDL_ConvertSurface(i,s->format,0);qsort(i->pixels,i->w*i->h,4,C);SDL_BlitSurface(i,0,s,0);for(;;SDL_Flip(s));}

```

compile and run : `gcc -I/usr/include/SDL snippet.c -lSDL -lSDL_image && ./a.out`

I don't usually golf in C, but I've just answered this challenge yesterday and I just wanted to continue playing with that new toy :)

thanks to @pseudonym117 for helping me saving 5 bytes
> 5 votes

# Answer
# Javascript ES6, 334 bytes

```
f=s=>{with((d=document).body.appendChild(c=d.createElement`canvas`).getContext`2d`)(i=new Image).src=s,drawImage(i,0,0,w=c.width=i.width,h=c.height=i.height),t=(g=getImageData(0,0,w,h)).data,t.set([...t].map(i=>(0+i.toString(16)).slice(-2)).join``.match(/.{8}/g).sort().join``.match(/../g).map(i=>parseInt(i,16))),putImageData(g,0,0)}

```

Ungolfed:

```
f=s=>{                                   // create function that accepts image name
 with((d=document).body.appendChild(     // use "with" to exclude having to prepend "<context2d>." to drawImage, getImageData and putImageData
   c=d.createElement`canvas`).getContext`2d`) // create canvas to get pixels from and draw output to
  (i=new Image).src=s,                   // create image, define source filename
  drawImage(i,0,0,w=c.width=i.width,     // draw image to canvas
                  h=c.height=i.height),
  t=(g=getImageData(0,0,w,h)).data,      // get image data from canvas in form of Uint8Array
  t.set([...t]                           // convert image data from Uint8Array to standard array
   .map(i=>(0+i.toString(16)).slice(-2)) // convert R,G,B,A bytes to base16 strings with leading zeros
   .join``.match(/.{8}/g)                // convert array of [R,G,B,A,R,G,B,A,...] to [RGBA,RGBA,...]
   .sort()                               // sort pixel values
   .join``.match(/../g)                  // convert array of [RGBA,RGBA,...] to [R,G,B,A,R,G,B,A,...]
   .map(i=>parseInt(i,16))),             // convert hex strings back to integers, reassign to image data
  putImageData(g,0,0)                    // dump image data onto canvas
}

```
> 5 votes

# Answer
## JavaScript (ES6), 452 ~~480~~ ~~484~~ ~~487~~ ~~511~~ bytes

Wow, this got longer than expected:

```
f=u=>{i=new Image;i.src=u;i.onload=_=>{c=(d=document).createElement`canvas`;c.width=w=i.width;c.height=h=i.height;x=c.getContext`2d`;x.drawImage(i,0,0,w,h);p=x.getImageData(0,0,w,h).data;t=[];f=[];for(j=0;j<p.length;++j)t.push([p[j],p[++j],p[++j],p[++j]]);t.sort((a,b)=>a[0]>b[0]||a[0]==b[0]&&a[1]>b[1]||a[0]==b[0]&&a[1]==b[1]&&a[2]>b[2]).map(u=>f.push.apply(f,u));x.putImageData(new ImageData(new Uint8ClampedArray(f),w,h),0,0);d.body.appendChild(c)}}

```

The function takes an URL as its input `f('test.jpg');` and draws the result into a `canvas`-element which is appended to the `body`.

Note that the source must be on the same domain or the script will halt with a security issue.

---

**Limitations**

I've tested it in Firefox 42 on OS X (10.10) on a machine with 2.5 GHz i7 and 16 GB RAM. The maximum image size I could process without Firefox asking to continue the script execution was *1600 x 1932 px*.

---

**Ungolfed**

```
f = u => {
    i = new Image;
    i.src = u;
    i.onload = _ => {
        c = (d = document).createElement`canvas`;
        c.width = w = i.width;
        c.height = h = i.height;

        x = c.getContext`2d`;
        x.drawImage(i, 0, 0, w, h);

        p = x.getImageData(0, 0, w, h).data;

        t = [];
        f = [];

        for (j = 0; j < p.length;++j)
            t.push([p[j], p[++j], p[++j], p[++j]]);

        t.sort( (a,b) => a[0] > b[0] || a[0] == b[0] && a[1] > b[1] || a[0] == b[0] && a[1] == b[1] && a[2] > b[2] )
         .map(u => f.push.apply(f,u));

        x.putImageData( new ImageData( new Uint8ClampedArray(f), w, h), 0, 0);
        d.body.appendChild(c)
    }
}

```

---

**Output**

For a better comparison I also took the "*American Gothic*" as the example source:

---

**Edits**

* *Saved* **24 bytes** by using `for (a in b)` instead of `for(;;)`. Thanks to ar34z
* *Saved* **3 bytes** by storing `document` in a variable.
* *Saved* **4 bytes** by dropping some of these `()`.
* *Saved* **10 bytes** by using tagged template strings, omitting the `()` on object creation and removing another pair of redundant `()`. Thanks to apsillers.
* *Saved* **14 bytes** by massive refactoring of the code that flattens the color array after sorting. Thanks a lot to apsillers and Ypnypn for undercutting each other.
* *Saved* **1 byte** by refactoring the `for`-loop that gets the colors of each pixel.
> 4 votes

# Answer
# Bash + GNU utils, 80

```
s()(sed 1d $1|cut -d\  -f$2)
sed 1q $1
s $1 2-|sort -t_ -k1.16|paste <(s $1 1) -

```

This assumes input/output format is in ImageMagick pixel enumeration .txt format. The input is passed as a filename and the output goes to STDOUT.

---

If the above is not considered a well-known image format, then we can add in the necessary conversions:

# Bash + GNU utils + ImageMagick, 108

```
s()(sed 1d t|cut -d\  -f$1)
convert $1 txt:t
(sed 1q t
s 2-|sort -t_ -k1.16|paste <(s 1) -)|convert txt:- $2

```

Input and output are specified as filenames. ImageMagick determines what file formats to use by the passed file extensions, so we can use any common ones:

```
$ ./sortpixels.sh 398px-Grant_Wood_-_American_Gothic_-_Google_Art_Project.jpg o.png
$

```

The resulting o.png looks like:
> 4 votes

# Answer
# Java, 316 bytes

```
import javax.imageio.*;class C{public static void main(String[]a)throws Exception{java.awt.image.BufferedImage i=ImageIO.read(new java.io.File(a[0]));int w=i.getWidth(),h=i.getHeight(),v[]=i.getRGB(0,0,w,h,null,0,w);java.util.Arrays.sort(v);i.setRGB(0,0,w,h,v,0,w);ImageIO.write(i,"png",new java.io.File("a.png"));}}

```

Places the hexadecimal values of the pixel colors in an array. The array is sorted, and the colors are remapped to the pixels in the image. The name of the resulting image is `a.png`.
> 3 votes

# Answer
# Python 2, 128 bytes

```
from PIL import*
a=Image.open('a')
b=a.load()
c,d=a.size
a.putdata(sorted(b[e,f]for f in range(d)for e in range(c)))
a.save('b')

```

Provided the image is a file named `a` with no extension, the output will be a file named `b` with no extension.
> 3 votes

# Answer
# SmileBASIC, ~~39~~ 35 bytes

Assuming the image is loaded onto the 512\*512 graphics page:

```
DIM A[0]GSAVE A,0SORT A
GLOAD A,0,1

```

Explained:

```
DIM IMG[0] 'create array
GSAVE IMG,0 'save graphics to array
SORT IMG 'sort
GLOAD IMG,0,1 'load graphics from array

```

It's that simple! ~~Unfortunately we have to use integers, which adds 4 bytes to the program size due to the type suffixes.~~
> 3 votes

# Answer
# Java, ~~424~~ ~~417~~ 404 bytes

Well, this is not a language you want to golf in...

```
import java.awt.image.*;import java.io.*;import javax.imageio.*;class F{public static void main(String[]x)throws Exception{BufferedImage i,o;i=ImageIO.read(new File(x[0]));o=new BufferedImage(i.getWidth(),i.getHeight(),BufferedImage.TYPE_INT_RGB);o.setData(i.getRaster());int[]p=((DataBufferInt)o.getRaster().getDataBuffer()).getData();java.util.Arrays.sort(p);ImageIO.write(o,"png",new File("o.png"));}}

```
> 2 votes

# Answer
# C#, 497 bytes

First time post, first golf. Clearly not the best for golfing

Not really respecting the piping. Takes an image path as an input and outputs it with letter "o" prepend to the name.

Works better with bitmaps, odds results with others

```
using System.Linq;using System.Drawing;using System.Runtime.InteropServices;class Program{static void Main(string[]args){using(var im=(Bitmap)Image.FromFile(args[0])){int h=im.Height;int w=im.Width;var b=im.LockBits(new Rectangle(0,0,w,h),System.Drawing.Imaging.ImageLockMode.ReadWrite,System.Drawing.Imaging.PixelFormat.Format32bppRgb);var p=new int[h*w];Marshal.Copy(b.Scan0,p,0,h*w);var q=p.ToList();q.Sort();p=q.ToArray();Marshal.Copy(p,0,b.Scan0,h*w);im.UnlockBits(b);im.Save("o"+args[0]);}}}

```
> 2 votes

# Answer
## Haskell, 195 bytes

```
import Data.List
import Graphics.GD
f p=do
 a<-loadPngFile p;(x,y)<-imageSize a;let l=[(i,j)|j<-[0..y],i<-[0..x]]
 mapM(flip getPixel a)l>>=mapM(\(d,c)->setPixel d c a).zip l.sort;savePngFile"o"a

```

This uses the `GD` library. Usage `f <filename>`. The input file must be in `png` format. The output file is named `o`.

How it works: straightforward, i.e. read the picture, walk over all coordinates and get the pixels, sort the pixels, walk over the coordinates again, but this time set the pixels in the order they appear in the sorted list, write the file to disk.
> 1 votes

# Answer
# Java + Processing, 198 bytes

```
PImage p;void settings(){p=loadImage("peter-olexa-Ejgx4k3iryg-unsplash.jpg");size(p.width,p.height);}void setup(){p.loadPixels();p.pixels=sort(p.pixels);p.updatePixels();image(p,0,0);save("2.jpg");}

```

Outputs image as `1.jpg`.
> 1 votes

# Answer
# Funky2, 156 Bytes

```
I=draw.newImageData(arg[1)C=[]j=i=0I.mapPixels((x,y,r)=>C[i++]=r)C=C:sort((a,b)=>(a.r*256+a.g+a.b/256)-(b.r*256+b.g+b.b/256))I.mapPixels(@C[j++)I.export"o"

```

Non-Competing due to the fact that several features were implemented in response to this challenge, including image exporting and list sorting. Also several features needed fixing.

## Ungolfed

```
I = draw.newImageData(arg[1]); $// Load the specified image as a bitmap.
c = []; $// Generate a table for us to pour colours into.
j=i=0;  $// Initialize two counters, I and J. Saves us resetting i, and c[#c] is prohibitively memory intensive
I.mapPixels((x,y,r)=>c[i++]=r); $// Iterate over the image and assign each colour (Contained as a list in r) into c.
c = c:sort((a,b)=>(a.r*256+a.g+a.b/256)-(b.r*256+b.g+b.b/256)); $// Sort C by R, then G, then B.
I.mapPixels(@c[j++]); $// Map back over I, replacing each colour with their index in C
I.export"o"; $// Export as "o"

```

There are a couple of opportunities for improvement; such as using `c[#c]` rather than `c[i++]` or Sorting `c` every time I try to index it in the second mapping; but these methods cause the execution time to skyrocket.

**Output, American Gothic**
> 1 votes

---
Tags: code-golf, graphical-output, image-processing
---


# Question
Title: Can I exit the maze with this many bombs?
Inspired by this Puzzling challenge, and easier version of my previous challenge.

## Challenge

A 2D rectangular grid is given, where each cell is either an empty space or a wall. You start at the top left cell, and you need to exit through the bottom right cell. You can move to one of four adjacent cells in one step.

You have some bombs, so that using one bomb will let you break exactly one cell-sized wall and go through it. Can you exit the maze using just what you have?

## Input and output

The input is the maze and the initial number of bombs. The maze can be taken as a matrix (or any equivalent) containing two distinct values to represent empty spaces and walls. The top left and bottom right cells are guaranteed to be empty. The number of bombs `n` is always a non-negative integer.

The output should be truthy if you can exit the maze using `n` or fewer bombs, falsy otherwise. You can output truthy/falsy using your language's convention (swapping is allowed), or use two distinct values to represent true or false respectively.

Standard rules apply. The shortest code in bytes wins.

## Test cases

Uses `.#` for spaces and walls.

```
Input maze
..#..#..
Output: false (for 0 or 1 bomb), true (≥2 bombs)

Input maze
.....
####.
.....
.####
.....
Output: true (for any bombs)

Input maze
.
Output: true (for any bombs)

Input maze
.#.#.
##.##
.###.
Output: false (for ≤2 bombs), true (≥3 bombs)

Input maze
.####
#####
##.##
#####
####.
Output: false (for ≤5 bombs), true (≥6 bombs)

```

# Answer
> 5 votes
# JavaScript (ES7),  ~~126 125~~  119 bytes

Expects `(matrix)(bombs)`, where the `matrix` is filled with `-1` for an empty cell and `-2` for a wall.

Returns *false* if we can exit the maze, or *true* if we can't.

```
m=>g=(b,X=0,Y=0)=>m.every((r,y)=>m[Y+1]||r[X+1]?r.every((v,x)=>r[x]*=v>0|(X-x)**2+(Y-y)**2!=1||g(b-~v,x,y,r[x]=1)):b<0)

```

Try it online!

### Commented

```
m =>                        // m[] = matrix
g = (                       // g is a recursive function taking:
  b,                        //   b = number of bombs
  X = 0, Y = 0              //   (X, Y) = current position, starting at (0, 0)
) =>                        //
  m.every((r, y) =>         // for each row r[] at position y in m[]:
    m[Y + 1] ||             //   if there's a row below the current cell
    r[X + 1] ?              //   or there's a column on the right:
      r.every((v, x) =>     //     for each value v at position x in r[]:
        r[x] *=             //       restore r[x] if any of these tests is true:
          v > 0 |           //         - v is greater than 0 (this cell was visited)
          (X - x) ** 2 +    //         - the squared distance between
          (Y - y) ** 2 != 1 //           (x, y) and (X, Y) is not equal to 1
          ||                //
          g(                //         - this recursive call is truthy:
            b - ~v,         //             decrement b if v = -2
            x, y,           //             use the new position (x, y)
            r[x] = 1        //             mark r[x] as visited by setting it to 1
          )                 //           end of recursive call
      )                     //     end of inner every()
    :                       //   else (bottom-right cell):
      b < 0                 //     return true if we've used too many bombs
  )                         // end of outer every()

```

# Answer
> 3 votes
# Charcoal, 81 bytes

```
≔⟦⟧θＷＳ⊞θι⊞υ⟦⊕Ｎω⟧≔⁰ηＦυ«⪫θ¶←Ｆ§ι¹✳κ+¿∨ⅈⅉＦruld«≔⌕….#§ι⁰∨⊟ＫＤ²✳κ+ζ¿⊕ζ⊞υＥι⁺λ⎇μκ±ζ»≔¹η⎚»η

```

Try it online! Link is to verbose version of code. Based on my answer to the previous challenge. Works better on grids with lots of walls. Bomb count is separated from the grid by a blank line. Outputs a Charcoal boolean, i.e. `-` for a path, nothing if not. Explanation:

```
≔⟦⟧θＷＳ⊞θι

```

Input the grid.

```
⊞υ⟦⊕Ｎω⟧

```

Start with an initial state of `n+1` bombs and no moves. (This is because the algorithm stops when you run out of bombs, rather than when you need a bomb to move.)

```
≔⁰η

```

We haven't found a path yet.

```
Ｆυ«

```

Perform a breadth-first search of the states.

```
⪫θ¶←

```

Draw the input to the canvas, leaving the cursor at the end point.

```
Ｆ§ι¹✳κ+

```

Draw the path so far.

```
¿∨ⅈⅉ

```

If the start hasn't been reached, then:

```
Ｆruld«

```

Loop over the orthogonal directions.

```
≔⌕….#§ι⁰∨⊟ＫＤ²✳κ+ζ

```

Look at the next character in that direction to see how many bombs we need (`-1` for an illegal move, including running out of bombs).

```
¿⊕ζ⊞υＥι⁺λ⎇μκ±ζ

```

If the move is legal then create a new state by subtracting the number of bombs and adding the current direction.

```
»≔¹η

```

But if the start was reached, then record that we found a path.

```
⎚»

```

Clear the canvas ready for the next state (or the final output).

```
η

```

Output the flag for whether we found a path.

# Answer
> 3 votes
# APL (Dyalog Unicode), 34 or 37 bytes

```
⎕≥⊃⌽,(0@0@0⊢⌊⎕+(⍉g∘⍉)⌊g←3⌊/,,⊣)⍣≡⍨9e9

```

Try it online!

this shorter version works in dyalog v18 but not on tio:

```
⎕≥⊃⌽,(0@0@0⊢⌊⎕+g⍤1⌊g←3⌊⌿⍪⍪⊣)⍣≡⍨9e9

```

inputs

`9e9` a very large number, used as a substitute for infinity

`(` `)⍣≡⍨9e9` apply the function train in parens until convergence, using `9e9` both as a constant always passed on the left, and a starting value initially passed on the right

`g←3⌊/,,⊣` helper function to compute the minimum of each cell and its two horizontal neighbours, using `9e9` for the boundary around the matrix

`(⍉g∘⍉)` same for vertical - this is `g` under transposition

`⎕+`.... min of horizontals and verticals, and add the original matrix (this accounts for the cost of 1 bomb when we encounter a wall)

`⊢⌊`.. update the matrix of best known path costs

`0@0@0` put a 0 in the top left cell

on the first iteration of `(` `)⍣≡`, the scalar `9e9` is extended to a matrix (the matrix of best costs) because of `⎕+`, and then it remains a matrix until the end.

`⊃⌽,` lower right cell

`⎕≥` compare with the number of bombs available

# Answer
> 2 votes
# Java (JDK), ~~235~~ ~~233~~ ~~227~~ ~~222~~ 211 bytes

```
int c(int[][]m,int x,int y,int b){int a=0,v;try{m[x][y]=(b-=v=m[x][y])*v<0?v/0:-1;a+=(x==m.length-1&y==m[0].length-1?1:0)+c(m,x+1,y,b)+c(m,x-1,y,b)+c(m,x,y+1,b)+c(m,x,y-1,b);if(a<1)m[x][y]=v;}finally{return a;}}

```

Try it online!

Requires an `int[][]` with `0` as field and `1` as wall.
Returns `0` on failure and `1` on success. I'm unsure if this is a vaild truthy/falsy value for Java though.

A rather simple approach: Walk around and bomb walls until you reach the exit or run out of bombs.

I removed the explaination, it got too messy for me to update because of line length.

### EDIT:

-2 thanks to ceilingcat!
-4, again thanks to ceilingcat!
-2 by optimising the goal check
-5, again thanks to ceilingcat! They also fixed the awful formatting in the TIO link.
-11 thanks to Kevin Cruijssen!

# Answer
> 2 votes
# Python 3, 130 bytes

```
def f(g,b,x=0,c=0):w=len(g[0])+1;l=w*len(g);return~x%w*(b>-1<x<l>c)and any(f(g,b-g[x//w][x%w],x+a,c+1)for a in(1,-1,w,-w))or-~x==l

```

Try it online!

Recursive function to find all paths. Takes a 2D matrix as input, with `0` for empty spaces and `1` for walls. The number of bombs `b` is reduced by 1 each time it encounters a wall. Recursion stops immediately when the edge of the grid `g` is detected, more steps `c` have been taken than the size `l` of the grid, or the number of bombs remaining falls below zero. Returns `True` when any of the paths reaches the final space and `False` otherwise.

Adaptation from my answer to Find the shortest route on an ASCII road.

# Answer
> 1 votes
# Javascript, 186 bytes

Is a function `(maze, width, height, bombs) => boolean` returning whether or not the maze can be solved with the given number of bombs. The `maze` should be supplied as a flatten list of booleans, `true` for walls and `false` for empty spaces.

```
(m,w,h,b)=>{s=Array(w*h).fill(1/0);i=d=s[0]=0;l:for(;;){for(i=0;i<w*h;i++)for(d of[-w,-1*!!(i%w),1*!!((i+1)%w),w])if(s[i+d]+m[i]<s[i]){s[i]=s[i+d]+m[i];continue l}return s[w*h-1]<=b;}}

```

Try it online!

Sadly, I wasn't able to get this below the other JS answer. I tip my hat to @Arnauld and look forward to reading how his works.

Degolfed and annotated:

```
S = (m, w, h, b) => {

s = Array(w*h).fill(1/0);  // initialize the scoreboard to infinity the scoreboard
                           // .. which holds the running minimum for number of
                           // .. bombs required to reach a certain grid cell
i = d = s[0] = 0;          // declare variables i and d and note on the scoreboard
                           // .. that we can reach the top-left cell with 0 bombs
l: for(;;) {               // repeat infinitely
 for (i = 0; i < w*h; i++) // loop over all grid cells
 for (d of [-w,            // for direction of [up,
           -1*!!(i%w),     // left, (note: if the cell is at the start of a row
                           // .. then -1 could wrap; handle this with `*!!(i%w)`)
           1*!!((i+1)%w),  // right, (likewise here for the end of a row)
           w])             // down].
 if (s[i+d] + m[i]<s[i]) { // if moving from the given direction onto this cell
                           // .. would take less bombs than what's currently in
                           // the scoreboard,
   s[i] = s[i + d] + m[i]; // then update the scoreboard
   continue l              // we've made a change to the scoreboard, so ensure we
                           // .. don't reach the below `return`
 }
 return s[w * h - 1] <= b; // return the score value for the bottom-right cell.
                           // .. due to the above `continue`, this statement will
                           // .. only be reached once no more changes to the
                           // .. scoreboard can be made
}

}

```

# Answer
> 1 votes
# Wolfram Language (Mathematica), 68 bytes

```
Last@GraphDistance[GridGraph[#2,EdgeWeight->{_b_:>#[[b]]}],1]>#3&

```

Try it online!

Returns `True` if there are not enough bombs, and `False` otherwise. Takes `[maze, {w,h}, bombs]`, where `maze` is a 1d list of `0`s (no wall) and `1`s (wall).

---
Tags: code-golf, decision-problem, grid, path-finding, maze
---


title: Stroke Count of a Chinese Numeral

text:
# Question
Title: Stroke Count of a Chinese Numeral
## Task

Given an integer \\$n\in\[0,10^{12})\\$ in any convenient format, return the number of strokes needed to write that character in simplified Chinese.

### Background

Chinese numerals are expressed in base 10 with a system of digits and places, with an important distinction that digits are in groups of four, rather than three.

The individual characters used to write Chinese can be described at the lowest level as a collection of strokes, laid out in a certain order and manner. The number of strokes required to write a character is that character's stroke count.

The (simplified) characters used to write numbers in Chinese are:

```
num   char　 strokes
0     零*    13
1     一     1
2     二**   2
3     三     3
4     四     5
5     五     4
6     六     4
7     七     2
8     八     2
9     九     2
10    十     2
100   百     6
1000  千     3
10^4  万     3
10^8  亿     3
 * 0 can also be written 〇, but we won't use that here.
** 两 is largely interchangeable with 二, apart from never
   appearing before 十. We won't consider it here for simplicity,
   but 两 is very common in actual usage.

```

For example, 9 8765 4321 is 九亿八千七百六十五万四千三百二十一: nine hundred-million (九 亿), eight thousand seven hundred sixty-five ten-thousand (八千七百六十五 万), four thousand three hundred twenty-one (四千三百二十一). In all, 53 strokes are needed to write this out.

There are additionally some special rules involving the digits `0` and `1`. These can vary slightly between dialects, but we'll choose these:

* When there are non-trailing `0`s in a 4-digit group, they are combined into a single . No place marker is used.
(This is because e.g. 一百二 is a common way to say 120. We won't consider that form.)

+ 1020 is 一千零二十.
    + 6 0708 is 六万零七百零八.
    + 3 0000 4005 is 三亿四千零五.
    + 0 is 零.
* If the number would begin with `一十`, the  is omitted.

+ Powers of 10 are 一, 十, 一百, 一千, 一万, 十万, 一百万, etc.
    + 111 is 一百一十一.

---

## Test Cases

```
           n  strokes  chinese
           0  13       零
          10  2        十
         236  17       二百三十六
        7041  26       七千零四十一
       50010  23       五万零一十
      100000  5        十万
      860483  42       八十六万零四百八十三
     4941507  52       四百九十四万一千五百零七
    51001924  38       五千一百万一千九百二十四
   105064519  70       一亿零五百零六万四千五百一十九
   300004005  31       三亿四千零五
   987654321  53       九亿八千七百六十五万四千三百二十一
  1240601851  56       十二亿四千零六十万一千八百五十一
608726402463  79       六千零八十七亿二千六百四十万二千四百六十三

```

@user202729 provided a script in the sandbox to help with reading the Chinese numbers: Try it online!

# Answer
## Python 2, ~~163~~ ~~162~~ ~~155~~ 151 bytes

<sup>-1 thanks to Arnauld, -7 thanks to ovs</sup>

Hexdump:

```
00000000: efbb bf65 7865 6322 789c 1dcb c10e c220  ...exec"x......
00000010: 0cc6 f157 e1c2 d2c2 88b4 6359 34c2 bb20  ...W......cY4..
00000020: b893 7689 274f 3ebb 65a7 36bf 7cff 3dbf  ..v.'O>.e.6.|.=.
00000030: eafb d1ab 91b9 e578 1353 a543 fc32 734a  .......x.S.C.2sJ
00000040: ebc2 d44b 114b d1a5 8956 f430 7e5d 9866  ...K.K...V.0~].f
00000050: 93bb 5a62 bba1 5fdc e067 1ace 1b97 d226  ..Zb.._..g.....&
00000060: c240 ca40 67ec 75ad 2de8 b92f 7852 440c  .@.@g.u.-../xRD.
00000070: 2039 0f62 f43b c885 e21c 7e5c 728f 8fd1   9.b.;....~\r...
00000080: 343a d7fe d154 254b 222e 6465 636f 6465  4:...T%K".decode
00000090: 2827 7a69 7027 29                        ('zip')

```

Unpacked:

```
f=lambda n,c=0:n and(0x222445321d>>n%10*4&15)+(n%10and c%4*9%12%7)+3*(n%1e4and 272>>c&1)-13*((1>n%10+c%4)+((c%4<3)>n%100))-(n==c%4<2)+f(n/10,-~c)or 13*0**c

```

Try it online!

### Commented version

```
f=lambda n,c=0:n and        # If n is not zero, sum the following:
                            # The strokes of the current digit
  (0x222445321d>>n%10*4&15) #   [13,1,2,3,5,4,4,2,2,2][n%10]
                            # The 10/100/1000 marker stroke count, if
                            # the currrent digit is zero
  + (n%10and c%4*9%12%7) #   [0, 2, 6, 3][n%10 and c%4]
                            # The 10^4 / 10^8 marker stroke count
  + 3*(n%1e4and 272>>c&1)   #   3*(n%10000 and c in [4,8])
  - 13*(                    # Subtract 13 because we overcounted the zeros if
                            #   It's a trailing zero in a group of four
        (1>n%10+c%4)        #     (n%10 == 0 and c%4 == 0)
                            #   and/or if it's part of a group of zeros
        +((c%4<3)>n%100)    #     (n%100 == 0 and c%4 != 3)
    )                       # Subtract 1 if we lead a group of four with 10
  - (n==c%4<2)              #   (n==1 and c%4==1)
  + f(n/10,-~c)             # Recurse on next digit
                            # Else if n=0, return 13 if we are at the first
or 13*0**c                  # digit, else 0

```

We also use the `exec"...".decode('zip')` to compress our code to save a few bytes.
> 10 votes

# Answer
# JavaScript (ES6),  ~~148 145~~  139 bytes

```
f=(n,i=0,z=n/1e4|0)=>n?f(z,3)-!(n/10^1)+(g=d=>~~d?n%1e4&&(q=n/d%10,x=296/~~q%5^20>>q&1,x?z=x+=3064/d&7:z*q&&(z=0,13))+g(d/10):i)(1e3):!i*13

```

Try it online!

### Black magic

* The number of strokes in the Chinese character for a single decimal digit \\$1\le q\le 9\\$ can be obtained with:

```
296 / q % 5 ^ 20 >> q & 1

```

Try it online!

This formula returns \\$0\\$ for \\$q=0\\$, which is actually what we want because zeros have to be processed separately anyway.
* The number of strokes in the Chinese character for a multiplier \\$d\in\\{1,10,100,1000\\}\\$ can be obtained with:

```
3064 / d & 7

```

Try it online!

### Commented

```
f = (                      // f is a recursive function taking:
  n,                       //   n = input
  i = 0,                   //   i = number of strokes for either 10^4 or 10^8
                           //       (0 for the first iteration, then 3)
  z = n / 1e4 | 0          //   z = next 4-digit block, also used as a flag to
) =>                       //       tell whether zeros must be inserted
n ?                        // if n is not equal to 0:
  f(z, 3)                  //   do a recursive call with the next 4-digit block
  - !(n / 10 ^ 1)          //   subtract 1 if n is the range [10...19] to account
  +                        //   for a leading '10' where the '1' must be omitted
  ( g = d =>               //   g is a recursive function taking a divisor d
    ~~d ?                  //     if d greater than or equal to 1:
      n % 1e4 &&           //       abort if the entire block is 0
      (                    //       otherwise:
        q = n / d % 10,    //         q = (n / d) mod 10
        x = 296 / ~~q % 5  //         compute the number of strokes x for this
            ^ 20 >> q & 1, //         decimal digit (i.e. the integer part of q)
        x ?                //         if it's not a zero:
          z =              //           set z to a non-zero value
            x +=           //           add to x ...
              3064 / d & 7 //           ... the number of strokes for d
        :                  //         else:
          z * q &&         //           if both z and q are non-zero:
            (z = 0, 13)    //             add 13 strokes and set z to 0 so that
      )                    //             other contiguous zeros are ignored
      + g(d / 10)          //       add the result of a recursive call with d / 10
    :                      //     else:
      i                    //       add i to account for either 10^4 or 10^8
  )(1e3)                   //   initial call to g with d = 1000
:                          // else:
  !i * 13                  //   stop the recursion; add 13 if the input was zero

```
> 8 votes

# Answer
# 05AB1E, ~~60~~ 59 bytes

```
Rv•мΛ~?•ÁyèyĀiŽAΘÁNè]IR4ô¤¤sg*2Q(sεÔ0Û}©J0¢I_+13*ŽE₃Á®¦ĀOèO

```

Try it online or verify all test cases.

**Explanation:**

Here a run-down of what each part of the code accomplishes, which we'll sum with `O` at the very end to get the result:

* Convert each digit to the amount of strokes used, except for edge case 0: `Rv•мΛ~?•Áyè]`
* Get the strokes of 10, 100, 1000 per part of four digits if the current digit is not a 0: `RvyĀiŽAΘÁNè]`
* Subtract 1 in case there is a leading `一十` conversion, by separating the number into parts of four digits, and checking if the first part starts with a `1` and contains just two digits: `IR4ô¤¤sg*2Q(`
* Add 13 for each 0, excluding those that are adjacent to other 0s, or trailing in every part of four digits: `IR4ôεÔ0Û}©J0¢13*`
* Take care of edge case \\$n=0\\$: `I_+13*`
* Get the strokes of the \\$10^4\\$ and \\$10^8\\$, excluding parts that consist solely of 0s: `ŽE₃Á®¦ĀOè`

```
R                # Reverse the (implicit) input-integer
 v               # Loop over each of its digits `y`:
  •мΛ~?•         #  Push compressed integer 1235442220
        Á        #  Rotate it once to the right: 0123544222
         yè      #  Index the current digit `y` into 0123544222
  yĀi            #  If `y` is NOT 0:
     ŽAΘ         #   Push compressed integer 2630
        Á        #   Rotate it once to the right: 0263
         Nè      #   Index the loop-index `N` into 0263 (modulair 0-based)
 ]               # Close both the if-statement and loop
IR               # Push the reversed input-integer again
  4ô             # Split it into parts of size 4, where the last part might be smaller
    ¤            # Push the last part (without popping the list itself)
     ¤           # Push the last digit (without popping the number itself)
      s          # Swap so the number is at the top of the stack again
       g         # Pop and push its length
        *        # Multiply it to the digit
         2Q      # Check if length * digit is equal to 2 (1 if truthy; 0 if falsey)
           (     # Negate that
    s            # Swap so the list of parts of 4 is at the top of the stack again
     ε           # Map each part to:
      Ô          #  Connected uniquify each digit (i.e. 0030 becomes 030)
       0Û        #  Remove any leading 0s
     }©          # After the map: store this list in variable `®` (without popping)
       J         # Join it together to a single string
        0¢       # Count the amount of 0s left in the string
          I_     # Check if the input-integer is 0 (1 if 0; 0 otherwise)
            +    # Add that to the count of 0s
             13* # Multiply it by 13
ŽE₃              # Push compressed integer 3660
   Á             # Rotate it once to the right: 0366
    ®            # Push the list from variable `®`
     ¦           # Remove the first item
      Ā          # Truthify each part ("" becomes 0; everything else becomes 1)
       O         # Sum this list
        è        # Index it into the 0366
O                # And finally sum all values on the stack
                 # (after which this sum is output implicitly as result)

```

See this 05AB1E tip of mine (section *How to compress large integers?*) to understand why `•мΛ~?•` is `1235442220`; `ŽAΘ` is `2630`; and `ŽE₃` is `3660`.
> 7 votes

# Answer
## Setanta, ~~243~~ ~~236~~ ~~235~~ 234 bytes

```
gniomh(n){gniomh S(n,c){gniomh Q(t){s:=0ma t s=1toradh s}toradh(n&[13,1,2,3,5,4,4,2,2,2][n%10]+[0,2,6,3][n%10&c%4]+3*Q((c==4|c==8)&n%10000)-13*(Q(n%10==0&c%4==0)+Q(n%100==0&c%4!=3))-Q(n==1&c%4==1)+S(n//10,c+1))|13*Q(!c)}toradh S(n,0)}

```

Try it here!

Notes:

* This is a port of the Python 2 solution.
* The correct keywords for 'function' and 'if' are `gníomh` and , respectively, but Setanta allows spelling them without the accents so I did that to save some bytes.
* I had to define the `Q` function above because Setanta doesn't coerce booleans to integers. As far as I know, there's nothing equivalent to ternary expressions in this language.
* There are no bitwise operations as far as I know, so I had to stick to array lookups.
> 4 votes

---
Tags: code-golf
---


# Question
Title: Easier communication means… ○ circular edition!
Based on this, but I was asked\[1\],\[2\] to do a graphical edition.

# Task

**Generate a circular display of the text `EASIER COMMUNICATION MEANS FASTER CODING MEANS FEWER CODERS MEANS` .** Any uppercase font goes. The letter tops must face out of the circle. Either the space before `FEWER` must be at the bottom (6 o'clock) or (for brownie points) your text must rotate.

# Answer
> 7 votes
# Octave, 168 bytes

```
for i=1:(n=nnz(s='CATION MEANS FASTER CODING MEANS FEWER CODERS MEANS EASIER COMMUNI')),text(sin(i/n*2*pi),cos(i/n*2*pi),s(i),'rotation',-i/n*360);axis([-2,2,-2,2]);end

```

**Try it here**

# Answer
> 7 votes
## Mathematica, 153 Bytes

```
x=Pi/33;Graphics@Map[StringSplit["CATION MEANS FASTER CODING MEANS FEWER CODERS MEANS EASIER COMMUNI",""][[#/x]]~Text~{Sin@#,Cos@#}~Rotate~-#&,Range@66x]

```

Increasing it to 178 bytes lets it rotate:

```
x=Pi/33;Graphics@Map[StringSplit["CATION MEANS FASTER CODING MEANS FEWER CODERS MEANS EASIER COMMUNI",""][[#/x-a]]~Text~{Sin@#,Cos@#}~Rotate~-#&,(a+Range@66)x]~Animate~{a,1,66,1}

```

# Answer
> 6 votes
# HTML and JS, 13 + 170 = ~~284~~ ~~227~~ ~~192~~ ~~186~~ 183 bytes

*41 bytes saved thanks to @Shaggy*

```
(x=c.getContext("2d")).translate(r=c.height/2,r)
for(i in t="CATION MEANS FASTER CODING MEANS FEWER CODERS MEANS EASIER COMMUNI"){x.rotate(2/21);x.fillText(t[i],0,-r*.8)}
```

```
<canvas id=c height=250 width=250>
```

The width and height set to 250 to allow better visibility were not included in the bytes count.

# Answer
> 5 votes
# Bash + ImageMagick, 168 bytes

```
convert -font $(convert -list font|grep Font:|head -n1|cut -d" " -f4) label:"FEWER CODERS MEANS EASIER COMMUNICATION MEANS FASTER CODING MEANS " -distort Arc 360 i.png

```

**Example output:**

# Answer
> 5 votes
# BBC BASIC, 248

Download interpreter at http://www.bbcbasic.co.uk/bbcwin/download.html

```
F.i=6TO767a=i DIV6*.024-.3r=(80+i MOD6)MOD83*9-369IFASC(M."?[O@Wo}@b}O@O|O@y?l@xHH@GGE@zmo@yo|@}oU@o?m@?L|@phh@GEE@?u?@?{o@Wk?@?I?@g|O@EGE@?]M@z}o@|_|@}?m@?MO@GDG@xhh@?k?@?oU@z}o@?LO@a?L",i DIV6))>>i MOD6A.1CIRCLEFILL500-r*COS(a),500+r*SIN(a),9
N.

```

**Ungolfed**

```
  FORi=6TO767
    a=i DIV 6 * .024 - .3
    r=(80 + i MOD 6) MOD 83 * 9 - 369            :REM iterate through ({80,81,82,0,1,2} - 41 = {39,40,41,-41,-40,-39}) * 9
    IF ASC(MID$("?[O@Wo}@b}O@O|O@y?l@xHH@GGE@zmo@yo|@}oU@o?m@?L|@phh@GEE@?u?@?{o@Wk?@?I?@g|O@EGE@?]M@z}o@|_|@}?m@?MO@GDG@xhh@?k?@?oU@z}o@?LO@a?L",i DIV6))>>i MOD6AND1    THEN
      CIRCLEFILL 500-r*COS(a), 500+r*SIN(a), 9   :REM dots are circle radius 9 according to documentation (but BBC BASIC actually draws them smaller.)
    ENDIF
  NEXT

```

Doing this in a language that doesn't support rotated text - for fun. Each character is a 3x3 dot matrix array. There are a total of 66 characters but there are 2 diametrically opposite spaces that we do not print (after CODERS and before FASTER.) Each byte of the magic string encodes for a single vertical row of one letter and a single vertical row of the opposite letter. In this way we only have to sweep through 180 degrees.

# Answer
> 3 votes
## SVG(HTML5), 212 bytes

```
<svg width=500 height=500><defs><path id=p d=M250,451a201,201,0,0,1,0,-402a201,201,0,0,1,0,402></defs><text font-size="32"><textPath xlink:href=#p>FEWER CODERS MEANS EASIER COMMUNICATION MEANS FASTER CODING MEANS
```

Static version of my answer to Easier communication means faster coding means fewer coders means … tweaked so that `FEWER` is at the bottom. Adjustable version in case your font doesn't match mine:

```
<p><input type=number value=0 min=0 max=9 oninput=p.setAttribute('d','M250,250m0,20_a20_,20_,0,1,1,20_,-20_a20_,20_,0,1,1,-20_,-20_a20_,20_,0,1,1,-20_,20_a20_,20_,0,1,1,20_,20_'.replace(/_/g,this.value))></p>
<svg width=500 height=500><defs><path id=p d=M250,250m0,200a200,200,0,1,1,200,-200a200,200,0,1,1,-200,-200a200,200,0,1,1,-200,200a200,200,0,1,1,200,200></defs><text font-size="32"><textPath xlink:href=#p>FEWER CODERS MEANS EASIER COMMUNICATION MEANS FASTER CODING MEANS
```

# Answer
> 2 votes
# Postscript (127 bytes)

Hex dump:

```
00000000: 2f53 2031 92a5 9233 2f4d 6f6e 6f20 3992  /S 1...3/Mono 9.
00000010: 8e88 6432 3030 926b 3292 9a28 4943 4154  ..d200.k2..(ICAT
00000020: 494f 4e20 4d45 414e 5320 4641 5354 4552  ION MEANS FASTER
00000030: 2043 4f44 494e 4720 4d45 414e 5320 4645   CODING MEANS FE
00000040: 5745 5220 434f 4445 5253 204d 4541 4e53  WER CODERS MEANS
00000050: 2045 4153 4945 5220 434f 4d4d 554e 297b   EASIER COMMUN){
00000060: 5388 0033 2032 9287 9278 5320 6661 6c73  S..3 2...xS fals
00000070: 6592 112d 352e 3435 9288 7d92 4992 a7    e..-5.45..}.I..

```

Ungolfed (text) version:

```
/S 1 string def
/Mono 9 selectfont
100 200 moveto
2 setlinejoin
(ICATION MEANS FASTER CODING MEANS FEWER CODERS MEANS EASIER COMMUN)
{S 0 3 2 roll put S false charpath -5.45 rotate} forall
stroke

```

Sample output:

# Answer
> 2 votes
# Java 8, 1,087 bytes

```
import javafx.animation.*;import javafx.application.*;import javafx.scene.*;import javafx.scene.layout.*;import javafx.scene.shape.*;import javafx.scene.text.*;import javafx.stage.*;import javafx.util.*;import java.util.*;import static javafx.animation.PathTransition.*;public class P extends Application{@Override public void start(Stage stage)throws Error{char[] t="EASIER COMMUNICATION MEANS FASTER CODING MEANS FEWER CODERS MEANS".toCharArray();Shape q=new Circle(400,400,120);List<Text>h=new ArrayList<>();List<Transition>z=new ArrayList<>();for(char c:t){Text o=new Text(c+"");h.add(o);z.add(x(q, o));}Pane p=new Pane();p.getChildren().addAll(h);int s=h.size();for(int i=0;i<s;i++){Transition w=z.get(i);w.jumpTo(Duration.seconds(10).multiply((i+.5)*1/s));w.play();}stage.setScene(new Scene(p,800,800));stage.show();}Transition x(Shape e,Text t){PathTransition v=new PathTransition(Duration.seconds(10),e,t);v.setCycleCount(INDEFINITE);v.setOrientation(OrientationType.ORTHOGONAL_TO_TANGENT);v.setInterpolator(Interpolator.LINEAR);return v;}static void main(String[]a){launch(a);}}

```

# Answer
> 1 votes
# LaTeX + TikZ, 380 bytes

```
\documentclass[tikz,border=20pt]{standalone}\usetikzlibrary{decorations,decorations.text,}\begin{document}\begin{tikzpicture}\draw[color=white,rotate=-90,postaction={decorate,decoration={text along path,raise=4pt,text align={align=center},text={FEWER CODERS MEANS EASIER COMMUNICATION MEANS FASTER CODING MEANS },reverse path}}](0,0)circle(2.425cm);\end{tikzpicture}\end{document}

```

Example output

# Answer
> 0 votes
# Tcl/Tk, 222

```
grid [canvas .c]
set i 90;lmap c [split "FEWER CODERS MEANS EASIER COMMUNICATION MEANS FASTER CODING MEANS" ""] {.c cr t [expr {93*cos([set i [expr $i+5.5]]/57.3)+99}] [expr {93*sin($i/57.3)+99}] -te $c -ang [expr -$i-90]}

```

Will golf it more later!

# Answer
> 0 votes
# Java 10, 690 bytes

```
import java.awt.*;v->new Frame(){{add(new Panel(){RenderingHints R;Math M;public void paint(Graphics g){Graphics2D G=(Graphics2D)g;G.setRenderingHint(R.KEY_ANTIALIASING,R.VALUE_ANTIALIAS_ON);var v=new Font("Hei",0,8).createGlyphVector(G.getFontRenderContext(),"ICATION MEANS FASTER CODING MEANS FEWER CODERS MEANS EASIER COMMUN");int n=99,i=0;for(G.translate(n,0);i<66;){double p=M.PI/2,t=-6/1.1*i*p/90+p;v.setGlyphPosition(i,new Point(0,0));var a=java.awt.geom.AffineTransform.getTranslateInstance((int)(n*M.cos(t)+n),(int)(-n*M.sin(t)+n));a.rotate(t=4*p*i/66);a.translate(n*M.cos(p-t),n*M.sin(p-t));a.translate(-4,0);G.fill(a.createTransformedShape(v.getGlyphOutline(i++)));}}});show();}}

```

Result:

**Explanation:**

```
import java.awt.*;              // Required import for almost everything
v->                             // Method with empty unused parameter and Frame return-type
  new Frame(){                  //  Create the Frame
   {                            //   In an inner code-block:
     add(new Panel(){           //    Add a Panel we can draw on:
       RenderingHints R;        //    RenderingHints null-instance on class-level, for static calls
       Math M;                  //    Math null-instance on class-level, for static calls
       public void paint(Graphics g){
                                //     Overwrite its paint method:
         Graphics2D G=(Graphics2D)g;
                                //      Cast Graphics to Graphics2D, to open up more methods
         G.setRenderingHint(R.KEY_ANTIALIASING,R.VALUE_ANTIALIAS_ON);
                                //      Apply a focus rendering, otherwise the text would be too vague
         var v=new Font("Hei",0,//      Create a "Hei"-Font (could alternatively be Kai)
                                //      (note that a centered monospace font would look a lot better
                                //       for letters like the 'I's)
               8)               //      of font-size 8
               .createGlyphVector(G.getFontRenderContext(),
                                //      Create a GlyphVector using that Font
           "ICATION MEANS FASTER CODING MEANS FEWER CODERS MEANS EASIER COMMUN");
                                //      With this String as text (first character will be at the top)
         int n=99,              //      Integer we'll use a bunch of times, radius of the circle 99
             i=0;               //      Index integer, starting at 0
         for(G.translate(n,0);  //      Set the starting coordinate at 99,0, otherwise the circle
                                //      would be partially outside of the Frame boundaries
             i<66;){            //      Loop `i` in the range [0, 66):
                                //      (66 is the amount of characters of the text)
           double p=M.PI/2,     //       Temp-double `p`, set to PI/2
                  t=            //       Temp-double, which will start at:
                    -6/1.1      //        -360 / 66
                    *i          //        multiplied by the current `i`
                    *p/90       //        multiplied by PI/180
                    +p;         //        And add PI/2 to that
           v.setGlyphPosition(i,new Point(0,0));
                                //        Start the `i`'th Glyph position at 0,0
           var a=java.awt.geom.AffineTransform.getTranslateInstance(
                                //        Create an AffineTransform Translation instance, starting at:
                  (int)(        //         x-coordinate as casted integer:
                   n            //          99
                   *M.cos(t)    //          multiplied by the cosine of `t`
                   +n),         //          Add 99 to that
                  (int)(        //         y-coordinate as casted integer:
                   -n           //          -99
                   *M.sin(t)    //          multiplied by the sinus of `t`
                   +n));        //          Add 99 to that
           a.rotate(            //         Then rotate the current letter by,
                    t=          //         and set `t` to it at the same time:
                      4*p       //          2*PI
                      *i        //          multiplied by the current `i`
                      /66);     //          divided by the length 66
           a.translate(         //         Then move the letter to the correct position:
                                //          x-coordinate:
              n                 //           99
              *M.cos(p-t),      //           multiplied by the cosine of PI/2 - the new `t`
                                //          y-coordinate:
              n                 //           99
              *M.sin(p-t));     //           multiplied by the sinus of PI/2 - the new `t`
           a.translate(-4,0);   //          Slightly adjust based on the font-size
           G.fill(              //          Then add to the Panel through the Graphics2D:
             a.createTransformedShape(
                                //           The shape of:
               v.getGlyphOutline(i++)));}}});
                                //            The current `i`'th Glyph we just created
     show();}}                  //    And afterwards show the Frame

```

---
Tags: code-golf, kolmogorov-complexity, graphical-output, animation, typography
---


title: Fastest semiprime factorization

text:
# Question
Title: Fastest semiprime factorization
Write a program to factorize a semi-prime number in the shortest amount of time.

For testing purposes, use this: \\$38!+1\\$ (\\$523022617466601111760007224100074291200000001\\$)

It is equal to: \\$14029308060317546154181 × 37280713718589679646221\\$

# Answer
## Python (w/ pypy2 v7.3.1) ~0.9s

Using a Multiple Polynomial Quadratic Sieve. I took this to be a code challenge, so I opted not to use any external libraries (other than the standard `log` function, I suppose). When timing, the PyPy JIT should be used, as it results in timings 4-5 times faster than that of cPython.

**Update (2013-07-29):** Since originally posting, I've made several minor, but significant changes which increase the overall speed by a factor of about 2.5x.

**Update (2014-08-27):** As this post is still receiving attention, I've updated `my_math.py` correcting two errors, for anyone who may be using it:

* `isqrt` was faulty, sometimes producing incorrect output for values very close to a perfect square. This has been corrected, and the performance increased by using a much better seed.
* `is_prime` has been updated. My previous attempt to remove perfect square 2-sprps was half-hearted, at best. I've added a 3-sprp check - a technique used by Mathmatica - to ensure that the tested value is square-free.

**Update (2014-11-24):** If at the end of the calculation no non-trivial congruencies are found, the progam now sieves additional polynomials. This was previously marked in the code as `TODO`.

**Update (2020-08-29):** A few improvements resulting in minor performance gains, including:

* Calculating the average contributation of the skipped primes, rather than roughly approximating.
* Actually sieving powers of the prime factor base. The modular square root of each power can be computed efficiently through use of Hensel Lifting.
* Reusing the sieve array in subsequent iterations, rather than creating a new one.

---

**mpqs.py**

```
import math
import my_math
import time

# Multiple Polynomial Quadratic Sieve
# assumes n is composite
def mpqs(n, verbose=False):
  if verbose:
    time1 = time.time()

  root_n = my_math.isqrt(n)
  root_2n = my_math.isqrt(n+n)

  # formula chosen by experimentation
  # seems to be close to optimal for n < 10^50
  bound = int(7.5 * math.log(n, 10)**2)

  prime = []
  mod_root = []
  log_p = []
  num_prime = 0

  # size of the sieve
  x_max = bound * 5
  x_max_2 = x_max+x_max

  # maximum value on the sieved range
  m_val = (x_max * root_2n) >> 1

  # find a number of small primes for which n is a quadratic residue
  pg = my_math.primes()
  p = pg.next()
  while p < bound or num_prime < 3:

    # legendre (n|p) is only defined for odd p
    if p > 2:
      leg = my_math.legendre(n, p)
    else:
      leg = n & 1

    if leg == 1:
      prime += [p]
      log_p += [math.log(p, 10)]
      r = int(my_math.mod_sqrt(n, p))
      roots = [r]
      q = p
      while q < x_max:
        # find all square roots mod p^n via Hensel Lifting
        r = int((r + (n - r*r)*my_math.mod_inv(r+r, q))%q)
        #assert r*r%q == n%q
        roots += [r]
        q *= p
      mod_root += [roots]
      num_prime += 1
    elif leg == 0:
      if verbose:
        print 'trial division found factors:'
        print p, 'x', n/p
      return p

    p = pg.next()

  # fudging the threshold down a bit makes it easier to find powers of partial-partial
  # relationships, but it also makes the smoothness check slower. reducing by twice the log
  # of the largest prime in the factor base results in cofactors less than that value squared
  thresh = math.log(m_val, 10) - log_p[-1]*2

  # skip small primes. they contribute very little to the log sum
  # and add a lot of unnecessary entries to the table
  # instead, fudge the threshold down a bit, according to expected number of factors
  min_prime = int(thresh*2)
  sp_idx = my_math.binary_search(prime, min_prime)
  sieve_primes = prime[sp_idx:]

  fudge = sum([log_p[i]/(prime[i]-1) for i in xrange(sp_idx)])

  sums = [fudge]*x_max_2

  if verbose:
    print 'smoothness bound:', bound
    print 'sieve size:', x_max
    print 'log threshold:', thresh
    print 'skipping primes less than:', min_prime

  smooth = []
  used_prime = set()
  partial = {}
  num_smooth = 0
  prev_num_smooth = 0
  num_used_prime = 0
  num_partial = 0
  num_poly = 0
  root_A = my_math.isqrt(root_2n / x_max)

  if verbose:
    print 'sieving for smooths...'
  while True:
    # find an integer value A such that:
    # A is =~ sqrt(2*n) / x_max
    # A is a perfect square
    # sqrt(A) is prime, and n is a quadratic residue mod sqrt(A)
    while True:
      root_A = my_math.next_prime(root_A)
      leg = my_math.legendre(n, root_A)
      if leg == 1:
        break
      elif leg == 0:
        if verbose:
          print 'dumb luck found factors:'
          print root_A, 'x', n/root_A
        return root_A

    A = root_A * root_A

    # solve for an adequate B
    # B*B is a quadratic residue mod n, such that B*B-A*C = n
    # this is unsolvable if n is not a quadratic residue mod sqrt(A)
    b = my_math.mod_sqrt(n, root_A)
    B = (b + (n - b*b) * my_math.mod_inv(b + b, root_A))%A

    # B*B-A*C = n <=> C = (B*B-n)/A
    C = (B*B - n) / A

    num_poly += 1

    # sieve for prime factors
    i = sp_idx
    for p in sieve_primes:
      logp = log_p[i]

      e = 0
      q = p
      while q < x_max:
        inv_A = my_math.mod_inv(A, q)
        # modular root of the quadratic
        a = int(((mod_root[i][e] - B) * inv_A)%q)
        b = int(((q - mod_root[i][e] - B) * inv_A)%q)

        amx = a+x_max
        bmx = b+x_max

        apx = amx-q
        bpx = bmx-q

        k = q
        while k < x_max:
          sums[apx+k] += logp
          sums[bpx+k] += logp
          sums[amx-k] += logp
          sums[bmx-k] += logp
          k += q

        q *= p
        e += 1

      i += 1

    # check for smooths
    x = -x_max
    i = 0
    while i < x_max_2:
      v = sums[i]
      if v > thresh:
        vec = set()
        sqr = []
        # because B*B-n = A*C
        # (A*x+B)^2 - n = A*A*x*x+2*A*B*x + B*B - n
        #               = A*(A*x*x+2*B*x+C)
        # gives the congruency
        # (A*x+B)^2 = A*(A*x*x+2*B*x+C) (mod n)
        # because A is chosen to be square, it doesn't need to be sieved
        sieve_val = (A*x + B+B)*x + C

        if sieve_val < 0:
          vec = {-1}
          sieve_val = -sieve_val

        for p in prime:
          while sieve_val%p == 0:
            if p in vec:
              # keep track of perfect square factors
              # to avoid taking the sqrt of a gigantic number at the end
              sqr += [p]
            vec ^= {p}
            sieve_val = int(sieve_val / p)

        if sieve_val == 1:
          # smooth
          smooth += [(vec, (sqr, (A*x+B), root_A))]
          used_prime |= vec
        elif sieve_val in partial:
          # combine two partials to make a (xor) smooth
          # that is, every prime factor with an odd power is in our factor base
          pair_vec, pair_vals = partial[sieve_val]
          sqr += list(vec & pair_vec) + [sieve_val]
          vec ^= pair_vec
          smooth += [(vec, (sqr + pair_vals[0], (A*x+B)*pair_vals[1], root_A*pair_vals[2]))]
          used_prime |= vec
          num_partial += 1
        else:
          # save partial for later pairing
          partial[sieve_val] = (vec, (sqr, A*x+B, root_A))
      x += 1

      # reset the value for the next go
      sums[i] = fudge
      i += 1

    prev_num_smooth = num_smooth
    num_smooth = len(smooth)
    num_used_prime = len(used_prime)

    if verbose:
      print 100 * num_smooth / num_prime, 'percent complete\r',

    if num_smooth > num_used_prime and num_smooth > prev_num_smooth:
      if verbose:
        print '%d polynomials sieved (%d values)'%(num_poly, num_poly*x_max_2)
        print 'found %d smooths (%d from partials) in %.3f seconds'%(num_smooth, num_partial, time.time()-time1)
        print 'solving for non-trivial congruencies...'

      # set up bit fields for gaussian elimination
      masks = []
      mask = 1
      bit_fields = [0]*num_used_prime
      for vec, vals in smooth:
        masks += [mask]
        i = 0
        for p in used_prime:
          if p in vec: bit_fields[i] |= mask
          i += 1
        mask += mask

      # row echelon form
      col_offset = 0
      null_cols = []
      for col in xrange(num_smooth):
        pivot = col-col_offset == num_used_prime or bit_fields[col-col_offset] & masks[col] == 0
        for row in xrange(col+1-col_offset, num_used_prime):
          if bit_fields[row] & masks[col]:
            if pivot:
              bit_fields[col-col_offset], bit_fields[row] = bit_fields[row], bit_fields[col-col_offset]
              pivot = False
            else:
              bit_fields[row] ^= bit_fields[col-col_offset]
        if pivot:
          null_cols += [col]
          col_offset += 1

      # reduced row echelon form
      for row in xrange(num_used_prime):
        # lowest set bit
        mask = bit_fields[row] & -bit_fields[row]
        for up_row in xrange(row):
          if bit_fields[up_row] & mask:
            bit_fields[up_row] ^= bit_fields[row]

      # check for non-trivial congruencies
      for col in null_cols:
        all_vec, (lh, rh, rA) = smooth[col]
        lhs = lh   # sieved values (left hand side)
        rhs = [rh] # sieved values - n (right hand side)
        rAs = [rA] # root_As (cofactor of lhs)
        i = 0
        for field in bit_fields:
          if field & masks[col]:
            vec, (lh, rh, rA) = smooth[i]
            lhs += list(all_vec & vec) + lh
            all_vec ^= vec
            rhs += [rh]
            rAs += [rA]
          i += 1

        factor = my_math.gcd(my_math.list_prod(rAs)*my_math.list_prod(lhs) - my_math.list_prod(rhs), n)
        if 1 < factor < n:
          break
      else:
        if verbose:
          print 'none found.'
        continue
      break

  if verbose:
    print 'factors found:'
    print factor, 'x', n/factor
    print 'time elapsed: %.3f seconds'%(time.time()-time1)
  return factor

if __name__ == "__main__":
  import argparse
  parser = argparse.ArgumentParser(description='Uses a MPQS to factor a composite number')
  parser.add_argument('composite', metavar='number_to_factor', type=long, help='the composite number to factor')
  parser.add_argument('--verbose', dest='verbose', action='store_true', help="enable verbose output")
  args = parser.parse_args()

  if args.verbose:
    mpqs(args.composite, args.verbose)
  else:
    time1 = time.time()
    print mpqs(args.composite)
    print 'time elapsed: %.3f seconds'%(time.time()-time1)

```

**my\_math.py**

```
import math
import mpqs

# primes less than 212
small_primes = [
    2,  3,  5,  7, 11, 13, 17, 19, 23, 29, 31, 37,
   41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89,
   97,101,103,107,109,113,127,131,137,139,149,151,
  157,163,167,173,179,181,191,193,197,199,211]

# pre-calced sieve of eratosthenes for n = 2, 3, 5, 7
indices = [
    1, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47,
   53, 59, 61, 67, 71, 73, 79, 83, 89, 97,101,103,
  107,109,113,121,127,131,137,139,143,149,151,157,
  163,167,169,173,179,181,187,191,193,197,199,209]

# distances between sieve values
offsets = [
  10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6,
   6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4,
   2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6,
   4, 2, 4, 6, 2, 6, 4, 2, 4, 2,10, 2]

# tabulated, mod 105
dindices =[
  0,10, 2, 0, 4, 0, 0, 0, 8, 0, 0, 2, 0, 4, 0,
  0, 6, 2, 0, 4, 0, 0, 4, 6, 0, 0, 6, 0, 0, 2,
  0, 6, 2, 0, 4, 0, 0, 4, 6, 0, 0, 2, 0, 4, 2,
  0, 6, 6, 0, 0, 0, 0, 6, 6, 0, 0, 0, 0, 4, 2,
  0, 6, 2, 0, 4, 0, 0, 4, 6, 0, 0, 2, 0, 6, 2,
  0, 6, 0, 0, 4, 0, 0, 4, 6, 0, 0, 2, 0, 4, 8,
  0, 0, 2, 0,10, 0, 0, 4, 0, 0, 0, 2, 0, 4, 2]

max_int = 2147483647

# returns the index of x in a sorted list a
# or the index of the next larger item if x is not present
# i.e. the proper insertion point for x in a
def binary_search(a, x):
  s = 0
  e = len(a)
  m = e >> 1
  while m != e:
    if a[m] < x:
      s = m
      m = (s + e + 1) >> 1
    else:
      e = m
      m = (s + e) >> 1
  return m

# divide and conquer list product
def list_prod(a):
  while len(a) > 1:
    a = [m*n for m,n in zip(a[::2], a[1::2] + [1])]
  return a[0]

# greatest common divisor of a and b
def gcd(a, b):
  while b:
    a, b = b, a%b
  return a

# extended gcd
def ext_gcd(a, m):
  a = int(a%m)
  x, u = 0, 1
  while a:
    x, u = u, x - (m/a)*u
    m, a = a, m%a
  return (m, x, u)

# legendre symbol (a|m)
# note: returns m-1 if a is a non-residue, instead of -1
def legendre(a, m):
  return pow(a, (m-1) >> 1, m)

# modular inverse of a mod m
def mod_inv(a, m):
  return ext_gcd(a, m)[1]

# modular sqrt(n) mod p
# p must be prime
def mod_sqrt(n, p):
  a = n%p
  if p%4 == 3:
    return pow(a, (p+1) >> 2, p)
  elif p%8 == 5:
    v = pow(a << 1, (p-5) >> 3, p)
    i = ((a*v*v << 1) % p) - 1
    return (a*v*i)%p
  elif p%8 == 1:
    # Shanks' method
    q = p-1
    e = 0
    while q&1 == 0:
      e += 1
      q >>= 1

    n = 2
    while legendre(n, p) != p-1:
      n += 1

    w = pow(a, q, p)
    x = pow(a, (q+1) >> 1, p)
    y = pow(n, q, p)
    r = e
    while True:
      if w == 1:
        return x

      v = w
      k = 0
      while v != 1 and k+1 < r:
        v = (v*v)%p
        k += 1

      if k == 0:
        return x

      d = pow(y, 1 << (r-k-1), p)
      x = (x*d)%p
      y = (d*d)%p
      w = (w*y)%p
      r = k
  else: # p == 2
    return a

#integer sqrt of n
def isqrt(n):
  c = (n << 2)/3
  d = c.bit_length()

  a = d>>1
  if d&1:
    x = 1 << a
    y = (x + (n >> a)) >> 1
  else:
    x = (3 << a) >> 2
    y = (x + (c >> a)) >> 1

  if x != y:
    x = y
    y = (x + n/x) >> 1
    while y < x:
      x = y
      y = (x + n/x) >> 1
  return x

# integer cbrt of n
def icbrt(n):
  d = n.bit_length()

  if d%3 == 2:
    x = 3 << d/3-1
  else:
    x = 1 << d/3

  y = (2*x + n/(x*x))/3
  if x != y:
    x = y
    y = (2*x + n/(x*x))/3
    while y < x:
      x = y
      y = (2*x + n/(x*x))/3
  return x

# strong probable prime
def is_sprp(n, b=2):
  if n < 2: return False
  d = n-1
  s = 0
  while d&1 == 0:
    s += 1
    d >>= 1

  x = pow(b, d, n)
  if x == 1 or x == n-1:
    return True

  for r in xrange(1, s):
    x = (x * x)%n
    if x == 1:
      return False
    elif x == n-1:
      return True

  return False

# lucas probable prime
# assumes D = 1 (mod 4), (D|n) = -1
def is_lucas_prp(n, D):
  P = 1
  Q = (1-D) >> 2

  # n+1 = 2**r*s where s is odd
  s = n+1
  r = 0
  while s&1 == 0:
    r += 1
    s >>= 1

  # calculate the bit reversal of (odd) s
  # e.g. 19 (10011) <=> 25 (11001)
  t = 0
  while s:
    if s&1:
      t += 1
      s -= 1
    else:
      t <<= 1
      s >>= 1

  # use the same bit reversal process to calculate the sth Lucas number
  # keep track of q = Q**n as we go
  U = 0
  V = 2
  q = 1
  # mod_inv(2, n)
  inv_2 = (n+1) >> 1
  while t:
    if t&1:
      # U, V of n+1
      U, V = ((U + V) * inv_2)%n, ((D*U + V) * inv_2)%n
      q = (q * Q)%n
      t -= 1
    else:
      # U, V of n*2
      U, V = (U * V)%n, (V * V - 2 * q)%n
      q = (q * q)%n
      t >>= 1

  # double s until we have the 2**r*sth Lucas number
  while r:
    U, V = (U * V)%n, (V * V - 2 * q)%n
    q = (q * q)%n
    r -= 1

  # primality check
  # if n is prime, n divides the n+1st Lucas number, given the assumptions
  return U == 0

## Baillie-PSW ##
# this is technically a probabalistic test, but there are no known pseudoprimes
def is_bpsw(n):
  if not is_sprp(n, 2): return False

  # idea shamelessly stolen from Mathmatica's PrimeQ
  # if n is a 2-sprp and a 3-sprp, n is necessarily square-free
  if not is_sprp(n, 3): return False

  a = 5
  s = 2
  # if n is a perfect square, this will never terminate
  while legendre(a, n) != n-1:
    s = -s
    a = s-a
  return is_lucas_prp(n, a)

# an 'almost certain' primality check
def is_prime(n):
  if n < 212:
    m = binary_search(small_primes, n)
    return n == small_primes[m]

  for p in small_primes:
    if n%p == 0:
      return False

  # if n is a 32-bit integer, perform full trial division
  if n <= max_int:
    p = 211
    while p*p < n:
      for o in offsets:
        p += o
        if n%p == 0:
          return False
    return True

  return is_bpsw(n)

# next prime strictly larger than n
def next_prime(n):
  if n < 2:
    return 2

  n += 1
  if n < 212:
    m = binary_search(small_primes, n)
    return small_primes[m]

  # find our position in the sieve rotation via binary search
  x = int(n%210)
  m = binary_search(indices, x)
  i = int(n + (indices[m] - x))

  # adjust offsets
  offs = offsets[m:] + offsets[:m]
  while True:
    for o in offs:
      if is_prime(i):
        return i
      i += o

# an infinite prime number generator
def primes(start = 0):
  for n in small_primes[start:]: yield n
  pg = primes(6)
  p = pg.next()
  q = p*p
  sieve = {221: 13, 253: 11}
  n = 211
  while True:
    for o in offsets:
      n += o
      stp = sieve.pop(n, 0)
      if stp:
        nxt = n/stp
        nxt += dindices[nxt%105]
        while nxt*stp in sieve:
          nxt += dindices[nxt%105]
        sieve[nxt*stp] = stp
      elif n < q:
        yield n
      else:
        sieve[q + dindices[p%105]*p] = p
        p = pg.next()
        q = p*p

# true if n is a prime power > 0
def is_prime_power(n):
  if n > 1:
    for p in small_primes:
      if n%p == 0:
        n /= p
        while n%p == 0: n /= p
        return n == 1

    r = isqrt(n)
    if r*r == n:
      return is_prime_power(r)

    s = icbrt(n)
    if s*s*s == n:
      return is_prime_power(s)

    p = 211
    while p*p < r:
      for o in offsets:
        p += o
        if n%p == 0:
          n /= p
          while n%p == 0: n /= p
          return n == 1

    if n <= max_int:
      while p*p < n:
        for o in offsets:
          p += o
          if n%p == 0:
            return False
      return True

    return is_bpsw(n)
  return False

max_trial = 1e10
max_pollard = 1e22

def factor(n):
  if n < max_trial:
    return factor_trial(n)
  for p in small_primes:
    if n%p == 0:
      return [p] + factor(n/p)
  if is_prime(n):
    return [n]
  if n < max_pollard:
    p = pollard_rho(n)
  else:
    p = lenstra_ecf(n) or mpqs.mpqs(n)
  return factor(p) + factor(n/p)

def factor_trial(n):
  a = []
  for p in small_primes:
    while n%p == 0:
      a += [p]
      n /= p
  i = 211
  while i*i < n:
    for o in offsets:
      i += o
      while n%i == 0:
        a += [i]
        n /= i
  if n > 1:
    a += [n]
  return a

def pollard_rho(n):
  # Brent's variant
  y, r, q = 0, 1, 1
  c, m = 9, 40
  g = 1
  while g == 1:
    x = y
    for i in range(r):
      y = (y*y + c) % n
    k = 0
    while k < r and g == 1:
      ys = y
      for j in range(min(m, r-k)):
        y = (y*y + c) % n
        q = q*abs(x-y) % n
      g = gcd(q, n)
      k += m
    r *= 2
  if g == n:
    ys = (ys*ys + c) % n
    g = gcd(n, abs(x-ys))
    while g == 1:
      ys = (ys*ys + c) % n
      g = gcd(n, abs(x-ys))
  return g

def ec_add((x1, z1), (x2, z2), (x0, z0), n):
  t1, t2 = (x1-z1)*(x2+z2), (x1+z1)*(x2-z2)
  x, z = t1+t2, t1-t2
  return (z0*x*x % n, x0*z*z % n)

def ec_double((x, z), (a, b), n):
  t1 = x+z; t1 *= t1
  t2 = x-z; t2 *= t2
  t3 = t1 - t2
  t4 = 4*b*t2
  return (t1*t4 % n, t3*(t4 + a*t3) % n)

def ec_multiply(k, p, C, n):
  # Montgomery ladder algorithm
  p0 = p
  q, p = p, ec_double(p, C, n)
  b = k >> 1
  while b > (b & -b):
    b ^= b & -b
  while b:
    if k&b:
      q, p = ec_add(p, q, p0, n), ec_double(p, C, n)
    else:
      q, p = ec_double(q, C, n), ec_add(p, q, p0, n),
    b >>= 1
  return q

def lenstra_ecf(n, m = 5):
  # Montgomery curves w/ Suyama parameterization.
  # Based on pseudocode found in:
  # "Implementing the Elliptic Curve Method of Factoring in Reconfigurable Hardware"
  # Gaj, Kris et. al
  # http://www.hyperelliptic.org/tanja/SHARCS/talks06/Gaj.pdf
  # Phase 2 is not implemented.
  B1, B2 = 8, 13
  for i in range(m):
    pg = primes()
    p = pg.next()
    k = 1
    while p < B1:
      k *= p**int(math.log(B1, p))
      p = pg.next()
    for s in range(B1, B2):
      u, v = s*s-5, 4*s
      x = u*u*u
      z = v*v*v
      t = pow(v-u, 3, n)
      P = (x, z)
      C = (t*(3*u+v) % n, 4*x*v % n)
      Q = ec_multiply(k, P, C, n)
      g = gcd(Q[1], n)
      if 1 < g < n: return g
    B1, B2 = B2, B1 + B2

if __name__ == "__main__":
  import argparse
  import time
  parser = argparse.ArgumentParser(description='Uses a various methods to factor a composite number')
  parser.add_argument('composite', metavar='number_to_factor', type=long, help='the composite number to factor')
  args = parser.parse_args()

  time1 = time.time()
  print factor(args.composite)
  print 'time elapsed: %.3f seconds'%(time.time()-time1)

```

---

**Sample I/O:**

```
$ pypy mpqs.py --verbose 523022617466601111760007224100074291200000001
smoothness bound: 14998
sieve size: 74990
log threshold: 18.7325510316
skipping primes less than: 37
sieving for smooths...
230 polynomials sieved (34495400 values)
found 849 smooths (242 from partials) in 0.786 seconds
solving for non-trivial congruencies...
factors found:
37280713718589679646221 x 14029308060317546154181
time elapsed: 1.011 seconds

```

Note: not using the `--verbose` option will give slightly better timings:

```
$ pypy mpqs.py 523022617466601111760007224100074291200000001
37280713718589679646221
time elapsed: 0.902 seconds

```

---

### Basic Concepts

In general, a quadratic sieve is based on the following observation: any odd composite \\$n\\$ may be represented as:

\\$n=(x+d)(x-d)=x^2-d^2\Rightarrow d^2=x^2-n\\$

This is not very difficult to confirm. Since \\$n\\$ is odd, the distance between any two cofactors of \\$n\\$ must be even \\$2d\\$, where \\$x\\$ is the mid point between them. Moreover, the same relation holds for any multiple of \\$n\\$

\\$abn=(ax+ad)(bx-bd)=abx^2-abd^2\Rightarrow abd^2=abx^2-abn\\$

Note that if any such \\$x\\$ and \\$d\\$ can be found, it will immediately result in a (not necessarily prime) factor of \\$n\\$, since \\$x+d\\$ and \\$x-d\\$ both divide \\$n\\$ by definition. This relation can be further weakened - at the consequence of allowing potential trivial congruencies - to the following form:

\\$d^2\equiv x^2\mod n\\$

So in general, if we can find two perfect squares which are equivalent mod \\$n\\$, then it's fairly likely that we can directly produce a factor of \\$n\\$ a la \\$\gcd(x±d,n)\\$. Seems pretty simple, right?

Except it's not. If we intended to conduct an exhaustive search over all possible \\$x\\$, we would need to search the entire range from \\$\sqrt{n}\\$ to \\$\sqrt{2n}\\$, which is marginally smaller than full trial division, but also requires an expensive `is_square` operation each iteration to confirm the value of \\$d\\$. Unless it is known beforehand that \\$n\\$ has factors very near \\$\sqrt{n}\\$, trial division is likely to be faster.

Perhaps we can weaken this relation even more. Suppose we chose an \\$x\\$, such that for

\\$y\equiv x^2\mod n\\$

a full prime factorization of \\$y\\$ is readily known. If we had enough such relations, we should be able to *construct* an adequate \\$d\\$, if we choose a number of \\$y\\$ such that their product is a perfect square; that is, all prime factors are used an even number of times. In fact, if we have more such \\$y\\$ than the total number of unique prime factors they contain, a solution is guaranteed to exist; It becomes a system of linear equations. The question now becomes, how do we chose such \\$x\\$? That's where sieving comes into play.

### The Sieve

Consider the polynomial:

\\$y(x)=x^2-n\\$

Then for any prime \\$p\\$ and integer \\$k\\$, the following is true:

\\$y(x+kp)=(x+kp)^2-n\\y(x+kp)=x^2+2xkp+(kp)^2-n\\y(x+kp)=y(x)+2xkp+(kp)^2\equiv y(x)\mod p\\$

This means that after solving for the roots of the polynomial mod \\$p\\$ \- that is, you've found an \\$x\\$ such that \\$y(x)\equiv 0\mod p\\$, ergo \\$y\\$ is divisible by \\$p\\$ \- then you have found an infinite number of such \\$x\\$. In this way, you can sieve over a range of \\$x\\$, identifying small prime factors of \\$y\\$, hopefully finding some for which all prime factors are small. Such numbers known as \\$k-smooth\\$, where \\$k\\$ is the largest prime factor used.

There's a few problems with this approach, though. Not all values of \\$x\\$ are adequate, in fact, there's only very few of them, centered around \\$\sqrt{n}\\$. Smaller values will become largely negative (due to the \\$-n\\$ term), and larger values will become too large, such that it is unlikely that their prime factorization consists only of small primes. There will be a number of such \\$x\\$, but unless the composite you're factoring is very small, it's highly unlikely that you'll find enough smooths to result in a factorization. And so, for larger \\$n\\$, it becomes necessary to sieve over *multiple* polynomials of a given form.

### Multiple Polynomials

So we need more polynomials to sieve? How about this:

\\$y(x)=(Ax+B)^2-n\\$

That'll work. Note that \\$A\\$ and \\$B\\$ could literally be any integer value, and the math still holds. All we need to do is choose a few random values, solve for the root of the polynomial, and sieve the values close to zero. At this point we could just call it good enough: if you throw enough stones in random directions, you're bound to break a window sooner or later.

Except, there's a problem with that too. If the slope of the polynomial is large at the x-intercept, there'll only be a few suitable values to sieve per polynomial. It'll work, but you'll end up sieving a whole lot of polynomials before you get what you need. Can we do better?

We can do better. An observation, as a result of Montgomery is as follows: if \\$A\\$ and \\$B\\$ are chosen such that there exists some \\$C\\$ satisfying

\\$B^2-n=AC\\$

then the entire polynomial can be rewritten as

\\$y(x)=(Ax+B)^2-n=(Ax)^2+2ABx+B^2-n=A(Ax^2+2Bx+C)\\$

Furthermore, if \\$A\\$ is chosen to be a perfect square, the leading \\$A\\$ term can be neglected while sieving, resulting in much smaller values, and a much flatter curve. For such a solution to exist, \\$n\\$ must be a quadratic residue mod \\$\sqrt{A}\\$, which can be known immediately by computing the Legendre symbol: \\$(n|\sqrt{A})=1\\$. Note that in order to solve for \\$B\\$, a complete prime factorization of \\$\sqrt{A}\\$ needs to be known (in order to take the modular square root \\$\sqrt{n}\mod\sqrt{A}\\$), which is why \\$\sqrt{A}\\$ is typically chosen to be prime.

It can then be shown that if \\$A\approx\frac{\sqrt{2n}}{M}\\$, then for all values of \\$x\in\[-M,M\]\\$:

\\$|y(x)|\le\frac{M\sqrt{2n}}{2}\\$

And now, finally, we have all the components necessary to implement our sieve. Or do we?

### Powers of Primes as Factors

Our sieve, as described above, has one major flaw. It can identify which values of \\$x\\$ will result in a \\$y\\$ divisible by \\$p\\$, but it cannot identify whether or not this \\$y\\$ is divisible by a *power* of \\$p\\$. In order to determine that, we would need to perform trial division on the value to be sieved, until it is no longer divisible by \\$p\\$.

> **Edit:** This is incorrect. The roots of \\$n\\$ mod \\$p^k\\$ can be computed directly from the roots mod \\$p\\$ through the use of Hensel Lifting. This current implementation does precisely this.

We seemed to have reached an impasse: the whole point of the sieve was so that we *didn't* have to do that. Time to check the playbook.

\\$\ln(a\cdot b\cdot c\cdot d\cdot\ldots)=\ln(a)+\ln(b)+\ln(c)+\ln(d)+\ln(\ldots)\\$

That looks pretty useful. If the sum of the \\$\ln\\$ of all of the small prime factors of \\$y\\$ is close to the expected value of \\$\ln(y)\\$, then it's almost a given that \\$y\\$ has no other factors. In addition, if we adjust the expected value down a little bit, we can also identify values as smooth which have several powers of primes as factors. In this way, we can use the sieve as a 'pre-screening' process, and only factor those values which are likely to be smooth.

This has a few other advantages as well. Note that small primes contribute very little to the \\$\ln\\$ sum, but yet they require the most sieve time. Sieving the value 3 requires more time than 11, 13, 17, 19, and 23 *combined*. Instead, we can just skip the first few primes, and adjust the threshold down accordingly, assuming a certain percentage of them would have passed.

Another result, is that a number of values will be allowed to 'slip through', which are mostly smooth, but contain a single large cofactor. We could just discard these values, but suppose we found another mostly smooth value, with exactly the same cofactor. We can then use these two values to construct a usable \\$y\\$; since their product will contain this large cofactor squared, it no longer needs to be considered.

### Putting it all together

The last thing we need to do is to use these values of \\$y\\$ construct an adequate \\$x\\$ and \\$d\\$. Suppose we only consider the non-square factors of \\$y\\$, that is, the prime factors of an odd power. Then, each \\$y\\$ can be expressed in the following manner:

\\$y\_0=p\_0^0\cdot p\_1^1\cdot p\_2^1\cdots p\_n^0\\y\_1=p\_0^1\cdot p\_1^0\cdot p\_2^1\cdots p\_n^1\\y\_2=p\_0^0\cdot p\_1^0\cdot p\_2^0\cdots p\_n^1\\y\_3=p\_0^1\cdot p\_1^1\cdot p\_2^0\cdots p\_n^0\\\vdots\\$

which can be expressed in the matrix form:

\\$M=\begin{bmatrix}0&1&1&\cdots&0\\1&0&1&\cdots&1\\0&0&0&\cdots&1\\1&1&0&\cdots&0\\\vdots\end{bmatrix}\\$

The problem then becomes to find a vector \\$v\\$ such that \\$vM=\vec{0}\mod 2\\$, where \\$\vec{0}\\$ is the null vector. That is, to solve for the left null space of \\$M\\$. This can be done in a number of ways, the simplest of which is to perform Gaussian Elimination on \\$M^T\\$, replacing the row addition operation with a row `xor`. This will result in a number of null space basis vectors, any combination of which will produce a valid solution.

The construction of \\$x\\$ is fairly straight-forward. It is simply the product of \\$Ax+B\\$ for each of the \\$y\\$ used. The construction of \\$d\\$ is slightly more complicated. If we were to take the product of all \\$y\\$, we will end up with a value with 10s of thousands, if not 100s of thousands of digits, for which we need to find the square root. This calcuation is impractically expensive. Instead, we can keep track of the even powers of primes during the sieving process, and then use `and` and `xor` operations on the vectors of non-square factors to reconstruct the square root.

~~I seem to have reached the 30000 character limit. Ahh well, I suppose that's good enough.~~ Saved a bunch of bytes by switching to \\$MathJax\\$.
> 69 votes

# Answer
# Wolfram Language (Mathematica), 0.34 sec

```
FactorInteger

```

Try it online!
> 3 votes

# Answer
Well, your 38!+1 broke my php script, not sure why. In fact, any semi-prime over 16 digits long breaks my script.

However, using 8980935344490257 ( 86028157 * 104395301 ) my script managed a time of **25.963 seconds** on my home computer (2.61GHz AMD Phenom 9950). A lot faster than my work computer which was nearly 31 seconds @ 2.93GHz Core 2 Duo.

**php - 757 chars** incl. new lines

```
<?php
function getTime() {
    $t = explode( ' ', microtime() );
    $t = $t[1] + $t[0];
    return $t;
}
function isDecimal($val){ return is_numeric($val) && floor($val) != $val;}
$start = getTime();
$semi_prime = 8980935344490257;
$slice      = round(strlen($semi_prime)/2);
$max        = (pow(10, ($slice))-1);
$i          = 3;
echo "\nFactoring the semi-prime:\n$semi_prime\n\n";

while ($i < $max) {
    $sec_factor = ($semi_prime/$i);
    if (isDecimal($sec_factor) != 1) {
        $mod_f = bcmod($i, 1);
        $mod_s = bcmod($sec_factor, 1);
        if ($mod_f == 0 && $mod_s == 0) {
            echo "First factor = $i\n";
            echo "Second factor = $sec_factor\n";
            $end=getTime();
            $xtime=round($end-$start,4).' seconds';
            echo "\n$xtime\n";
            exit();
        }
    }
    $i += 2;
}
?>

```

I'd be interested to see this same algorithm in c or some other compiled language.
> 2 votes

# Answer
# Pari/GP, 0.21s

```
factor(input())

```

Pari/GP seems to be slightly more performant in this task than Mathmatica.

Try it online!
> 2 votes

---
Tags: fastest-code, primes
---


title: Convert LifeOnTheEdge to LifeOnTheSlope

text:
# Question
Title: Convert LifeOnTheEdge to LifeOnTheSlope
Your task here is to take a LifeOnTheEdge pattern and convert it to LifeOnTheSlope.

A LifeOnTheEdge pattern is composed of these four characters: `|_L` . A pattern corresponds to a certain arrangement of "on" edges in a square grid. The pattern is placed in the grid first with the characters in the cells, and each of the four letters specifies the state of the edges on the left and the bottom of that cell. `|` means the edge on the left is on, `_` means the bottom edge is on, `L` means both of them are on and  means neither of them are on.

For example the following LifeOnTheEdge:

```
|_L
 |

```

translates to:

```
. . . . .
|   |
. ._._. .
  |
. . . . .

```

Your task is however convert it to LifeOnTheSlope. LifeOnTheSlope is a LifeOnTheEdge equivalent but only uses three symbols: `/\` . You should rotate the pattern 45-degree clockwise, for example the above example translates to:

```
/

/\/
  \

```

This is , so shortest answer wins.

EDIT: You can have leading and/or trailing spaces/newlines as long as the output itself (somewhere on the screen) is correct.

## Test Cases

Input:
|_L
 |
Output:
 /

 /\/
   \

Input:
 |
_L
L_|
L_|
L_|
L_|
 |
 |
Output:
        /
     /\/
    /\ \
   /\ \/
  /\ \/
  \ \/
  /\/
 /

Input:
__ __ __ __   __ __   __
|  | | ||L_   |_ | || L_
L_ L_|_L|L_   L_|L_|L |

Output:
 /\
/  \
\
 \  /\
   /  \
   \  /
    \/  \
        /\
      \/ /
       \/ /\
         /\ \
         \ \
          \

               /\
              /  \
              \ \
               \/ /\
                 /  \
                 \  /
                  \/ /
                    /
                    \  /\
                      /\ \
                        \
```
</code></pre>

# Answer
> 5 votes
# JavaScript (ES6),  ~~135 120~~  119 bytes

*Saved 16 bytes thanks to @Shaggy!*

Expects either an array of strings or a matrix of characters. Returns a matrix of characters.

```
a=>[...a,...a[0]].map((_,y,b)=>b.map((_,x)=>'\\/ '[(c=(a[~(x-=y+a.length)>>1]||0)[x+y-~y>>1])=='_|'[x&=1]|c=='L'?x:2]))

```

Try it online!

# Answer
> 4 votes
# Charcoal, 27 bytes

```
ＷＳ«Ｐ↘⭆ι§ /№|Lκ↓Ｐ↘⭆ι§ \№_Lκ←

```

Try it online! Link is to verbose version of code. Takes input as a newline-terminated list. Explanation:

```
ＷＳ«

```

For each line in the list...

```
Ｐ↘⭆ι§ /№|Lκ

```

Diagonally print `/`s for `|`s and `L`s otherwise print spaces.

Move to the next line.

```
Ｐ↘⭆ι§ \№_Lκ

```

Diagonally print `\`s for `_`s and `L`s otherwise print spaces.

Move to the previous column.

# Answer
> 2 votes
# C (gcc), ~~195~~ 186 bytes

```
T,N,a,A,b,x,y,z;f(m,X,Y)char*m;{for(x=T=y=0;x<Y*3;x++,y++,T+=~Y)for(z=X-~Y,a=x,b=y;z--;T++,b--)A=a++/2-Y,N=m[A+b/2*X],putchar(z?(A|b)>=0&A<X&b<Y*2?T%2?N%19?32:92:(N-76)%48?32:47:32:10);}

```

Try it online!

* Saved 11 thanks to @ceilingcat

Example size X=3, Y=2

```
00 10 20  (x,y) indexes
01 11 21  of LOnH
      mapped to below

-2 0|-2 0|-1-1|-1-1| 0-2|
-2 0|-1 0|-1 0| 0-1| 0-1|
-1 1|-1 0| 0 0| 0 0| 1-1|
-1 1| 0 1| 0 0| 1 0| 1 0|
 0 2| 0 1| 1 1| 1 0| 2 0|
 0 2| 1 2| 1 1| 2 1| 2 0|
 1 3| 1 2| 2 2| 2 1| 3 1|
 1 3| 2 3| 2 2| 3 2| 3 1|

```

We start with x=-Y,y=0
We decrement x / increment y every 2 on first column and we do the inverse on each row.

Then we mask taking only valid indexes.

```
    |    |    |    |    |
    |    |    |    |    |
    |    | 0 0|    |    |
    | 0 1| 0 0| 1 0|    |
    | 0 1| 1 1| 1 0| 2 0|
    |    | 1 1| 2 1| 2 0|
    |    |    | 2 1|    |
    |    |    |    |    |

```

Now we have doubled (x,y) pairs one for `|`,`L` and the second for `_`,`L`

---
Tags: code-golf, ascii-art
---