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-2023 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-genie.h"
  28  #include "a68g-prelude.h"
  29  #include "a68g-double.h"
  30  #include "a68g-mp.h"
  31  
  32  //! @brief Return "pi" with "digs" precision, using Borwein & Borwein AGM.
  33  
  34  MP_T *mp_pi (NODE_T * p, MP_T * api, int mod, int digs)
  35  {
  36    int gdigs = FUN_DIGITS (digs);
  37    if (gdigs > A68_MP (mp_pi_size)) {
  38  //
  39  // No luck with the kept value, hence we generate a longer "pi".
  40  // Calculate "pi" using the Borwein & Borwein AGM algorithm.
  41  // This AGM doubles the numbers of digits every iteration.
  42  //
  43      a68_free (A68_MP (mp_pi));
  44      a68_free (A68_MP (mp_half_pi));
  45      a68_free (A68_MP (mp_two_pi));
  46      a68_free (A68_MP (mp_sqrt_two_pi));
  47      a68_free (A68_MP (mp_sqrt_pi));
  48      a68_free (A68_MP (mp_ln_pi));
  49      a68_free (A68_MP (mp_180_over_pi));
  50      a68_free (A68_MP (mp_pi_over_180));
  51  //
  52      ADDR_T pop_sp = A68_SP;
  53      MP_T *pi_g = nil_mp (p, gdigs);
  54      MP_T *two = lit_mp (p, 2, 0, gdigs);
  55      MP_T *x_g = lit_mp (p, 2, 0, gdigs);
  56      MP_T *y_g = nil_mp (p, gdigs);
  57      MP_T *u_g = nil_mp (p, gdigs);
  58      MP_T *v_g = nil_mp (p, gdigs);
  59      (void) sqrt_mp (p, x_g, x_g, gdigs);
  60      (void) add_mp (p, pi_g, x_g, two, gdigs);
  61      (void) sqrt_mp (p, y_g, x_g, gdigs);
  62      BOOL_T iterate = A68_TRUE;
  63      while (iterate) {
  64  // New x.
  65        (void) sqrt_mp (p, u_g, x_g, gdigs);
  66        (void) rec_mp (p, v_g, u_g, gdigs);
  67        (void) add_mp (p, u_g, u_g, v_g, gdigs);
  68        (void) half_mp (p, x_g, u_g, gdigs);
  69  // New pi.
  70        (void) plus_one_mp (p, u_g, x_g, gdigs);
  71        (void) plus_one_mp (p, v_g, y_g, gdigs);
  72        (void) div_mp (p, u_g, u_g, v_g, gdigs);
  73        (void) mul_mp (p, v_g, pi_g, u_g, gdigs);
  74  // Done yet?.
  75        if (same_mp (p, v_g, pi_g, gdigs)) {
  76          iterate = A68_FALSE;
  77        } else {
  78          (void) move_mp (pi_g, v_g, gdigs);
  79  // New y.
  80          (void) sqrt_mp (p, u_g, x_g, gdigs);
  81          (void) rec_mp (p, v_g, u_g, gdigs);
  82          (void) mul_mp (p, u_g, y_g, u_g, gdigs);
  83          (void) add_mp (p, u_g, u_g, v_g, gdigs);
  84          (void) plus_one_mp (p, v_g, y_g, gdigs);
  85          (void) div_mp (p, y_g, u_g, v_g, gdigs);
  86        }
  87      }
  88  // Keep the result for future restore.
  89      (void) shorten_mp (p, api, digs, pi_g, gdigs);
  90      A68_MP (mp_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
  91      (void) move_mp (A68_MP (mp_pi), api, digs);
  92      A68_MP (mp_half_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
  93      (void) half_mp (p, A68_MP (mp_half_pi), api, digs);
  94      A68_MP (mp_sqrt_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
  95      (void) sqrt_mp (p, A68_MP (mp_sqrt_pi), api, digs);
  96      A68_MP (mp_ln_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
  97      (void) ln_mp (p, A68_MP (mp_ln_pi), api, digs);
  98      A68_MP (mp_two_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
  99      (void) mul_mp_digit (p, A68_MP (mp_two_pi), api, (MP_T) 2, digs);
 100      A68_MP (mp_sqrt_two_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
 101      (void) sqrt_mp (p, A68_MP (mp_sqrt_two_pi), A68_MP (mp_two_pi), digs);
 102      A68_MP (mp_pi_over_180) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
 103      (void) div_mp_digit (p, A68_MP (mp_pi_over_180), api, 180, digs);
 104      A68_MP (mp_180_over_pi) = (MP_T *) get_heap_space ((unt) SIZE_MP (digs));
 105      (void) rec_mp (p, A68_MP (mp_180_over_pi), A68_MP (mp_pi_over_180), digs);
 106      A68_MP (mp_pi_size) = gdigs;
 107      A68_SP = pop_sp;
 108    }
 109    switch (mod) {
 110    case MP_PI:          return move_mp (api, A68_MP (mp_pi), digs);
 111    case MP_HALF_PI:     return move_mp (api, A68_MP (mp_half_pi), digs);
 112    case MP_TWO_PI:      return move_mp (api, A68_MP (mp_two_pi), digs);
 113    case MP_SQRT_TWO_PI: return move_mp (api, A68_MP (mp_sqrt_two_pi), digs);
 114    case MP_SQRT_PI:     return move_mp (api, A68_MP (mp_sqrt_pi), digs);
 115    case MP_LN_PI:       return move_mp (api, A68_MP (mp_ln_pi), digs);
 116    case MP_180_OVER_PI: return move_mp (api, A68_MP (mp_180_over_pi), digs);
 117    case MP_PI_OVER_180: return move_mp (api, A68_MP (mp_pi_over_180), digs);
 118    default:             return NaN_MP; // Should not be here.
 119    }
 120  }
 121