[Solved] Prime Numbers Calculation

[UEZ]

Hi,

I tried to convert a PHP script to calculate prime numbers but something is very odd:

PHP code:

expandcollapse popup

<?php
$timeA = microtime(true); 
$limit = 1000; //1299710; 
$is_prime = array(); 
$is_prime[2] = true;
$is_prime[3] = true;
$is_prime[4] = false;
for ($i = 5; $i <= $limit; $i++) 
{ 
 $is_prime[$i] = false; 
} 
for ($x=1; $x <= sqrt($limit); $x++) 
{ 
 for ($y=1; $y <= sqrt($limit); $y++) 
 { 
 $n = 4 * $x*$x + $y*$y; 
 if ($n <= $limit && ($n%12 == 1 || $n%12 == 5)) 
 { 
 if ($is_prime[$n]) 
 { 
 $is_prime[$n] = false; 
 } 
 else 
 { 
 $is_prime[$n] = true; 
 } 
 } 
 $n = 3 * $x*$x + $y*$y; 
 if ($n <= $limit && $n%12 == 7) 
 { 
 if ($is_prime[$n]) 
 { 
 $is_prime[$n] = false; 
 } 
 else 
 { 
 $is_prime[$n] = true; 
 } 
 } 
 $n = 3 * $x*$x - $y*$y; 
 if ($x > $y && $n <= $limit && $n%12 == 11) 
 { 
 if ($is_prime[$n]) 
 {
 $is_prime[$n] = false; 
 } 
 else 
 { 
 $is_prime[$n] = true; 
 } 
 } 
 } 
} 
/* for ($n = 2; $n <= $limit; $n++) 
{ 
 if ($is_prime[$n]) 
 { 
 echo $n.";";
 } 
} */
?><br><?php
for ($n = 5; $n <= sqrt($limit); $n++) 
{ 
 if ($is_prime[$n]) 
 { 
 for ($k=1; $k * $n*$n <= $limit; $k++) 
 { //echo $k.";";
 //echo $k * $n*$n.";";
 $is_prime[$k * $n*$n] = false; 
 } 
 } 
} 
?><br><?php
$c=0;
for ($n = 2; $n <= $limit; $n++) 
{ 
 if ($is_prime[$n]) 
 { 
 //echo $n." ";
 $c++; 
 } 
} 
?><br><?php
?><br><?php
$timeB = microtime(true); 
 echo "Prime numbers found: $c";
 ?><br><?php
 echo 'Time: '; 
 echo $timeB - $timeA; 
?>

AutoIt code:

expandcollapse popup

#include <Array.au3>
Local $limit = 1000
Dim $is_prime[$limit]

$is_prime[2] = 1
$is_prime[3] = 1
$is_prime[4] = 0

For $i = 5 To $limit - 1
    $is_prime[$i] = 0
Next

For $x = 1 To Sqrt($limit)
    For $y = 1 To Sqrt($limit)
        $n = 4 * $x * $x + $y * $y
        If ($n <= $limit) And (Mod($n, 12) = 1 Or Mod($n, 12) = 5) Then
            If $is_prime[$n] Then
                $is_prime[$n] = 0
            Else
                $is_prime[$n] = 1
            EndIf
        EndIf
        $n = 3 * $x * $x + $y * $y
        If ($n <= $limit And Mod($n, 12) = 7) Then
            If $is_prime[$n] Then
                $is_prime[$n] = 0
            Else
                $is_prime[$n] = 1
            EndIf
        EndIf
        $n = 3 * $x * $x - $y * $y
        If ($x > $y) And ($n <= $limit) And (Mod($n, 12) = 11) Then
            If $is_prime[$n] Then
                $is_prime[$n] = 0
            Else
                $is_prime[$n] = 1
            EndIf
        EndIf
    Next
Next

;~ For $n = 2 To $limit - 1
;~  If $is_prime[$n] Then ConsoleWrite($n & ";")
;~ Next
;~ ConsoleWrite(@CRLF)

