v0.14.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
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,
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
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
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
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
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
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
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
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
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
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
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
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
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
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
994
995  Licensing:
996
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
1020 */
1021 {
1022  double exact = 1.0;
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 */
1047
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
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
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
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
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
a2
constexpr double a2
Definition: hcurl_check_approx_in_2d.cpp:39
r8_factorial
double r8_factorial(int n)
Definition: gm_rule.c:730
i4_power
int i4_power(int i, int j)
Definition: gm_rule.c:541
r8vec_uniform_01_new
double * r8vec_uniform_01_new(int n, int *seed)
Definition: gm_rule.c:817
i4_max
int i4_max(int i1, int i2)
Definition: gm_rule.c:459
sdf.r
int r
Definition: sdf.py:8
simplex_unit_sample
double * simplex_unit_sample(int dim_num, int n, int *seed)
Definition: gm_rule.c:1052
monomial_value
double * monomial_value(int dim_num, int point_num, double x[], int expon[])
Definition: gm_rule.c:619
comp_next
void comp_next(int n, int k, int a[], int *more, int *h, int *t)
Definition: gm_rule.c:10
timestamp
void timestamp(void)
Definition: gm_rule.c:1268
gm_rule_size
int gm_rule_size(int rule, int dim_num)
Definition: gm_rule.c:294
ContactOps::scale
double scale
Definition: EshelbianContact.hpp:23
gm_rule.h
a1
constexpr double a1
Definition: hcurl_check_approx_in_2d.cpp:38
a
constexpr double a
Definition: approx_sphere.cpp:30
r8_abs
double r8_abs(double x)
Definition: gm_rule.c:689
double
i4_choose
int i4_choose(int n, int k)
Definition: gm_rule.c:353
r8vec_dot_product
double r8vec_dot_product(int n, double a1[], double a2[])
Definition: gm_rule.c:776
h
double h
Definition: photon_diffusion.cpp:60
t
constexpr double t
plate stiffness
Definition: plate.cpp:59
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
convert.n
n
Definition: convert.py:82
simplex_unit_to_general
double * simplex_unit_to_general(int dim_num, int point_num, double t[], double ref[])
Definition: gm_rule.c:1137
TIME_SIZE
#define TIME_SIZE
i4_min
int i4_min(int i1, int i2)
Definition: gm_rule.c:500
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
gm_rule_set
void gm_rule_set(int rule, int dim_num, int point_num, double w[], double x[])
Definition: gm_rule.c:152
simplex_unit_volume
double simplex_unit_volume(int dim_num)
Definition: gm_rule.c:1224
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
i4_huge
int i4_huge(void)
Definition: gm_rule.c:428
sdf_wavy_2d.w
int w
Definition: sdf_wavy_2d.py:6