view release on metacpan or  search on metacpan

devel/rationals-tree.pl  view on Meta::CPAN

  my $seq = Math::NumSeq::PlanePathTurn->new (planepath_object => $path,
                                              turn_type => 'Right');
  for (my $n = $seq->i_start; $n <= 16384; $n+=1) {
    # next if $n % 2;
    if (is_pow2($n)) {
      printf "\n%5d ", $n;
    }

    # my $turn = $seq->ith($n);
    my ($x,$y) = $path->n_to_xy($n);
    # my $turn = ($x ^ $y) & 1;
    my $turn = ($x&1) + 2*($y&1);

    # if ($n % 8 == 0) { print " "; }
    print "$turn";
  }
  print "\n";
  exit 0;
}

{
  # HCS vs Bird
  require Math::NumSeq::PlanePathTurn;
  my $hcs  = Math::PlanePath::RationalsTree->new(tree_type => 'HCS');
  my $bird = Math::PlanePath::RationalsTree->new(tree_type => 'Bird');
  my $n = 0b1000000010000000000;
  my ($x,$y) = $hcs->n_to_xy($n);
  my $nb = $bird->xy_to_n($x,$y);
  printf "%10b\n", $n;
  printf "%10b\n", $nb;
  exit 0;
}

{
  # Minkowski question mark
  #
  # cf = [0,a1,a2,...] range 0to1
  #             (-1)^(k-1)
  # ? = sum  -----------
  #      k    2^(a1+...+ak-1)
  # (-1)^(1-1)/2^a1 = 1/2^a1 = 0.000..001 binary

  # + (-1)^(1-2)/2^(a1+a2) = -1/2^(a1+a2)
  #   = 0.0001 - 0.000000001
  #   = 0.000011111
  #
  # 0to1 cf = [0,a0,a1,...]
  # ? = 2*(1 - 2^-a0 + 2^-(a0+a1) - 2^-(a0+a1+a2) + ...)
  #
  # ? =
  #
  # ?(1/k^n) = 1/2^(k^n-1)
  # ?(0) = 0
  # ?(1/3) = 1/4
  require Math::BaseCnv;
  require Math::BigRat;
  my $path = Math::PlanePath::RationalsTree->new (tree_type => 'SB');

  # ?(1/3)=1/4  ?(1/2)=1/2  ?(2/3)=3/4
  foreach my $xy ('1/3', '1/2', '2/3') {
    my ($x,$y) = split m{/}, $xy;
    try ($x,$y);
  }

  foreach my $n ($path->n_start .. 64) {
    my ($x,$y) = $path->n_to_xy($n);
    try ($x,$y);
  }

  foreach my $xy ('1/3', '1/2', '2/3') {
    my ($x,$y) = split m{/}, $xy;
    try ($x,$y);
  }

  sub try {
    my ($x,$y) = @_;
    require Math::ContinuedFraction;
    my $cfrac = Math::ContinuedFraction->from_ratio($x,$y);
    my $cfrac_str = $cfrac->to_ascii;
    my $n = $path->xy_to_n($x,$y);
    my $nbits = Math::BaseCnv::cnv($n,10,2);
    my $mp = minkowski_by_path($x,$y);
    my $mc = minkowski_by_cfrac($x,$y);
    my $mpstr = to_binary($mp);
    my $mcstr = to_binary($mc);
    print "$x/$y  $nbits  p=$mp c=$mc   $cfrac_str\n";
  }

  # pow=2^level <= N
  # ? = (2*(N-pow) + 1) / pow
  #   = (2N - 2pow + 1) / pow
  #   = (2N+1)/pow - 2pow/pow
  #   = (2N+1)/pow - 2
  #   = 2*((N+1/2)/pow - 1)
  sub minkowski_by_path {
    my ($x,$y) = @_;
    my $n = $path->xy_to_n($x,$y);
    my ($pow,$exp) = round_down_pow($n,2);
    return Math::BigRat->new(2*$n+1) / $pow - 2;
    return Math::BigRat->new(2*($n-$pow) + 1) / $pow;

    return (2*($n-$pow) + 1) / $pow;
    return (2*$pow-1 - $n) / $pow;
    return $n / (2*$pow);
  }

  # q0, q1, ...
  #               1          1           1
  # ? = 2 * (1 - --- * (1 - ---- * (1 - ---- * (... 
  #             2^q0        2^q1        2^q2
  #
  sub minkowski_by_cfrac {
    my ($x,$y) = @_;
    require Math::ContinuedFraction;
    my $cfrac = Math::ContinuedFraction->from_ratio($x,$y);
    my $aref = $cfrac->to_array;  # first to last
    ### $aref
    my $ret = 1;
    foreach my $q (reverse @$aref) {
      $ret = 1 - 1/Math::BigRat->new(2)**$q * $ret;
    }
    return 2*$ret;
  }

  # q0, q1, ...
  #                (-1)^k
  # ? = sum -------------------
  #      k  2^(q0+q1+...qk - 1)
  sub minkowski_by_cfrac_cumul {
    my ($x,$y) = @_;
    require Math::ContinuedFraction;