mp-pi.c

     
   1  //! @file mp-pi.c
   2  //! @author J. Marcel van der Veer
   3  
   4  //! @section Copyright
   5  //!
   6  //! This file is part of Algol68G - an Algol 68 compiler-interpreter.
   7  //! Copyright 2001-2024 J. Marcel van der Veer [algol68g@xs4all.nl].
   8  
   9  //! @section License
  10  //!
  11  //! This program is free software; you can redistribute it and/or modify it 
  12  //! under the terms of the GNU General Public License as published by the 
  13  //! Free Software Foundation; either version 3 of the License, or 
  14  //! (at your option) any later version.
  15  //!
  16  //! This program is distributed in the hope that it will be useful, but 
  17  //! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 
  18  //! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for 
  19  //! more details. You should have received a copy of the GNU General Public 
  20  //! License along with this program. If not, see [http://www.gnu.org/licenses/].
  21  
  22  //! @section Synopsis
  23  //!
  24  //! [LONG] LONG REAL value of pi by AGM.
  25  
  26  #include "a68g.h"
  27  #include "a68g-mp.h"
  28  
  29  //! @brief Return "pi" with "digs" precision, using a Borwein & Borwein AGM.
  30  
  31  MP_T *mp_pi (NODE_T * p, MP_T * api, int mod, int digs)
  32  {
  33    int gdigs = FUN_DIGITS (digs);
  34    if (gdigs > A68_MP (mp_pi_size)) {
  35  // No luck with the cached value, so we compute a longer "pi".
  36  // Calculate "pi" using a Borwein & Borwein AGM algorithm that doubles 
  37  // the numbers of digits on each iteration.
  38  //
  39  // J. M. Borwein, P. B. Borwein. "Pi and the AGM – A Study in Analytic Number
  40  // Theory and Computational Complexity", Wiley, New York, 1987.
  41      a68_free (A68_MP (mp_pi));
  42      a68_free (A68_MP (mp_half_pi));
  43      a68_free (A68_MP (mp_two_pi));
  44      a68_free (A68_MP (mp_sqrt_two_pi));
  45      a68_free (A68_MP (mp_sqrt_pi));
  46      a68_free (A68_MP (mp_ln_pi));
  47      a68_free (A68_MP (mp_180_over_pi));
  48      a68_free (A68_MP (mp_pi_over_180));
  49      ADDR_T pop_sp = A68_SP;
  50      MP_T *pi_g = nil_mp (p, gdigs);
  51      MP_T *two = lit_mp (p, 2, 0, gdigs);
  52      MP_T *x_g = lit_mp (p, 2, 0, gdigs);
  53      MP_T *y_g = nil_mp (p, gdigs);
  54      MP_T *u_g = nil_mp (p, gdigs);
  55      MP_T *v_g = nil_mp (p, gdigs);
  56      (void) sqrt_mp (p, x_g, x_g, gdigs);
  57      (void) add_mp (p, pi_g, x_g, two, gdigs);
  58      (void) sqrt_mp (p, y_g, x_g, gdigs);
  59      BOOL_T iterate = A68_TRUE;
  60      while (iterate) {
  61  // New x.
  62        (void) sqrt_mp (p, u_g, x_g, gdigs);
  63        (void) rec_mp (p, v_g, u_g, gdigs);
  64        (void) add_mp (p, u_g, u_g, v_g, gdigs);
  65        (void) half_mp (p, x_g, u_g, gdigs);
  66  // New pi.
  67        (void) plus_one_mp (p, u_g, x_g, gdigs);
  68        (void) plus_one_mp (p, v_g, y_g, gdigs);
  69        (void) div_mp (p, u_g, u_g, v_g, gdigs);
  70        (void) mul_mp (p, v_g, pi_g, u_g, gdigs);
  71  // Done yet?.
  72        if (same_mp (p, v_g, pi_g, gdigs)) {
  73          iterate = A68_FALSE;
  74        } else {
  75          (void) move_mp (pi_g, v_g, gdigs);
  76  // New y.
  77          (void) sqrt_mp (p, u_g, x_g, gdigs);
  78          (void) rec_mp (p, v_g, u_g, gdigs);
  79          (void) mul_mp (p, u_g, y_g, u_g, gdigs);
  80          (void) add_mp (p, u_g, u_g, v_g, gdigs);
  81          (void) plus_one_mp (p, v_g, y_g, gdigs);
  82          (void) div_mp (p, y_g, u_g, v_g, gdigs);
  83        }
  84      }
  85  // Keep the result for future restore.
  86      (void) shorten_mp (p, api, digs, pi_g, gdigs);
  87      A68_MP (mp_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
  88      (void) move_mp (A68_MP (mp_pi), api, digs);
  89      A68_MP (mp_half_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
  90      (void) half_mp (p, A68_MP (mp_half_pi), api, digs);
  91      A68_MP (mp_sqrt_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
  92      (void) sqrt_mp (p, A68_MP (mp_sqrt_pi), api, digs);
  93      A68_MP (mp_ln_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
  94      (void) ln_mp (p, A68_MP (mp_ln_pi), api, digs);
  95      A68_MP (mp_two_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
  96      (void) mul_mp_digit (p, A68_MP (mp_two_pi), api, (MP_T) 2, digs);
  97      A68_MP (mp_sqrt_two_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
  98      (void) sqrt_mp (p, A68_MP (mp_sqrt_two_pi), A68_MP (mp_two_pi), digs);
  99      A68_MP (mp_pi_over_180) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
 100      (void) div_mp_digit (p, A68_MP (mp_pi_over_180), api, 180, digs);
 101      A68_MP (mp_180_over_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
 102      (void) rec_mp (p, A68_MP (mp_180_over_pi), A68_MP (mp_pi_over_180), digs);
 103      A68_MP (mp_pi_size) = gdigs;
 104      A68_SP = pop_sp;
 105    }
 106    switch (mod) {
 107    case MP_PI: {
 108        return move_mp (api, A68_MP (mp_pi), digs);
 109      }
 110    case MP_HALF_PI: {
 111        return move_mp (api, A68_MP (mp_half_pi), digs);
 112      }
 113    case MP_TWO_PI: {
 114        return move_mp (api, A68_MP (mp_two_pi), digs);
 115      }
 116    case MP_SQRT_TWO_PI: {
 117        return move_mp (api, A68_MP (mp_sqrt_two_pi), digs);
 118      }
 119    case MP_SQRT_PI: {
 120        return move_mp (api, A68_MP (mp_sqrt_pi), digs);
 121      }
 122    case MP_LN_PI: {
 123        return move_mp (api, A68_MP (mp_ln_pi), digs);
 124      }
 125    case MP_180_OVER_PI: {
 126        return move_mp (api, A68_MP (mp_180_over_pi), digs);
 127      }
 128    case MP_PI_OVER_180: {
 129        return move_mp (api, A68_MP (mp_pi_over_180), digs);
 130      }
 131    default: {
 132        return NaN_MP; // Should not be here.
 133      }
 134    }
 135  }
 136  
     


© 2002-2024 J.M. van der Veer (jmvdveer@xs4all.nl)