v0.9.0
gm_rule.c
Go to the documentation of this file.
1 # include <stdlib.h>
2 # include <stdio.h>
3 # include <math.h>
4 # include <time.h>
5 
6 # include "gm_rule.h"
7 
8 /******************************************************************************/
9 
10 void comp_next ( int n, int k, int a[], int *more, int *h, int *t )
11 
12 /******************************************************************************/
13 /*
14  Purpose:
15 
16  COMP_NEXT computes the compositions of the integer N into K parts.
17 
18  Discussion:
19 
20  A composition of the integer N into K parts is an ordered sequence
21  of K nonnegative integers which sum to N. The compositions (1,2,1)
22  and (1,1,2) are considered to be distinct.
23 
24  The routine computes one composition on each call until there are no more.
25  For instance, one composition of 6 into 3 parts is
26  3+2+1, another would be 6+0+0.
27 
28  On the first call to this routine, set MORE = FALSE. The routine
29  will compute the first element in the sequence of compositions, and
30  return it, as well as setting MORE = TRUE. If more compositions
31  are desired, call again, and again. Each time, the routine will
32  return with a new composition.
33 
34  However, when the LAST composition in the sequence is computed
35  and returned, the routine will reset MORE to FALSE, signaling that
36  the end of the sequence has been reached.
37 
38  This routine originally used a STATICE statement to maintain the
39  variables H and T. I have decided (based on an wasting an
40  entire morning trying to track down a problem) that it is safer
41  to pass these variables as arguments, even though the user should
42  never alter them. This allows this routine to safely shuffle
43  between several ongoing calculations.
44 
45 
46  There are 28 compositions of 6 into three parts. This routine will
47  produce those compositions in the following order:
48 
49  I A
50  - ---------
51  1 6 0 0
52  2 5 1 0
53  3 4 2 0
54  4 3 3 0
55  5 2 4 0
56  6 1 5 0
57  7 0 6 0
58  8 5 0 1
59  9 4 1 1
60  10 3 2 1
61  11 2 3 1
62  12 1 4 1
63  13 0 5 1
64  14 4 0 2
65  15 3 1 2
66  16 2 2 2
67  17 1 3 2
68  18 0 4 2
69  19 3 0 3
70  20 2 1 3
71  21 1 2 3
72  22 0 3 3
73  23 2 0 4
74  24 1 1 4
75  25 0 2 4
76  26 1 0 5
77  27 0 1 5
78  28 0 0 6
79 
80  Licensing:
81 
82  This code is distributed under the GNU LGPL license.
83 
84  Modified:
85 
86  02 July 2008
87 
88  Author:
89 
90  Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf.
91  C version by John Burkardt.
92 
93  Reference:
94 
95  Albert Nijenhuis, Herbert Wilf,
96  Combinatorial Algorithms for Computers and Calculators,
97  Second Edition,
98  Academic Press, 1978,
99  ISBN: 0-12-519260-6,
100  LC: QA164.N54.
101 
102  Parameters:
103 
104  Input, int N, the integer whose compositions are desired.
105 
106  Input, int K, the number of parts in the composition.
107 
108  Input/output, int A[K], the parts of the composition.
109 
110  Input/output, int *MORE.
111  Set MORE = FALSE on first call. It will be reset to TRUE on return
112  with a new composition. Each new call returns another composition until
113  MORE is set to FALSE when the last composition has been computed
114  and returned.
115 
116  Input/output, int *H, *T, two internal parameters needed for the
117  computation. The user should allocate space for these in the calling
118  program, include them in the calling sequence, but never alter them!
119 */
120 {
121  int i;
122 
123  if ( !( *more ) )
124  {
125  *t = n;
126  *h = 0;
127  a[0] = n;
128  for ( i = 1; i < k; i++ )
129  {
130  a[i] = 0;
131  }
132  }
133  else
134  {
135  if ( 1 < *t )
136  {
137  *h = 0;
138  }
139  *h = *h + 1;
140  *t = a[*h-1];
141  a[*h-1] = 0;
142  a[0] = *t - 1;
143  a[*h] = a[*h] + 1;
144  }
145 
146  *more = ( a[k-1] != n );
147 
148  return;
149 }
150 /******************************************************************************/
151 
152 void gm_rule_set ( int rule, int dim_num, int point_num, double w[],
153  double x[] )
154 
155 /******************************************************************************/
156 /*
157  Purpose:
158 
159  GM_RULE_SET sets a Grundmann-Moeller rule.
160 
161  Discussion:
162 
163  This is a revised version of the calculation which seeks to compute
164  the value of the weight in a cautious way that avoids intermediate
165  overflow. Thanks to John Peterson for pointing out the problem on
166  26 June 2008.
167 
168  This rule returns weights and abscissas of a Grundmann-Moeller
169  quadrature rule for the DIM_NUM-dimensional unit simplex.
170 
171  The dimension POINT_NUM can be determined by calling GM_RULE_SIZE.
172 
173  Licensing:
174 
175  This code is distributed under the GNU LGPL license.
176 
177  Modified:
178 
179  30 August 2012
180 
181  Author:
182 
183  John Burkardt
184 
185  Reference:
186 
187  Axel Grundmann, Michael Moeller,
188  Invariant Integration Formulas for the N-Simplex
189  by Combinatorial Methods,
190  SIAM Journal on Numerical Analysis,
191  Volume 15, Number 2, April 1978, pages 282-290.
192 
193  Parameters:
194 
195  Input, int RULE, the index of the rule.
196  0 <= RULE.
197 
198  Input, int DIM_NUM, the spatial dimension.
199  1 <= DIM_NUM.
200 
201  Input, int POINT_NUM, the number of points in the rule.
202 
203  Output, double W[POINT_NUM], the weights.
204 
205  Output, double X[DIM_NUM*POINT_NUM], the abscissas.
206 */
207 {
208  int *beta;
209  int beta_sum;
210  int d;
211  int dim;
212  int h;
213  int i;
214  int j;
215  int j_hi;
216  int k;
217  int more;
218  int n;
219  double one_pm;
220  int s;
221  int t;
222  double weight;
223 
224  s = rule;
225  d = 2 * s + 1;
226  k = 0;
227  n = dim_num;
228  one_pm = 1.0;
229 
230  beta = ( int * ) malloc ( ( dim_num + 1 ) * sizeof ( int ) );
231 
232  for ( i = 0; i <= s; i++ )
233  {
234  weight = ( double ) one_pm;
235 
236  j_hi = i4_max ( n, i4_max ( d, d + n - i ) );
237 
238  for ( j = 1; j <= j_hi; j++ )
239  {
240  if ( j <= n )
241  {
242  weight = weight * ( double ) ( j );
243  }
244  if ( j <= d )
245  {
246  weight = weight * ( double ) ( d + n - 2 * i );
247  }
248  if ( j <= 2 * s )
249  {
250  weight = weight / 2.0;
251  }
252  if ( j <= i )
253  {
254  weight = weight / ( double ) ( j );
255  }
256  if ( j <= d + n - i )
257  {
258  weight = weight / ( double ) ( j );
259  }
260  }
261 
262  one_pm = - one_pm;
263 
264  beta_sum = s - i;
265  more = 0;
266  h = 0;
267  t = 0;
268 
269  for ( ; ; )
270  {
271  comp_next ( beta_sum, dim_num + 1, beta, &more, &h, &t );
272 
273  w[k] = weight;
274  for ( dim = 0; dim < dim_num; dim++ )
275  {
276  x[dim+k*dim_num] = ( double ) ( 2 * beta[dim+1] + 1 )
277  / ( double ) ( d + n - 2 * i );
278  }
279  k = k + 1;
280 
281  if ( !more )
282  {
283  break;
284  }
285  }
286  }
287 
288  free ( beta );
289 
290  return;
291 }
292 /******************************************************************************/
293 
294 int gm_rule_size ( int rule, int dim_num )
295 
296 /******************************************************************************/
297 /*
298  Purpose:
299 
300  GM_RULE_SIZE determines the size of a Grundmann-Moeller rule.
301 
302  Discussion:
303 
304  This rule returns the value of POINT_NUM, the number of points associated
305  with a GM rule of given index.
306 
307  After calling this rule, the user can use the value of POINT_NUM to
308  allocate space for the weight vector as W(POINT_NUM) and the abscissa
309  vector as X(DIM_NUM,POINT_NUM), and then call GM_RULE_SET.
310 
311  Licensing:
312 
313  This code is distributed under the GNU LGPL license.
314 
315  Modified:
316 
317  30 August 2012
318 
319  Author:
320 
321  John Burkardt
322 
323  Reference:
324 
325  Axel Grundmann, Michael Moeller,
326  Invariant Integration Formulas for the N-Simplex
327  by Combinatorial Methods,
328  SIAM Journal on Numerical Analysis,
329  Volume 15, Number 2, April 1978, pages 282-290.
330 
331  Parameters:
332 
333  Input, int RULE, the index of the rule.
334  0 <= RULE.
335 
336  Input, int DIM_NUM, the spatial dimension.
337  1 <= DIM_NUM.
338 
339  Output, int GM_RULE_SIZE, the number of points in the rule.
340 */
341 {
342  int arg1;
343  int point_num;
344 
345  arg1 = dim_num + rule + 1;
346 
347  point_num = i4_choose ( arg1, rule );
348 
349  return point_num;
350 }
351 /******************************************************************************/
352 
353 int i4_choose ( int n, int k )
354 
355 /******************************************************************************/
356 /*
357  Purpose:
358 
359  I4_CHOOSE computes the binomial coefficient C(N,K).
360 
361  Discussion:
362 
363  The value is calculated in such a way as to avoid overflow and
364  roundoff. The calculation is done in integer arithmetic.
365 
366  The formula used is:
367 
368  C(N,K) = N! / ( K! * (N-K)! )
369 
370  Licensing:
371 
372  This code is distributed under the GNU LGPL license.
373 
374  Modified:
375 
376  03 May 2008
377 
378  Author:
379 
380  John Burkardt
381 
382  Reference:
383 
384  ML Wolfson, HV Wright,
385  Algorithm 160:
386  Combinatorial of M Things Taken N at a Time,
387  Communications of the ACM,
388  Volume 6, Number 4, April 1963, page 161.
389 
390  Parameters:
391 
392  Input, int N, K, are the values of N and K.
393 
394  Output, int I4_CHOOSE, the number of combinations of N
395  things taken K at a time.
396 */
397 {
398  int i;
399  int mn;
400  int mx;
401  int value;
402 
403  mn = i4_min ( k, n - k );
404 
405  if ( mn < 0 )
406  {
407  value = 0;
408  }
409  else if ( mn == 0 )
410  {
411  value = 1;
412  }
413  else
414  {
415  mx = i4_max ( k, n - k );
416  value = mx + 1;
417 
418  for ( i = 2; i <= mn; i++ )
419  {
420  value = ( value * ( mx + i ) ) / i;
421  }
422  }
423 
424  return value;
425 }
426 /******************************************************************************/
427 
428 int i4_huge ( void )
429 
430 /******************************************************************************/
431 /*
432  Purpose:
433 
434  I4_HUGE returns a "huge" I4.
435 
436  Licensing:
437 
438  This code is distributed under the GNU LGPL license.
439 
440  Modified:
441 
442  29 August 2006
443 
444  Author:
445 
446  John Burkardt
447 
448  Parameters:
449 
450  Output, int I4_HUGE, a "huge" integer.
451 */
452 {
453  static int value = 2147483647;
454 
455  return value;
456 }
457 /******************************************************************************/
458 
459 int i4_max ( int i1, int i2 )
460 
461 /******************************************************************************/
462 /*
463  Purpose:
464 
465  I4_MAX returns the maximum of two I4's.
466 
467  Licensing:
468 
469  This code is distributed under the GNU LGPL license.
470 
471  Modified:
472 
473  29 August 2006
474 
475  Author:
476 
477  John Burkardt
478 
479  Parameters:
480 
481  Input, int I1, I2, are two integers to be compared.
482 
483  Output, int I4_MAX, the larger of I1 and I2.
484 */
485 {
486  int value;
487 
488  if ( i2 < i1 )
489  {
490  value = i1;
491  }
492  else
493  {
494  value = i2;
495  }
496  return value;
497 }
498 /******************************************************************************/
499 
500 int i4_min ( int i1, int i2 )
501 
502 /******************************************************************************/
503 /*
504  Purpose:
505 
506  I4_MIN returns the smaller of two I4's.
507 
508  Licensing:
509 
510  This code is distributed under the GNU LGPL license.
511 
512  Modified:
513 
514  29 August 2006
515 
516  Author:
517 
518  John Burkardt
519 
520  Parameters:
521 
522  Input, int I1, I2, two integers to be compared.
523 
524  Output, int I4_MIN, the smaller of I1 and I2.
525 */
526 {
527  int value;
528 
529  if ( i1 < i2 )
530  {
531  value = i1;
532  }
533  else
534  {
535  value = i2;
536  }
537  return value;
538 }
539 /******************************************************************************/
540 
541 int i4_power ( int i, int j )
542 
543 /******************************************************************************/
544 /*
545  Purpose:
546 
547  I4_POWER returns the value of I^J.
548 
549  Licensing:
550 
551  This code is distributed under the GNU LGPL license.
552 
553  Modified:
554 
555  23 October 2007
556 
557  Author:
558 
559  John Burkardt
560 
561  Parameters:
562 
563  Input, int I, J, the base and the power. J should be nonnegative.
564 
565  Output, int I4_POWER, the value of I^J.
566 */
567 {
568  int k;
569  int value;
570 
571  if ( j < 0 )
572  {
573  if ( i == 1 )
574  {
575  value = 1;
576  }
577  else if ( i == 0 )
578  {
579  fprintf ( stderr, "\n" );
580  fprintf ( stderr, "I4_POWER - Fatal error!\n" );
581  fprintf ( stderr, " I^J requested, with I = 0 and J negative.\n" );
582  exit ( 1 );
583  }
584  else
585  {
586  value = 0;
587  }
588  }
589  else if ( j == 0 )
590  {
591  if ( i == 0 )
592  {
593  fprintf ( stderr, "\n" );
594  fprintf ( stderr, "I4_POWER - Fatal error!\n" );
595  fprintf ( stderr, " I^J requested, with I = 0 and J = 0.\n" );
596  exit ( 1 );
597  }
598  else
599  {
600  value = 1;
601  }
602  }
603  else if ( j == 1 )
604  {
605  value = i;
606  }
607  else
608  {
609  value = 1;
610  for ( k = 1; k <= j; k++ )
611  {
612  value = value * i;
613  }
614  }
615  return value;
616 }
617 /******************************************************************************/
618 
619 double *monomial_value ( int dim_num, int point_num, double x[], int expon[] )
620 
621 /******************************************************************************/
622 /*
623  Purpose:
624 
625  MONOMIAL_VALUE evaluates a monomial.
626 
627  Discussion:
628 
629  This routine evaluates a monomial of the form
630 
631  product ( 1 <= dim <= dim_num ) x(dim)^expon(dim)
632 
633  where the exponents are nonnegative integers. Note that
634  if the combination 0^0 is encountered, it should be treated
635  as 1.
636 
637  Licensing:
638 
639  This code is distributed under the GNU LGPL license.
640 
641  Modified:
642 
643  30 August 2012
644 
645  Author:
646 
647  John Burkardt
648 
649  Parameters:
650 
651  Input, int DIM_NUM, the spatial dimension.
652 
653  Input, int POINT_NUM, the number of points at which the
654  monomial is to be evaluated.
655 
656  Input, double X[DIM_NUM*POINT_NUM], the point coordinates.
657 
658  Input, int EXPON[DIM_NUM], the exponents.
659 
660  Output, double MONOMIAL_VALUE[POINT_NUM], the value of the monomial.
661 */
662 {
663  int dim;
664  int point;
665  double *value;
666 
667  value = ( double * ) malloc ( point_num * sizeof ( double ) );
668 
669  for ( point = 0; point < point_num; point++ )
670  {
671  value[point] = 1.0;
672  }
673 
674  for ( dim = 0; dim < dim_num; dim++ )
675  {
676  if ( 0 != expon[dim] )
677  {
678  for ( point = 0; point < point_num; point++ )
679  {
680  value[point] = value[point] * pow ( x[dim+point*dim_num], expon[dim] );
681  }
682  }
683  }
684 
685  return value;
686 }
687 /******************************************************************************/
688 
689 double r8_abs ( double x )
690 
691 /******************************************************************************/
692 /*
693  Purpose:
694 
695  R8_ABS returns the absolute value of an R8.
696 
697  Licensing:
698 
699  This code is distributed under the GNU LGPL license.
700 
701  Modified:
702 
703  07 May 2006
704 
705  Author:
706 
707  John Burkardt
708 
709  Parameters:
710 
711  Input, double X, the quantity whose absolute value is desired.
712 
713  Output, double R8_ABS, the absolute value of X.
714 */
715 {
716  double value;
717 
718  if ( 0.0 <= x )
719  {
720  value = + x;
721  }
722  else
723  {
724  value = - x;
725  }
726  return value;
727 }
728 /******************************************************************************/
729 
730 double r8_factorial ( int n )
731 
732 /******************************************************************************/
733 /*
734  Purpose:
735 
736  R8_FACTORIAL computes the factorial of N.
737 
738  Discussion:
739 
740  factorial ( N ) = product ( 1 <= I <= N ) I
741 
742  Licensing:
743 
744  This code is distributed under the GNU LGPL license.
745 
746  Modified:
747 
748  26 June 2008
749 
750  Author:
751 
752  John Burkardt
753 
754  Parameters:
755 
756  Input, int N, the argument of the factorial function.
757  If N is less than 1, the function value is returned as 1.
758 
759  Output, double R8_FACTORIAL, the factorial of N.
760 */
761 {
762  int i;
763  double value;
764 
765  value = 1.0;
766 
767  for ( i = 1; i <= n; i++ )
768  {
769  value = value * ( double ) ( i );
770  }
771 
772  return value;
773 }
774 /******************************************************************************/
775 
776 double r8vec_dot_product ( int n, double a1[], double a2[] )
777 
778 /******************************************************************************/
779 /*
780  Purpose:
781 
782  R8VEC_DOT_PRODUCT computes the dot product of a pair of R8VEC's.
783 
784  Licensing:
785 
786  This code is distributed under the GNU LGPL license.
787 
788  Modified:
789 
790  26 July 2007
791 
792  Author:
793 
794  John Burkardt
795 
796  Parameters:
797 
798  Input, int N, the number of entries in the vectors.
799 
800  Input, double A1[N], A2[N], the two vectors to be considered.
801 
802  Output, double R8VEC_DOT_PRODUCT, the dot product of the vectors.
803 */
804 {
805  int i;
806  double value;
807 
808  value = 0.0;
809  for ( i = 0; i < n; i++ )
810  {
811  value = value + a1[i] * a2[i];
812  }
813  return value;
814 }
815 /******************************************************************************/
816 
817 double *r8vec_uniform_01_new ( int n, int *seed )
818 
819 /******************************************************************************/
820 /*
821  Purpose:
822 
823  R8VEC_UNIFORM_01_NEW returns a unit pseudorandom R8VEC.
824 
825  Discussion:
826 
827  This routine implements the recursion
828 
829  seed = 16807 * seed mod ( 2^31 - 1 )
830  unif = seed / ( 2^31 - 1 )
831 
832  The integer arithmetic never requires more than 32 bits,
833  including a sign bit.
834 
835  Licensing:
836 
837  This code is distributed under the GNU LGPL license.
838 
839  Modified:
840 
841  19 August 2004
842 
843  Author:
844 
845  John Burkardt
846 
847  Reference:
848 
849  Paul Bratley, Bennett Fox, Linus Schrage,
850  A Guide to Simulation,
851  Second Edition,
852  Springer, 1987,
853  ISBN: 0387964673,
854  LC: QA76.9.C65.B73.
855 
856  Bennett Fox,
857  Algorithm 647:
858  Implementation and Relative Efficiency of Quasirandom
859  Sequence Generators,
860  ACM Transactions on Mathematical Software,
861  Volume 12, Number 4, December 1986, pages 362-376.
862 
863  Pierre L'Ecuyer,
864  Random Number Generation,
865  in Handbook of Simulation,
866  edited by Jerry Banks,
867  Wiley, 1998,
868  ISBN: 0471134031,
869  LC: T57.62.H37.
870 
871  Peter Lewis, Allen Goodman, James Miller,
872  A Pseudo-Random Number Generator for the System/360,
873  IBM Systems Journal,
874  Volume 8, Number 2, 1969, pages 136-143.
875 
876  Parameters:
877 
878  Input, int N, the number of entries in the vector.
879 
880  Input/output, int *SEED, a seed for the random number generator.
881 
882  Output, double R8VEC_UNIFORM_01_NEW[N], the vector of pseudorandom values.
883 */
884 {
885  int i;
886  int i4_huge = 2147483647;
887  int k;
888  double *r;
889 
890  if ( *seed == 0 )
891  {
892  fprintf ( stderr, "\n" );
893  fprintf ( stderr, "R8VEC_UNIFORM_01_NEW - Fatal error!\n" );
894  fprintf ( stderr, " Input value of SEED = 0.\n" );
895  exit ( 1 );
896  }
897 
898  r = ( double * ) malloc ( n * sizeof ( double ) );
899 
900  for ( i = 0; i < n; i++ )
901  {
902  k = *seed / 127773;
903 
904  *seed = 16807 * ( *seed - k * 127773 ) - k * 2836;
905 
906  if ( *seed < 0 )
907  {
908  *seed = *seed + i4_huge;
909  }
910 
911  r[i] = ( double ) ( *seed ) * 4.656612875E-10;
912  }
913 
914  return r;
915 }
916 /******************************************************************************/
917 
918 double simplex_unit_monomial_int ( int dim_num, int expon[] )
919 
920 /******************************************************************************/
921 /*
922  Purpose:
923 
924  SIMPLEX_UNIT_MONOMIAL_INT integrates a monomial over a simplex.
925 
926  Discussion:
927 
928  This routine evaluates a monomial of the form
929 
930  product ( 1 <= dim <= dim_num ) x(dim)^expon(dim)
931 
932  where the exponents are nonnegative integers. Note that
933  if the combination 0^0 is encountered, it should be treated
934  as 1.
935 
936  Licensing:
937 
938  This code is distributed under the GNU LGPL license.
939 
940  Modified:
941 
942  30 August 2012
943 
944  Author:
945 
946  John Burkardt
947 
948  Parameters:
949 
950  Input, int DIM_NUM, the spatial dimension.
951 
952  Input, int EXPON[DIM_NUM], the exponents.
953 
954  Output, double SIMPLEX_UNIT_MONOMIAL_INT, the value of the integral
955  of the monomial.
956 */
957 {
958  int dim;
959  int i;
960  int k;
961  double value;
962 
963  value = 1.0;
964 
965  k = 0;
966 
967  for ( dim = 0; dim < dim_num; dim++ )
968  {
969  for ( i = 1; i <= expon[dim]; i++ )
970  {
971  k = k + 1;
972  value = value * ( double ) ( i ) / ( double ) ( k );
973  }
974  }
975 
976  for ( dim = 0; dim < dim_num; dim++ )
977  {
978  k = k + 1;
979  value = value / ( double ) ( k );
980  }
981 
982  return value;
983 }
984 /******************************************************************************/
985 
986 double simplex_unit_monomial_quadrature ( int dim_num, int expon[],
987  int point_num, double x[], double w[] )
988 
989 /******************************************************************************/
990 /*
991  Purpose:
992 
993  SIMPLEX_UNIT_MONOMIAL_QUADRATURE: quadrature of monomials in a unit simplex.
994 
995  Licensing:
996 
997  This code is distributed under the GNU LGPL license.
998 
999  Modified:
1000 
1001  30 August 2012
1002 
1003  Author:
1004 
1005  John Burkardt
1006 
1007  Parameters:
1008 
1009  Input, int DIM_NUM, the spatial dimension.
1010 
1011  Input, int EXPON[DIM_NUM], the exponents.
1012 
1013  Input, int POINT_NUM, the number of points in the rule.
1014 
1015  Input, double X[DIM_NUM*POINT_NUM], the quadrature points.
1016 
1017  Input, double W[POINT_NUM], the quadrature weights.
1018 
1019  Output, double SIMPLEX_UNIT_MONOMIAL_QUADRATURE, the quadrature error.
1020 */
1021 {
1022  double exact = 1.0;
1023  double quad;
1024  double quad_error;
1025  double scale;
1026  double *value;
1027  double volume;
1028 /*
1029  Get the exact value of the integral of the unscaled monomial.
1030 */
1031  scale = simplex_unit_monomial_int ( dim_num, expon );
1032 /*
1033  Evaluate the monomial at the quadrature points.
1034 */
1035  value = monomial_value ( dim_num, point_num, x, expon );
1036 /*
1037  Compute the weighted sum and divide by the exact value.
1038 */
1039  volume = simplex_unit_volume ( dim_num );
1040  quad = volume * r8vec_dot_product ( point_num, w, value ) / scale;
1041 
1042  free ( value );
1043 /*
1044  Error:
1045 */
1046  quad_error = r8_abs ( quad - exact );
1047 
1048  return quad_error;
1049 }
1050 /******************************************************************************/
1051 
1052 double *simplex_unit_sample ( int dim_num, int n, int *seed )
1053 
1054 /******************************************************************************/
1055 /*
1056  Purpose:
1057 
1058  SIMPLEX_UNIT_SAMPLE returns uniformly random points from a general simplex.
1059 
1060  Discussion:
1061 
1062  The interior of the unit DIM_NUM dimensional simplex is the set of
1063  points X(1:DIM_NUM) such that each X(I) is nonnegative, and
1064  sum(X(1:DIM_NUM)) <= 1.
1065 
1066  This routine is valid for any spatial dimension DIM_NUM.
1067 
1068  Licensing:
1069 
1070  This code is distributed under the GNU LGPL license.
1071 
1072  Modified:
1073 
1074  30 August 2012
1075 
1076  Author:
1077 
1078  John Burkardt
1079 
1080  Reference:
1081 
1082  Reuven Rubinstein,
1083  Monte Carlo Optimization, Simulation, and Sensitivity
1084  of Queueing Networks,
1085  Krieger, 1992,
1086  ISBN: 0894647644,
1087  LC: QA298.R79.
1088 
1089  Parameters:
1090 
1091  Input, int DIM_NUM, the dimension of the space.
1092 
1093  Input, int N, the number of points.
1094 
1095  Input/output, int *SEED, a seed for the random number generator.
1096 
1097  Output, double UNIFORM_IN_SIMPLEX01_MAP[DIM_NUM*N], the points.
1098 */
1099 {
1100  double *e;
1101  int i;
1102  int j;
1103  double total;
1104  double *x;
1105 /*
1106  The construction begins by sampling DIM_NUM+1 points from the
1107  exponential distribution with parameter 1.
1108 */
1109  x = ( double * ) malloc ( dim_num * n * sizeof ( double ) );
1110 
1111  for ( j = 0; j < n; j++ )
1112  {
1113  e = r8vec_uniform_01_new ( dim_num+1, seed );
1114 
1115  for ( i = 0; i <= dim_num; i++ )
1116  {
1117  e[i] = -log ( e[i] );
1118  }
1119 
1120  total = 0.0;
1121  for ( i = 0; i <= dim_num; i++ )
1122  {
1123  total = total + e[i];
1124  }
1125 
1126  for ( i = 0; i < dim_num; i++ )
1127  {
1128  x[i+dim_num*j] = e[i] / total;
1129  }
1130  free ( e );
1131  }
1132 
1133  return x;
1134 }
1135 /******************************************************************************/
1136 
1137 double *simplex_unit_to_general ( int dim_num, int point_num, double t[],
1138  double ref[] )
1139 
1140 /******************************************************************************/
1141 /*
1142  Purpose:
1143 
1144  SIMPLEX_UNIT_TO_GENERAL maps the unit simplex to a general simplex.
1145 
1146  Discussion:
1147 
1148  Given that the unit simplex has been mapped to a general simplex
1149  with vertices T, compute the images in T, under the same linear
1150  mapping, of points whose coordinates in the unit simplex are REF.
1151 
1152  The vertices of the unit simplex are listed as suggested in the
1153  following:
1154 
1155  (0,0,0,...,0)
1156  (1,0,0,...,0)
1157  (0,1,0,...,0)
1158  (0,0,1,...,0)
1159  (...........)
1160  (0,0,0,...,1)
1161 
1162  Thanks to Andrei ("spiritualworlds") for pointing out a mistake in the
1163  previous implementation of this routine, 02 March 2008.
1164 
1165  Licensing:
1166 
1167  This code is distributed under the GNU LGPL license.
1168 
1169  Modified:
1170 
1171  30 August 2012
1172 
1173  Author:
1174 
1175  John Burkardt
1176 
1177  Parameters:
1178 
1179  Input, int DIM_NUM, the spatial dimension.
1180 
1181  Input, int POINT_NUM, the number of points to transform.
1182 
1183  Input, double T[DIM_NUM*(DIM_NUM+1)], the vertices of the
1184  general simplex.
1185 
1186  Input, double REF[DIM_NUM*POINT_NUM], points in the
1187  reference triangle.
1188 
1189  Output, double SIMPLEX_UNIT_TO_GENERAL[DIM_NUM*POINT_NUM],
1190  corresponding points in the physical triangle.
1191 */
1192 {
1193  int dim;
1194  double *phy;
1195  int point;
1196  int vertex;
1197 
1198  phy = ( double * ) malloc ( dim_num * point_num * sizeof ( double ) );
1199 //
1200 // The image of each point is initially the image of the origin.
1201 //
1202 // Insofar as the pre-image differs from the origin in a given vertex
1203 // direction, add that proportion of the difference between the images
1204 // of the origin and the vertex.
1205 //
1206  for ( point = 0; point < point_num; point++ )
1207  {
1208  for ( dim = 0; dim < dim_num; dim++ )
1209  {
1210  phy[dim+point*dim_num] = t[dim+0*dim_num];
1211 
1212  for ( vertex = 1; vertex < dim_num + 1; vertex++ )
1213  {
1214  phy[dim+point*dim_num] = phy[dim+point*dim_num]
1215  + ( t[dim+vertex*dim_num] - t[dim+0*dim_num] ) * ref[vertex-1+point*dim_num];
1216  }
1217  }
1218  }
1219 
1220  return phy;
1221 }
1222 /******************************************************************************/
1223 
1224 double simplex_unit_volume ( int dim_num )
1225 
1226 /******************************************************************************/
1227 /*
1228  Purpose:
1229 
1230  SIMPLEX_UNIT_VOLUME computes the volume of the unit simplex.
1231 
1232  Discussion:
1233 
1234  The formula is simple: volume = 1/DIM_NUM!.
1235 
1236  Licensing:
1237 
1238  This code is distributed under the GNU LGPL license.
1239 
1240  Modified:
1241 
1242  04 September 2003
1243 
1244  Author:
1245 
1246  John Burkardt
1247 
1248  Parameters:
1249 
1250  Input, int DIM_NUM, the dimension of the space.
1251 
1252  Output, double SIMPLEX_UNIT_VOLUME, the volume of the cone.
1253 */
1254 {
1255  int i;
1256  double volume;
1257 
1258  volume = 1.0;
1259  for ( i = 1; i <= dim_num; i++ )
1260  {
1261  volume = volume / ( ( double ) i );
1262  }
1263 
1264  return volume;
1265 }
1266 /******************************************************************************/
1267 
1268 void timestamp ( void )
1269 
1270 /******************************************************************************/
1271 /*
1272  Purpose:
1273 
1274  TIMESTAMP prints the current YMDHMS date as a time stamp.
1275 
1276  Example:
1277 
1278  31 May 2001 09:45:54 AM
1279 
1280  Licensing:
1281 
1282  This code is distributed under the GNU LGPL license.
1283 
1284  Modified:
1285 
1286  24 September 2003
1287 
1288  Author:
1289 
1290  John Burkardt
1291 
1292  Parameters:
1293 
1294  None
1295 */
1296 {
1297 # define TIME_SIZE 40
1298 
1299  static char time_buffer[TIME_SIZE];
1300  const struct tm *tm;
1301  size_t len;
1302  time_t now;
1303 
1304  (void)(len);
1305 
1306  now = time ( NULL );
1307  tm = localtime ( &now );
1308 
1309  len = strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
1310 
1311  fprintf ( stdout, "%s\n", time_buffer );
1312 
1313  return;
1314 # undef TIME_SIZE
1315 }
1316 
int gm_rule_size(int rule, int dim_num)
Definition: gm_rule.c:294
int i4_min(int i1, int i2)
Definition: gm_rule.c:500
void gm_rule_set(int rule, int dim_num, int point_num, double w[], double x[])
Definition: gm_rule.c:152
double r8_abs(double x)
Definition: gm_rule.c:689
double * r8vec_uniform_01_new(int n, int *seed)
Definition: gm_rule.c:817
double r8vec_dot_product(int n, double a1[], double a2[])
Definition: gm_rule.c:776
double simplex_unit_monomial_int(int dim_num, int expon[])
Definition: gm_rule.c:918
constexpr double a2
int i4_choose(int n, int k)
Definition: gm_rule.c:353
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
double r8_factorial(int n)
Definition: gm_rule.c:730
double * simplex_unit_to_general(int dim_num, int point_num, double t[], double ref[])
Definition: gm_rule.c:1137
int i4_max(int i1, int i2)
Definition: gm_rule.c:459
double * monomial_value(int dim_num, int point_num, double x[], int expon[])
Definition: gm_rule.c:619
constexpr double a1
double simplex_unit_volume(int dim_num)
Definition: gm_rule.c:1224
#define TIME_SIZE
double * simplex_unit_sample(int dim_num, int n, int *seed)
Definition: gm_rule.c:1052
double simplex_unit_monomial_quadrature(int dim_num, int expon[], int point_num, double x[], double w[])
Definition: gm_rule.c:986
void comp_next(int n, int k, int a[], int *more, int *h, int *t)
Definition: gm_rule.c:10
void timestamp(void)
Definition: gm_rule.c:1268
int i4_power(int i, int j)
Definition: gm_rule.c:541
int i4_huge(void)
Definition: gm_rule.c:428