For $n = 5 To Sqrt($limit)
    If $is_prime[$n] Then
;~      ConsoleWrite($n & ";")
        For $k = 1 To ($k * $n * $n <= $limit)
;~          ConsoleWrite($k & ";")
;~          ConsoleWrite($k * $n * $n & ";")
            $is_prime[$k * $n * $n] = 0
        Next
    EndIf
Next

;~ _ArrayDisplay($is_prime)
$c = 0
For $n = 2 To $limit - 1
    If $is_prime[$n] Then $c += 1
Next
;~ ConsoleWrite(@CRLF)
ConsoleWrite("Prime numbers found: " & $c & @CRLF)

This example calculates the prime number till 1000! PHP returns 168 numbers which is correct and AutoIt returns 185 numbers!

I cannot see the problem here!

Can somebody help?

Thanks,

UEZ

Edited November 8, 2010 by UEZ

[JRowe]

From the results, these are obviously wrong:

41: 175

55: 245

61: 275

70: 325

87: 425

97: 475

113: 575

119: 605

139: 725

149: 775

160: 845

173: 925

There are probably a couple squares in there, too. Not sure what the relationship is yet, this should help.

[UEZ]

Yes, I know and probably the code

For $n = 5 To Sqrt($limit)
    If $is_prime[$n] Then
;~  ConsoleWrite($n & ";")
    For $k = 1 To ($k * $n * $n <= $limit)
;~  ConsoleWrite($k & ";")
;~  ConsoleWrite($k * $n * $n & ";")
    $is_prime[$k * $n * $n] = 0
    Next
    EndIf
Next

is the problem but when I compare PHP with AutoIt I cannot see the difference why these numbers are calculated!

UEZ

Edited March 3, 2010 by UEZ

[PsaltyDS]

Two things I see:

1. This is likely to produce a non-integer value: Sqrt($limit)

Maybe change to: Floor(Sqrt($limit))

2. This evaluates a Boolean compare then uses it to control a loop? For $k = 1 To ($k * $n * $n <= $limit)

So that evaluates to either: For $k = 1 To True, or For $k = 1 To False.

[UEZ]

The problem is when you implement it 1:1 then the array dimension will be exceeded:

For $n = 5 To Sqrt($limit)
 If $is_prime[$n] Then
;~ ConsoleWrite($n & ";")
 For $k = 1 To $k * $n * $n
;~ ConsoleWrite($k & ";")
;~ ConsoleWrite($k * $n * $n & ";")
 $is_prime[$k * $n * $n] = 0
 Next
 EndIf
Next

Array variable has incorrect number of subscripts or subscript dimension range exceeded.:

$is_prime[($k * $n * $n)] = 0

^ ERROR

although it is the same implemenation as in PHP which is working properly!

A little workaround (which should not needed because in PHP all calculated values are in the range of the array) will reduce the amount of prime number to 171 but it has still not prime numbers as prime marked!

For $n = 5 To Sqrt($limit)
    If $is_prime[$n] Then
;~      ConsoleWrite($n & ";")
        For $k = 1 To ($k * $n * $n)
;~          ConsoleWrite($k & ";")
;~          ConsoleWrite(($k * $n * $n) & ";")
            If ($k * $n * $n) <= $limit Then $is_prime[($k * $n * $n)] = 0
        Next
    EndIf
Next

Floor(Sqrt($limit)) produces the same result.

UEZ

Edited March 3, 2010 by UEZ

[PsaltyDS]

You should have declared the array as:

Dim $is_prime[$limit + 1]

I translate that other part this way:

#cs
    for ($n = 5; $n <= sqrt($limit); $n++)
    {
    
        if ($is_prime[$n])
        {
            for ($k=1; $k * $n*$n <= $limit; $k++)
            { //echo $k.";";
                //echo $k * $n*$n.";";
                $is_prime[$k * $n*$n] = false;
            }
        }
    }
