poniedziałek, 14 lutego 2022

Palindromiczne liczby pierwsze (Palindromic prime numbers)

Naiwnie
 select( i -> i == fromdigits( Vecrev( digits( i ) )), primes( [ 2, 1000 ] ) )

Lepiej

  1. p10( n ) = 10^n;
  2.  
  3. rew( m, c ) = {
  4.     local( t, n );
  5.     t = 0; n = m;
  6.     for(i=1, c,
  7.         t = t*10 + n%10;
  8.         n \= 10 );
  9.     return( t ) }
  10.  
  11. range( p, w, disp = 0 ) = {
  12.     local( w10, mi, mj, z, pal, q ,= -1);
  13.     w10 = p * p10( w ) + p;
  14.     mi  = p10( w \ 2 + 1 );
  15.     mj  = p10( w \ 2     );
  16.     z   = p10( w \ 2 - 1 ) - 1;
  17.     for( i = 0, z,
  18.         pal = rew( i, w\2 );
  19.         q = w10  + i * mi + pal;
  20.         for( j = 0, 9,
  21.             if( isprime(+ j * mj ),
  22.                 k++;
  23.                 if( disp,
  24.                    if((% 8)==0,print());
  25.                    print1( q + j * mj, "\t") ) ) ) );
  26.     return( [ k+1, q + 9*mj ]); }

  27. gener( disp=0 ) = {
  28.     local( t=[ 1, 3, 7, 9], s=5, x,start );
  29.     start = getabstime();
  30.     for( w = 1, 8,
  31.         for( i = 1, 20 - 2*w, print1(" "));
  32.         print1( p10(w*2));
  33.         for( i = 1, 4,
  34.             print1(".");
  35.             x=range(t[i], w*2, disp);
  36.             s+=x[1]; );
  37.     printf( "\t %20d  # %8d   %8.3g [sec]\n", 
  38.                 x[2] s,  (getabstime()-start)/1000.0) )
  39. }

czwartek, 7 października 2021

Combination list

 Crv ( k, v, d ) = {
    if( d == k,
        print ( vecextract( v , "2..-2" ) )
    ,   
        for( i = v[ d + 1 ] + 1, #v,
            v[ d + 2 ] = i;
            Crv( k, v, d + 1 ) ));

combRV( n, k ) = Crv ( k, vector( n, X, X-1), 0 );
\\--------------------------------------------------------------
Cr ( c, z, b, n, k ) = {
    if( z < b,  print1( c, " " );
        if( n>0, Cr( c+1, z  , b*  k  \n, n-1, k - 1 ))
    ,
        if( n>0, Cr( c+1, z-b, b*(n-k)\n, n-1, k     ))
    );
}

combR( n, k ) = { 
    local( 
        bnk = binomial( n,   k ),
        b11 = bnk * k \ n );          \\binomial( n-1, k-1 ) );
    for( z = 0,  bnk - 1,
        Cr( 1, z, b11, n-1, k-1 );
        print
    );
}
\\--------------------------------------------------------------
Ci( z, b, n, k ) = { local( c = 1 );
    n--; k--;
    while( k >= 0 ,
        if( z < b,
            print1(c, "-");
            c++;
            if( n > 0, 
                b = b*k \ n);
            n--; k--;
        ,
            c++;
            z -= b;
            b = b*(n-k)\n;
            n--
        )
    );
}

combI( n, k ) = { 
    local(  bnk = binomial( n, k ), 
            b11 = bnk * k \ n );            \\ binomial( n-1, k-1 )
    for( z = 0, bnk - 1,
        Ci(z,   b11,   n, k );
        print
    );
}

poniedziałek, 19 kwietnia 2021

Dzielniki liczby (naiwnie)

 sumaD=0; liczbaD = 0; 

dzielniki( n ) = { local( j, d=2, g=n );
sumaD = 1; liczbaD = 1;
while( d < g, 
j = n\d; 
if( d * j == n,
sumaD   += d + j; 
liczbaD += 2 );
d++ ;  g = j
);
if( g*j == n,  \\ n == x^2
liczbaD--; sumaD -= j )
};

sobota, 23 lutego 2019

Piękna myśl

Jeśli chcesz się nauczyć - czytaj i słuchaj,
Jeśli chcesz coś zrozumieć - pisz o tym,
Jeśli chcesz osiągnąć mistrzostwo - nauczaj o tym.

piękna myśl na dziś.

ludzie podchodzą do tego od wszystkich czterech możliwych stron
tylko czterech?

