Diskuze: Ověření největšího nalezeného prvočísla

Caster

Člen

Caster:29.10.2024 22:50

Pokouším se ověřit, zda je aktuálně největší nalezené prvočíslo skutečně prvočíslem pomocí vlastního programu, který byl generován ChatGPT a mnou upraven.

Zkusil jsem: Vygeneroval jsem program na oveření metodou Lucas-Lehmer Primality Testing (viz konec článku), který používá "jen" 3 operace s velkými čísly:

mocnina čísla (pro urychlení násobení se používá FFT a iFFT)
odečtení čísla 2
modulo

Číslo 2^136279841 - 1 jsem uložil na konec alokované paměti 3 GB RAM jako _m256i epi64 (můj notobook umí AVX2). Při ladění ve Visual Studio lze jeho začátek v okně memory zabrazit po zadání "ptr + First₂₅₆_offset", jeho konec po zadání "prt + (RAM_BYTES - 32 ) / 8". Alokovaná paměr RAM ja definována jako uint64_t * proto je nutné dělit 8. Číslo je tedy uloženo od pozice First₂₅₆_offset (nejnižší bity) až po konec paměti (nejvyšší bity) tj. standardní formát little-endian po 64bitech ve formátu _m256i (256 bitů).

Jakékoliv číslo 2ⁿ - 1 představuje binárně jen samé "1". Např. 2⁴ = 0b10000 po odečtení 1 pak 0b1111 tj, 2ⁿ je n binárních jedniček za sebou.

Chápu, že jde o extrémní výpočet na běžném notebooku. Zjistil jsem, že program "nekonečně" dlouho pracuje ve funkci bit_reversal_avx2.

Zkoušel jsem také použít i Precompiled FFTW 3.3.5 Windows DLLs, ale po zadání adresy knihoven v nastavení programu (C++ a Linker) mi to hlásí, že nelze nalézt fftw3.h Po zadání celé cesty chyba zmizí, ale pak mi to hlásí 2 nové chyby, že nelze nalézt funkce.

Vlastní program:

// AVX2_LL.cpp : Tento soubor obsahuje funkci main. Provádění programu se tam zahajuje a ukončuje.
//

#define _USE_MATH_DEFINES

#include <iostream>
#include <iostream>
#include <immintrin.h>
#include <complex>
#include <cmath>

#include <chrono>      // Pridal jsem
#include "windows.h"   // Pridal jsem pro BYTE* ptr

using namespace std;
using namespace std::chrono;


int64_t* ptr = NULL;

// Číslo má 41 024 320 číslic
constexpr size_t GB = 3;
constexpr size_t TOTAL_BITS_N = 136279841;       // Total bits in the big number
constexpr size_t RAM_BYTES = GB * 1ULL << 30;    // Total bytes of allocated RAM
constexpr int FFT_SIZE = 134217728; // Suggested FFT size (2^27) for handling 17 MB numbers
constexpr int TOTAL_256 = (TOTAL_BITS_N + 255) / 256;  // 532 344 Total 256i in the big number, výraz obsahuje int zaokrouhlení nahoru
constexpr size_t First_256_offset  = (RAM_BYTES - TOTAL_256 * 256) / 8; // Offset for the large number in POLE; / 8 viz uint64_t *


unsigned char* POLE;

// Load the large number from memory at the specified offset as an array of _m256i
void load_large_number(__m256i* buffer) {
    // Calculate the starting address of the large number in POLE
    unsigned char* number_start = POLE + First_256_offset;

    // Copy the large number into the buffer
    for (int i = 0; i < TOTAL_256; ++i) {
        buffer[i] = _mm256_loadu_si256(reinterpret_cast<const __m256i*>(number_start + i * sizeof(__m256i)));
    }
}

// Bit-reversal for in-place FFT using POLE
void bit_reversal_avx2(double* data, int n) {
    int j = 0;
    for (int i = 0; i < n; i++) {
        if (i < j) {
            std::swap(data[2 * i], data[2 * j]);         // Swap real parts
            std::swap(data[2 * i + 1], data[2 * j + 1]); // Swap imaginary parts
        }
        int bit = n >> 1;
        while (bit <= j) {
            j -= bit;
            bit >>= 1;
        }
        j += bit;
    }
}

// In-place Cooley-Tukey FFT using AVX2
void FFT_avx2(double* data, int n) {
    bit_reversal_avx2(data, n);

    for (int len = 2; len <= n; len <<= 1) {
        double angle = -2 * M_PI / len;
        __m256d wlen_real = _mm256_set1_pd(cos(angle));
        __m256d wlen_imag = _mm256_set1_pd(sin(angle));

        for (int i = 0; i < n; i += len) {
            __m256d w_real = _mm256_set1_pd(1.0);
            __m256d w_imag = _mm256_set1_pd(0.0);

            for (int j = 0; j < len / 2; j += 4) { // Process 4 elements at a time
                __m256d u_real = _mm256_loadu_pd(&data[2 * (i + j)]);
                __m256d u_imag = _mm256_loadu_pd(&data[2 * (i + j) + 4]);
                __m256d v_real = _mm256_loadu_pd(&data[2 * (i + j + len / 2)]);
                __m256d v_imag = _mm256_loadu_pd(&data[2 * (i + j + len / 2) + 4]);

                // Multiply complex numbers (v_real + i * v_imag) * (w_real + i * w_imag)
                __m256d real_part = _mm256_sub_pd(_mm256_mul_pd(v_real, w_real),
                    _mm256_mul_pd(v_imag, w_imag));
                __m256d imag_part = _mm256_add_pd(_mm256_mul_pd(v_real, w_imag),
                    _mm256_mul_pd(v_imag, w_real));

                // Update data[i + j]
                _mm256_storeu_pd(&data[2 * (i + j)], _mm256_add_pd(u_real, real_part));
                _mm256_storeu_pd(&data[2 * (i + j) + 4], _mm256_add_pd(u_imag, imag_part));
                _mm256_storeu_pd(&data[2 * (i + j + len / 2)], _mm256_sub_pd(u_real, real_part));
                _mm256_storeu_pd(&data[2 * (i + j + len / 2) + 4], _mm256_sub_pd(u_imag, imag_part));

                // Update w using complex multiplication
                __m256d temp_real = _mm256_sub_pd(_mm256_mul_pd(w_real, wlen_real),
                    _mm256_mul_pd(w_imag, wlen_imag));
                w_imag = _mm256_add_pd(_mm256_mul_pd(w_real, wlen_imag),
                    _mm256_mul_pd(w_imag, wlen_real));
                w_real = temp_real;
            }
        }
    }
}

// Inverse FFT (IFFT) with AVX2
void IFFT_avx2(double* data, int n) {
    for (int i = 0; i < n; ++i) {
        data[2 * i + 1] = -data[2 * i + 1]; // Conjugate imaginary part
    }
    FFT_avx2(data, n);

    __m256d inv_n = _mm256_set1_pd(1.0 / n);
    for (int i = 0; i < n / 4; ++i) {
        __m256d real = _mm256_loadu_pd(&data[2 * i]);
        __m256d imag = _mm256_loadu_pd(&data[2 * i + 4]);
        _mm256_storeu_pd(&data[2 * i], _mm256_mul_pd(real, inv_n));
        _mm256_storeu_pd(&data[2 * i + 4], _mm256_mul_pd(imag, inv_n));
    }
}

// Modular reduction function using POLE
void mod_reduce_avx2(double* data, size_t length, __m256i* number) {
    for (size_t i = 0; i < length / 4; ++i) {
        __m256d element = _mm256_loadu_pd(&data[2 * i]);

        // Modular reduction operation for each _m256d element against `number`
        for (int j = 0; j < TOTAL_256; ++j) {
            __m256i modulus_block = number[j];
            int64_t* mod_vals = reinterpret_cast<int64_t*>(&modulus_block);
            double* elem_vals = reinterpret_cast<double*>(&element);

            // Modular reduction on each 64-bit chunk in the 256-bit block
            for (int k = 0; k < 4; ++k) {
                while (elem_vals[k] >= static_cast<double>(mod_vals[k])) {
                    elem_vals[k] -= static_cast<double>(mod_vals[k]);
                }
            }
        }

        _mm256_storeu_pd(&data[2 * i], element);
    }
}

// Square function using FFT and modular reduction, using POLE for storage
void square_using_fft_avx2(double* data, int n, __m256i* number) {
    FFT_avx2(data, n);
    for (int i = 0; i < n; ++i) {
        double real = data[2 * i];
        double imag = data[2 * i + 1];
        data[2 * i] = real * real - imag * imag; // Real part of square
        data[2 * i + 1] = 2 * real * imag;       // Imaginary part of square
    }
    IFFT_avx2(data, n);
    mod_reduce_avx2(data, n, number); // Reduce modulo the large number
}

// Lucas-Lehmer test using FFT and modular arithmetic with AVX2
bool lucas_lehmer_primality_test() {
    int n = FFT_SIZE;
    double* S = reinterpret_cast<double*>(ptr); // bylo POLE

    S[0] = 4.0; // Real part of initial S
    S[1] = 0.0; // Imaginary part of initial S
    for (int i = 1; i < n; ++i) {
        S[2 * i] = 0.0;
        S[2 * i + 1] = 0.0;
    }

    __m256i* number = reinterpret_cast<__m256i*>(ptr + First_256_offset);   // bylo POLE

    for (int i = 0; i < 5; ++i) {
        square_using_fft_avx2(S, n, number);
        S[0] -= 2.0; // Subtract 2 from the real part of S[0]
        mod_reduce_avx2(S, n, number);
    }

    // Check if result is zero
    for (int i = 0; i < n; ++i) {
        if (S[2 * i] != 0.0 || S[2 * i + 1] != 0.0) {
            return false;
        }
    }
    return true;
}

