API tools faq
paste
Advertisement
Le_BuG63

Pi

Le_BuG63
Sep 28th, 2013
122
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 11.55 KB | None | 0 0
raw download clone embed print report
  1. /*  CE PROGRAMME CALCULE JUSQU'A 150000 (OU 644690 SUIVANT TYPE DE MACHINE)
  2.     DECIMALES DE PI A L'AIDE DE LA FORMULE DE MACHIN:
  3.  
  4.             PI=16*ARCTG(1/5)-4*ARCTG(1/239)
  5.  
  6.     CHAQUE ARCTG EST LA SOMME DES Vn AVEC
  7.  
  8.             Vn=(-1)^n/((2*n+1)*X^(2*n+1))
  9.  
  10.     X VALANT RESPECTIVEMENT 5 ET 239. CHAQUE TERME Vn EST OBTENU PAR
  11.     RECURRENCE SUIVANT LES FORMULES:
  12.  
  13.             V0   = 1/X
  14.             Vn+1 = [(-1)*(2*n+1)/((2*n+3)*X*X)]*Vn.
  15.  
  16.     LES CALCULS SONT FAITS EN MULTIPRECISION SUR 32 BITS (RESP. 64 BITS)
  17.     LA CONTRAINTE DUE A LA LIMITATION DES ENTIERS SUR 32 BITS (RESP. 64
  18.     BITS) FAIT QUE LES CALCULS DES DEUX ARCTG SONT FAITS DIFFEREMMENT
  19.     DANS UN SOUCI D'OPTIMISATION.
  20.  
  21.     LES CALCULS SONT FAITS SUR DES ENTIERS SIGNES, DONC ON UTILISE
  22.     LES 31 BITS (RESP. 63 BITS) DE CES NOMBRES. LA BASE DE NUMERATION
  23.     UTILISEE EST 1 000 000 (RESP. 10^15) JUSQU'A 1500 (RESP. 6649)
  24.     DECIMALES, 100 000 (RESP. 10^14) JUSQU'A 15 000 (RESP. 66471) ET
  25.     10 000 (RESP. 10^13) JUSQU'A 150 000 (RESP. 644690) DECIMALES CECI
  26.     POUR EVITER DE DEPASSER LES 31 BITS (RESP. 63 BITS) IMPOSES.
  27.     IL Y A DONC UNE BAISSE DE PERFORMANCE AU DELA DE 1500 (RESP. 6649)
  28.     DECIMALES.
  29.  
  30.  
  31.     SUR UNE MACHINE A MOTS DE 64 BITS L'ALGORITHME EST PLUS PERFORMANT
  32.     CAR ON UTILISE UNE BASE DE NUMERATION BEAUCOUP PLUS GRANDE.
  33.     A LA COMPILATION IL FAUT DONC CHOISIR SON TYPE DE MACHINE:
  34.    
  35.     cc -DS64BITS -o pim2 pim2.c POUR MACHINE A MOTS DE 64 BITS
  36.     cc -DS32BITS -o pim2 pim2.c POUR MACHINE A MOTS DE 32 BITS
  37.  
  38.     SI L'ON VEUT TRAVAILLER AVEC DES MOTS DE 32 BITS SUR UNE MACHINE
  39.     A MOTS DE 64 BITS IL FAUT UTILISER L'OPTION DE COMPILATION:
  40.    
  41.     cc -DS32BITSFORCE -o pim2 pim2.c SUR LA MACHINE A MOTS DE 64 BITS.
  42.  
  43.     CELA PERMET DE COMPARER LES PERFORMANCES DES MACHINES AVEC LE
  44.     MEME ALGORITHME.
  45.  
  46.  
  47.         Usage: pim2 [-2] [entier]
  48.  
  49.             -2: sortie des decimales en base 2.
  50.             entier: nombre de decimales maxi.
  51.  
  52.     DANS LE CAS D'UNE SORTIE DES DECIMALES EN BASE 2, LA PERFORMANCE
  53.     DU PROGRAMME EST MOINS BONNE CAR LE CALCUL EST TOUJOURS FAIT EN
  54.     BASE 10 MAIS ON FAIT ENSUITE UNE CONVERSION EN BASE 2 QUI EST PLUS
  55.     LOURDE QU'UNE CONVERSION EN BASE 10. AVOIR PI EN BASE 2 N'A D'INTERET
  56.     QUE POUR OBTENIR UNE SUITE DE 0 ET DE 1 CONSIDERES COMME ALEATOIRE.
  57.  
  58.     fichier pim2.c  Andre' Brouty (CNRS / ENST de Bretagne) decembre 1987
  59.  
  60.         portage pour machine a mots de 64 bits: decembre 1995
  61.         portage pour sortie des decimales en base 2: fevrier 1996
  62.  
  63.     [ http://www.brouty.fr/Maths/pi.html ]      */
  64.  
  65.  
  66.  
  67. #ifdef S64BITS  /* MACHINE CALCULANT SUR 64 BITS */
  68. #define TypInt  long
  69. #define NBOCT   8           /* NOMBRE D'OCTETS SUR 64 BITS  */
  70. #define BASE1   1000000000000000L   /* BASE DE CALCUL: 10^15    */
  71. #define BASE2   100000000000000L    /* BASE DE CALCUL: 10^14    */
  72. #define BASE3   10000000000000L     /* BASE DE CALCUL: 10^13    */
  73. #define NBZ1    15          /* NOMBRE DE ZEROS DE BASE1 */
  74. #define NBZ2    14          /* NOMBRE DE ZEROS DE BASE2 */
  75. #define NBZ3    13          /* NOMBRE DE ZEROS DE BASE3 */
  76. #define NBIT2   13  /* 2^13 MULTPLICATEUR MINI POUR PASSER EN BASE 2*/
  77. #define GBASE1  4/* NOMBRE DE BLOCKS DE NBZ1 CHIFFRES AFFICHES PAR LIGNE*/
  78. #define GBASE2  4/* NOMBRE DE BLOCKS DE NBZ2 CHIFFRES AFFICHES PAR LIGNE*/
  79. #define GBASE3  5/* NOMBRE DE BLOCKS DE NBZ3 CHIFFRES AFFICHES PAR LIGNE*/
  80. #define LIMP1   6450    /* NOMBRE MAX DE DECIMALES CALCULEES POUR BASE1 */
  81. #define LIMP2   64471   /* NOMBRE MAX DE DECIMALES CALCULEES POUR BASE2 */
  82. #define LIMP3   644691  /* NOMBRE MAX DE DECIMALES CALCULEES POUR BASE3 */
  83. #endif
  84.  
  85. #ifdef S32BITS  /* MACHINE CALCULANT SUR 32 BITS */
  86. #define TypInt  long
  87. #define NBOCT   4           /* NOMBRE D'OCTETS SUR 32 BITS  */
  88. #define BASE1   1000000L
  89. #define BASE2   100000L
  90. #define BASE3   10000L
  91. #define NBZ1    6
  92. #define NBZ2    5
  93. #define NBZ3    4
  94. #define NBIT2   11  /* 2^11 MULTPLICATEUR MINI POUR PASSER EN BASE 2*/
  95. #define GBASE1  10/* NOMBRE DE BLOCKS DE NBZ1 CHIFFRES AFFICHES PAR LIGNE*/
  96. #define GBASE2  12/* NOMBRE DE BLOCKS DE NBZ2 CHIFFRES AFFICHES PAR LIGNE*/
  97. #define GBASE3  14/* NOMBRE DE BLOCKS DE NBZ3 CHIFFRES AFFICHES PAR LIGNE*/
  98. #define LIMP1   1501    /* NOMBRE MAX DE DECIMALES CALCULEES POUR BASE1 */
  99. #define LIMP2   15001   /* NOMBRE MAX DE DECIMALES CALCULEES POUR BASE2 */
  100. #define LIMP3   150001  /* NOMBRE MAX DE DECIMALES CALCULEES POUR BASE3 */
  101. #endif
  102.  
  103. /* ON PEUT FORCER LE CALCUL SUR 32 BITS SUR UNE MACHINE A MOTS DE 64 BITS */
  104. /* AVEC L'OTION DE COMPILATION SUIVANTE. CELA PERMET DE COMPARER LE GAIN  */
  105. /* DE PERFORMANCE QUAND ON UTILISE REELLEMENT LES 64 BITS DE LA MACHINE.  */
  106. /*  cc -DS32BITSFORCE -o pim2 pim2.c                  */
  107.  
  108. #ifdef S32BITSFORCE
  109. #define TypInt  int
  110. #define NBOCT   4           /* NOMBRE D'OCTETS SUR 32 BITS  */
  111. #define BASE1   1000000
  112. #define BASE2   100000
  113. #define BASE3   10000
  114. #define NBZ1    6
  115. #define NBZ2    5
  116. #define NBZ3    4
  117. #define NBIT2   11  /* 2^11 MULTPLICATEUR MINI POUR PASSER EN BASE 2*/
  118. #define GBASE1  10/* NOMBRE DE BLOCKS DE NBZ1 CHIFFRES AFFICHES PAR LIGNE*/
  119. #define GBASE2  12/* NOMBRE DE BLOCKS DE NBZ2 CHIFFRES AFFICHES PAR LIGNE*/
  120. #define GBASE3  14/* NOMBRE DE BLOCKS DE NBZ3 CHIFFRES AFFICHES PAR LIGNE*/
  121. #define LIMP1   1501    /* NOMBRE MAX DE DECIMALES CALCULEES POUR BASE1 */
  122. #define LIMP2   15001   /* NOMBRE MAX DE DECIMALES CALCULEES POUR BASE2 */
  123. #define LIMP3   150001  /* NOMBRE MAX DE DECIMALES CALCULEES POUR BASE3 */
  124. #endif
  125.  
  126. #define LOGb10_2 0.30103/* LOG BASE 10 DE 2 POUR CONVERSION EN BINAIRE */
  127.  
  128. TypInt *resul;  /* TROIS TABLEAUX DYNAMIQUES CONTENANT RESPECTIVEMENT */
  129. TypInt *resul1; /* LE RESULTAT FINAL ET LES RESULTATS INTERMEDIAIRES   */
  130. TypInt *resul2;
  131. int K;      /* GESTION DES CASES NULLES: INDICE DE LA 1ERE CASE NON NULLE */
  132. int MAX, MAXB;  /* DIMENSION DES TABLEAUX */
  133. TypInt BASE;    /* BASE DE NUMERATION CORRESPOND AU GROUPEMENT DES CHIFFRES */
  134. int GBASE;  /* NOMBRE DE GROUPEMENT DES CHIFFRES POUR LA SORTIE */
  135. int I, J;   /* NOMBRE D'ITERATIONS POUR CHAQUE ARCTG */
  136.  
  137. main(argc,argv)
  138. int argc;
  139. char *argv[];
  140.  
  141. {
  142.     int k, i, B2;
  143.     int lim, limb;  /* NOMBRE DE DECIMALES DEMANDEES */
  144.     int NB;
  145.     int Dec;    /* 2^Dec = 1 << Dec */
  146.     char *nam;
  147.     nam = argv[0];
  148.     B2 = 0;
  149.     --argc; ++argv;
  150.     lim = 0;
  151.     while(argc) {
  152.         if(*argv[0] == '-') {
  153.             switch(*(argv[0] + 1)) {
  154.                 case '2':
  155.                     B2 = 1;
  156.                     break;
  157.                 default:
  158.                     Us(nam);
  159.                     break;
  160.             }
  161.         }
  162.         else
  163.             lim=atoi(argv[0]);
  164.         --argc; ++argv;
  165.     }
  166.     if(lim == 0) {
  167.         printf("Nombre de decimales souhaitees ( <= %d): ",LIMP3-1);
  168.         scanf("%d",&lim);
  169.     }
  170.     if(B2) {
  171.         limb = lim;
  172.         lim = (int)((double)lim * LOGb10_2);
  173.     }
  174.     if(lim < LIMP1) {
  175.         BASE=BASE1;
  176.         NB=NBZ1;
  177.         GBASE = B2 ? 1 + GBASE1/2 : GBASE1;
  178.         Dec = NBIT2;
  179.     }
  180.     else {
  181.         if(lim > LIMP3 - 1) {
  182.             printf("Pas plus de %d decimales merci.\n",LIMP3-1);
  183.             exit(0);
  184.         }
  185.         if(lim < LIMP2) {
  186.             BASE=BASE2;
  187.             NB=NBZ2;
  188.             GBASE = B2 ? GBASE1/2 - 1 : GBASE2;
  189.             Dec = NBIT2 + 3;
  190.         }
  191.         else {
  192.             BASE=BASE3;
  193.             NB=NBZ3;
  194.             GBASE = B2 ? GBASE1/2 - 1 : GBASE3;
  195.             Dec = NBIT2 + 6;
  196.         }
  197.        
  198.     }
  199.     if(B2)
  200.         MAXB = 5+limb/Dec;
  201.     MAX = 5+B2+lim/NB;
  202.     resul=(TypInt *)malloc(NBOCT*MAX);
  203.     resul1=(TypInt *)malloc(NBOCT*MAX);
  204.     resul2=(TypInt *)malloc(NBOCT*MAX);
  205.     memset((char *)resul1,0, NBOCT*MAX);
  206.     memset((char *)resul2,0, NBOCT*MAX);
  207.     memset((char *)resul,0, NBOCT*MAX);
  208.     I=arctg_5(resul1,5,16);
  209.     memcpy((char *)resul2,(char *)resul, NBOCT*MAX);
  210.     memset((char *)resul,0, NBOCT*MAX);
  211. #ifdef S64BITS
  212.     if(lim < LIMP1)
  213.         J=arctg_239(resul1,239,4);
  214.     else
  215.         J=arctg_5(resul1,239,4);    /* ON OPTIMISE:239^2 CA PASSE */
  216. #else
  217.     if(lim < LIMP2)
  218.         J=arctg_239(resul1,239,4);
  219.     else
  220.         J=arctg_5(resul1,239,4);    /* ON OPTIMISE:239^2 CA PASSE */
  221. #endif
  222.     K=2;
  223.     sous(resul2,resul);
  224.     if(B2) {
  225.         printf("\n Decimales affichees: %d  nombre d'iterations: %d , %d (%d)\n\n",Dec*(MAXB-4),I,J,I+J);
  226.         tobin(resul2,Dec);
  227.     }
  228.     else {
  229.         printf("\n Decimales affichees: %d  nombre d'iterations: %d , %d (%d)\n\n",NB*(MAX-4),I,J,I+J);
  230.         todec(resul2);
  231.     }
  232. }
  233.  
  234. todec(tab)  /* CONVERSION EN DECIMAL */
  235. TypInt tab[];
  236.  
  237. {
  238.     register i;
  239.     TypInt a;
  240.     TypInt BS10 = BASE/10;
  241.     printf("pi = %d,",tab[2]);
  242.     for(i=3;i<MAX-1;++i) {
  243.         a=tab[i];
  244.         a = a == 0 ? 1L : a;
  245.         while(a < BS10) {
  246.             putchar('0');
  247.             a=10*a;
  248.         }
  249.         printf("%lu%s",tab[i],(i-2)%GBASE ? " " : "\n       ");
  250.     }
  251.     putchar('\n');
  252. }
  253.  
  254. tobin(tab,nb)   /* CONVERSION EN BINAIRE */
  255. TypInt tab[];
  256. int nb;
  257.  
  258. {
  259.     register i;
  260.     TypInt a;
  261.     char *bas2();
  262.     a = (TypInt)(1 << (nb - 1));
  263.     printf("pi= 11,");
  264.     for(i=3;i<MAXB-1;++i) {
  265.         tab[0] = 0;
  266.         tab[1] = 0;
  267.         tab[2] = 0;
  268.         mul2(tab,nb);
  269.         if(tab[1])
  270.             tab[2] = BASE*BASE*tab[0] + BASE*tab[1] + tab[2];
  271.         printf("%s%s",bas2(tab[2],a,nb),(i-2)%GBASE ? " " : "\n       ");
  272.     }
  273.     putchar('\n');
  274. }
  275.  
  276. char *bas2(b,p,n)/* AFFICHE LES CHIFFRES BINAIRES DE b PAR BLOC DE n CHIFFRES */
  277. TypInt b, p;
  278. int n;
  279.  
  280. {
  281.     short i;
  282.     static char Base[256];
  283.     memset(Base,0,256);
  284.     for(i = 0; i < n ; ++i) {
  285.         Base[i] = (b & (p >> i)) > 0 ? '1' : '0';
  286.     }
  287.     return(Base);
  288. }
  289.  
  290. add(a,b)    /* ADDITION DE DEUX TABLEAUX DE NOMBRES */
  291. TypInt a[], b[];
  292.  
  293. {
  294.     register i;
  295.     TypInt register r=0;        /* r EST LA RETENUE */
  296.     for(i=MAX-1;i >= K;--i) {
  297.         r=a[i]+r+b[i];
  298.         if(r >= BASE) {
  299.             a[i]=r-BASE;
  300.             r=1;
  301.         }
  302.         else {
  303.             a[i]=r;
  304.             r=0;
  305.         }
  306.     }
  307.     while(r) {  /* PROPAGATION DE LA RETENUE */
  308.         r=a[i]+r;
  309.         if(r >= BASE) {
  310.             a[i]=r-BASE;
  311.             r=1;
  312.         }
  313.         else {
  314.             a[i]=r;
  315.             r=0;
  316.         }
  317.         --i;
  318.     }
  319. }
  320.  
  321. sous(a,b)   /* SOUSTRACTION DE DEUX TABLEAUX DE NOMBRES. LE SIGNE DU */
  322. TypInt a[], b[];    /* RESULTAT EST SUPPOSE POSITIF              */
  323.  
  324. {
  325.     register i;
  326.     TypInt register r=0;
  327.     for(i=MAX-1;i >= K;--i) {
  328.         r=a[i]-r-b[i];
  329.         if(r < 0) {
  330.             a[i]=r+BASE;
  331.             r=1;
  332.         }
  333.         else {
  334.             a[i]=r;
  335.             r=0;
  336.         }
  337.     }
  338.     while(r) {  /* PROPAGATION DE LA RETENUE */
  339.         r=a[i]-r;
  340.         if(r < 0) {
  341.             a[i]=r+BASE;
  342.             r=1;
  343.         }
  344.         else {
  345.             a[i]=r;
  346.             r=0;
  347.         }
  348.         --i;
  349.     }
  350. }
  351.  
  352. mul(tab,n)     /* MULTIPLICATION D'UN TABLEAU DE NOMBRES PAR UN ENTIER POSITIF*/
  353. TypInt tab[], n;
  354.  
  355. {
  356.     register i;
  357.     TypInt register S=0;
  358.     for(i=MAX-1;i>=K;--i) {
  359.         tab[i]=tab[i]*n+S;
  360.         S=tab[i]/BASE;
  361.         tab[i] = tab[i]%BASE;
  362.     }
  363.     tab[i]=tab[i]*n+S;
  364.     while((S=tab[i])>=BASE) {
  365.         tab[i-1] = tab[i-1] + S/BASE;
  366.         tab[i] = S%BASE;
  367.         --i; --K;
  368.     }
  369.     --K;
  370.     while(tab[K] == 0) ++K;
  371. }
  372.  
  373. mul2(tab,n)     /* MULTIPLICATION D'UN TABLEAU DE NOMBRES PAR 2^n */
  374. TypInt tab[], n;
  375.  
  376. {
  377.     register i;
  378.     TypInt register S=0;
  379.     for(i=MAX-1;i>=K;--i) {
  380.         tab[i] = (tab[i] << n) + S;
  381.         S=tab[i]/BASE;
  382.         tab[i] = tab[i]%BASE;
  383.     }
  384.     tab[i] = (tab[i] << n) + S;
  385.     while((S=tab[i])>=BASE) {
  386.         tab[i-1] = tab[i-1] + S/BASE;
  387.         tab[i] = S%BASE;
  388.         --i; --K;
  389.     }
  390.     --K;
  391.     while(tab[K] == 0) ++K;
  392. }
  393.  
  394. divisi(t,n) /* DIVISION D'UN TABLEAU DE NOMBRES PAR UN ENTIER */
  395. TypInt t[], n;
  396.  
  397. {
  398.     register i;
  399.     TypInt register a;
  400.     a=t[K];
  401.     for(i=K;i<MAX;++i) {
  402.         t[i]=a/n;
  403.         a=BASE*(a%n)+t[i+1];
  404.     }
  405.     while(t[K] == 0) ++K;   /* GESTION DES CASES NULLES */
  406. }
  407.  
  408. arctg_239(tab,p,b)  /* b*arctg(1/p) POUR p DE L'ORDRE DE GRANDEUR DE 239 */
  409. TypInt tab[], p, b;
  410.  
  411. {
  412.     register i;
  413.     K=2;
  414.     tab[2]=b;
  415.     divisi(tab,p);
  416.     add(resul,tab);
  417.     for(i=0;K < MAX;++i) {
  418.         divisi(tab,(2*i+3));
  419.         divisi(tab,p);/*ON FAIT 2 DIVISIONS SUCCESSIVES AU LIEU D'UNE*/
  420.         divisi(tab,p);/*POUR EVITER DE DEPASSER LES 31 (OU 63) BITS  */
  421.         mul(tab,(2*i + 1));
  422.         (i+1)%2 ? sous(resul,tab) : add(resul,tab);
  423.     }
  424.     return(i);
  425. }
  426.  
  427. arctg_5(tab,p,b)    /* b*arctg(1/p) POUR p DE L'ORDRE DE GRANDEUR DE 5 */
  428. TypInt tab[], p, b;
  429.  
  430. {
  431.     register i;
  432.     K=2;
  433.     tab[2]=b;
  434.     divisi(tab,p);
  435.     add(resul,tab);
  436.     p=p*p;
  437.     for(i=0;K < MAX;++i) {
  438.         divisi(tab,(2*i+3));
  439.         divisi(tab,p);/*ICI p EST TROP PETIT POUR RISQUER DE DEPASSER */
  440.         mul(tab,(2*i + 1));     /*       LES 31 (OU 63) BITS.         */
  441.         (i+1)%2 ? sous(resul,tab) : add(resul,tab);
  442.     }
  443.     return(i);
  444. }
  445.  
  446. Us(nom) /* USAGE DE LA COMMANDE */
  447. char *nom;
  448.  
  449. {
  450.     printf("Usage: %s [-2] [entier]\n\n\t-2: sortie en base 2\n\tentier: nombre de decimales desirees.\n",nom);
  451.     exit(1);
  452. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!