#ce

For $n = 5 To Floor(Sqrt($limit)) Step +2 ; Only odd numbers can be prime, no point in testing even numbers
    If $is_prime[$n] Then
        $k = 1
        While $k * $n ^ 2 <= $limit
            $is_prime[$k * $n ^ 2] = False
            $k += 1
        WEnd
    EndIf
Next

[UEZ]

Thanks PsaltyDS for a solution but what I don't understand is why the original code is not working properly!

I made several tests and all values are equal with PHP one until that for next loop!

Why it is not possible to translate it from

for ($n = 5; $n <= sqrt($limit); $n++) {    if ($is_prime[$n])  {       for ($k=1; $k * $n*$n <= $limit; $k++)      {           $is_prime[$k * $n*$n] = false;      }   } }[/PHP]

to

For $n = 5 To Sqrt($limit)
    If $is_prime[$n] Then
        For $k = 1 To ($k * $n * $n)
            $is_prime[($k * $n * $n)] = 0
        Next
    EndIf
Next

easily? That's what I don't understand!

Thanks,

UEZ

Edited April 9, 2010 by UEZ

[PsaltyDS]

The problem with that is here: for ($k=1; $k * $n*$n <= $limit; $k++)

That says loop until ($k * $n*$n) is greater than $limit.

NOT until $k = ($k * $n*$n).

That's why I recast it as a while loop.

[UEZ]

Yes, that is the odd thing!

E.g.

In PHP:

for ($n = 5; $n <= sqrt($limit); $n++) 
{ 
    
    if ($is_prime[$n]) 
    {   
        for ($k=1; ($k * $n*$n) <= $limit; $k++) 
        { //echo $k.";";
            echo ($k * $n*$n).";";
            $is_prime[$k * $n*$n] = false; 
        } 
    } 
}

The echo result until 1000 is:

25;50;75;100;125;150;175;200;225;250;275;300;325;350;375;400;425;450;475;500;525;550;575;600;625;650;675;700;725;750;775;800;825;850;875;900;925;950;975;1000;49;98;147;196;245;294;343;392;441;490;539;588;637;686;735;784;833;882;931;980;121;242;363;484;605;726;847;968;169;338;507;676;845;289;578;867;361;722;529;841;961;

As you can see the array will never be exceeded!

Why the array in AutoIt is exceeded although is same code?

UEZ

Edited March 3, 2010 by UEZ

[PsaltyDS]

I just told you, it's not the same code. This:

for ($k=1; $k * $n*$n <= $limit; $k++)

Is NOT the same thing as this:

For $k = 1 To ($k * $n * $n)

[UEZ]

Yes, indeed, now I understand it!

It is something like

For $k = 1 To ($k * $n * $n <= $limit) but in AutoIt it is true or false!

expandcollapse popup

#AutoIt3Wrapper_UseX64=n
#include <Memory.au3>

Local $until = 1299709
Local $start = TimerInit()
;~ Local $pz = Eratosthenes_Calc($until)
Local $pz = Eratosthenes_Calc_ASM($until)
;~ Local $pz = Atkin_Prim_Calc($until)
;~ Local $pz = Simple_Prime_Calc($until)
;~ Local $pz = Eratosthenes_Calc_Variant($until)
Local $stop = TimerDiff($start)
Local $rowl = Ceiling(($until ^ 0.5) ^ 0.5)
Local $aTemp = StringSplit(StringMid($pz, 1, StringLen($pz) - 1), " ")

#region GUI ini
Local $iMemo, $hGUI
Local $font_size = 9
Local $display = 0

If $display Then
    $w = StringLen($until) * $rowl * $font_size
    If $w < 400 Then
        $w = 400
    ElseIf $w > @DesktopWidth * 0.95 Then
        $w = @DesktopWidth * 0.95
    EndIf
    $h = @DesktopHeight * 0.9
Else
    $w = 400
    $h = 150
EndIf
#endregion GUI ini