int main() {
    //// Allocate 3 GB memory for POLE
    //POLE = new unsigned char[RAM_BYTES];
    //if (!POLE) {
    //    std::cerr << "Memory allocation failed." << std::endl;
    //    return 1;
    //}

      //* Alokujeme pamět pro celý soubor *//
    static const SIZE_T giga = 1024 * 1024 * 1024;
    static const SIZE_T size = GB * giga;
    ptr = static_cast<int64_t*>(VirtualAlloc(NULL, size, MEM_COMMIT, PAGE_READWRITE));
    POLE = reinterpret_cast<unsigned char*>(ptr);

    // Bacha + konstanta k pointru posunu pointer o 8x konstanta bytů
    memset((ptr + First_256_offset), 0x11, 0x81F77E0LL); // Initialize POLE with "1"; 136279841
    _mm256_maskstore_epi64(ptr + ((RAM_BYTES - 32) / 8), _mm256_setr_epi64x(-1, -1, -1, -1), _mm256_setr_epi64x(0x1111111111111111, 0x0000000000000001LL, 0x0, 0x0));

    // Load the large number into an array of _m256i
    __m256i* number = reinterpret_cast<__m256i*>(ptr + First_256_offset);
    //load_large_number(number);

    auto start = high_resolution_clock::now();

    if (lucas_lehmer_primality_test()) {
        std::cout << "The number is prime." << std::endl;
    }
    else {
        std::cout << "The number is not prime." << std::endl;
    }

    auto stop = high_resolution_clock::now();
    auto duration = duration_cast<microseconds>(stop - start);
    printf("Time  = %lli (us DEC)\n", duration.count());


    //delete[] POLE;
    return 0;
}

Chci docílit: Jde mi o to, zda by nešlo program nebo jeho jednotlivé funkce lépe naprogramovat, aby běžel podstatně rychleji.

Odpovědět

29.10.2024 22:50

Peter Mlich

Člen

Peter Mlich:30.10.2024 9:36

Cecko nepouzivam, ani AI, takze moc neporadim. Ale, mozna ta funkce nefunguje spravne.
Jeste by slo pouzit pole prvocisel nebo je vygenerovat do 50% velikosti toho cisla a zkouset to jimi podelit.

Ale, kdyz mas obri prvocislo, a chces pouzit dnes bezne metody overeni, tak musis pocitat, ze to nejaky cas trva. Ikdyz, nedavno byl teda clanek, ze objevili nejakou novou metodu.
V zasade je s metodami problem, ze nejsou 100% Lze treba 100% rici, pro nektera cisla, ze to prvocislo neni. Ale pro jina ten vzorec nefunguje. Takze, pro 100% jistotu stejne musis pouzit deleni prvocisly.

Jinak bych teda prvocisla resil na binarni urovni. Tam by meli platit trochu jine metody nez pro desitkovou soustavu.
Podle meho, prvocisla nemaji vyznam a nema smysl se jimi zabyvat. Pro sifrovani a-b-c stejne dobre jde pouzit jiny vzorecek. Neni treba hledat 3 unikatni cisla. Zvlast, kdyz tech cisel je pomalu proti ostatnim kombinacim.
Stejne dobre funguje xor.
a xor b = c
a xor c = b
Kde a-b nejsou vubec spesl unikatni cisla, prvocisla.

Nahoru Odpovědět

30.10.2024 9:36

Peter Mlich

Člen

Peter Mlich:30.10.2024 9:42

Xor se pouziva pro samoopravne kody, treba na cd, dvd.
Navic, muzes pouzit neomezeny pocet klicu.
a xor b xor c xor d xor e xor f = G
a xor b xor G xor d xor e xor f = c

Nahoru Odpovědět

30.10.2024 9:42

Caster

Člen

Caster:30.10.2024 16:23

Díky za komentář. Podle článku Implementing fast carryless multiplication a grafu na konci článku by mělo násobení mého čísla trvat cca 0,5 sekundy. Vyzkouším to.

Nahoru Odpovědět

30.10.2024 16:23

Caster

Člen

Caster:3.11.2024 23:37

Na doporučení expertů z Francie jsem si na Windows 10 notebooku pod WSL 2 (Ubuntu 20.04.6 LTS) nainstaloval jejich knihovny pro práci s velkými čísly a poznámky:

svn checkout https://subversion.renater.fr/…scm/svn/mmx/

cd mmx
./configure --enable-justinline
make -j8
source set-devel-paths
mmc mmx/justinline/bench/polynomial_f2_60_bench.mmx

a zkusil spustit testovací program:
caster@John-AS:/mnt/c/Windows/system32/mmx$ ./justinline/bench/polynomial_f2_60_bench
Bench mul_pol_f2_60 in size 16384 2.00000000000000000
Bench mul_pol_f2_60 in size 32768 4.00000000000000000
Bench mul_pol_f2_60 in size 65536 12.0000000000000000
Bench mul_pol_f2_60 in size 131072 26.0000000000000000
Bench mul_pol_f2_60 in size 262144 50.0000000000000000
Bench mul_pol_f2_60 in size 524288 148.000000000000000
Bench mul_pol_f2_60 in size 1048576 258.000000000000000
Bench mul_pol_f2_60 in size 2097152 571.000000000000000
Bench mul_pol_f2_60 in size 4194304 1145.00000000000000
Bench mul_pol_f2_60 in size 8388608 2808.00000000000000
Bench mul_pol_f2_60 in size 16777216 7247.00000000000000

Požádal jsem jen také, aby mi poradili jak udělat krátký Linux *.mmx program a otestovat zda je to nalezené číslo skutečně prvočíslo.

Pomocí ChatGPT jsem si nechal vygenerovat program C++ Karatsuba pro násobení velkých uint čísel a trochu jsem ho upravil. Násobení čísla (2^136279841)-1 na druhou mi trvá cca 3 minuty viz výstup programu Time = 231968096 (us DEC).

Nejsem si ale jistý, že to program počítá správně, hlavně pokud jde o rekurzivní volání funkce a práci s _m256i (integer), když má data uložena v POLE (RAM) jako uint64_t i když o to vygenerovaný program ví a měl by s tím počítat.

Alokoval jsem si 3 GB RAM. Vlastní číslo je uloženo od pozice POLE + First₂₅₆_offset, konec čísla pak POLE + First₂₅₆_offset + TOTAL₂₅₆ * 4 - 4. Číslo má všech 136279841 bitů nastaveno na "1". V okně paměti se na něj lze při ladění podívat po zadání POLE + First₂₅₆_offset (LSB začátek čísla) a POLE + First₂₅₆_offset + TOTAL₂₅₆ * 4 - 4. Číslo je uloženo ve tvaru _m256 uint64_t jako little-endian.

Výsledek je uložen od začátku POLE a má dvojnásobnou velikost. Blbý je, že druhou mocninu čísla, která má 41,024,320 cifer nikde neověřím na online kalkulačce ;-) .

// AVX2_Karatsuba_mult.cpp : Tento soubor obsahuje funkci main. Provádění programu se tam zahajuje a ukončuje.
//

#include <immintrin.h>
#include <iostream>
#include <cstdint>
#include <chrono>      // Pridal jsem
#include "windows.h"   // Pridal jsem pro BYTE* ptr

using namespace std;
using namespace std::chrono;

// Číslo má 41 024 320 číslic
constexpr size_t GB = 3;
constexpr size_t RAM_BYTES = GB * 1ULL << 30;    // Total bytes of allocated RAM
constexpr size_t TOTAL_BITS_N = 136279841;       // Total bits in the big number
constexpr int TOTAL_256 = (TOTAL_BITS_N + 255) / 256;  // 532 344 Total 256i in the big number, výraz obsahuje int zaokrouhlení nahoru
constexpr size_t First_256_offset = (RAM_BYTES - TOTAL_256 * 256) / 8; // Offset for the large number in POLE; / 8 viz uint64_t *

// Function to perform the Karatsuba multiplication on 256-bit chunks
void karatsuba256(__m256i* A, __m256i* C, int n) {
    if (n == 1) {
        // Use appropriate multiplication for AVX2 with unsigned integers
        __m256i a_lo = _mm256_and_si256(A[0], _mm256_set1_epi64x(0xFFFFFFFF));
        __m256i a_hi = _mm256_srli_epi64(A[0], 32);

        __m256i a_lo_a_lo = _mm256_mul_epu32(a_lo, a_lo);
        __m256i a_hi_a_lo = _mm256_mul_epu32(a_hi, a_lo);
        __m256i a_lo_a_hi = _mm256_mul_epu32(a_lo, a_hi);
        __m256i a_hi_a_hi = _mm256_mul_epu32(a_hi, a_hi);

        __m256i lo = _mm256_add_epi64(a_lo_a_lo, _mm256_slli_epi64(a_hi_a_lo, 32));
        __m256i hi = _mm256_add_epi64(_mm256_slli_epi64(a_lo_a_hi, 32), _mm256_slli_epi64(a_hi_a_hi, 64));

        C[0] = _mm256_add_epi64(lo, hi);
        return;
    }

    int m = n / 2;
    __m256i* A1 = A;
    __m256i* A2 = A + m;

    __m256i* P1 = C;         // P1 = A1 * A1
    __m256i* P2 = C + n;     // P2 = A2 * A2
    __m256i* P3 = C + 2 * n;   // P3 = (A1 + A2) * (A1 + A2) - P1 - P2

    karatsuba256(A1, P1, m);
    karatsuba256(A2, P2, m);

    // Calculate (A1 + A2)
    __m256i* A_sum = new __m256i[m];
    for (int i = 0; i < m; ++i) {
        A_sum[i] = _mm256_add_epi64(A1[i], A2[i]);
    }

    karatsuba256(A_sum, P3, m);

    for (int i = 0; i < n; ++i) {
        P3[i] = _mm256_sub_epi64(P3[i], P1[i]);
        P3[i] = _mm256_sub_epi64(P3[i], P2[i]);
    }

    // Combine results
    for (int i = 0; i < n; ++i) {
        C[m + i] = _mm256_add_epi64(C[m + i], P3[i]);
    }

    delete[] A_sum;
}

