Math-NumSeq
view release on metacpan or search on metacpan
devel/balanced-binary.pl view on Meta::CPAN
if ($value < $prev) {
die $i;
}
$prev = $value;
}
print "$prev\n";
exit 0;
}
{
# formula
require Math::NumSeq::Catalan;
my $seq = Math::NumSeq::Catalan->new;
my $cumul = 0;
foreach (1 .. 20) {
my ($i, $value) = $seq->next;
my $formula = 0;
foreach my $k (1 .. $i-1) {
$formula += $seq->ith($i-$k)*$k + $seq->ith($k)
}
print "$i value=$value formula=$formula\n";
$cumul += $value;
}
exit 0;
}
{
# cumulative
require Math::NumSeq::Catalan;
my $seq = Math::NumSeq::Catalan->new;
devel/factorials.pl view on Meta::CPAN
}
exit 0;
}
{
{
package Math::Symbolic::Custom::MySimplification;
use base 'Math::Symbolic::Custom::Simplification';
# use Math::Symbolic::Custom::Pattern;
# my $formula = Math::Symbolic->parse_from_string("TREE_a * (TREE_b / TREE_c)");
# my $pattern = Math::Symbolic::Custom::Pattern->new($formula);
use Math::Symbolic::Custom::Transformation;
my $trafo = Math::Symbolic::Custom::Transformation::Group->new
(
',',
# 'TREE_a * (TREE_b / TREE_c)' => '(TREE_a * TREE_b) / TREE_c',
# 'TREE_a * (TREE_b + TREE_c)' => 'TREE_a * TREE_b + TREE_a * TREE_c',
# '(TREE_b + TREE_c) * TREE_a' => 'TREE_b * TREE_a + TREE_c * TREE_a',
#
# # '(TREE_a / TREE_b) / TREE_c' => 'TREE_a / (TREE_b * TREE_c)',
devel/haferman-carpet.pl view on Meta::CPAN
'digit_split_lowtohigh',
'digit_join_lowtohigh';
use Math::NumSeq::HafermanCarpet;
# uncomment this to run the ### lines
# use Smart::Comments;
{
foreach my $n (0 .. 5) {
my $Seq_1s_init0 = Seq_1s_init0($n);
my $Seq_1s_init0_by_formula = Seq_1s_init0_by_formula($n);
my $Seq_1s_init1 = Seq_1s_init1($n);
my $Seq_1s_init1_by_formula = Seq_1s_init1_by_formula($n);
printf("%d %6d %6d%s %6d %6d%s %6d %6d\n",
$n,
$Seq_1s_init0,
$Seq_1s_init0_by_formula,
$Seq_1s_init0 == $Seq_1s_init0_by_formula ? '' : '******',
$Seq_1s_init1,
$Seq_1s_init1_by_formula,
$Seq_1s_init1 == $Seq_1s_init1_by_formula ? '' : '******',
Array_1s_init0($n),
Array_1s_init1($n),
);
}
exit 0;
# num black cells after n iterations (9^(k+1) - (-5)^(k+1))/14 = 1,4,61,424
# Array1s(k+1) = 9^(k+1) - 5*Array1s(k)
# Array1s(0) = 1
devel/haferman-carpet.pl view on Meta::CPAN
# Sones = (9^(n+1) - (-5)^(n+1)) / 14 - (1 + (-1)**$n)/2 * 5^n
# = (9^(n+1) - (-5)^(n+1) - 7*(1 + (-1)^n) * 5^n) / 14
# = (9^(n+1) - (-1)^(n+1)*5*5^n - 7*(1 + (-1)^n) * 5^n) / 14
# = (9^(n+1) + (- (-1)^(n+1)*5 - 7*(1 + (-1)^n)) * 5^n) / 14
# = (9^(n+1) + (- (-1)^(n+1)*5 - 7 - 7*(-1)^n) * 5^n) / 14
# = (9^(n+1) + (5*(-1)^n - 7 + 7*(-1)^n) * 5^n) / 14
# = (9^(n+1) + (-2*(-1)^n - 7) * 5^n) / 14
# = (9^(n+1) - (2*(-1)^n + 7) * 5^n) / 14 # 7+2=9 or 7-2=5
sub Seq_1s_init0_by_formula {
my ($n) = @_;
return (9**($n+1) - (2*(-1)**$n + 7) * 5**$n) / 14;
return Array_1s_init1($n) - (1 + (-1)**$n)/2 * Hbox($n);
}
# S1ones = (9^(n+1) - (2*(-1)^n + 7) * 5^n) / 14 + 5^n
# = (9^(n+1) - (2*(-1)^n + 7) * 5^n + 14 * 5^n) / 14
# = (9^(n+1) - (2*(-1)^n + 7 - 14) * 5^n) / 14
# = (9^(n+1) - (2*(-1)^n - 7) * 5^n) / 14
sub Seq_1s_init1_by_formula {
my ($n) = @_;
return (9**($n+1) - (2*(-1)**$n - 7) * 5**$n) / 14;
return Array_1s_init1($n) - (1 + (-1)**$n)/2 * Hbox($n);
}
BEGIN {
my @count;
my $seq = Math::NumSeq::HafermanCarpet->new (initial_value => 0);
sub Seq_1s_init0 {
devel/lib/Math/NumSeq/Padovan.pm view on Meta::CPAN
The pattern starts from an equilateral triangle of side length 1, the lower
left one shown above. Then put a new triangle successively to the right,
left and lower sides, each time the side length of the triangle is the width
of the figure so far.
The effect is to add the immediately preceding triangle side and the side of
the 5th previous,
P[i] = P[i-1] + P[i-5]
This is the same as the i-2,i-3 formula shown above since
P[i] = P[i-2] + P[i-3] # formula above
= P[i-4]+P[i-5] + P[i-3]
= P[i-3]+P[i-4] + P[i-5]
= P[i-1] + P[i-5] # geometric form
=head1 FUNCTIONS
See L<Math::NumSeq/FUNCTIONS> for behaviour common to all sequence classes.
=over 4
lib/Math/NumSeq/Catalan.pm view on Meta::CPAN
values_type => "odd"
1, 1, 1, 5, 7, 21, 33, 429, 715, 2431, 4199, ... (A098597)
starting i=0
The number of 2s in C(i) is
num2s = (count-1-bits of i+1) - 1
The odd part is always monotonically increasing. When i increments num2s
increases by at most 1, ie. a single factor of 2. In the formula above
C(i) = C(i-1) * 2*(2i-1)/(i+1)
it can be seen that C(i) gains at least 1 factor of 2, so after dividing out
2^num2s it's still greater than C(i-1).
=head1 FUNCTIONS
See L<Math::NumSeq/FUNCTIONS> for behaviour common to all sequence classes.
lib/Math/NumSeq/Fibonacci.pm view on Meta::CPAN
bit = least significant bit of i
if bit == 1
F[2k+1] = (2F[k]+F[k-1])*(2F[k]-F[k-1]) + add
else
F[2k] = F[k]*(F[k]+2F[k-1])
The "add" amount is 2*(-1)^k which means +2 or -2 according to k odd or
even, which means whether the previous bit taken from i was 1 or 0. That
can be easily noted from each bit, to be used in the following loop
iteration or the final step F[2k+1] formula.
For small i it's usually faster to just successively add F[k+1]=F[k]+F[k-1],
but when in bignums the doubling k-E<gt>2k by two squares is faster than
doing k many individual additions for the same thing.
=head2 Value to i Estimate
F[i] increases as a power of phi, the golden ratio. The exact value is
F[i] = (phi^i - beta^i) / (phi - beta) # exactly
lib/Math/NumSeq/GolayRudinShapiroCumulative.pm view on Meta::CPAN
index)
=over
L<http://user42.tuxfamily.org/alternate/index.html>
=back
=cut
# ENHANCE-ME: Cross-ref to AlternatePaper formulas section when its X,Y
# calculation is described.
=pod
=head1 SEE ALSO
L<Math::NumSeq>,
L<Math::NumSeq::GolayRudinShapiro>
L<Math::PlanePath::AlternatePaper>
lib/Math/NumSeq/HafermanCarpet.pm view on Meta::CPAN
The difference between the two is 5^k,
Seq1s_init1 = Seq1s_init0 + 5^k
This difference is the box fractal 1s described above which are gained in
the initial_value=1 form. They're at positions where i has only even digits
0,2,4,6,8 in base 9, so 5 digit possibilities at each of k many digit
positions giving 5^k.
X<Weinstein, Eric>The formulas can be calculated by considering how the 0s
and 1s expand. The array morphism form with initial_value=1 is given by
Eric Weinstein,
=over
L<http://mathworld.wolfram.com/HafermanCarpet.html>,
L<http://oeis.org/A118005>
=back
lib/Math/NumSeq/HafermanCarpet.pm view on Meta::CPAN
Array1s_init0(k) = ------------------
9 - (-5)
For the sequence forms here the initial value does not change, unlike the
array alternating 0E<lt>-E<gt>1. The sequence takes the array starting 0 or
1 according as k is even or odd, thereby ensuring the first value is 0. So,
Seq1s_init0(k) = / Array1s_init0(k) if k even
\ Array1s_init1(k) if k odd
The term S<(2*(-1)^k - 7)> seen above in the Seq1s_init0() formula selects
between 9 or 5 as the multiplier for (-5)^k. That 9 or 5 multiplier is the
difference between the two Array1s forms.
=head1 SEE ALSO
L<Math::NumSeq>
L<Math::PlanePath::ZOrderCurve>,
L<Math::PlanePath::PeanoCurve>,
L<Math::PlanePath::WunderlichSerpentine>,
lib/Math/NumSeq/LucasNumbers.pm view on Meta::CPAN
log(L[i]) ~= i*log(phi)
i ~= log(L[i]) * 1/log(phi)
Or the same using log base 2 which can be estimated from the highest bit
position of a bignum,
log2(L[i]) ~= i*log2(phi)
i ~= log2(L[i]) * 1/log2(phi)
This is very close to the Fibonacci formula (see
L<Math::NumSeq::Fibonacci/Value to i Estimate>), being bigger by
Lestimate(value) - Festimate(value)
= log(value) / log(phi) - (log(value) + log(phi-beta)) / log(phi)
= -log(phi-beta) / log(phi)
= -1.67
On that basis it could even be close enough to take Lestimate = Festimate-1,
or vice-versa.
lib/Math/NumSeq/Polygonal.pm view on Meta::CPAN
The letters "a", "b" "c" show the extra added onto the previous figure to
grow its points. Each side except two are extended. In general the
k-gonals increment by k-2 sides of i points, plus 1 at the end of the last
side, so
P(i+1) = P(i) + (k-2)*i + 1
=head2 Second Kind
Option C<pairs =E<gt> 'second'> gives the polygonals of the second kind,
which are the same formula but with a negative i.
S(i) = P(-i) = (k-2)/2 * i*(i-1) + (k-3)*i
The result is still positive values, bigger than the plain P(i). For
example the pentagonals are 0,1,5,12,22,etc and the second pentagonals are
0,2,7,15,26,etc.
=head2 Both Kinds
C<pairs =E<gt> 'both'> gives the firsts and seconds interleaved. P(0) and
lib/Math/NumSeq/RadixWithoutDigit.pm view on Meta::CPAN
if ($digit == 0) {
my $len = 1;
while (($r1**$len - $r1) / $r2 <= $i) {
$len++;
}
$len--;
### $len
### base: ($r1**$len - $r1) / $r2
### sentinel: $r1**$len
### adj add: $r1**$len - ($r1**$len - $r1) / $r2
### adj formula: (($r2-1)*$r1**$len + $r1) / $r2
### assert: $i >= ($r1 ** $len - $r1) / $r2
### assert: $i < ($r1 ** $len - $r1) / $r2 + $r1 ** $len
$i += (($r2-1)*$r1**$len + $r1) / $r2;
$value = -2 * $radix**$len;
### i remainder: $i
### $value
### assert: $i >= $r1 ** $len
lib/Math/NumSeq/ReRound.pm view on Meta::CPAN
# i = (-1 + sqrt(1 + 4*2v/m)) / 2
# = -1/2 + sqrt(1 + 4*2v/m)/2
# = -1/2 + sqrt(4 + 2v/m)
# i -> sqrt(2v/m)
#
# as extra_multiples dominates the estimate changes from
# extra_multiples==0 est = sqrt(pi*value)
# up to
# extra_multiples==m large est = sqrt(2*value/m)
#
# What formula that progressively morphs pi to 2/m ?
#
sub value_to_i_estimate {
my ($self, $value) = @_;
if ($value < 0) { return 0; }
if ($self->{'extra_multiples'} == 0) {
# use sqrt(pi)~=296/167 to preserve Math::BigInt
return int( (sqrt(int($value)) * 296) / 167);
} else {
return int(sqrt(int (2 * $value / $self->{'extra_multiples'})));
lib/Math/NumSeq/ReRound.pm view on Meta::CPAN
value =~ m * i*(i+1)/2 for large m
i =~ sqrt(value*2/m)
As m increases the factor sqrt(pi) progressively morphs towards sqrt(2/m).
For m even it might be sqrt(pi)/2, 3/8*sqrt(pi), 5/16*sqrt(pi),
35/128*sqrt(pi), etc, and for m odd it might be 2*sqrt(1/pi),
4/3*sqrt(1/pi), 16/15*sqrt(1/pi), 32/35*sqrt(1/pi), 256/315*sqrt(1/pi), etc.
What's the pattern?
The current code uses the "large m" formula for any mE<gt>0, which is no
more than roughly a factor 1.25 too big.
=head1 SEE ALSO
L<Math::NumSeq>,
L<Math::NumSeq::ReReplace>
=head1 HOME PAGE
L<http://user42.tuxfamily.org/math-numseq/index.html>
lib/Math/NumSeq/SpiroFibonacci.pm view on Meta::CPAN
Value 36 on the right is 31+5, being the immediately preceding 31 and the
value on the next inward loop closest to that new 36 position.
At the corners the same inner value is used three times, so for example
42=36+6, then 48=42+6 and 54=48+6, all using the corner "6". For the
innermost loop SF[2] through SF[7] the "0" at the origin is the inner value,
hence the run of seven 1s at the start.
=head2 Absolute Differences
Optional C<recurrence_type =E<gt> 'absdiff'> changes the recurrence formula
to an absolute difference
SF[i] = abs (SF[i-1] - SF[i-k])
With the default initial values SF[0]=0 and SF[1]=1 this behaves as an XOR,
always giving 0 or 1.
0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, ...
The result plotted around the square spiral is similar to some of the
lib/Math/NumSeq/Tetrahedral.pm view on Meta::CPAN
sub pred {
my ($self, $value) = @_;
### Tetrahedral pred(): $value
$value *= 6;
my $i = _cbrt_floor($value);
return ($i*($i+1)*($i+2) == $value);
}
# Cubic root formula
# 6*T(i) = i^3 + 3i^2 + 2i
# i^3 + 3i^2 + 2i - value = 0
# subst j=i+1, i=j-1
# (j-1)^3 + 3(j-1)^2 + 2(j-1) - value = 0
# j^3-3j^2+3j-1 + 3j^2-6j+3 + 2j-2 - value = 0
# j^3 + (-3+3)^2 + (3-6+2)j + (-1+3+-2) - value = 0
# j^3 - j - value = 0 p=-1 q=-v
#
# x^3+px+q=0
# x=a-b
xt/oeis/Fibbinary-oeis.t view on Meta::CPAN
my ($count) = @_;
my $seq = Math::NumSeq::Fibbinary->new;
my @got = (0);
for (my $n = 0; @got < $count; $n++) {
my $value = $seq->ith($n+1);
push @got, ($value >> 1) & 1;
}
return \@got;
});
# A188009 -- [nr]-[nr-kr]-[kr]
# second lowest bit of Zeck per Wolfdieter Lang formula A123740
MyOEIS::compare_values
(anum => 'A188009',
func => sub {
my ($count) = @_;
my $seq = Math::NumSeq::Fibbinary->new;
my @got = (0,0,0);
for (my $n = 0; @got < $count; $n++) {
my $value = $seq->ith($n+1);
push @got, ($value >> 1) & 1;
}
xt/oeis/Fibbinary-oeis.t view on Meta::CPAN
# GP-Test got=concat(got,vector(n,k,A346434_T(n,k))); \
# GP-Test if(#got>=#v, print(" which is rows 1..",n); break)); \
# GP-Test got==v
#
# GP-DEFINE vector_reps(v,n) = if(n==0,[], concat(vector(n,i,v)));
# GP-Test /* comment 10s and 01s */ \
# GP-Test vector(5,n, vector(n,k, A346434_T(n,k))) == \
# GP-Test vector(5,n, vector(n,k, \
# GP-Test fromdigits(concat(vector_reps([1,0],k), vector_reps([0,1],n-k)))))
#
# GP-Test /* formula */ \
# GP-Test vector(50,n, vector(n,k, A346434_T(n,k))) == \
# GP-Test vector(50,n, vector(n,k, (10*100^n - 9*100^(n-k) - 1)/99 ))
#
# GP-Test /* formula */ \
# GP-Test vector(50,n, vector(n,k, A346434_T(n,k))) == \
# GP-Test vector(50,n, vector(n,k, A014417_to_Zeckendorf(A210619_Zeckendorf_euqal_0s1s_T(n,k)) ))
# GP-Test /* example table */ \
# GP-Test my(n=1); vector(1,k, A346434_T(n,k)) == [ 10 ]
# GP-Test my(n=2); vector(2,k, A346434_T(n,k)) == [ 1001, 1010 ]
# GP-Test my(n=3); vector(3,k, A346434_T(n,k)) == [ 100101, 101001, 101010 ]
# GP-Test my(n=4); vector(4,k, A346434_T(n,k)) == [ 10010101, 10100101, 10101001, 10101010 ]
# GP-Test A346434_T(5,3) == 1010100101
#
# diagonal
xt/oeis/Fibbinary-oeis.t view on Meta::CPAN
# GP-DEFINE A255773(n) = {
# GP-DEFINE my(x=0,y=0);
# GP-DEFINE for(i=0,logint(n,2),
# GP-DEFINE [x,y]=[y,x+y+1];
# GP-DEFINE if(bittest(n,i), [x,y]=[y,x+y]));
# GP-DEFINE y;
# GP-DEFINE }
# GP-Test my(v=OEIS_data("A255773")); /* OFFSET=1 */ \
# GP-Test vector(#v,n, A255773(n)) == v
# GP-Test /* formula */ \
# GP-Test vector(1024,n, A255773(n)) == \
# GP-Test vector(1024,n, A095903(A092754(n)))
#
# GP-Test vector(1024,n, A255773(n)) == \
# GP-Test vector(1024,n, A095903(insert_highbit(n,0)-1))
# GP-Test vector(1024,n, A255773(n)) == \
# GP-Test vector(1024,n, A095903(A004754_insert_high_0(n)-1))
# GP-Test A095903(1) == 1
# GP-Test A095903(2) == 2
# GP-Test A095903(3) == 3
# GP-DEFINE A255774(n) = {
# GP-DEFINE my(x=0,y=0);
# GP-DEFINE for(i=0,logint(n,2),
# GP-DEFINE y++; [x,y]=[y,x+y];
# GP-DEFINE if(bittest(n,i), [x,y]=[y,x+y]));
# GP-DEFINE y;
# GP-DEFINE }
# GP-Test my(v=OEIS_data("A255774")); /* OFFSET=1 */ \
# GP-Test vector(#v,n, A255774(n)) == v
# GP-Test /* formula */ \
# GP-Test vector(1024,n, A255774(n)) == \
# GP-Test vector(1024,n, A095903(A206332(n)))
#
# GP-Test vector(1024,n, A255774(n)) == \
# GP-Test vector(1024,n, A095903(insert_highbit(n,1)-1))
# GP-Test vector(1024,n, A255774(n)) == \
# GP-Test vector(1024,n, A095903(A004755_insert_high_0(n)-1))
#------
xt/oeis/SternDiatomic-oeis.t view on Meta::CPAN
# GP-Test vector(8,n, A337277_row(n))
# GP-DEFINE \\ code by MFH in A002487
# GP-DEFINE A002487(n, a=1, b=0) = {
# GP-DEFINE if(n==0,return(0));
# GP-DEFINE for(i=0, logint(n, 2), if(bittest(n, i), b+=a, a+=b)); b;
# GP-DEFINE }
# GP-Test my(v=OEIS_data("A002487")); /* OFFSET=0 */ \
# GP-Test vector(#v,n,n--; A002487(n)) == v
#
# GP-Test /* formula by Alois P. Heinz in A337277 */ \
# GP-Test vector(100,n,n--; A337277_T(n,n)) == \
# GP-Test vector(100,n,n--; A002487(n+1))
# GP-Test A002487(0) == 0
# GP-Test /* each row the diatomic sequence forward and backward */ \
# GP-Test vector(8,n, A337277_row(n)) == \
# GP-Test vector(8,n, vector(2^(n+1)-1,k,k--; \
# GP-Test if(k<2^n, A002487(k+1), A002487(2^(n+1)-1-k))))
( run in 0.391 second using v1.01-cache-2.11-cpan-26ccb49234f )