#region GUI
$hGUI = GUICreate("Prime Numbers Calculator by UEZ 2010", $w, $h, 0, 0)
$iMemo = GUICtrlCreateEdit("", 2, 2, $w - 4, $h - 6, 0x00200000 + 0x00100000 + 0x0800)
GUICtrlSetLimit(-1, 0x7FFFFFFF)
GUICtrlSetFont($iMemo, $font_size, 400, 0, "Courier New")
GUICtrlSetBkColor($iMemo, 0xFFFFFF)

MemoWrite("Sieving prime numbers till " & $until & ": " & @CRLF)
MemoWrite("Found " & $aTemp[0] & " prime numbers!" & @CRLF)
MemoWrite("Benchmark: " & Round($stop / 1000, 4) & " seconds to sieve prime numbers." & @CRLF)
If $display Then
    MemoWrite(@CRLF & "Numbers:" & @CRLF)
    For $z = 1 To $aTemp[0]
        MemoWrite(StringFormat("%" & StringLen($until) + 1 & "s", $aTemp[$z] & " "))
        If Mod($z, $rowl) = 0 Then MemoWrite(@CRLF)
    Next
EndIf
MemoWrite(@CRLF)
GUISetState()
#endregion GUI

MouseClick("left", 100, 85, 1, 1)

While GUIGetMsg() <> -3
WEnd

Exit

Func MemoWrite($sMessage = "")
    GUICtrlSetData($iMemo, $sMessage, 1)
EndFunc ;==>MemoWrite

Func Eratosthenes_Calc($l)
    Local $p[$l + 1], $i, $j, $z, $n
    For $z = 0x3 To $l Step 0x2
        $p[$z] = 0x1
    Next
    For $i = 0x3 To Sqrt($l) Step 0x2
        If $p[$i] And Mod($p[$i], 0x2) Then
            For $j = $i * $i To $l Step $i
                $p[$j] = 0
            Next
            $n &= $i & " "
        EndIf
    Next
    If Not Mod($i, 0x2) Then $i += 0x1
    For $z = $i To $l Step 0x2
        If $p[$z] Then $n &= $z & " "
    Next
    Return "2 " & $n
EndFunc ;==>Eratosthenes_Calc