int main() {
    const int num_bits = 136279841;
    const int num_uint64 = (num_bits + 63) / 64; // Number of 64-bit integers required
    const int num_256_chunks = (num_uint64 + 3) / 4; // Number of 256-bit chunks required

    //uint64_t* POLE = new uint64_t[num_uint64];

    //* Alokujeme pamět pro celý soubor *//
    static const SIZE_T giga = 1024 * 1024 * 1024;
    static const SIZE_T size = GB * giga;
    uint64_t* POLE = static_cast<uint64_t*>(VirtualAlloc(NULL, size, MEM_COMMIT, PAGE_READWRITE));

    __m256i mask = _mm256_setr_epi64x(-1, -1, -1, -1);
    __m256i values = _mm256_setr_epi64x(UINT64_MAX, UINT64_MAX, UINT64_MAX, UINT64_MAX);

    uint64_t i = TOTAL_256;
    while (i-- > 0) {
        _mm256_maskstore_epi64((long long int*) & POLE[First_256_offset + i * 4], mask, values);
    }
    _mm256_maskstore_epi64((long long int*) & POLE[First_256_offset + TOTAL_256 * 4 - 4], _mm256_setr_epi64x(-1, -1, -1, -1), _mm256_setr_epi64x(0x1111111111111111, 0x0000000000000001LL, 0x0, 0x0));


    __m256i* A = reinterpret_cast<__m256i*>(& POLE[First_256_offset]);

    //// Initialize POLE with all bits set
    //for (int i = 0; i < num_uint64; ++i) {
    //    POLE[i] = ~0ULL; // All bits set
    //}

    // Allocate space for the result
    //__m256i* C = new __m256i[4 * num_256_chunks];
    __m256i* C = reinterpret_cast<__m256i*>(POLE);

    auto start = high_resolution_clock::now();

    karatsuba256(A, C, num_256_chunks);

    auto stop = high_resolution_clock::now();
    auto duration = duration_cast<microseconds>(stop - start);
    printf("Time  = %lli (us DEC)\n", duration.count());

    //// Output the result
    //for (int i = 0; i < 4 * num_256_chunks; ++i) {
    //    uint64_t* res = reinterpret_cast<uint64_t*>(&C[i]);
    //    std::cout << "C[" << i << "]: " << std::hex << res[0] << " " << res[1] << " " << res[2] << " " << res[3] << std::endl;
    //}

    /*delete[] C;
    delete[] POLE;*/
    return 0;
}

Nahoru Odpovědět

3.11.2024 23:37

Peter Mlich

Člen

Peter Mlich:4.11.2024 7:58

Jj, pro praci s velkymi cisly je treba extra knihovny ve vsech jazycich.
Jeste mne napada, proc to delas v C++, kdyz jsou specializovane programy jako math-lab, ktere maji optimalizovane knihovny pro vypocty a umi vypocet rozkladat na vsechna cpu/gpu jadra v pocitaci. Kolegyne na tom pocitala nejake fyzikalni integraly z urychlovace a vychazelo ji to na 3 dny Takze tech 0.5s je naprosto zanedbatelne

Nahoru Odpovědět

4.11.2024 7:58

Caster

Člen

Caster:18.11.2024 16:22

Vypadá to, že mi program na násobní dvou čísel (2^136279841)-1 funguje. Nejsem si ale jistý správností výsledku. Program by měl vypsat u prvních *5ti hodnot: 1, 0, 0, 0 a 0. Vypíše ale:

Time = 68148010 (us DEC)
Pole[0]: 0x1
Pole[1]: 0
Pole[2]: 0
Pole[3]: 0xffffffffffffffff
Pole[4]: 0xfffffffffffffffd

*Poznámka: ověřuji pomocí funkce BitAnd[BitShiftRight[(2^136279841 - 1)^2, 0], 2⁶⁴ - 1] na wolframalpha.com Pro POLE[x] se ve vzorci místo 0 napíše x * 64 tj. pro Pole[3] se místo 0 zadá 192.

V programu si alokuji 3 GB RAM. Číslo (2^136279841) - 1 je uloženo ke konci paměti od pozice First₂₅₆_offset, do okna při ladění zadat POLE + First₂₅₆_offset/ 8, konec čísla je pak pod adresou POLE + First₂₅₆_offset/ 8 + n (zadat do okna paměti). Výsledek je uložen od začátku alokovaného paměti tj. POLE (zadat do okna paměti). Prostor paměti nad výsledkem a pod číslem používá program pro pomocná pole výpočtů.

Vlastní funkční program:

// AVX2_openai-01-preview_Karatsuba_test.cpp

#include <iostream>

#include <immintrin.h>
#include <cstdint>
#include <cstring>

#include <chrono>      // Pridal jsem
#include "windows.h"   // Pridal jsem pro BYTE* ptr

using namespace std;
using namespace std::chrono;

constexpr size_t GB = 3;

// Function prototypes
void big_add(uint64_t* result, const uint64_t* a, const uint64_t* b, size_t n);
void big_sub(uint64_t* result, const uint64_t* a, const uint64_t* b, size_t n);
void karatsuba_mul(uint64_t* result, const uint64_t* x, const uint64_t* y, size_t n, uint64_t* temp);

// The main function performing Karatsuba multiplication
void karatsuba_multiplication(uint64_t* POLE, size_t First_256_offset, size_t n_words) {
    // The two numbers are the same and start at POLE + First_256_offset
    uint64_t* x = POLE + First_256_offset / sizeof(uint64_t);
    uint64_t* y = x; // Since we are multiplying the number by itself

    // The result will be stored in POLE[0...2*n_words - 1]
    uint64_t* result = POLE; // Use the beginning of POLE for the result

    // Allocate temporary space from POLE[2*n_words...First_256_offset - 1]
    uint64_t* temp = POLE + 2ULL * n_words;
    //size_t temp_size = First_256_offset - 2ULL * n_words;
    size_t temp_size = (size_t)(GB * 0x40000000ULL) - 3ULL * n_words;

    // Ensure we have enough temporary space
    if (temp_size < 6 * n_words) {
        // Not enough temporary space
        // Handle error appropriately (e.g., return, print error)
        return;
    }

    karatsuba_mul(result, x, y, n_words, temp);
}

// Function to add two big numbers: result = a + b
void big_add(uint64_t* result, const uint64_t* a, const uint64_t* b, size_t n) {
    unsigned char carry = 0;
    for (size_t i = 0; i < n; ++i) {
        carry = _addcarry_u64(carry, a[i], b[i], &result[i]);
    }
}


// Function to subtract two big numbers: result = a - b
void big_sub(uint64_t* result, const uint64_t* a, const uint64_t* b, size_t n) {
    unsigned char borrow = 0;
    for (size_t i = 0; i < n; ++i) {
        borrow = _subborrow_u64(borrow, a[i], b[i], &result[i]);
    }
}

// Karatsuba multiplication function
void karatsuba_mul(uint64_t* result, const uint64_t* x, const uint64_t* y, size_t n, uint64_t* temp) {
    // Base case: if the size is small, use traditional multiplication
    if (n <= 32) { // Threshold can be adjusted
        // Zero the result
        memset(result, 0, 2 * n * sizeof(uint64_t));

        // Perform _umul128 multiplication
        for (size_t i = 0; i < n; ++i) {
            uint64_t carry = 0;
            for (size_t j = 0; j < n; ++j) {
                uint64_t low, high;
                low = _umul128(x[i], y[j], &high);

                // Add the low part and the current result value
                uint64_t sum = result[i + j] + low + carry;

                // Update the result with the sum's lower part
                result[i + j] = sum;

                // Calculate the new carry
                carry = high + (sum < low) + ((sum == low) && carry);
            }
            result[i + n] = carry;
        }
        return;
    }

    size_t k = n / 2;

    // Allocate space for temporary variables
    uint64_t* x_low = const_cast<uint64_t*>(x);
    uint64_t* x_high = const_cast<uint64_t*>(x + k);
    uint64_t* y_low = const_cast<uint64_t*>(y);
    uint64_t* y_high = const_cast<uint64_t*>(y + k);

    // Z0 = x_low * y_low
    uint64_t* z0 = result;
    karatsuba_mul(z0, x_low, y_low, k, temp);

    // Z2 = x_high * y_high
    uint64_t* z2 = result + n;
    karatsuba_mul(z2, x_high, y_high, n - k, temp);

    // Allocate temp arrays for (x_low + x_high) and (y_low + y_high)
    uint64_t* x_sum = temp;
    uint64_t* y_sum = temp + k;
    temp += 2 * k; // Advance temp pointer

    // x_sum = x_low + x_high
    big_add(x_sum, x_low, x_high, k);
    if (n % 2) {
        x_sum[k] = x_high[k]; // Copy the extra word if n is odd
    }

    // y_sum = y_low + y_high
    big_add(y_sum, y_low, y_high, k);
    if (n % 2) {
        y_sum[k] = y_high[k];
    }

    // Z1 = (x_low + x_high) * (y_low + y_high)
    uint64_t* z1 = temp; // Use remaining temp space
    karatsuba_mul(z1, x_sum, y_sum, k + (n % 2), temp + (n - k) * 4);

    // Compute Z1 = Z1 - Z0 - Z2
    big_sub(z1, z1, z0, n);
    big_sub(z1, z1, z2, n);

    // Assemble the final result:
    // result = Z2 * 2^{2k} + Z1 * 2^{k} + Z0

    // Add Z1 * 2^{k} to result
    for (size_t i = 0; i < n; ++i) {
        uint64_t carry = 0;
        size_t idx = i + k;
        uint64_t sum = result[idx] + z1[i] + carry;
        carry = (sum < result[idx]) || (carry && sum == result[idx]);
        result[idx] = sum;
        // Propagate carry
        size_t j = idx + 1;
        while (carry && j < 2 * n) {
            sum = result[j] + carry;
            carry = (sum == 0);
            result[j] = sum;
            ++j;
        }
    }

    // Add Z2 * 2^{2k} to result
    for (size_t i = 0; i < n; ++i) {
        uint64_t carry = 0;
        size_t idx = i + 2 * k;
        if (idx >= 2 * n) break;
        uint64_t sum = result[idx] + z2[i] + carry;
        carry = (sum < result[idx]) || (carry && sum == result[idx]);
        result[idx] = sum;
        // Propagate carry
        size_t j = idx + 1;
        while (carry && j < 2 * n) {
            sum = result[j] + carry;
            carry = (sum == 0);
            result[j] = sum;
            ++j;
        }
    }
}