jedno chciałem zrozumieć, a tylko słuchałem i czytałem
czegoś chciałem się nauczyć, tylko pisałem o tym i w tym
chciałem być mistrzem, taaa nauczałem,

nie silnia ale klika (grafy)
można by zapytać OEIS ile jest dróg, ,,,,,,,

niedziela, 20 maja 2012

najmłodsza cyfra silni


int ld(int n){
        int q, t, x, z, ai;
        q=0; t=0; x=0;
        if( ((n%5)&1)==0 ) t=n%5;
        n /= 5;
        while(n){
                ai=n%5; n/=5;
                q+=ai; x+=q;
                if((ai&1)==0)
                        t += ai; }
        z = (x+t/2)%4; 
        if( z==0 ) return 6;
        else return 1<<z;}

niedziela, 13 maja 2012

binomial


unsigned long binomial (int n, int m) { 
        unsigned long cnm = 1UL; 
        int i, f; 
        if (m*2 >n) 
                m = n-m; 
        for (i=1 ; i <= m; n--, i++) { 
                if ((f=n) % i == 0) 
                        f /= i; 
                else 
                        cnm /= i; 
                cnm *= f; } 
        return cnm; }

niedziela, 8 stycznia 2012

Sito Eratostenesa

#include <iostream>
#include <bitset>
using namespace std;
 
#define N 590000000
bitset <N/2+1> t;
 
int main(){
        int i,j,k;
        printf("N=%d\n", N);
        for( i=1,k=3; k*k<N; k+=2,i++) 
                if( !t.test(i) )
                        for( j=(i*i+i)*2; j+j<N; j += 2*i+1 ) 
                                t.set(j);
        cout << t.size()-t.count()-1;
        return 0;
}

sobota, 7 stycznia 2012

Podzbiory

main() {
        int i, k=1, n, x[100];
        scanf("%d", &n);
        x[1] = 1;
        while (k) {
                for( i=1; i <= k; i++ ) 
                        printf("%d ", x[i]);
                puts("");
                if ( x[k] == n ) {
                        k--;
                        x[k]++;
                } else {
                        k++;
                        x[k] = x[k-1] + 1;}}}

Zamiana

a=(a+b)-(b=a);

piątek, 6 stycznia 2012

Partycje

#include <stdio.h>
 
void denew(int B, int * N, int * K){
        int n,k,i,b;
        for( n=1; n<B; n++ ){
                b=i=n;
                for( k=2; 2*k<=n;  k++ ) {
                        b = b*(--i)/k;
                        if( b==B ) {
                                *N= n; *K= k; return; }
                        if( b>B ) break;
                }
        }
        *N=B; *K=1;
}
#define MAX 20
 
int  n,k,licznik=0,wsie=0;
int  p[MAX+2];
 
int min( int x, int y) { return x<y ? x : y;}
 
long long int silnia(int n){
   long long int s=1;
   while(n) s*=n--; return s;}
   
void PrintIt(int t, int z) {
   if( z>0 && t!=z) return;
   int i,pop; long long int perm;
   printf("%3d, %3d  | ", ++licznik, t);
   p[t+1]=0; pop=1; perm=silnia(t); 
   for(i=1; i<=t; i++){
       if(p[i]!=p[pop]){
           perm/=silnia(i-pop);
           pop=i;} 
       printf("%2d ", p[i]);}
   perm/=silnia(i-pop);
   while(i++<23) printf("   "); printf(" %20llu\n", perm);
   wsie+=perm;}
 
void P (int n, int k, int t, int z) {
   int j;
   p[t] = k;
   if (n==k) 
       PrintIt(t, z);
   for (j=min(k,n-k); j>=1; j--) 
       P(n-k,j,t+1,z);}
 
int main() {
   int i, z ,d;
   scanf("%d %d",&i, &z);
   P(i*2, i, 0, z);
   for(i=1;i<23;i++) printf("   ");
   denew(wsie, &i, &z);
   printf("Wszystkich do kupy %14d=newton(%d,%d)\n", wsie, i, z);
   z=wsie & (-wsie);
   if( z== wsie){
      for(i=1;i<23;i++) printf("   ");
      i=1; z=2;
      while( z != wsie) {
         i++; z*=2; }
   printf("                                 = 2 ^ %d\n", i);
   }
   return 0;}

Szybkie potęgowanie

double bpow(double x, long long int n){
        double r=1;
        double xf, xi;
        xf=modf(x, &xi);
        while( 1 ){
                if( n & 1 ){
                        r = r*xi + r*xf; }
                if( n>>=1 ) {
                        xf  = xf*xf + 2*xi*xf;
                        xi *= xi; }
                else 
                        return r; }}