Func Eratosthenes_Calc_ASM($limit);coded by UEZ - my 1st ASM code (thanks to AndyG for helping me :-)!
    #region bytecode
    Local $bytecode = "0x558B6C2408C7450C020000008B5D006BDB04B80800000089C1C1C90241894C050883C00839D876EF89EF83C70883C70101DFBB030000008B4D04B8080000008B74050839CB7F3383FE00750883C00883C302EBEB53505189D90FAFDB3B5D007F1489D86BC00483E804C74405080000000001CBEBE759585BEBD28B75006BF604B80800000089C1C1C902418B5C050883FB00750983C00839F076EA5DC3E802000000EBF05052535189D8BB0A00000031C931D2F7F3665280C10109C075F366580C30AAE2F9C6072083C701595B5A58C3"
    #cs
        :00000000 55    push ebp
        :00000001 8B6C2408  mov ebp, dword ptr [esp+08]
        :00000005 C7450C02000000    mov [ebp+0C], 00000002
        :0000000C 8B5D00    mov ebx, dword ptr [ebp+00]
        :0000000F 6BDB04    imul ebx, 00000004
        :00000012 B808000000    mov eax, 00000008
        :00000017 89C1  mov ecx, eax
        :00000019 C1C902    ror ecx, 02
        :0000001C 41    inc ecx
        :0000001D 894C0508  mov dword ptr [ebp+eax+08], ecx
        :00000021 83C008    add eax, 00000008
        :00000024 39D8  cmp eax, ebx
        :00000026 76EF  jbe 00000017
        :00000028 89EF  mov edi, ebp
        :0000002A 83C708    add edi, 00000008
        :0000002D 83C701    add edi, 00000001
        :00000030 01DF  add edi, ebx
        :00000032 BB03000000    mov ebx, 00000003
        :00000037 8B4D04    mov ecx, dword ptr [ebp+04]
        :0000003A B808000000    mov eax, 00000008
        :0000003F 8B740508  mov esi, dword ptr [ebp+eax+08]
        :00000043 39CB  cmp ebx, ecx
        :00000045 7F33  jg 0000007A
        :00000047 83FE00    cmp esi, 00000000
        :0000004A 7508  jne 00000054
        :0000004C 83C008    add eax, 00000008
        :0000004F 83C302    add ebx, 00000002
        :00000052 EBEB  jmp 0000003F
        :00000054 53    push ebx
        :00000055 50    push eax
        :00000056 51    push ecx
        :00000057 89D9  mov ecx, ebx
        :00000059 0FAFDB    imul ebx, ebx
        :0000005C 3B5D00    cmp ebx, dword ptr [ebp+00]
        :0000005F 7F14  jg 00000075
        :00000061 89D8  mov eax, ebx
        :00000063 6BC004    imul eax, 00000004
        :00000066 83E804    sub eax, 00000004
        :00000069 C744050800000000  mov [ebp+eax+08], 00000000
        :00000071 01CB  add ebx, ecx
        :00000073 EBE7  jmp 0000005C
        :00000075 59    pop ecx
        :00000076 58    pop eax
        :00000077 5B    pop ebx
        :00000078 EBD2  jmp 0000004C
        :0000007A 8B7500    mov esi, dword ptr [ebp+00]
        :0000007D 6BF604    imul esi, 00000004
        :00000080 B808000000    mov eax, 00000008
        :00000085 89C1  mov ecx, eax
        :00000087 C1C902    ror ecx, 02
        :0000008A 41    inc ecx
        :0000008B 8B5C0508  mov ebx, dword ptr [ebp+eax+08]
        :0000008F 83FB00    cmp ebx, 00000000
        :00000092 7509  jne 0000009D
        :00000094 83C008    add eax, 00000008
        :00000097 39F0  cmp eax, esi
        :00000099 76EA  jbe 00000085
        :0000009B 5D    pop ebp
        :0000009C C3    ret

        :0000009D E802000000    call 000000A4
        :000000A2 EBF0  jmp 00000094

        :000000A4 50    push eax
        :000000A5 52    push edx
        :000000A6 53    push ebx
        :000000A7 51    push ecx
        :000000A8 89D8  mov eax, ebx
        :000000AA BB0A000000    mov ebx, 0000000A
        :000000AF 31C9  xor ecx, ecx
        :000000B1 31D2  xor edx, edx
        :000000B3 F7F3  div ebx
        :000000B5 6652  push dx
        :000000B7 80C101    add cl, 01
        :000000BA 09C0  or eax, eax
        :000000BC 75F3  jne 000000B1
        :000000BE 6658  pop ax
        :000000C0 0C30  or al, 30
        :000000C2 AA    stosb
        :000000C3 E2F9  loop 000000BE
        :000000C5 C60720    mov byte ptr [edi], 20
        :000000C8 83C701    add edi, 00000001
        :000000CB 59    pop ecx
        :000000CC 5B    pop ebx
        :000000CD 5A    pop edx
        :000000CE 58    pop eax
        :000000CF C3    ret
    #ce
    #endregion bytecode
    Local $limes = Floor(Sqrt($limit))
    Local $bin_len = Ceiling(StringLen($bytecode) * 0.5)
    Local $pRemoteCode = _MemVirtualAlloc(0, $bin_len, $MEM_COMMIT, $PAGE_EXECUTE_READWRITE)
    Local $tCodeBuffer = DllStructCreate("byte[" & $bin_len & "]", $pRemoteCode)
    DllStructSetData($tCodeBuffer, 1, $bytecode)
    Local $tStruct = DllStructCreate("uint limit;uint limes;uint sieve[" & $limit + 1 & "];char prime[" & $limit & "]")
    DllStructSetData($tStruct, "limit", $limit)
    DllStructSetData($tStruct, "limes", $limes)
    Local $ret = DllCall("user32.dll", "int", "CallWindowProcW", "ptr", $pRemoteCode, "ptr", DllStructGetPtr($tStruct), "int", 0, "int", 0, "int", 0)
    Local $p = "2 3 5" & DllStructGetData($tStruct, "prime")
    _MemVirtualFree($pRemoteCode, $bin_len, $MEM_DECOMMIT)
    $tStruct = ""
    $tCodeBuffer = ""
    Return $p