// Main function demonstrating usage
int main() {
    size_t num_bits = 136279841;
    size_t num_uint64 = (num_bits + 255) / 256 * 4;
    size_t n = num_uint64;

    // Correct calculation for First_256_offset
    size_t First_256_offset = (GB * 0x40000000ULL) - ((2ULL + 1ULL) * num_uint64 * sizeof(uint64_t));

    //* Alokujeme pamět pro celý soubor *//
    static const SIZE_T giga = 1024 * 1024 * 1024;
    static const SIZE_T size = GB * giga;
    uint64_t* POLE = static_cast<uint64_t*>(VirtualAlloc(NULL, size, MEM_COMMIT, PAGE_READWRITE));
    uint64_t* number = POLE + First_256_offset / sizeof(uint64_t);

     // Store the number using _mm256_maskstore_epi64 in a loop
    __m256i ones = _mm256_set1_epi64x(-1);
    size_t i = 0;
    for (; i < (num_uint64)-4; i += 4) {
        _mm256_store_si256((__m256i*) & number[i], ones);
    }

    // Handle remaining bits
    _mm256_maskstore_epi64((long long int*) & number[i], _mm256_setr_epi64x(-1, -1, -1, -1), _mm256_setr_epi64x(0x1FFFFFFFF, 0x0, 0x0, 0x0));
    uint64_t ahoj = (uint64_t)number[i];

    auto start = high_resolution_clock::now();

    karatsuba_multiplication(POLE, First_256_offset, n);

    auto stop = high_resolution_clock::now();
    auto duration = duration_cast<microseconds>(stop - start);
    printf("Time  = %lli (us DEC)\n", duration.count());

    printf("Pole[0]: %#I64x\nPole[1]: %#I64x\nPole[2]: %#I64x\nPole[3]: %#I64x\nPole[4]: %#I64x\n", POLE[0], POLE[1], POLE[2], POLE[3], POLE[4]);

    return 0;
}

Nahoru Odpovědět

18.11.2024 16:22

Caster

Člen

Caster:21.11.2024 0:51

Pro ověření správného výsledku násobení čísla (2^136279841 - 1)^2 pomocí Karatsuba algoritmu jsem AI nechal udělat testovací program v RUSTu. Program jsek pak spustil na svém NB a výsledek souhlasí. Po vynásobení čísel je prvních 5 LSB uint64_t: 1, 0, 0, 0 a 0. Výsledek ověřen na wolframalpha(tečka)com Program v RUSTu ale na mém NB běžel více než hodinu. V C++ trvá násobení cca 1 minutu. Pracuji na finální verzi C++ programu, který pomocí Lucas-Lehmerova testu viz začátek ověří zda je výše uvedené číslo prvočíslo.

use std::cmp::min;

const BASE: u128 = 1 << 64; // 2^64 for splitting the numbers into chunks of u64
const CHUNK_BITS: usize = 64;

fn add(a: &mut Vec<u64>, b: &Vec<u64>, start: usize) {
    let mut carry: u128 = 0;
    let end = min(a.len(), b.len() + start);  // Ensure we're within bounds
    for i in start..end {
        let sum = a[i] as u128 + b[i - start] as u128 + carry;
        a[i] = (sum % BASE) as u64;
        carry = sum / BASE;
    }
    let mut i = end;
    while carry > 0 {
        if i >= a.len() {
            a.push(0); // Extend the vector if necessary
        }
        let sum = a[i] as u128 + carry;
        a[i] = (sum % BASE) as u64;
        carry = sum / BASE;
        i += 1;
    }
}

fn sub(a: &mut Vec<u64>, b: &Vec<u64>) {
    let mut borrow: i128 = 0;
    for i in 0..b.len() {
        let diff = a[i] as i128 - b[i] as i128 - borrow;
        if diff < 0 {
            a[i] = (diff + BASE as i128) as u64;
            borrow = 1;
        } else {
            a[i] = diff as u64;
            borrow = 0;
        }
    }
    for i in b.len()..a.len() {
        if borrow == 0 {
            break;
        }
        let diff = a[i] as i128 - borrow;
        if diff < 0 {
            a[i] = (diff + BASE as i128) as u64;
            borrow = 1;
        } else {
            a[i] = diff as u64;
            borrow = 0;
        }
    }
}

fn karatsuba(x: &Vec<u64>, y: &Vec<u64>) -> Vec<u64> {
    if x.len() <= 32 || y.len() <= 32 {
        // Base case: multiply directly for small sizes
        let mut result = vec![0; x.len() + y.len()];
        for i in 0..x.len() {
            let mut carry: u128 = 0;
            for j in 0..y.len() {
                let prod = x[i] as u128 * y[j] as u128 + result[i + j] as u128 + carry;
                result[i + j] = (prod % BASE) as u64;
                carry = prod / BASE;
            }
            result[i + y.len()] = carry as u64;
        }
        return result;
    }

    let m = min(x.len(), y.len()) / 2;
    let (x1, x0) = x.split_at(m);
    let (y1, y0) = y.split_at(m);

    let x1 = x1.to_vec();
    let x0 = x0.to_vec();
    let y1 = y1.to_vec();
    let y0 = y0.to_vec();

    let z2 = karatsuba(&x1, &y1);
    let z0 = karatsuba(&x0, &y0);

    let mut x0_plus_x1 = x0.clone();
    add(&mut x0_plus_x1, &x1, 0);
    let mut y0_plus_y1 = y0.clone();
    add(&mut y0_plus_y1, &y1, 0);

    let mut z1 = karatsuba(&x0_plus_x1, &y0_plus_y1);
    sub(&mut z1, &z2);
    sub(&mut z1, &z0);

    let mut result = vec![0; x.len() + y.len()];
    add(&mut result, &z0, 0);
    add(&mut result, &z1, m);
    add(&mut result, &z2, 2 * m);

    result
}

fn main() {
    let total_bits = 136279841;
    let num_chunks = (total_bits + CHUNK_BITS - 1) / CHUNK_BITS;
    let mut number = vec![0u64; num_chunks];

    // Construct the number 2^136279841 - 1: all bits set to 1
    for i in 0..num_chunks {
        if i < num_chunks - 1 {
            number[i] = !0u64; // Set all bits to 1
        } else {
            number[i] = (1 << (total_bits % CHUNK_BITS)) - 1; // Set remaining bits to 1
        }
    }

    let squared = karatsuba(&number, &number);

    // Print the first five least significant 64-bit chunks in hex
    for i in 0..5 {
        if i < squared.len() {
            println!("{:016x}", squared[i]);
        } else {
            println!("0000000000000000");
        }
    }
}

Nahoru Odpovědět

21.11.2024 0:51

DarkCoder

Člen

DarkCoder:21.11.2024 9:36

Pokud chceš ještě vyšší výkon pro násobení dvou velmi velkých čísel než Karatsuba, pak se podívej na Schönhage-Strassen algoritmus. Má časovou složitost O(nlognloglogn), což je velmi efektivní pro čísla o velikosti desítek milionů cifer. Popřípadě Fürer's algoritmus, ten teoreticky může být ještě efektivnější, nicméně jeho implementace je velmi složitá.

Dále může být užitečné, pokud se podíváš na pojem Mersennova čísla a práci s nimi. Což jsou čísla blízká mocnině 2 a existují optimalizované algoritmy pro jejich násobení.

Nahoru Odpovědět

21.11.2024 9:36

"I ta nejlepší poučka postrádá na významu, není-li patřičně předána." - DarkCoder

Caster

Člen

Caster:21.11.2024 10:28

Díky za tipy. Včera jsem zkusil spustit již celý program v C++ na test, zda je dané číslo prvočíslo. Po více jak 9 hodinách to sice spočítal chybně viz text dole ale vím, že je v programu drobná chybka. Urychlím také operaci MOD, která je nyní zdlouhavě počítána prostým odečítáním čísla.

Napadlo mě, že bude rychlejší dělitele vynásobit (Karatsuba) kvalifikovaným násobkem čísla 2 (jen bitové posuny), odečíst od prvního operandu a teprve pak použít prosté odečítání. Např. 25² - 2 MOD 25, 25² potřebuje pro uložení výsledku 9 bitů, 25 pak 4. Rozdíl je 9 - 5 = 4. Od 25² - 2 = 623 odečtu 2⁴ * 25 tj. 623 - 400 = 223 a pak už jen odečítat 25 až dostanu výsledek 23.

Time  = 33451692331 (us DEC) t.j 9:17:31
The number is NOT prime

Nahoru Odpovědět

21.11.2024 10:28

Peter Mlich

Člen

Peter Mlich:21.11.2024 15:58

"by mělo násobení mého čísla trvat cca 0,5 sekundy" - jsi tam nekde psal. To nic neznamena, zalezi na tom, jaky mas pc a pres jaka jadra to spousti. Jak jsem psal, mathlab a jine programy nachystaji ulohu a pak ji spusti pres vsechna jadra, co najdou, cpu, gpu, podle nastaveni a umi snad i omezit cinnost win, zrusit cinnosti, kdy win rozpozna, ze s pocitacem nepracujes, neovladas mys, klavesnici. tez zalezi na hw toho pc. Vis jak, ja ma v praci treba nic moc pocitac a soukromy notebook je uplne nekde jinde. Treba, kdyz maji k dispozici serverovy pocitac a na nem to spusti, tak to muze dopadnout uplne jinak

Ale jste mne nahlodali, ted budu premyslet, jakym zpusobem bych pocital prvocisla, abych to nedelal veky Musi byt nejaka jednoducha finta, jak zjistit, zda je cislo nasobkem jineho mensiho cisla.
U serazovacich algoritmu je treba nekolik fint. Treba pro n_log_n se pouziva spojovani serazenych poli 1-1 do 2, 2-2 do 4... az nakonec spojis 50% + 50% do 100%. A pritom, na prvni urovni, kde mas 1-1 muzes vsechny resit paralelne a tim cas stahnes treba na 30%. Treba neco podobneho jde udelat i u prvocisel.

Jeste by mohlo byt zajimave, zkusit na prvocisla furrierovu transformaci. Zda by to nehodilo zajimave koeficienty pro cely interval a pro pulku intervalu. Krivku bud sin+cos neno zkusit 2 na n. Nebo zkusit AI, zda by dokazala odhadnout pristi prvocislo
A az ted koukam na to pdf a o tom se prave bavite

Editováno 21.11.2024 16:00

Nahoru Odpovědět

21.11.2024 15:58

Caster

Člen

Caster:21.11.2024 16:33

Někde jsem četl, že pro velká čísla je rychlejší Karatsuba.

Teď jsem ale přišel na to, že nechápu jak můžete ten Lucas-Lehmer test pro provočíslo (2^136279841) - 1 fungovat. V každém z 5ti kroků je to číslo tj. B u A MOD B vždy větší než A a výsledkem MOD je:

krok 4² - 2 MOD prvočíslo = 14
krok 14² - 2 MOD prvočíslo = 194
krok 194² - 2 MOD provočíslo = 37634
krok 37634² -2 MOD provočíslo = 1 416 317 954
krok 1416317954² - 2 MOD provočíslo = 2005956546822746114

Prvočíslo (2^136279841) - 1 má více jak 40 mil. desítkových číslic... Pro A MOD B kde A není nikdy větší než B, je výsledek 5ti kroků vždy jen A. Počítal jsem to zbytečně. Je má úvaha správná ?

Nahoru Odpovědět

21.11.2024 16:33

Caster

Člen

Caster:21.11.2024 21:08

Díval jsem se na popis Lucas–Lehmer primality test na Wikipedii a zjistil jsem, že to prvočíslo se netestuje v 5. krocích, ale p - 2, kde p je exponent u 2ky. V našem případě bych chtěl testovat prvočíslo 2^136279841 - 1 tj. p = 136279841 - 2 = 136 279 839 kroků.

I pouhý prázdný cyklus 136 279 839 v x64 asm bude chvíli trvat. Výpis Hello World! mi trvá 137 µs

V popisu se dále píše, že v testu jsou dvě časově náročné operace: mocnina čísla a MOD. Místo dělení u MOD lze použit snadné p-bit sčítání. Zůstane nám tak časově náročné pouze násobení.

I tak ověření daného prvočísla na obyčejném notebooku asi nemá šanci vzhledem k časové náročnosti. Asi jen zkusím zefektivnit Karatsuba násobení a změřit čas, kam se mohu dostat. První možností je omezit práci s bity je na nezbytně nutnou velikost (nyní vždy počítám s 2 129 376 uint64_t daty tj. 17 035 000 bytů což je 136 280 064 bitů) např. při výpočtu mocniny čísla 4, u Karatsuba algoritmu zvážit rekorisvní volání a současnou práci na více jádrech CPU najednou. Mohu také případně zkusit naprogramovat Karatsuba násobení na mé slabé grafické kartě NVIDIA GeForce GTX 960M (4GB RAM, 5 SMs, CUDA capability 5.0).

Nahoru Odpovědět

21.11.2024 21:08

Caster

Člen

Caster:22.11.2024 18:07

Přeci jen mi to nedalo a napadlo mě, že by šel možný výpočet čtverce velkého čísla urychlit ve dvou krocích. Nejdříve ho vynásobit mocninou 2ⁿ, cože je jen bitový posun a někdy ani to ne viz dále, kde n je počet bitů čísla "x" v binárním tvaru. V některých případech je ale nutné pro posun použít hodnotu jen n - 1.

Ve druhém kroku by se v "předpřipraveném" výsledku v cyklu přičítalo (případně udečítalo) číslo "x: tak aby bylo výsledkem x². Jak ale zjistit hodnotu "i" pro cyklus ?

Tabulka, která to více objasní

Bitový posun v C++

size_t num_bits = 136279841;
size_t num_uint64 = ((num_bits + 255) / 256) * 4; // Number of uint64_t elements
__m256i tmp;
for (size_t i = 0; i < (num_uint64); i += 4) {
    tmp = _mm256_load_si256((__m256i*) & src[i]);
    _mm256_store_si256((__m256i*) & dst[i], tmp);
}

Na mém notebooku, mi výše uvedený cyklus pro num_uint = 532 344 trvá jen cca 11 ms. Číslo x mám tím "předvynásobeno" 2ⁿ, kde n = 136 280 064

Nahoru Odpovědět

22.11.2024 18:07

Caster

Člen

Caster:24.11.2024 12:05

Čtverec čísla 2^136279841 - 1 lze vypočítat jen posunem 136279841 bitů doleva a odečtením 2^136279841 - 1 od výsledku. Už mi to funguje správně. Pro jiná čísla bude nutné po posunu bitů v cyklu přičíst nebo odečíst dané číslo. Vím, jak spočítat kolikrát (číslo - 2ⁿ), záporná hodnota znamená odečítat. Čtverec velkého čísla jsem tím zredukoval na bitový posun a sčítání, případně odečítání.

Zbývá poslední problém a to MOD, ale to i by mělo jít vyřešit bez násobení a dělení. Dám vědět za pár dnů.

P.S. Pro uložení čísla 2^136279841 - 1 v paměti stačí cca 17 MB. Desítkově přitom obsahuje přes 40 mil. cifer. Čtverec čísla pak zabere 2x tolik.

Nahoru Odpovědět

24.11.2024 12:05

Peter Mlich

Člen

Peter Mlich:28.11.2024 15:39

1*1 + 1 = 2
1*1 + 2 = 3
2*2 + 1 = 5
2*2 + 3 = 7
3*3 + 2 = 11
3*3 + 4 = 13
4*4 + 1 = 17
4*4 + 3 = 19
4*4 + 7 = 23
5*5 + 4 = 29
5*5 + 6 = 31
6*6 + 1 = 37 // 5*5 + 12
6*6 + 5 = 41
6*6 + 7 = 43 (ctverec + plus prvocislo nebo 2 * nasobek prvocisla, ale ne vzdy to prvocislo je)
6*6 + 11 = 47
7*7 + 4 = 53
7*7 + 8(1+7) = 57 (tady to neni dvojnasobel prvocisla)
7*7 + 10(3+7) = 59
7*7 + 12(5+7) = 61
8*8 + 3 (8/2 - 1) = 67 // strida se vzdycky 1* prvocislo nebo 2*prvocislo, pro cely ctverec
8*8 + 7 (8-1) = 71
8*8 + 9 (8+1) = 73 // za pluskem uz to neni dvojnasobek prvocisla!
8*8 + 15 (8*2-1) = 79 (4 bity pro 8, tak tak pocet 8*8 budou 4)
9*9 + 2 = 83
9*9 + 8 = 89
9*9 + 16 = 97
9*9 + 20 = 101 // 10*10 + 1
9*9 + 22 = 103 // 10*10 + 3

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,

Cili, odectenim ctverce si muzes urcite zkratit vypocet. Ale nedokazi rici pri velkych cislech, do jake miry to jeste pujde. Hadam, ze by to mohlo jit az pro 2/3 cifer. Ale nedokazi rici, jakym presne algoritmem se ridi pocet prvocisel na danem stupni, nebo, kolik jich bude nutne testovat. U techto malych cisel to vypada, ze to pujde, ale netusim, zda to plati u velkych take. Jestli staci testovat jen do

n*n + A
n*n + B ...
n*n + X < (n+1)*(n+1) + 1

Editováno 28.11.2024 15:41

Nahoru Odpovědět

28.11.2024 15:39

Caster

Člen

Caster:29.11.2024 14:36

Díky za tipy. Pro informaci posílám tabulku s posledními nalezenými prvočísly viz č. 14. 15 jsem pomocí polynomu 3 stupně odhadl, kde by mělo asi být. Nyní pracuji na dokončení programu pro ověření té 14ky.

Nahoru Odpovědět

29.11.2024 14:36

Caster

Člen

Caster:29.11.2024 19:08

Program, který provádí jeden krok Lucas-Lehmerova testu na prvočíslo.
Znamená to vypočítat čtverec čísla tj. (2^136279841 - 1)^2, odečíst 2 a vypočítat zbytek po dělení tj. MOD(výsledek, 2^136279841 - 1).

Pro urychlení výpočtu čverce čísla ho nejdříve vynásobím mocninou čísla 2 tj. 2^136279841 případně o 1 méně (zjistím výpočtem), jde jen o bitový posun a k výsledku pak připočtu či odečtu vypočtený násobek (pomocí Karatsuba algoritmu) čísla 2^136279841 abych dostal správný čtverec.

Výpočet bez funkce MOD trvá po kompilaci Release necelé 2 minuty, s MOD to musím ještě vyzkoušet. Náročnost ověření daného čísla, zda jde o prvočíslo vystihuje skutečnost, že tento jeden krok je nutné opakovat 136279841x. Pokud je pak výsledek 0, dané číslo je prvočíslo.

Program generován s využitím OpenAI's o1 Preview FREE, ručně upraven a po částech vyzkoušen.

V režimu Debug lze dát v programu stopku např. na řádek 411 (MOD) a vypsat si zadané prvočíslo v okně paměti zadáním src (source) - v paměti jsou jen samé "1" tj. 0xFF na konci je pak sekvence 0x01FFFFFFFF, kde ta "1" je poslední, nejvyšší bit z 136279841 - 1 bitů a výsledek pomocí dst (destination), před operací MOD je 2x větší než src.

Program jinak nic nevypisuje, mimo času v µs, jak dlouho trvá výpočet.

// AVX2_square.cpp : Tento soubor obsahuje funkci main. Provádění programu se tam zahajuje a ukončuje.
//

#include <iostream>

#include <immintrin.h>
#include <cstdint>

#include <windows.h>   // For VirtualAlloc and VirtualFree
#include <chrono>      // Pridal jsem

using namespace std;
using namespace std::chrono;

// Constants and definitions
constexpr size_t GB = 3;
constexpr size_t giga = 1024 * 1024 * 1024;
size_t num_bits = 136279841;
//size_t num_bits = 3;
size_t num_uint64 = ((num_bits + 255) / 256) * 4; // Number of uint64_t elements

// Function prototypes
size_t num_bits_in_large_uint64(uint64_t* x, size_t len);
void big_integer_abs_difference(uint64_t* src, size_t len, size_t n, uint64_t* diff);
int big_integer_compare(uint64_t* A, uint64_t* B, size_t len);
void big_integer_subtract_in_place(uint64_t* A, uint64_t* B, size_t len);
void big_integer_add_in_place(uint64_t* result, uint64_t* A, size_t len);
void subtract_small_value(uint64_t* result, int64_t value, size_t num_uint64);
void karatsuba_multiply(uint64_t* x, uint64_t* y, uint64_t* result, size_t len, uint64_t* temp);
void bitwise_square(uint64_t* src, uint64_t* dst, uint64_t* tmp, size_t len);


// Function to determine the number of bits in a large integer
size_t num_bits_in_large_uint64(uint64_t* x, size_t len) {
    size_t bits = len * 64;
    while (bits > 0 && ((x[(bits - 1) / 64] & (1ULL << ((bits - 1) % 64))) == 0)) {
        bits--;
    }
    return bits;
}

// Function to compute absolute difference between src and 2 ^ n: abs(2 ^ n - src)
void big_integer_abs_difference(uint64_t* src, size_t len, size_t n, uint64_t* diff) {
    // Initialize diff to zero
    for (size_t i = 0; i < len; ++i) {
        diff[i] = 0;
    }

    size_t index = n / 64;
    uint64_t mask = 1ULL << (n % 64);

    // Set bit n in temp (representing 2^n)
    diff[index] = mask;

    // Compare diff (2^n) with src to determine which is larger
    int cmp = 0; // -1 if src > 2^n, 1 if src < 2^n, 0 if equal
    for (size_t i = len; i-- > 0;) {
        if (src[i] > diff[i]) {
            cmp = -1;
            break;
        }
        else if (src[i] < diff[i]) {
            cmp = 1;
            break;
        }
    }

    if (cmp >= 0) {
        // diff = diff - src
        uint8_t borrow = 0;
        for (size_t i = 0; i < len; ++i) {
            uint64_t temp;
            borrow = _subborrow_u64(borrow, diff[i], src[i], &temp);
            diff[i] = temp;
        }
    }
    else {
        // diff = src - diff
        uint8_t borrow = 0;
        for (size_t i = 0; i < len; ++i) {
            uint64_t temp;
            borrow = _subborrow_u64(borrow, src[i], diff[i], &temp);
            diff[i] = temp;
        }
    }
}

// Function to compare two big integers
int big_integer_compare(uint64_t* A, uint64_t* B, size_t len) {
    for (size_t i = len; i-- > 0;) {
        if (A[i] > B[i]) return 1;
        if (A[i] < B[i]) return -1;
    }
    return 0;
}

// Function to subtract big integer B from A: A = A - B
void big_integer_subtract_in_place(uint64_t* A, uint64_t* B, size_t len) {
    uint8_t borrow = 0;
    for (size_t i = 0; i < len; ++i) {
        uint64_t diff;
        borrow = _subborrow_u64(borrow, A[i], B[i], &diff);
        A[i] = diff;
    }
}

// Function to add big integer A to result: result = result + A
void big_integer_add_in_place(uint64_t* result, uint64_t* A, size_t len) {
    uint8_t carry = 0;
    for (size_t i = 0; i < len; ++i) {
        uint64_t sum;
        carry = _addcarry_u64(carry, result[i], A[i], &sum);
        result[i] = sum;
    }
}

// Function to subtract a small integer value from a large integer
void subtract_small_value(uint64_t* result, int64_t value, size_t num_uint64) {
    uint8_t borrow = 0;
    borrow = _subborrow_u64(borrow, result[0], value, &result[0]);

    for (size_t i = 1; i < num_uint64 && borrow; ++i) {
        borrow = _subborrow_u64(borrow, result[i], 0, &result[i]);
    }
}

// Function to perform modulo reduction: dst = dst mod src
void big_integer_mod(uint64_t* dst, size_t dst_len, uint64_t* src, size_t src_len) {
    // Ensure dst_len >= src_len
    if (dst_len < src_len) {
        return;
    }

    // Get the number of bits in dst and src
    size_t dst_bits = num_bits_in_large_uint64(dst, dst_len);
    size_t src_bits = num_bits_in_large_uint64(src, src_len);

    if (dst_bits < src_bits) {
        // dst is already less than src
        return;
    }

    // Temporary storage for shifted src and multiplication results
    // Since we cannot allocate new memory, we must use preallocated memory
    uint64_t* temp = dst + dst_len; // Ensure this is within allocated memory
    size_t temp_len = dst_len;

    // Zero out temp
    memset(temp, 0, temp_len * sizeof(uint64_t));

    while (dst_bits >= src_bits) {
        size_t shift = dst_bits - src_bits;

        // Shift src left by 'shift' bits and store in temp
        memset(temp, 0, temp_len * sizeof(uint64_t));
        size_t word_shift = shift / 64;
        size_t bit_shift = shift % 64;

        for (size_t i = 0; i < src_len; ++i) {
            size_t index = i + word_shift;
            if (index < temp_len) {
                temp[index] |= src[i] << bit_shift;
                if (bit_shift != 0 && index + 1 < temp_len) {
                    temp[index + 1] |= src[i] >> (64 - bit_shift);
                }
            }
        }

        // Estimate quotient q by dividing the highest word(s)
        uint64_t dst_high = dst[dst_len - 1];
        uint64_t temp_high = temp[dst_len - 1];

        // Handle cases where the highest word might be zero
        size_t dst_high_index = dst_len - 1;
        while (dst_high_index > 0 && dst[dst_high_index] == 0) {
            dst_high_index--;
        }
        dst_high = dst[dst_high_index];

        size_t temp_high_index = dst_len - 1;
        while (temp_high_index > 0 && temp[temp_high_index] == 0) {
            temp_high_index--;
        }
        temp_high = temp[temp_high_index];

        uint64_t q = dst_high / (temp_high + 1); // Add 1 to avoid division by zero
        if (q == 0) q = 1;

        // Multiply temp by q and store in temp2
        uint64_t* temp2 = temp + temp_len; // Use another part of preallocated memory
        memset(temp2, 0, temp_len * sizeof(uint64_t));

        uint8_t carry = 0;
        //for (size_t i = 0; i < temp_len; ++i) {
        //    __uint128_t product = (__uint128_t)temp[i] * q + carry;
        //    temp2[i] = (uint64_t)product;
        //    carry = (uint8_t)(product >> 64);
        //}

        for(size_t i = 0; i < temp_len; ++i) {
            uint64_t low;
            uint64_t high;
            uint64_t a = temp[i];
            uint64_t b = q;

            // Perform 64-bit multiplication
            low = _umul128(a, b, &high);

            // Add carry to the low part
            uint64_t new_low = low + carry;
            if (new_low < low) { // Check for carry out of the low part
                high++;
            }

            temp2[i] = new_low;
            carry = static_cast<uint8_t>(high);
        }

        // Subtract temp2 from dst
        big_integer_subtract_in_place(dst, temp2, dst_len);

        // Update dst_bits
        dst_bits = num_bits_in_large_uint64(dst, dst_len);

        // If dst < src, we're done
        if (big_integer_compare(dst, src, src_len) < 0) {
            break;
        }
    }

    // At this point, dst contains the remainder
}


// Karatsuba multiplication for big integers
void karatsuba_multiply(uint64_t* x, uint64_t* y, uint64_t* result, size_t len, uint64_t* temp) {
    // Base case for recursion
    if (len <= 32) {
        // For small sizes, use schoolbook multiplication
        memset(result, 0, 2 * len * sizeof(uint64_t));
        for (size_t i = 0; i < len; ++i) {
            uint64_t carry = 0;
            for (size_t j = 0; j < len; ++j) {
                uint64_t low, high;
                low = _umul128(x[i], y[j], &high);

                // Add the low part and the current result value
                uint64_t sum = result[i + j] + low + carry;

                // Update the result with the sum's lower part
                result[i + j] = sum;

                // Calculate the new carry
                carry = high + (sum < low) + ((sum == low) && carry);
            }
            result[i + len] = carry;
        }
        return;
    }

    size_t half = len / 2;

    uint64_t* x_low = x;
    uint64_t* x_high = x + half;
    uint64_t* y_low = y;
    uint64_t* y_high = y + half;

    uint64_t* z0 = result;                  // Placeholder for z0
    uint64_t* z1 = result + len;            // Placeholder for z1 (needs 2 * half space)
    uint64_t* z2 = result + len * 2;        // Placeholder for z2

    uint64_t* temp1 = temp;                 // Temporary space for (x_low + x_high)
    uint64_t* temp2 = temp + half + 1;      // Temporary space for (y_low + y_high)
    // Ensure temp has enough space

    // z0 = x_low * y_low
    karatsuba_multiply(x_low, y_low, z0, half, temp + 2 * (half + 1));

    // z2 = x_high * y_high
    karatsuba_multiply(x_high, y_high, z2, len - half, temp + 2 * (half + 1));

    // temp1 = x_low + x_high
    memset(temp1, 0, (half + 1) * sizeof(uint64_t));
    memcpy(temp1, x_low, half * sizeof(uint64_t));
    big_integer_add_in_place(temp1, x_high, len - half);

    // temp2 = y_low + y_high
    memset(temp2, 0, (half + 1) * sizeof(uint64_t));
    memcpy(temp2, y_low, half * sizeof(uint64_t));
    big_integer_add_in_place(temp2, y_high, len - half);

    // z1 = (temp1) * (temp2)
    karatsuba_multiply(temp1, temp2, z1, half + 1, temp + 2 * (half + 1));

    // z1 = z1 - z0 - z2
    big_integer_subtract_in_place(z1, z0, len * 2);
    big_integer_subtract_in_place(z1, z2, len * 2);

    // Combine results
    // result += z0 + (z1 << (half * 64)) + (z2 << (2 * half * 64))
    // z0 is already in place
    // Add z1 shifted by half words
    for (size_t i = 0; i < len * 2 - half; ++i) {
        uint64_t sum;
        uint8_t carry = _addcarry_u64(0, result[i + half], z1[i], &sum);
        result[i + half] = sum;
        size_t k = i + half + 1;
        while (carry && k < len * 2) {
            carry = _addcarry_u64(carry, result[k], 0, &sum);
            result[k] = sum;
            ++k;
        }
    }
    // Add z2 shifted by len words
    for (size_t i = 0; i < len * 2 - 2 * half; ++i) {
        uint64_t sum;
        uint8_t carry = _addcarry_u64(0, result[i + 2 * half], z2[i], &sum);
        result[i + 2 * half] = sum;
        size_t k = i + 2 * half + 1;
        while (carry && k < len * 2) {
            carry = _addcarry_u64(carry, result[k], 0, &sum);
            result[k] = sum;
            ++k;
        }
    }
}