EndFunc ;==>Eratosthenes_Calc_ASM

Func Eratosthenes_Calc_Variant($limes)
    $limes += 1
    Local $prime[$limes], $i, $ii, $j, $n, $prime_numbers
    $prime[2] = 1
    $prime[3] = 1
    For $n = 5 To $limes - 1 Step 2
        $prime[$n] = 1
    Next

    For $i = 0x3 To $limes - 1
        If $prime[$i] Then
            $prime_numbers &= $i & " "
            $j = 0x2 * $i
            While $j < $limes
                $prime[$j] = 0x0
                $j += $i
            WEnd
        EndIf
    Next
    Return "2 " & $prime_numbers
EndFunc ;==>Eratosthenes_Calc_Variant

Func Atkin_Prim_Calc($limes) ;using Atkin's algorithm
    $limes += 1
    Local $prime[$limes]
    Local $i, $k, $n, $nn, $w, $x, $xx, $y, $yy, $prime_numbers
    $w = Sqrt($limes)
    For $x = 1 To $w
        $xx = $x * $x ;x^2
        For $y = 1 To $w

            $yy = $y * $y ;y^2
            $n = $xx + $xx + $xx + $xx + $yy ;4 * x^2 + y^2
            If ($n < $limes) And (Mod($n, 12) = 1 Or Mod($n, 12) = 5) Then
                If $prime[$n] Then
                    $prime[$n] = 0
                Else
                    $prime[$n] = 1
                EndIf
            EndIf

            $n = $xx + $xx + $xx + $yy ;3 * x^2 + y^2
            If $n < $limes And Mod($n, 12) = 7 Then
                If $prime[$n] Then
                    $prime[$n] = 0
                Else
                    $prime[$n] = 1
                EndIf
            EndIf

            $n = $xx + $xx + $xx - $yy ;4 * x^2 - y^2
            If $x > $y And $n < $limes And Mod($n, 12) = 11 Then
                If $prime[$n] Then
                    $prime[$n] = 0
                Else
                    $prime[$n] = 1
                EndIf
            EndIf

        Next
    Next

    For $n = 5 To $w Step 2 ; Only odd numbers can be prime, no point in testing even numbers
        If $prime[$n] Then
            $nn = $n * $n
            $k = 1
            While $k * $nn < $limes
                $prime[$k * $nn] = 0
                $k += 1
            WEnd
        EndIf
    Next

    For $n = 3 To $limes - 1 Step 2
        If $prime[$n] Then $prime_numbers &= $n & " "
    Next
    Return "2 " & $prime_numbers
EndFunc ;==>Atkin_Prim_Calc

Func Simple_Prime_Calc($limit) ;simple and slow
    Local $v, $x
    Local $primes
    Local $z = 1
    While $z < $limit - 2
        $z += 2
        $x = 3
        $v = 2
        While 1
            If Mod($z, $x) = 0 Then
                ExitLoop
            Else
                If $x <= Floor($z / $v) - 2 Then
                    $v = $x
                    $x += 2
                Else
                    $primes &= $z & " "
                    $z += 2
                    $x = 3
                    $v = 2
                EndIf
            EndIf
        WEnd
    WEnd
    Return "2 3 " & $primes
EndFunc ;==>Simple_Prime_Calc

BR,

UEZ

Edited April 11, 2010 by UEZ