// Modified bitwise_square function with Karatsuba multiplication
void bitwise_square(uint64_t* src, uint64_t* dst, uint64_t* tmp, size_t len) {
    // Initialize destination array
    memset(dst, 0, len * 2 * sizeof(uint64_t)); // The result can be twice the size of the input

    size_t total_bits = num_bits_in_large_uint64(src, len);

    // Determine whether to adjust total_bits
    // Compute abs(src - 2^total_bits)
    uint64_t* diff1 = tmp;              // Reuse tmp for diff1
    uint64_t* diff2 = tmp + len;        // Use additional space in tmp for diff2

    big_integer_abs_difference(src, len, total_bits, diff1);
    big_integer_abs_difference(src, len, total_bits - 1, diff2);

    // Compare diff1 and diff2
    int cmp = big_integer_compare(diff1, diff2, len);
    if (cmp > 0) {
        // diff1 > diff2, so difference with 2^(total_bits - 1) is smaller
        total_bits--;
    }
    // If cmp <= 0, we do not change total_bits (either diff1 < diff2 or they are equal)

    // Step 1: Shift src left by total_bits bits
    size_t shift_words = total_bits / 64;
    size_t shift_bits = total_bits % 64;

    // Perform the shift
    memset(dst, 0, (len * 2) * sizeof(uint64_t));
    for (size_t i = 0; i < len; ++i) {
        size_t dst_index = i + shift_words;
        if (dst_index < len * 2) {
            dst[dst_index] |= src[i] << shift_bits;
        }
        if (shift_bits != 0 && dst_index + 1 < len * 2) {
            dst[dst_index + 1] |= src[i] >> (64 - shift_bits);
        }
    }

    // Subtract src from dst
    //big_integer_subtract_in_place(dst, src, len);

    // Compute i_value = src - 2^total_bits
    uint64_t* i_value = tmp + len * 2; // Use additional space in tmp
    big_integer_abs_difference(src, len, total_bits, i_value);

    // Depending on the comparison, decide whether to add or subtract the product
    cmp = big_integer_compare(src, i_value, len);
    bool is_negative = false;
    if (cmp >= 0) {
        // i_value = src - 2^total_bits (positive)
        is_negative = false;
    }
    else {
        // i_value = 2^total_bits - src (negative)
        is_negative = true;
    }

    // Now, multiply src by i_value using Karatsuba multiplication
    uint64_t* multiply_result = tmp + len * 3; // Ensure we have enough space
    size_t multiply_len = len; // Lengths should be the same for Karatsuba multiplication

    // The Karatsuba multiplication requires additional temporary space
    size_t karatsuba_temp_space = multiply_len * 4;
    uint64_t* karatsuba_temp = tmp + len * 4;

    // Perform multiplication
    karatsuba_multiply(src, i_value, multiply_result, multiply_len, karatsuba_temp);

    // Add or subtract the multiplication result from dst
    if (!is_negative) {
        // dst = dst + multiply_result
        big_integer_add_in_place(dst, multiply_result, len * 2);
    }
    else {
        // dst = dst - multiply_result
        big_integer_subtract_in_place(dst, multiply_result, len * 2);
    }

    // Subtract 2 from dst
    subtract_small_value(dst, 2, len * 2);

    // MOD
    big_integer_mod(dst, len * 2, src, len);
}

// Main function
int main() {
    // Allocate 3GB memory
    static const size_t giga = 1024 * 1024 * 1024;
    static const size_t size = GB * giga;
    uint64_t* POLE = static_cast<uint64_t*>(VirtualAlloc(NULL, size, MEM_RESERVE | MEM_COMMIT, PAGE_READWRITE));

    if (!POLE) {
        std::cerr << "Memory allocation failed" << std::endl;
        return -1;
    }

    // Initialize the large number in POLE (example initialization)
    uint64_t* x = POLE;
    uint64_t* result = POLE + num_uint64;
    uint64_t* tmp = POLE + 3 * num_uint64;

    //// Example initialization of x
    //for (size_t i = 0; i < num_uint64; ++i) {
    //    x[i] = i + 1; // Just an example; replace with actual initialization
    //}

    // Store the number 2^136279841 - 1 using _mm256_maskstore_epi64 in a loop
    __m256i ones = _mm256_set1_epi64x(-1);
    size_t i = 0;
    for (; i < (num_uint64)-4; i += 4) {
        _mm256_store_si256((__m256i*) & x[i], ones);
    }

    // Handle remaining bits
    _mm256_maskstore_epi64((long long int*) & x[i], _mm256_setr_epi64x(-1, -1, -1, -1), _mm256_setr_epi64x(0x01FFFFFFFF, 0, 0, 0));

    //// 4
    //size_t i = 0;
    //// Handle remaining bits
    //_mm256_maskstore_epi64((long long int*) & x[i], _mm256_setr_epi64x(-1, -1, -1, -1), _mm256_setr_epi64x(1416317954, 0, 0, 0));


    auto start = high_resolution_clock::now();

    // Perform squaring
    bitwise_square(x, result, tmp, num_uint64);

    auto stop = high_resolution_clock::now();
    auto duration = duration_cast<microseconds>(stop - start);
    printf("Time  = %lli (us DEC)\n", duration.count());

    // Free allocated memory
    VirtualFree(POLE, 0, MEM_RELEASE);

    return 0;
}

Nahoru Odpovědět

29.11.2024 19:08

Caster

Člen

Caster:7.12.2024 23:00

Postupně program dolaďuji. Čtverec čísla do délky 64 bitů počítám klasicky viz výstup, pro větší čísla to pak bude vynásobení mocninou čísla 2 (bitový posun) a přičtení či odečtení rozdílu x * násobek (pomocí Karatsuba algoritmu) k výsledku tak, abych dostal správný čtverec čísla x. Krok je test prvočísla x = MOD(x² - 2, 2^136279841 - 1)

Krok: 1, x = 4
Krok: 2, x = 14
Krok: 3, x = 194
Krok: 4, x = 37634
Krok: 5, x = 1416317954
Krok: 6, x = 2005956546822746114

if (total_bits <= 64) { // We don't need to use Karatsuba for multiplication
    uint64_t low, high;
    low = _umul128(src[0], src[0], &high);
    dst[0] = low;
    dst[1] = high;
    printf("Krok: %llu, x = %llu\n", n, src[0]);
}
else {

Nahoru Odpovědět

7.12.2024 23:00

Caster

Člen

Caster:9.12.2024 22:55

Pro zajímavost jsem vyzkoušel násobení velkých uint čísel pomocí instrukce mulx (instrukce adox, adcx nebyly použity) v x64 asm funkci, kterou volám z C++ programu. Adresa alokované RAM paměti je předána v registru rcx (viz Example of argument passing 1 - all integers).

Pro čtverec (násobení) čísla o velikosti 2 129 376 * uint64_t je ale třeba dvou vnořených cyklů tj. (2 129 376)^2 = 4,53 E12 což je časově nesmysl.

P.S. Nedělá mi problém programovat MCU/CPU 8/32/64bit i v assembleru. Na konci kódu je nepoužitá funkce MultiplyBy10, který vynásobí malé uint číslo předané v rcx z C++ * 10 a výsledek vráti do C++ programu.

// AVX2_mulx_test.cpp : Tento soubor obsahuje funkci main. Provádění programu se tam zahajuje a ukončuje.
//

#include <iostream>

#include <immintrin.h>
#include <cstdint>

#include <windows.h>   // For VirtualAlloc and VirtualFree
#include <chrono>      // Pridal jsem

extern "C" uint64_t MultiplyBy10(uint64_t);
extern "C" void Multiply(uint64_t*);

using namespace std;
using namespace std::chrono;

// Constants and definitions
constexpr size_t GB = 3;
constexpr size_t giga = 1024 * 1024 * 1024;
size_t num_bits = 136279841;
//size_t num_bits = 3;
size_t num_uint64 = ((num_bits + 255) / 256) * 4; // Number of uint64_t elements

int main()
{
    // Allocate 3GB memory
    static const size_t giga = 1024 * 1024 * 1024;
    static const size_t size = GB * giga;
    uint64_t* POLE = static_cast<uint64_t*>(VirtualAlloc(NULL, size, MEM_RESERVE | MEM_COMMIT, PAGE_READWRITE));

    if (!POLE) {
        std::cerr << "Memory allocation failed" << std::endl;
        return -1;
    }

    // Initialize the large number in POLE (example initialization)
    uint64_t* x = POLE;
    uint64_t* result = POLE + num_uint64;
    uint64_t* prime = POLE + 3 * num_uint64;
    uint64_t* tmp = POLE + 4 * num_uint64;

    // Store the number 2^136279841 - 1 using _mm256_maskstore_epi64 in a loop
    __m256i ones = _mm256_set1_epi64x(-1);
    size_t i = 0;
    for (; i < (num_uint64)-4; i += 4) {
        _mm256_store_si256((__m256i*) & prime[i], ones);
    }

    // Handle remaining bits
    _mm256_maskstore_epi64((long long int*) & prime[i], _mm256_setr_epi64x(-1, -1, -1, -1), _mm256_setr_epi64x(0x01FFFFFFFF, 0, 0, 0));

    auto start = high_resolution_clock::now();

    // volat mulx
    //uint64_t value = 25;
    //uint64_t res = MultiplyBy10(value);
    //printf("25 * 10 = %lli\n", res);

    Multiply(x);

    auto stop = high_resolution_clock::now();
    auto duration = duration_cast<microseconds>(stop - start);
    printf("Time  = %lli (us DEC)\n", duration.count());

    // Free allocated memory
    VirtualFree(POLE, 0, MEM_RELEASE);

    return 0;
}

;mulx.asm

.data
x_ptr       DQ  0;
d_ptr       DQ  0;
p_ptr       DQ  0;

.code

Multiply      PROC
    push rbx
    push r12
    push rsi
    push rdi

    mov x_ptr, rcx      ; rcx = uint64_t* POLE
    mov r8, 103EF00h    ; uint64_t* dst = POLE + 1 * num_uint64;
    add rcx, r8
    mov d_ptr, rcx
    mov r8, 207DE00h    ; uint64_t* prm = POLE + 3 * num_uint64;
    add rcx, r8
    mov p_ptr, rcx

    ; Initialize pointers
    mov rsi, [p_ptr]                        ; rsi points to num array
    mov rdi, [d_ptr]                        ; rdi points to result array

    ; Initialize result to 0
    xor r8, r8
    xor r9, r9
    xor r10, r10
    xor r11, r11

    ; Zero initialize the result array
    mov rcx, 4258752
    mov rdi, d_ptr
zero_init:
    mov qword ptr [rdi], 0
    add rdi, 8
    loop zero_init

    ; Reset pointers
    mov rsi, [p_ptr]                      ; rsi points to num array
    mov rdi, [d_ptr]                      ; rdi points to result array

    ; Outer loop to handle each segment of num
    xor rbx, rbx                          ; outer loop counter

outer_loop:
    ; Reset rsi to the base address of num at the start of each outer loop
    mov rsi, [p_ptr]                      ; rsi points to the start of num
    mov rdx, [rsi + rbx*8]                ; load num[rbx]

    ; Inner loop to handle multiplication with all segments
    mov rcx, 2129376                      ; number of 64-bit segments
    mov r12, rbx
    shl r12, 3                            ; lea rdi, [result + rbx*8]
    add r12, [d_ptr]
    lea rdi, [r12]                        ; rdi points to result[2*rbx]
    mov rsi, [p_ptr]                      ; rsi points to num

    xor r8, r8                            ; reset accumulator for carries

inner_loop:
    mov rax, [rsi]                        ; load num[inner loop counter]
    mulx r10, r11, rdx                    ; (r11:r10) = rdx * rax

    ; Store the low and high parts of the result
    add qword ptr [rdi], r10              ; result[rbx] += r10 (low part)

    ; Accumulate carry using r12 as a temporary register
    mov r12, r11                          ; Move the high part (r11) to r12
    add r12, [rdi + 8]                    ; Add to the next result element (high part)
    mov [rdi + 8], r12                    ; Store updated high part in result[rbx+1]

    ; Handle carry for the low part using r8 (carry accumulator)
    add r8, r11                           ; Add the carry to the accumulator
    adc [rdi], r8                         ; result[rbx] += carry (r8)

    ; Move to the next segment
    add rsi, 8                            ; move to the next segment in num
    add rdi, 16                           ; move to the next segment in result

    ; Decrement inner loop counter and repeat if not zero
    loop inner_loop

    ; Store the final carry
    mov [rdi], r8

    ; Increment outer loop counter
    inc rbx
    cmp rbx, 2129376
    jl outer_loop

    pop rdi
    pop rsi
    pop r12
    pop rbx
    ret
Multiply      ENDP


MultiplyBy10  PROC

    shl           RCX, 1
    mov           RAX, RCX
    shl           RCX, 2
    add           RAX, RCX

    ret

MultiplyBy10  ENDP

END

Nahoru Odpovědět

9.12.2024 22:55

Caster

Člen

Caster:14.12.2024 10:36

Pro zajímavost funkce, pomocí které zjištuji počet bitů v uloženém čísle. Vlastní algoritmus s částečnou pomocí AI (test na první nejvyšší bit v uint64_t), která původně testovala v cyklu 136 279 841 každý jednotlivý bit...
Nově je cyklus jen 532 344, kdy načítám _m256i z uloženého čísla zprava.

// Function to determine the number of bits in a large integer
size_t num_bits_in_large_uint64(uint64_t* x, size_t len) {
    size_t k = num_uint64 / 4 - 1;
    while (k-- > 0) {
        __m256i v = _mm256_loadu_si256(reinterpret_cast<const __m256i*>(&x[k*4]));
        if (_mm256_testz_si256(v, v) == 0) {    // returns 1 if all bits in the result are zero, 0 otherwise
            // Extract each 64-bit segment from the __m256i variable
            uint64_t segments[4];
            _mm256_storeu_si256(reinterpret_cast<__m256i*>(segments), v);
            // Iterate over each 64-bit segment to find the first set bit
            for (int i = 3; i >= 0; --i) {
                if (segments[i] != 0) {
                    unsigned long index;
                    if (_BitScanReverse64(&index, segments[i])) {
                        return i * 64 + static_cast<int>(index) + 1 + k * 256;
                    }

                }
            }
        }
    }
}

Nahoru Odpovědět

14.12.2024 10:36

Caster

Člen

Caster:15.12.2024 13:53

Tak jsem tu poslední funkci ještě zjednodušil, nemá cenu načítat _m256i a pak rozložit na 4 * uint64_t

// Function to determine the number of bits in a large integer
size_t num_bits_in_large_uint64(uint64_t* x, size_t len) {
        size_t k = len - 1;
        while (k-- > 0) {
                if (x[k] != 0) {    // returns 1 if all bits in the result are zero, 0 otherwise
                        // Iterate over each 64-bit segment to find the first set bit
                        if (x[k] != 0) {
                                unsigned long index;
                                if (_BitScanReverse64(&index, x[k])) {
                                        return static_cast<int>(index) + 1 + k * 64;
                                }

                        }
                }
        }
}

Program doladím a upravím, aby variabilně pracoval se skutečnou délkou zpracovávaného čísla v daném kroku testu a ne vždy s jeho maximální velikostí t.j. 532 344 * uint64_t. Tím by se měl znatelně urychlit výpočet. Pokud délka aktuálního čísla v daném kroku nepřekročí 136279841 bitů t.j. 2^136279841 - 1 nemusím počítat MOD(x, 2^136279841 - 1) je to x.

Nahoru Odpovědět

15.12.2024 13:53

Jiří Kocián

Člen

Jiří Kocián:2. března 2:00

Jestli stačí pouze zjištění zda číslo je prvočíslo a počet bitů, pak je nejlepší použít big number knihovnu miracl zde MIRACL

Nahoru Odpovědět

2. března 2:00

When in doubt, use brute force.

Caster

Člen

Caster:2. března 21:48

Díku za tip MIRACL. Knihovna je primárně určena pro šifrování s čísly až do velikosti 4 096 bitů. Já ale pracuji s 136 279 841 bity ;-) . Na to je nutné použít speciální algoritmy. Nechodí mi správně jen funkce modulo (Knut D algoritmus). Pořeším s AI ;-) .

Nahoru Odpovědět

2. března 21:48

Jiří Kocián

Člen

Jiří Kocián:4. března 2:50

Není pravda, že MIRACL pracuje s čísly pouze do velikosti 4096 bitů.
na webu jsem našel toto:
MIRACL (Multiprecision Integer and Rational Arithmetic C Library) používá variabilní délku čísel a není pevně omezeno na konkrétní maximální velikost čísla, ale praktické omezení závisí na:

Dostupné paměti – Čísla jsou uložena jako pole 32bitových nebo 64bitových slov, takže největší číslo je omezeno dostupnou pamětí.
Nastavení při kompilaci – V MIRACL je možné nastavit maximální velikost čísla pomocí makra MR_BIG_RESERVE, což určuje počet slov, která budou použita k uložení velkého čísla.
Datový typ big – MIRACL pracuje s typem big, což je struktura obsahující ukazatel na pole s číslem.
Pokud není velikost omezena ručně, lze pracovat s čísly řádově tisíců až milionů bitů, pokud to dovoluje paměť systému. Pro konkrétní aplikaci je dobré zkontrolovat nastavení MIRACL a její konfiguraci při kompilaci.

Nahoru Odpovědět

4. března 2:50

When in doubt, use brute force.

Caster

Člen

Caster:5. března 23:53

AI Grog:

The MIRACL (Multiprecision Integer and Rational Arithmetic C/C++ Library) library is designed to handle very large numbers, specifically for cryptographic applications, and its capacity for big number operations is quite extensive. The exact size of numbers it can work with depends on how the library is configured and the underlying hardware, but here’s a breakdown based on its design and documentation:

Flexible Size Configuration:
MIRACL uses a multiprecision arithmetic approach where large integers are represented as arrays of smaller "digits" (typically based on the word size of the processor, e.g., 32-bit or 64-bit). The size of these big numbers is determined when the library is initialized via the mirsys() function, which takes a parameter specifying the number of digits and the base for each digit.

For example, calling miracl *mip = mirsys(500, 10) initializes the system to handle numbers up to 500 decimal digits. The base can be adjusted (e.g., 2, 10, 16, or even a power of 2 up to the maximum supported by the system), affecting the total bit size.
Bit Length Capacity:
The library can theoretically handle numbers with bit lengths far exceeding typical cryptographic requirements (e.g., 256-bit, 512-bit, or 4096-bit keys). In practice, the maximum size is limited by the memory available and the word size of the underlying architecture.

For a 32-bit system, if each "digit" is a 32-bit word, and you specify 500 digits, the library can handle numbers up to approximately 500×32=16,000500×32=16,000 bits. On a 64-bit system, this doubles to 32,000 bits for the same number of digits.

Documentation mentions optimizations for numbers up to 8192 bits in benchmark examples, indicating that MIRACL is well-tested for sizes commonly used in cryptography (e.g., RSA with 2048-bit or 4096-bit keys).
Practical Limits:
The actual limit depends on the memory allocated for each big variable. MIRACL allocates fixed-length arrays for its big numbers, and the size must be specified at initialization. If intermediate results (e.g., during multiplication) exceed this size, special routines like mad (Multiply, Add, Divide) are used to handle double-length products, which are then reduced.

For very large numbers (e.g., thousands of digits), the library can still operate efficiently, but performance may degrade if the numbers exceed the pre-allocated size, requiring careful configuration by the programmer.
Special Features:
MIRACL supports a big data type for large integers and a flash type for rational numbers (floating-slash arithmetic), both of which can scale to large sizes. The big type is optimized for integer operations critical to cryptography, such as modular exponentiation.

The library includes assembly language optimizations for time-critical operations on popular architectures (e.g., Intel 80x86), which can handle large numbers efficiently up to at least 8192 bits, as shown in its benchmarks.
Theoretical Maximum:
There’s no strict upper bound hardcoded into the library beyond what memory constraints impose. By adjusting the number of digits in mirsys() and ensuring sufficient memory, MIRACL can handle numbers with tens of thousands of bits or more—far larger than needed for most cryptographic applications.

In summary, the MIRACL library can work with numbers ranging from a few hundred bits to tens of thousands of bits (or thousands of decimal digits), depending on configuration and available memory. For cryptographic purposes, it’s typically used with numbers in the 256-bit to 8192-bit range, but its design allows for much larger numbers if properly initialized—making it a powerful tool for big number arithmetic beyond just standard cryptography needs.

Nahoru Odpovědět

5. března 23:53

Nejčastěji vyhledáváné

Diskuze: Ověření největšího nalezeného prvočísla - 41 024 320 cifer